在To P or not to P: on the evidence nature of P-values and their place in scientic inference中，Michael Lew 表明，至少对于 t 检验，可以解释单边 p 值和样本量作为给定似然函数的“地址”（我的术语）。我在下面重复了他的一些数字，稍作修改。左列显示了由于不同效应大小（均值/合并标准差之间的差异）和样本大小的理论而预期的 p 值分布。水平线标记了“切片”，我们从中得到右图所示的 p=0.50 和 p=0.025 的似然函数。
这些结果与蒙特卡罗模拟一致。对于这个图，我通过 t 检验比较了 n=10 的两组，在许多不同的效应大小下，将 10,000 个 p 值合并到每个效应大小的 0.01 区间。具体来说，有一组平均值 = 0，sd = 1，第二组平均值从 -4 变化到 4，sd = 1。
（上面的数字可以直接与上面链接的论文中的图 7/8 进行比较，并且非常相似，我发现热图比那篇论文中使用的“云”更能提供信息，并且还希望独立地复制他的结果。）
如果我们检查由 p 值“索引”的似然函数，拒绝/接受假设或忽略基于截止值（在任何地方使用的任意 0.05 或由成本确定的 p 值）大于 0.05 的结果的行为-benefit) 似乎是荒谬的。为什么我不应该从 n=100, p=0.5 的情况下得出“当前证据表明任何影响（如果存在）都很小”的结论？当前的做法是要么“接受”没有效果（假设检验），要么说“需要更多数据”（显着性检验）。我不明白为什么我应该做这两件事。
也许当一个理论预测了一个精确的点值时，拒绝一个假设可能是有道理的。但是，当假设的形式为“mean1=mean2 或 mean1!=mean2”时，我认为没有任何价值。在经常使用这些测试的条件下，随机化并不能保证所有混杂因素在各组之间都是平衡的，并且应该始终担心潜伏变量，因此据我所知，拒绝 mean1 完全等于 mean2 的假设没有科学价值。
是否有超出 t 检验的案例不适用此论点？我是否错过了拒绝先验概率低的假设为研究人员提供的有价值的东西？忽略任意截止值以上的结果似乎导致了广泛的发表偏倚。忽略结果对科学家有什么有用的作用？
Michael Lew 的 R 代码来计算 p 值分布
LikeFromStudentsTP<-function(n,x,Pobs,test.type, alt='one.sided'){
# test.type can be 'one.sample', 'two.sample' or 'paired'
# n is the sample size (per group for test.type = 'two.sample')
# Pobs is the observed P-value
# h is a small number used in the trivial differentiation
h<-10^-7
PowerDn<-power.t.test('n'=n, 'delta'=x, 'sd'=1,
'sig.level' = Pobs-h, 'type'= test.type, 'alternative'=alt)
PowerUp<-power.t.test('n'=n, 'delta'=deltaOnSigma, 'sd'=1,
'sig.level' = Pobs+h, 'type'= test.type, 'alternative'=alt)
PowerSlope<-(PowerUp$power-PowerDn$power)/(h*2)
L<-PowerSlope
}
图 1 的 R 代码
deltaOnSigma <- 0.01*c(-400:400)
type<-'two.sample'
alt='one.sided'
p.vals<-seq(0.001,.999,by=.001)
#dev.new()
par(mfrow=c(4,2))
for(n in c(3,5,10,100)){
m<-matrix(nrow=length(deltaOnSigma), ncol=length(p.vals))
cnt<-1
for(P in p.vals){
m[,cnt]<-LikeFromStudentsTP(n,deltaOnSigma,P,type, alt)
cnt<-cnt+1
}
#remove very small values
m[which(m/max(m,na.rm=T)<10^-5)]<-NA
m2<-log(m)
par(mar=c(4.1,5.1,2.1,2.1))
image.plot(m2, axes=F,
breaks=seq(min(m2, na.rm=T),max(m2, na.rm=T),length=1000), col=rainbow(999),
xlab="Effect Size", ylab="P-value"
)
title(main=paste("n=",n))
axis(side=1, at=seq(0,1,by=.25), labels=seq(-4,4,by=2))
axis(side=2, at=seq(0,1,by=.05), labels=seq(0,1,by=.05))
axis(side=4, at =.5, labels="Log-Likelihood", pos=.95, tick=F)
abline(v=0.5, lwd=1)
abline(h=.5, lwd=3, lty=1)
abline(h=.025, lwd=3, lty=2)
par(mar=c(5.1,4.1,4.1,2.1))
plot(deltaOnSigma,m[,which(p.vals==.025)], type="l", lwd=3, lty=2,
xlab="Effect Size", ylab="Likelihood", xlim=c(-4,4),
main=paste("Likelihood functions for","n=",n)
)
lines(deltaOnSigma,m[,which(p.vals==.5)], lwd=3, lty=1)
legend("topleft", legend=c("p=.5","p=.025"), lty=c(1,2),lwd=1, bty="n")
}
图 2 的 R 代码
p.vals<-seq(0,1,by=.01)
deltaOnSigma <- 0.01*c(-400:400)
n<-10
n2<-10
sd2<-1
num.sims<-10000
sp<-sqrt((9*1^2 +(n2-1)*sd2^2)/(n+n2-2))
p.out=matrix(nrow=num.sims*length(deltaOnSigma) ,ncol=2)
m<-matrix(0,nrow=length(deltaOnSigma),ncol=length(p.vals))
pb<-txtProgressBar(min = 0, max = length(deltaOnSigma) ,style = 3)
cnt<-1
cnt2<-1
for(i in deltaOnSigma ){
for(j in 1:num.sims){
m2<-i
a<-rnorm(n,0,1)
b<-rnorm(n,m2,sd2)
p<-t.test(a,b, alternative="less")$p.value
r<-end(which(deltaOnSigma<=m2/sp))[1]
m[r,end(which(p.vals<p))[1]]<-m[r,end(which(p.vals<p))[1]]+1
p.out[cnt,]<-cbind(m2/sp,p)
cnt<-cnt+1
}
cnt2<-cnt2+1
setTxtProgressBar(pb, cnt2)
}
close(pb)
m[which(m==0)]<-NA
m2<-log(m)
dev.new()
par(mfrow=c(2,1))
par(mar=c(4.1,5.1,2.1,2.1))
image.plot(m2, axes=F,
breaks=seq(min(m2, na.rm=T),max(m2, na.rm=T),length=1000), col=rainbow(999),
xlab="Effect Size", ylab="P-value"
)
title(main=paste("n=",n))
axis(side=1, at=seq(0,1,by=.25), labels=seq(-4,4,by=2))
axis(side=2, at=seq(0,1,by=.05), labels=seq(0,1,by=.05))
axis(side=4, at =.5, labels="Log-Count", pos=.95, tick=F)
abline(h=.5, lwd=3, lty=1)
abline(h=.025, lwd=3, lty=2)
abline(v=.5, lwd=2, lty=1)
par(mar=c(5.1,4.1,4.1,2.1))
hist(p.out[which(p.out[,2]>.024 & p.out[,2]<.026),1],
xlim=c(-4,4), xlab="Effect Size", col=rgb(1,0,0,.5),
main=paste("Effect Sizes for","n=",n)
)
hist(p.out[which(p.out[,2]>(.499) & p.out[,2]<.501),1], add=T,
xlim=c(-4,4),col=rgb(0,0,1,.5)
)
legend("topleft", legend=c("0.499<p<0.501","0.024<p<0.026"),
col=c("Blue","Red"), lwd=3, bty="n")