14

catsデータセット（-MASS-パッケージから）からブートストラップサンプルを作成するためのスクリプトを作成しています。

Davidson and Hinkleyの教科書[1]に従って、単純な線形回帰を実行し、iid観測からのブートストラップ、つまりペアのリサンプリングに基本的なノンパラメトリック手法を採用しました。

元のサンプルは次の形式です。

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

単変量線形モデルを通じて、猫の炉床重量を脳重量で説明したいと思います。

コードは次のとおりです。

library(MASS)
library(boot)


##################
#   CATS MODEL   #
##################

cats.lm <- glm(Hwt ~ Bwt, data=cats)
cats.diag <- glm.diag.plots(cats.lm, ret=T)


#######################
#   CASE resampling   #
#######################

cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt)) 
statistic.coef <- function(data, i) cats.fit(data[i,]) 

bootl <- boot(data=cats, statistic=statistic.coef, R=999)

ここで、クラスタリング変数が存在するとしますcluster = 1, 2,..., 24（たとえば、各猫は特定の同腹子に属しています）。簡単にするために、データのバランスが取れていると仮定します。クラスターごとに6つの観測値があります。したがって、24匹の同腹子のそれぞれは、6匹の猫（すなわちn_cluster = 6、n = 144）で構成されています。

次の方法で偽のcluster変数を作成できます。

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)

関連する質問が2つあります。

（クラスター化された）データセット構造に従ってサンプルをシミュレートする方法は？つまり、クラスターレベルでリサンプリングする方法は？置換を使用してクラスターをサンプリングし、選択した各クラスター内の観測値を元のデータセットと同じように設定します（つまり、クラスターを置換してサンプリングし、各クラスター内の観測値を置換せずにサンプリングします）。

これは、Davidson（p。100）によって提案された戦略です。B = 100サンプルを描画するとします。それらのそれぞれは、24のおそらく再発するクラスター（たとえばcluster = 3, 3, 1, 4, 12, 11, 12, 5, 6, 8, 17, 19, 10, 9, 7, 7, 16, 18, 24, 23, 11, 15, 20, 1）で構成されている必要があり、各クラスターには、元のデータセットの同じ6つの観測値が含まれている必要があります。でそれを行う方法はR？（パッケージの有無にかかわらず-boot-。）続行するための代替案はありますか？

2番目の質問は、初期回帰モデルに関するものです。クラスターレベルの切片を使用した固定効果モデルを採用するとします。採用したリサンプリング手順は変わりますか？

[1] Davidson、AC、Hinkley、DV（1997）。ブートストラップ法とその応用。ケンブリッジ大学出版局。

4

2 に答える 2

7

私があなたを正しく理解していれば、これはあなたが入力としてやろうとしていることですc.data:

  • クラスターを置換してリサンプリングする
  • ランダム サンプル内の各クラスターと元のデータ セット (iecdata) からのそのポイントとの間の関連付けを維持します。
  • サンプリングされたクラスターでブートストラップを作成する

これを実現するスクリプトを次に示します。これを関数にラップして R 回繰り返すことができます。ここで、R はブートストラップの複製数です。

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)

# get a vector with all clusters
c <- sort(unique(c.data$cluster))

# group the data points per cluster
clust.group <- function(c) {
    c.data[c.data$cluster==c,]
}

clust.list <- lapply(c,clust.group)

# resample clusters with replacement
c.sample <- sample(c, replace=T)

clust.sample <- clust.list[c.sample]

clust.size <- 6

# combine the cluster list back to a single data matrix
clust.bind <- function(c) {
    matrix(unlist(c),nrow=clust.size)
}

c.boot <- do.call(rbind,lapply(clust.sample,clust.bind))

# Just to maintain columns name
colnames(c.boot) <- names(c.data)

# the new data set (single bootstrap replicate)
c.boot
于 2013-01-03T01:26:03.597 に答える
0

私は次の方法で問題を解決しようとしました。機能しますが、速度と「優雅さ」の点でおそらく改善される可能性があります。また、可能であれば、-boot-パッケージを使用する方法を見つけることを好みました。これにより、ブートストラップされた信頼区間の数を自動的に計算できるようになりますboot.ci...

簡単にするために、最初のデータセットは、6 つの研究所 (クラスタリング変数) にネストされた 18 匹の猫 (「下位レベル」の観察) で構成されます。データセットはバランスが取れています (n_cluster = 3クラスターごと)。xを説明するためのリグレッサーが 1 つありyます。

結果を保存する偽のデータセットとマトリックスは次のとおりです。

  # fake sample 
  dat <- expand.grid(cat=factor(1:3), lab=factor(1:6))
  dat <- cbind(dat, x=runif(18), y=runif(18, 2, 5))

  # empty matrix for storing coefficients estimates and standard errors of x
  B <- 50 # number of bootstrap samples
  b.sample <- matrix(nrow=B, ncol=3, dimnames=list(c(), c("sim", "b_x", "se_x")))
  b.sample[,1] <- rep(1:B)

各B反復で、次のループは、置換なしでサンプリングされた 3 匹の猫で構成される 6 つのクラスターを置換してサンプリングします (つまり、クラスターの内部構成は変更されずに維持されます)。リグレッサー係数とその標準誤差の推定値は、以前に作成した行列に格納されます。

  ####################################
  #   loop through "b.sample" rows   #
  ####################################

  for (i in seq(1:B)) {

  ###   sampling with replacement from the clustering variable   

    # sampling with replacement from "cluster" 
    cls <- sample(unique(dat$lab), replace=TRUE)
    cls.col <- data.frame(lab=cls)

    # reconstructing the overall simulated sample
    cls.resample <- merge(cls.col, dat, by="lab")


  ###   fitting linear model to simulated data    

    # model fit
    mod.fit <- function(data) glm(data$y ~ data$x)

    # estimated coefficients and standard errors
    b_x <- summary(mod.fit(data=cls.resample))$coefficients[2,1]
    se_x <- summary(mod.fit(data=cls.resample))$coefficients[2,2]

    b.sample[i,2] <- b_x
    b.sample[i,3] <- se_x

  }

最終的なブートストラップ標準誤差は次のとおりです。

 boot_se_x <- sum(b.sample[,3])/(B-1) 
 boot_se_x

-boot-これを行うためにパッケージを採用する方法はありますか?

また、単純な線形回帰の代わりにクラスターレベルでの固定効果の採用に関して、私の主な疑問は、シミュレートされたサンプルの一部で、初期クラスターの一部が選択されていない可能性があることです。特定のインターセプトは識別できません。私が投稿したコードを見ると、「機械的」な観点からは問題にならないはずです (反復ごとに、サンプリングされたクラスターのインターセプトだけで異なる FE モデルを適合させることができます)。

このすべてに統計上の問題があるかどうか疑問に思っていました

于 2013-01-03T01:34:09.223 に答える