摘要:
拟合度检验富于教学内涵,没有应用意义。文中图解
拟合度检验中的正态随机向量,示以数值实例。文末给出(与全文其实无关的)应用建议和代码示例。数据:随机抛32枚硬币,22枚正面朝上[1]。双尾虚无假设:正面朝上总体期望0.5。
Excel中的解:随机抛32枚均匀硬币,正面朝上服从二项分布。百分位数21的百分等级97.49%,本题的右尾p值为100%-97.49%。分布对称,双尾p 值为右尾p 值的两倍。
二项分布检验的Excel解答binom.test(x = 22,n = 32)
##
## Exact binomial test
##
## data: 22 and 32
## number of successes = 22, number of trials = 32, p-value = 0.0501
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.4999224 0.8388153
## sample estimates:
## probability of success
## 0.6875
在没有个人电脑、只能纸笔查表的年代,这类问题通过正态近似解决。题中的二项分布随机数X 近似服从的正态分布总体均值
![]()
=32×0.5=16、标准误总体参数
![]()
=
![]()
。下图解释了为什么要做±0.5的校正——Pr(离散的二项随机数X≥22)≈Pr(对应连续正态随机数Z>21.5),约等号左边的22恰好取到的概率近似于Pr(22.5≥Z>21.5)。
Pr(离散的二项随机数X≥22)≈Pr(对应连续正态随机数Z>21.5)本例21.5的标准化z 值1.9445,对应的双尾p 值为0.05183。
prop.test(x = 22,n = 32)
##
## 1-sample proportions test with continuity correction
##
## data: 22 out of 32, null probability 0.5
## X-squared = 3.7812, df = 1, p-value = 0.05183
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.4986377 0.8325051
## sample estimates:
## p
## 0.6875
但如果没有作连续性校正,本例22的标准化z 值就是2.1213,对应的双尾p 值为0.03389,
![]()
=4.5000。
![]()
的右尾
p 值即对应的
z 双尾
p 值
。prop.test(x = 22,n = 32,correct = F)
##
## 1-sample proportions test without continuity correction
##
## data: 22 out of 32, null probability 0.5
## X-squared = 4.5, df = 1, p-value = 0.03389
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.5143332 0.8204746
## sample estimates:
## p
## 0.6875
- 二水平拟合度之
![]()
检验
R (ver. 4.0.0) 与SPSS (ver. 23) 的
![]()
拟合度检验都不作连续性校正。
chisq.test(c(22,10))
##
## Chi-squared test for given probabilities
##
## data: c(22, 10)
## X-squared = 4.5, df = 1, p-value = 0.03389
SPSS的菜单是⌥A ▶ N (▶ L)▶ C ▶ ... ▶ ⌥X ▶ ⌥E ▶...。SPSS的结果可以同时给出二项分布检验的Exact Sig.双尾p 值0.0501。
二水平(Head变量)、三水平(Tri变量) χ2拟合度SPSS结果- 拟合度之
![]()
检验:从二水平到三水平
二水平的
![]()
拟合度检验
df =1,等价于不作±0.5连续性校正的
z 检验。这个
z 检验也就是百分比的点估计
prop=22/32相对虚无假设
![]()
=0.5的
z 检验,
![]()
。
从二水平向三水平的情形推广,将正面朝上的硬币分为两种情形:5枚硬币没有弹起、直接正面朝上(Tri=1),17枚硬币落地后弹起,最后正面朝上(Tri=-1)。此时虚无假设参数是一个三维空间中的二维度向量,Tri三种取值(-, 0, +)的百分比总体期望满足
![]()
=1。缺省的虚无假设是三个百分比总体期望都是
![]()
。SPSS报告的
![]()
值与R的结果一致,都等于
![]()
。
chisq.test(c(17,10,5))
##
## Chi-squared test for given probabilities
##
## data: c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317
受二水平的情形启发,这个三水平情形的
![]()
值应当与某一个二维正态向量的长度平方对应。为了使二水平情形与三水平情形更好类比,可以将二水平的问题表述为:Head的两种取值(0, 1)的百分比总体期望满足
![]()
=1。这正是下图斜率为-1的对角线上的点集。 一维随机向量以概率
![]()
出现蓝色箭头,以概率
![]()
出现粉色箭头。32枚个案有22次粉色箭头,10次蓝色箭头,平均值以黑点描出。黑点距离出发点
![]()
的标准误总体参数=无限枚箭头个案的总体标准差
![]()
。根据中心极限定理,计算出的从出发点到黑点向量以标准误作标准化,得到的向量长度平方值(在样本量不太小的时候)近似服从
![]()
分布。
μ_{Head1}=μ_{Head0}=1/2情形可以验证,即使
![]()
取一般的数值比如0.4,以上的几何解读同样成立。
对一般的
![]()
作图、验证的R代码如下。
mu1 <- 0.4; mu0 <- 1 - mu1;
plot(c(0,1),c(1,0),col='grey',type='l',
asp=1,xlab = "Head=1 Proportion",ylab="Head=0 Proportion")
abline(h=0,v=0,lty=3)
points(22/32,10/32,pch=19)
points(mu1,mu0)
arrows(mu1,mu0,0,1,col='blue')
arrows(mu1,mu0,1,0,col='pink')
chisq.test(xc <- c(22,10),p = mu<- c(mu1,mu0))
##
## Chi-squared test for given probabilities
##
## data: xc <- c(22, 10)
## X-squared = 11.021, df = 1, p-value = 0.0009009
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1)
## [1] 11.02083
(sigma2 <- mu1*sum((c(1,0)-mu)^2) + mu0*sum((c(0,1)-mu)^2))
## [1] 0.48
sum((xc/sum(xc) - mu)^2)/(sigma2/sum(xc))
## [1] 11.02083
μ_{Head1}=0.4, μ_{Head0}=0.6 情形这样就不难从二维推广到三维。下图三个坐标和=1的平面构成(
![]()
,
![]()
,
![]()
)的取值空间。蓝、红、绿三种向量各自出现的概率为
![]()
。如果这三个总体参数都是1/3,三种向量的平均数根据中心极限定理近似服从(三维空间中)二维的正态分布。根据对称性,这个二维正态分布各向方差一致。可以计算:彩色个案的
![]()
=各彩色向量长度平方用概率加权求和,再除以自由度2。黑向量的标准误总体参数平方
![]()
=
![]()
。在样本量不太小的时候,黑向量长度平方/
![]()
得到的数值近似服从
![]()
分布。
总体参数各向对称的情形作图和数值验证的R代码如下——
## 数据与对称的参数
mu_neg <- 1/3; mu_neu <- 1/3;
mu_pos <- 1- mu_neg - mu_neu;
mu <- c(mu_neg,mu_neu,mu_pos)
chisq.test(x = xc<-c(17,10,5),p = mu)
##
## Chi-squared test for given probabilities
##
## data: xc <- c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1) # 6.8125
## [1] 6.8125
## 三维作图
if(!require(plot3Drgl)){install.packages('plot3Drgl');require(plot3Drgl);}
## Loading required package: plot3Drgl
## Loading required package: rgl
## Loading required package: plot3D
open3d()
## wgl
## 1
plot3d(c(0,1,0,0),c(0,0,1,0),c(1,0,0,1),type='l',col='grey',
asp=1,xlab='Tri<0 Prop',ylab = 'Tri=0 Prop',zlab = 'Tri>0 Prop')
# points3d(mu_neg,mu_neu,mu_pos,pch=1)
arrow3d(mu,c(1,0,0),col='blue')
arrow3d(mu,c(0,1,0),col='red')
arrow3d(mu,c(0,0,1),col='green')
arrow3d(mu,xc/sum(xc),col='black',type='r')
## 对称情形的数值复核
(sigma2 <- sum(mu*c(
sum((c(1,0,0)-mu)^2),
sum((c(0,1,0)-mu)^2),
sum((c(0,0,1)-mu)^2)))/2) # df=2
## [1] 0.3333333
sum((xc/sum(xc) - mu)^2)/(sigma2/sum(xc))
## [1] 6.8125
- 拟合度之
![]()
检验:三水平不对称参数的一般情形
如果蓝、红、绿三种向量各自出现的概率
![]()
不对称,比如 (0.3, 0.5, 0.2),三种向量的平均数根据中心极限定理仍然近似服从(三维空间中)二维的正态分布。但这个二维正态分布各向方差不再一致,需要在旋转后的新坐标架下估算图示中三角形截面上二维正态分布的方差-协方差矩阵总体参数。参考本专栏同主题文章
[2]给出的旋转变换矩阵,可将三角截面上所有向量变换为第一坐标为0的二维向量——向量起点、终点变换后的第一坐标都是同一常数0.5773503,相减为0)。
总体参数各向不对称的情形## 坐标旋转
(P_mt <- cbind(c(1,1,1)/ 3^.5,
c(1,-1/2,-1/2)/ 1.5^.5,
c(0,1,-1)/ 2^.5))
## [,1] [,2] [,3]
## [1,] 0.5773503 0.8164966 0.0000000
## [2,] 0.5773503 -0.4082483 0.7071068
## [3,] 0.5773503 -0.4082483 -0.7071068
(prop_t <- ((xc <- c(17,10,5))/sum(xc)) %*% P_mt)
## [,1] [,2] [,3]
## [1,] 0.5773503 0.2423974 0.1104854
## 不对称的参数及其数值复核
mu_neg <- 0.3; mu_neu <- 0.5;
mu_pos <- 1- mu_neg - mu_neu;
mu <- c(mu_neg,mu_neu,mu_pos)
chisq.test(x = xc,p = mu)
##
## Chi-squared test for given probabilities
##
## data: xc
## X-squared = 8.2604, df = 2, p-value = 0.01608
sum((xc- sum(xc)*mu)^2 * (sum(xc)*mu)^-1)
## [1] 8.260417
(mu_t <- mu %*% P_mt)
## [,1] [,2] [,3]
## [1,] 0.5773503 -0.04082483 0.212132
(brg <- P_mt[,2:3]-cbind(rep(1,3))%*%mu_t[2:3])
## [,1] [,2]
## [1,] 0.8573214 -0.2121320
## [2,] -0.3674235 0.4949747
## [3,] -0.3674235 -0.9192388
(sigma2 <- t(brg)%*%diag(mu)%*%brg)
## [,1] [,2]
## [1,] 0.31500000 -0.07794229
## [2,] -0.07794229 0.30500000
rbind(prop_t[2:3] - mu_t[2:3])%*%solve(sigma2/sum(xc))%*%(prop_t[2:3] - mu_t[2:3])
## [,1]
## [1,] 8.260417
即使不是保持尺度不变的正交旋转,一般的线性变化也不改变样本均值(黑箭头)与个案总体分布(彩色箭头)的相对尺度。可以把三角截面简单投影的第1、2轴平面,或者第1、3轴平面,或者第2、3轴平面,都有同样的结果——
t(sapply(list(c(1,2),c(1,3),c(2,3)),
function(jy){
Sigma2 <-
t(brg2 <- (diag(3)-cbind(rep(1,3))%*%mu)[,jy]) %*%
diag(mu) %*% brg2;
z2 <- rbind(D2 <- (xc/sum(xc)-mu)[jy]) %*%
solve(Sigma2/sum(xc)) %*% D2;
c(jy, z2)
}
)
)
## [,1] [,2] [,3]
## [1,] 1 2 8.260417
## [2,] 1 3 8.260417
## [3,] 2 3 8.260417
- SPSS (ver. 23)的多水平拟合度Exact Sig.与多项分布检验结果不一致
在上图对Tri变量作的拟合度
![]()
检验结果中,倒数第二行报告Exact Sig.为0.034,与多项分布检验缺省的报告结果(0.04119207≈0.0412)不一致。R中的EMT包提供了多项分布检验功能,与R原生的多项分布密度函数编程得到的结果一致。
chisq.test(xc<-c(17,10,5))
##
## Chi-squared test for given probabilities
##
## data: xc <- c(17, 10, 5)
## X-squared = 6.8125, df = 2, p-value = 0.03317
if(!require(EMT)){ install.packages("EMT");require(EMT);}
## Loading required package: EMT
multinomial.test(xc ,prob =mu <- rep(1/3,3))
##
## Exact Multinomial Test, distance measure: p
##
## Events pObs p.value
## 561 9e-04 0.0412
## Verified (Exact Multinomial Test)
(d_cri <- dmultinom(xc,prob = mu))
## [1] 0.0009168089
dar <- array(dim=c(33,33))
for (i in 0:32) for (j in 0:32){k <- 32- i - j; dar[i+1,j+1]=ifelse(k<0,0,dmultinom(c(i,j,k),prob = mu))}
sum(dar) # all probabilities
## [1] 1
sum(dar[dar<=d_cri]) # p_value
## [1] 0.04119207
SPSS报告的Exact Sig.可以通过复制表格到Excel得到更精确的数值0.034497。注意到EMT包的多项分布检验提供了useChisq=T选项,在此设置下,虚无假设下各事件的极端程度等级由
![]()
值定义,与多项分布密度定义的极端程度等级略有出入,正如《F 检验与Tukey HSD检验的备择假设方向图解》的图解,超出单维度之后,两种备择假设梯度不同。考虑到
![]()
值在学术史上被引入这个问题只是因为特定时代计算不便被迫查表,在多项分布检验各事件极端等级的比较中可能没必要越俎代庖替代多项分布密度本身。
验算确证,SPSS报告的Exact Sig.是虚无假设下看到更大或相等
![]()
值的概率,但只是多项分布的概率而不是
![]()
分布的概率。
## Exact Sig. of Chisq
multinomial.test(xc ,prob =mu,useChisq = T)
##
## Exact Multinomial Test, distance measure: chisquare
##
## Events chi2Obs p.value
## 561 6.8125 0.0345
(c_cri <- sum((xc-(ne <- mu*sum(xc)))^2 / ne))
## [1] 6.8125
dar.c <- dar
for (i in 0:32) for (j in 0:32){
k <- 32- i - j;
dar.c[i+1,j+1]=ifelse(k<0,0,sum((c(i,j,k)-ne)^2 / ne));
dar[i+1,j+1]=ifelse(k<0,0,dmultinom(c(i,j,k),prob = mu));
}
sum(dar[dar.c>=c_cri]) # p_value
## [1] 0.03449689
- 给不教统计课同行的应用建议(以及给统计课同行的吐槽)
- 二水平拟合度检验问题,推荐在 R 中作二项分布检验,报告百分比总体参数的置信区间。SPSS没有给出置信区间,但仍在
![]()
拟合度检验⌥A ▶ N(▶ L)▶ C ... 结果表格报告了Exact Sig.,即二项分布检验的
p 值,应用中以此为准,不必参考![]()
的
p 值。可以这么说:对于二水平拟合度检验问题,![]()
拟合度检验只有教学价值没有应用价值,可以、也应当被二项分布检验全面替代。
- 三水平(以及较少水平数的)拟合度检验问题,应用中有两类情形应当被区别对待。其中一种情形虽然非常流行但几乎从来都不是研究者的真实意图:只关注虚无假设是否被拒绝,不报告各个水平具体的百分比总体参数置信区间。这种情形建议在 R 中通过EMT包作多项分布检验。同样可以这么说:对于多水平拟合度检验问题,
![]()
拟合度检验同样只有教学价值,如果与多项分布检验结果有出入,以多项分布检验
p 值为准。那么,实验学科的统计课为什么要教![]()
拟合度检验,只是因为它的公式比较简洁?可能学术史的逻辑真就是这么回事。在没有计算机的时代,
![]()
统计量方便手算然后查表对临界值。到了个人计算机普及的时代,统计课程如果仍然坚持教
![]()
拟合度检验,或许需要这样一个理由:
![]()
拟合度检验可以进一步讲解中心极限定理,这也正是本文的意图。
- 三水平(以及较少水平数的)拟合度检验问题,应用中常见情形应当报告各个百分比参数的置信区间,可以在R中作多个二项分布检验。在事后检验的多重比较场合,必须确保一系列(family-wise)置信区间整体置信水平,可通过Bonferroni方法校正单个区间的一类错误率。理论上也不难推算百分比总体参数向量的置信区域,但不便于文本方式传播,在实际研究中极少有同行报告置信区域。三水平拟合度Bonferroni方法事后检验示例——
alpha_familywise <- 0.05
xc<-c(17,10,5)
K <- length(xc)
alpha_corrected <- alpha_familywise/K
cl<-sapply(xc,FUN = function(x) binom.test(x,sum(xc),conf.level = 1-alpha_corrected)$conf)
rbind(Lower=cl[1,],Est=xc/sum(xc),Upper=cl[2,])
## [,1] [,2] [,3]
## Lower 0.3128146 0.1367970 0.03996758
## Est 0.5312500 0.3125000 0.15625000
## Upper 0.7413646 0.5383827 0.36617910
参考
- ^数据背景: Kahneman, D., Frederickson, B. L., Schreiber, C. A., & Redelmeier, D. A. (1993). When More Pain Is Preferred to Less: Adding a Better End, Psychological Science, 4, 401-405.
- ^F检验与Tukey HSD检验的备择假设方向图解 https://zhuanlan.zhihu.com/p/50996318