1

私はこの質問に関して結論に至らなかったので、言い換えてもう一度質問すると思いました。

データセットを 10,000 回サブサンプリングして、各応答の平均と 95% CI を生成したいと考えています。

以下は、データ セットがどのように構成されているかの例です。

x <- read.table(tc <- textConnection("
study      expt    variable  value1  value2
  1         1         A       1.0      1.1 
  1         2         B       1.1      2.1 
  1         3         B       1.2      2.9
  1         4         C       1.5      2.3 
  2         1         A       1.7      0.3 
  2         2         A       1.9      0.3 
  3         1         A       0.2      0.5"), header = TRUE); close(tc)

各研究/変数の組み合わせを 1 回だけサブサンプリングしたいと思います。たとえば、サブセット化されたデータセットは次のようになります。

study      expt    variable  value1  value2
  1         1         A       1.0      1.1 
  1         2         B       1.1      2.1 
  1         4         C       1.5      2.3 
  2         1         A       1.7      0.3 
  3         1         A       0.2      0.5

行 3 と 6 がなくなっていることに注意してください。どちらも変数を 2 回測定したためです (最初のケースでは B、2 番目のケースでは A)。

サブサンプリングされたデータセットを何度も描画したいので、各変数の 95% CI で value1 と value2 の全体的な平均を導き出すことができます。したがって、サブサンプリングルーチン全体の後の出力は次のようになります。

variable   mean_value1   lower_value1  upper_value1  mean_value2  etc....
   A            2.3           2.0          2.6           2.1
   B            2.5           2.0          3.0           2.5
   C            2.1           1.9          2.3           2.6

サブセットを取得するために必要なコードを次に示します。

 subsample<-function(x, B){
samps<-ddply(x, .(study,variable), nrow)[,3] #for each study/variable combination, 
                                                  #how many experiments are there
expIdx<-which(!duplicated(x$study)) #what is the first row of each study
n<-length(samps) #how many studies are there

sapply(1:B, function(a) { #use sapply for the looping, as it's more efficient than for
    idx<-floor(runif(n, rep(0,n), samps)) #get the experiment number-1 for each study
    x$value[idx+expIdx] #now get a vector of values
})

どんな助けでも大歓迎です。これは複雑だと認識していますので、説明が必要な場合はお知らせください。

4

2 に答える 2

3

研究、実験、変数ごとにデータを分割し、各サブセットにブートストラップを適用します。これを行うには、次のような多くの方法があります。

sdfr <- with(dfr, split(dfr, list(Study, Experiment, Variable)))
sdfr <- Filter(nrow, sdfr)   #to remove empty data frames

lapply(sdfr, function(x) 
{
  boot(x$Response1, statistic = mean, R = 10000, sim = "parametric")
})
于 2011-07-25T15:12:31.417 に答える
2

これが解決策ですが、公正な警告ですが、それはひどくうまくスケーリングすることはなく、私はこの種のスキームの統計的妥当性に気づいていません：

#Replicate your example data
set.seed(1)
dat <- expand.grid(Study = 1:4,Experiment = 1:3, Response = LETTERS[1:4])
dat$Value1 <- runif(48)
dat$Value2 <- runif(48)

#Function to apply to each Response level
#Note the rather inefficient use of ddply 
# in a for loop to do the 'stratified' 
# subsampling you describe
myFun <- function(x,B){
    rs <- matrix(NA,B,2)
    for (i in 1:B){
        temp <- ddply(x,.(Study), .fun = function(x) x[sample(1:nrow(x),1),])
        rs[i,] <- colMeans(temp[,4:5])
    }
    c(Value1 = mean(x$Value1), quantile(rs[,1],probs=c(0.025,0.975)),
            Value2 = mean(x$Value2), quantile(rs[,2],probs=c(0.025,0.975)))
}

ddply(dat,.(Response),.fun = myFun,B=50)

出力例

  Response    Value1      2.5%     97.5%    Value2      2.5%     97.5%
1        A 0.4914725 0.2721876 0.8311799 0.4600546 0.2596446 0.6909686
2        B 0.5941457 0.4018281 0.8047503 0.5241470 0.2865285 0.7099486
3        C 0.4596998 0.2752685 0.6340614 0.5761497 0.3546133 0.8115933
4        D 0.5550651 0.2717772 0.7298913 0.4645609 0.1868757 0.7985816
于 2011-07-26T14:53:57.760 に答える