我有如下问题：
1) 每个个体有六个测量值，具有较大的主体内差异
2）有两组（治疗和控制）
3) 每组由 5 个人组成
4）我想对两组进行显着性检验，以了解组均值是否彼此不同。
数据如下所示：
我已经使用此代码运行了一些模拟，该代码进行 t 测试以比较组均值。组均值是通过取单个均值的均值来计算的。这忽略了主体内的可变性：
n.simulations<-10000
pvals=matrix(nrow=n.simulations,ncol=1)
for(k in 1:n.simulations){
subject=NULL
for(i in 1:10){
subject<-rbind(subject,as.matrix(rep(i,6)))
}
#set.seed(42)
#Sample Subject Means
subject.means<-rnorm(10,100,2)
#Sample Individual Measurements
values=NULL
for(sm in subject.means){
values<-rbind(values,as.matrix(rnorm(6,sm,20)))
}
out<-cbind(subject,values)
#Split into GroupA and GroupB
GroupA<-out[1:30,]
GroupB<-out[31:60,]
#Add effect size to GroupA
GroupA[,2]<-GroupA[,2]+0
colnames(GroupA)<-c("Subject", "Value")
colnames(GroupB)<-c("Subject", "Value")
#Calculate Individual Means and SDS
GroupA.summary=matrix(nrow=length(unique(GroupA[,1])), ncol=2)
for(i in 1:length(unique(GroupA[,1]))){
GroupA.summary[i,1]<-mean(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])
GroupA.summary[i,2]<-sd(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])
}
colnames(GroupA.summary)<-c("Mean","SD")
GroupB.summary=matrix(nrow=length(unique(GroupB[,1])), ncol=2)
for(i in 1:length(unique(GroupB[,1]))){
GroupB.summary[i,1]<-mean(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])
GroupB.summary[i,2]<-sd(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])
}
colnames(GroupB.summary)<-c("Mean","SD")
Summary<-rbind(cbind(1,GroupA.summary),cbind(2,GroupB.summary))
colnames(Summary)[1]<-"Group"
pvals[k]<-t.test(GroupA.summary[,1],GroupB.summary[,1], var.equal=T)$p.value
}
这是绘图的代码：
#Plots
par(mfrow=c(2,2))
boxplot(GroupA[,2]~GroupA[,1], col="Red", main="Group A",
ylim=c(.9*min(out[,2]),1.1*max(out[,2])),
xlab="Subject", ylab="Value")
stripchart(GroupA[,2]~GroupA[,1], vert=T, pch=16, add=T)
#abline(h=mean(GroupA[,2]), lty=2, lwd=3)
for(i in 1:length(unique(GroupA[,1]))){
m<-mean(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])
ci<-t.test(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])$conf.int[1:2]
points(i-.2,m, pch=15,cex=1.5, col="Grey")
segments(i-.2,
ci[1],i-.2,
ci[2], lwd=4, col="Grey"
)
}
legend("topleft", legend=c("Individual Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")
boxplot(GroupB[,2]~GroupB[,1], col="Light Blue", main="Group B",
ylim=c(.9*min(out[,2]),1.1*max(out[,2])),
xlab="Subject", ylab="Value")
stripchart(GroupB[,2]~GroupB[,1], vert=T, pch=16, add=T)
#abline(h=mean(GroupB[,2]), lty=2, lwd=3)
for(i in 1:length(unique(GroupB[,1]))){
m<-mean(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])
ci<-t.test(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])$conf.int[1:2]
points(i-.2,m, pch=15,cex=1.5, col="Grey")
segments(i-.2,
ci[1],i-.2,
ci[2], lwd=4, col="Grey"
)
}
legend("topleft", legend=c("Individual Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")
boxplot(Summary[,2]~Summary[,1], col=c("Red","Light Blue"), xlab="Group", ylab="Average Value",
ylim=c(.9*min(Summary[,2]),1.1*max(Summary[,2])),
main="Individual Averages")
stripchart(Summary[,2]~Summary[,1], vert=T, pch=16, add=T)
points(.9, mean(GroupA.summary[,1]), pch=15,cex=1.5, col="Grey")
segments(.9,
t.test(GroupA.summary[,1])$conf.int[1],.9,
t.test(GroupA.summary[,1])$conf.int[2], lwd=4, col="Grey"
)
points(1.9, mean(GroupB.summary[,1]), pch=15,cex=1.5, col="Grey")
segments(1.9,
t.test(GroupB.summary[,1])$conf.int[1],1.9,
t.test(GroupB.summary[,1])$conf.int[2], lwd=4, col="Grey"
)
legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")
hist(pvals, breaks=seq(0,1,by=.05), col="Grey",
main=c(paste("# sims=", n.simulations),
paste("% Sig p-values=",100*length(which(pvals<0.05))/length(pvals)))
)
现在，在我看来，因为每个单独的均值本身就是一个估计值，所以我们对组均值的确定性应该低于上图中左下方面板所示的 95% 置信区间所示。因此，计算的 p 值低估了真实的可变性，如果我们希望推断未来的数据，应该会导致假阳性增加。
那么分析这些数据的正确方法是什么？
奖金：
上面的例子是一个简化。对于实际数据：
1) 受试者内方差与均值正相关。
2) 值只能是 2 的倍数。
3）个别结果不是大致正态分布。它们受到零地板效应的影响，并且在正端有长尾。
4) 每组的科目数不一定相等。
以前的文献使用 t 检验忽略了受试者内部的变异性和其他细微差别，就像对上述模拟所做的那样。这些结果可靠吗？如果我可以从图中提取一些平均值和标准误差，我将如何计算“正确”的 p 值。
编辑：
好的，这是实际数据的样子。还有三组而不是两组：
dput() 的数据：
structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3,
3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6,
6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10,
10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12,
12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15,
15, 15, 15, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18,
18, 18, 18, 18, 18, 2, 0, 16, 2, 16, 2, 8, 10, 8, 6, 4, 4, 8,
22, 12, 24, 16, 8, 24, 22, 6, 10, 10, 14, 8, 18, 8, 14, 8, 20,
6, 16, 6, 6, 16, 4, 2, 14, 12, 10, 4, 10, 10, 8, 4, 10, 16, 16,
2, 8, 4, 0, 0, 2, 16, 10, 16, 12, 14, 12, 8, 10, 12, 8, 14, 8,
12, 20, 8, 14, 2, 4, 8, 16, 10, 14, 8, 14, 12, 8, 14, 4, 8, 8,
10, 4, 8, 20, 8, 12, 12, 22, 14, 12, 26, 32, 22, 10, 16, 26,
20, 12, 16, 20, 18, 8, 10, 26), .Dim = c(108L, 3L), .Dimnames = list(
NULL, c("Group", "Subject", "Value")))
编辑2：
回应 Henrik 的回答：因此，如果我改为对单个平均值执行 anova，然后执行 TukeyHSD 程序，如下所示，我可以将其解释为将我的 p 值低估了大约 3-4 倍？
我对这部分问题的目标是了解我作为期刊文章的读者，如何根据他们选择的分析方法更好地解释以前的结果。例如，他们有那些“权威之星”向我显示 0.01>p>.001。因此，如果我接受 0.05 作为合理的截止值，我应该接受他们的解释吗？唯一的附加信息是平均值和 SEM。
#Get Invidual Means
summary=NULL
for(i in unique(dat[,2])){
sub<-which(dat[,2]==i)
summary<-rbind(summary,cbind(
dat[sub,1][3],
dat[sub,2][4],
mean(dat[sub,3]),
sd(dat[sub,3])
)
)
}
colnames(summary)<-c("Group","Subject","Mean","SD")
TukeyHSD(aov(summary[,3]~as.factor(summary[,1])+ (1|summary[,2])))
# Tukey multiple comparisons of means
# 95% family-wise confidence level
#
# Fit: aov(formula = summary[, 3] ~ as.factor(summary[, 1]) + (1 | summary[, 2]))
#
# $`as.factor(summary[, 1])`
# diff lwr upr p adj
# 2-1 -0.672619 -4.943205 3.597967 0.9124024
# 3-1 7.507937 1.813822 13.202051 0.0098935
# 3-2 8.180556 2.594226 13.766885 0.0046312
编辑3： 我认为我们正在接近我的理解。这是@Stephane 的评论中描述的模拟：
#Get Subject Means
means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)
#Initialize "dat2" dataframe
dat2<-dat
#Initialize within-Subject sd
s<-.001
pvals=matrix(nrow=10000,ncol=2)
for(j in 1:10000){
#Sample individual measurements for each subject
temp=NULL
for(i in 1:nrow(means)){
temp<-c(temp,rnorm(6,means[i,3], s))
}
#Set new values
dat2[,3]<-temp
#Take means of sampled values and fit to model
dd2 <- aggregate(Value~Group+Subject, data=dat2, FUN=mean)
fit2 <- lm(Value~Group, data=dd2)
#Save sd and pvalue
pvals[j,]<-cbind(s,anova(fit2)[[5]][5])
#Update sd
s<-s+.001
}
plot(pvals[,1],pvals[,2], xlab="Within-Subject SD", ylab="P-value")
