5

目标：从重复测量方差分析中计算 eta 值的 95% 置信区间。

该设计是一个双因子设计（因子 1 有 3 个水平，因子 2 有 7 个水平），它与 30 个受试者完全交叉。但是，我不确定我是否正确地进行了引导来计算 eta 的置信区间。我使用 boot() 函数按主题引导样本，但不考虑因素和因素中的水平。因此，我不确定我是否正确执行此操作-我是否需要进行更复杂的引导/重新采样，按因子/级别重新采样，或者仅在主题级别上进行是否可以。我的代码似乎给出了合理的结果......

library(ez)
library(boot)
library(reshape2)

###create a data.frame for a 2-factor (Factor1-3 levels, Factor2-7 levels) fully crossed design with 30 subjects & fill fake data values
subject.number<-factor(rep(1:30,each=21))
factor1.levels<-rep(rep(c("level1","level2","level3"),each=7),30)
factor2.levels<-rep(rep(c("level1","level2","level3","level4","level5","level6","level7"),3),30)
set.seed(1234)
fake.data<-rnorm(630,mean=3)
dframe<-data.frame(subject.number,factor1.levels,factor2.levels,fake.data)
names(dframe)<-c("Subject","Factor1","Factor2","OutcomeValue")


###to work with boot() convert from long to wide format
dframe.wide<-dcast(dframe,Subject~Factor1+Factor2,value.var="OutcomeValue")


###function to use with boot() to calculate generalized eta value for Factor1, Factor2, and Factor1xFactor2 interaction in a repeated measures ANOVA

generalized_eta<-function(data,indices){
  d.wide<-data[indices,] #use boot() indices to sample data

  #now that have used indices from boot(), convert data back to long with correct Factor labeling
  dframe.long<-melt(d.wide,value.name="OutcomeValue",id="Subject")
  dframe.long<-cbind(dframe.long,colsplit(dframe.long$variable,"_",c("Factor1","Factor2")))
  dframe.long$Factor1<-factor(dframe.long$Factor1)
  dframe.long$Factor2<-factor(dframe.long$Factor2)
  dframe.long$Subject<-factor(dframe.long$Subject)

  #do repeated measures ANOVA with ezANOVA() which calculates generalized eta
  aov.ez = ezANOVA(data = dframe.long, dv = .(OutcomeValue), wid = .(Subject), within = .(Factor1,Factor2), type = 1)
  #return the three generalized eta values - Factor1, Factor2, Factor1xFactor2
  return(aov.ez[[1]]$ges)
}


###call boot() to do the bootstrap - only 200 to make it fast
results<-boot(data=dframe.wide,statistic=generalized_eta,R=200)


###plot the bootstrap results
plot(results,index=1) #for Factor1
plot(results,index=2) #for Factor2
plot(results,index=3) #for Factor1xFactor2

###create 95%-CI from bootstrap results
boot.ci(results,type="bca",index=1) #for Factor1
boot.ci(results,type="bca",index=2) #for Factor2
boot.ci(results,type="bca",index=3) #for Factor3
4

0 回答 0