杰罗姆·康菲尔德写道：
费舍尔革命最好的成果之一是随机化的想法，而在其他事情上几乎没有同意的统计学家至少同意这一点。但是，尽管达成了这一协议，并且尽管随机分配程序在临床和其他形式的实验中得到了广泛使用，但其逻辑状态，即它所执行的确切功能，仍然是模糊的。
麦田，杰罗姆 (1976)。“最近对临床试验的方法论贡献”。美国流行病学杂志 104 (4): 408–421。
在整个网站和各种文献中，我始终看到关于随机化力量的自信主张。诸如“它消除了混杂变量的问题”之类的强术语很常见。例如，请参见此处。然而，出于实际/道德原因，很多时候实验都是用小样本（每组 3-10 个样本）进行的。这在使用动物和细胞培养物的临床前研究中非常常见，研究人员通常报告 p 值以支持他们的结论。
这让我想知道，随机化在平衡混淆方面有多好。对于这个图，我模拟了一个比较治疗组和对照组的情况，其中一个混淆可能以 50/50 的机会取两个值（例如 type1/type2，男性/女性）。它显示了各种小样本量研究的“不平衡百分比”（治疗和对照样本之间的类型 1 的差异除以样本量）的分布。红线和右侧轴显示 ecdf。
小样本随机化下不同程度平衡的概率：
从这个情节中有两件事很清楚（除非我在某个地方搞砸了）。
1) 获得完全平衡样本的概率随着样本量的增加而降低。
2) 获得非常不平衡的样本的概率随着样本量的增加而降低。
3) 在两组 n=3 的情况下，有 3% 的机会获得一组完全不平衡的组（对照组中的所有类型 1，治疗中的所有类型 2）。N=3 常见于分子生物学实验（例如用 PCR 测量 mRNA，或用蛋白质印迹测量蛋白质）
当我进一步检查 n=3 的情况时，我观察到 p 值在这些条件下的奇怪行为。左侧显示了在 type2 子组的不同均值条件下使用 t 检验计算的 p 值的总体分布。type1 的平均值为 0，两组的 sd=1。右侧面板显示了从 0.05 到 0.0001 的名义“显着性截止值”的相应误报率。
当通过 t 检验（10000 次蒙特卡罗运行）进行比较时，n=3 的 p 值分布与两个子组和第二个子组的不同均值：
以下是两组 n=4 的结果：
对于两组 n=5：
对于两组 n=10：
从上面的图表中可以看出，样本量和子组之间的差异之间似乎存在交互作用，导致在不均匀的原假设下出现各种 p 值分布。
那么我们能否得出结论，对于小样本量的适当随机和对照实验，p 值不可靠？
第一个绘图的 R 代码
require(gtools)
#pdf("sim.pdf")
par(mfrow=c(4,2))
for(n in c(3,4,5,6,7,8,9,10)){
#n<-3
p<-permutations(2, n, repeats.allowed=T)
#a<-p[-which(duplicated(rowSums(p))==T),]
#b<-p[-which(duplicated(rowSums(p))==T),]
a<-p
b<-p
cnts=matrix(nrow=nrow(a))
for(i in 1:nrow(a)){
cnts[i]<-length(which(a[i,]==1))
}
d=matrix(nrow=nrow(cnts)^2)
c<-1
for(j in 1:nrow(cnts)){
for(i in 1:nrow(cnts)){
d[c]<-cnts[j]-cnts[i]
c<-c+1
}
}
d<-100*abs(d)/n
perc<-round(100*length(which(d<=50))/length(d),2)
hist(d, freq=F, col="Grey", breaks=seq(0,100,by=1), xlab="% Unbalanced",
ylim=c(0,.4), main=c(paste("n=",n))
)
axis(side=4, at=seq(0,.4,by=.4*.25),labels=seq(0,1,,by=.25), pos=101)
segments(0,seq(0,.4,by=.1),100,seq(0,.4,by=.1))
lines(seq(1,100,by=1),.4*cumsum(hist(d, plot=F, breaks=seq(0,100,by=1))$density),
col="Red", lwd=2)
}
图 2-5 的 R 代码
for(samp.size in c(6,8,10,20)){
dev.new()
par(mfrow=c(4,2))
for(mean2 in c(2,3,10,100)){
p.out=matrix(nrow=10000)
for(i in 1:10000){
d=NULL
#samp.size<-20
for(n in 1:samp.size){
s<-rbinom(1,1,.5)
if(s==1){
d<-rbind(d,rnorm(1,0,1))
}else{
d<-rbind(d,rnorm(1,mean2,1))
}
}
p<-t.test(d[1:(samp.size/2)],d[(1+ samp.size/2):samp.size], var.equal=T)$p.value
p.out[i]<-p
}
hist(p.out, main=c(paste("Sample Size=",samp.size/2),
paste( "% <0.05 =", round(100*length(which(p.out<0.05))/length(p.out),2)),
paste("Mean2=",mean2)
), breaks=seq(0,1,by=.05), col="Grey", freq=F
)
out=NULL
alpha<-.05
while(alpha >.0001){
out<-rbind(out,cbind(alpha,length(which(p.out<alpha))/length(p.out)))
alpha<-alpha-.0001
}
par(mar=c(5.1,4.1,1.1,2.1))
plot(out, ylim=c(0,max(.05,out[,2])),
xlab="Nominal alpha", ylab="False Postive Rate"
)
par(mar=c(5.1,4.1,4.1,2.1))
}
}
#dev.off()