这一节的主要内容有：模型比较的似然比检验、拟合优度检验与偏差、未分组数据的检验方法、logit模型的残差、logistic回归的影响诊断。每一点都会用例子来辅助理解。

模型比较的似然比检验

似然比检验是我们的老朋友啦！它比较的是一个特定模型与更复杂的模型，如果更复杂的模型没有拟合得更好，那么就一定程度上保证了所选的模型是充分的。什么是更复杂的模型呢？——包含非线性效应，如二次项、交互项等等。

栗子 ：母鲎及其追随者

考虑只用宽度来预测具有追随者的概率的模型为： logit[π(x)]=α+βxlogit[\pi(x)]=\alpha+\beta x ，一种检验就是比较这个简单模型与具有二次项的模型： logit[π(x)]=α+β1x+β2x2logit[\pi(x)]=\alpha+\beta_1 x+\beta_2 x^2 。

> library(cdabookdb) > data(horseshoecrabs) > m1=glm(Satellites>0~Width,data=horseshoecrabs,family=binomial()) > m2=glm(Satellites>0~Width+I(Width^2),data=horseshoecrabs,family=binomial()) > summary(m1) Call: glm(formula = Satellites > 0 ~ Width, family = binomial(), data = horseshoecrabs) Deviance Residuals: Min 1Q Median 3Q Max -2.0281 -1.0458 0.5480 0.9066 1.6942 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -12.3508 2.6287 -4.698 2.62e-06 *** Width 0.4972 0.1017 4.887 1.02e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 225.76 on 172 degrees of freedom Residual deviance: 194.45 on 171 degrees of freedom AIC: 198.45 Number of Fisher Scoring iterations: 4 > summary(m2) Call: glm(formula = Satellites > 0 ~ Width + I(Width^2), family = binomial(), data = horseshoecrabs) Deviance Residuals: Min 1Q Median 3Q Max -2.1185 -1.0441 0.5067 0.9483 1.5406 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 14.59156 30.22371 0.483 0.629 Width -1.59572 2.35195 -0.678 0.497 I(Width^2) 0.04047 0.04566 0.886 0.375 (Dispersion parameter for binomial family taken to be 1) Null deviance: 225.76 on 172 degrees of freedom Residual deviance: 193.63 on 170 degrees of freedom AIC: 199.63 Number of Fisher Scoring iterations: 5

可以看出 β^2=0.040,SE=0.046,p−value=0.375\hat{\beta}_2=0.040,SE=0.046,p-value= 0.375 ，并没有很强的证据支持加入二次项。

似然比统计量为： LR=194.45−193.63=0.82,df=1,p−value=0.365LR=194.45-193.63=0.82,df=1,p-value=0.365 ，我们没有充分的证据去拒绝原假设 H0:β2=0H_0:\beta_2=0，因此二次项并不是必须的。

拟合优度检验与偏差

一种探测模型拟合不佳的方式是探索模型表现欠佳的所有方式。一种拟合优度检验是比较模型的拟合值与真实数据。这其实是把数据本身当做最复杂的模型，每一个观测都会有一个对应的参数。

记工作模型为MM ，为了检验 MM 的拟合，我们检验在饱和模型中但是不在 MM 中的所有参数是否都为 00 。这个检验的似然比统计量就是模型的偏差。某些情况下，似然比统计量服从大样本卡方分布。

当预测变量都是二元属性变量的时候，数据可以由列联表中的频数来概括。对于预测变量第ii 组的 nin_i 个个体，用两种属性结果的估计概括乘以 nin_i 就可以得到两个该组的两个拟合值。

偏差统计量的形式为： G2(M)=2∑observed[log⁡(observed/fitted)]G^2(M)=2\sum observed[\log (observed/fitted)] 皮尔逊统计量的形式为： X2(M)=∑(observed−fitted)2/fittedX^2(M)=\sum (observed-fitted)^2/fitted

对于固定个数的分组，当拟合频数都不小于 55的时候，上述两个统计量都在原假设成立的前提下服从卡方分布。自由度为饱和模型与工作模型参数个数的差值。当两个统计量较大的时候，说明工作模型拟合不及饱和模型。

当拟合值非常小的时候， G2G^2 与 X2X^2 在原假设情况下并不近似服从卡方分布。

栗子：AID症状，AZT使用以及种族

> library(cdabookdb) > library(tidyr) Warning message: 程辑包‘tidyr’是用R版本3.6.3 来建造的 > data(AZT) > AZT_df=spread(as.data.frame(AZT),Symptoms,Freq) > AZT_df Race AZTUse Yes No 1 White Yes 14 93 2 White No 32 81 3 Black Yes 11 52 4 Black No 12 43

数据中，如果立即使用AZT则令x=1x=1 ，否则 x=0x=0 。 z=1z=1 表示白人， z=0z=0 表示黑人。

> m=glm(cbind(Yes,No)~(Race=="White")+(AZTUse=="Yes"),data=AZT_df,family=binomial("logit")) > summary(m) Call: glm(formula = cbind(Yes, No) ~ (Race == "White") + (AZTUse == "Yes"), family = binomial("logit"), data = AZT_df) Deviance Residuals: 1 2 3 4 -0.5547 0.4253 0.7035 -0.6326 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.07357 0.26294 -4.083 4.45e-05 *** Race == "White"TRUE 0.05548 0.28861 0.192 0.84755 AZTUse == "Yes"TRUE -0.71946 0.27898 -2.579 0.00991 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 8.3499 on 3 degrees of freedom Residual deviance: 1.3835 on 1 degrees of freedom AIC: 24.86 Number of Fisher Scoring iterations: 4

利用AZT使用情况与种族对模型拟合的结果为：

logit[π^]=−0.174+0.055Race−0.719AZTuselogit[\hat{\pi}]=-0.174+0.055 Race-0.719AZTuse

由于没有交互项，模型默认了每个预测变量与响应变量之间的优势比在其它变量各个类别之间是相同的。但是其次关联的假设合理吗？

> m$fitted.values 1 2 3 4 0.1496245 0.2653998 0.1427012 0.2547241

可以看出立即接受AZT治疗的白人老兵在研究中AIDS症状发展的概率估计为0.15，那么未发展的概论估计就是0.85。类似地，我们可以计算出所有类型的老兵AIDS症状发展与未发展的概率，从而进行拟合优度检验。

> #X2 and G2s df > df=nrow(AZT_df)-length(coef(m)) > #X2 > X2=sum(resid(m,type="pearson")^2) > X2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=X2_pvalue) X2 pvalue 1.3910259 0.2382319 > #G2 > G2=sum(resid(m,type="deviance")^2) > G2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=G2_pvalue) G2 pvalue 1.3835301 0.2395008

两个检验统计量的结果都说明我们的模型拟合效果是不错的（原假设为工作模型的参数没有冗余，而我们没有充分的证据去拒绝原假设）。

分组数据、未分组数据以及连续预测变量

拟合优度检验只有对分组数据才是有意义的。当列联表绝大部分的拟合频数超过 55时，皮尔逊统计量与偏差统计量的大样本理论才是适用的。

如果数据是连续变量或者近似连续，我们可以对每个预测变量建立几个类别——例如按照分位数进行划分，然后再对这些分组数据进行拟合优度检验。但是，随着解释变量个数的增加，对每一个变量都分组会产生一个单元数极多的列联表，数据被划分的很单薄，拟合优度检验的效果也会随之降低。

针对非分组数据的另一种方法，是基于估计概率来对观测值和拟合值进行划分。比如说，我们想划分成10个组，那么被分到第一组的是估计概率最高的n/10n/10 个观测对应的观测频数与拟合频数，以此类推。分完组后，再利用Hosmer-Lemeshow检验。（我们对这个检验不做具体的介绍，知道就可以了）

栗子：鲎及其追随者

> library(ResourceSelection) > data("horseshoecrabs") > horseshoecrabs$psat=as.integer(horseshoecrabs$Satellites>0) > m=glm(psat~factor(Color)+Width,data=horseshoecrabs,family=binomial()) > hoslem.test(m$y,fitted(m)) Hosmer and Lemeshow goodness of fit (GOF) test data: m$y, fitted(m) X-squared = 4.4581, df = 8, p-value = 0.8136 > m_W=glm(psat~Width,data=horseshoecrabs,family=binomial()) > hoslem.test(m_W$y,fitted(m_W)) Hosmer and Lemeshow goodness of fit (GOF) test data: m_W$y, fitted(m_W) X-squared = 4.3855, df = 8, p-value = 0.8208

检验的结果表明，利用颜色与宽度来预测是否有追随者的模型，拟合是充分的。仅仅利用宽度来进行预测的拟合也是充分的。

logit模型的残差

一般来说，比较工作模型与更复杂模型的方法比整体的拟合优度检验更加有用。拟合优度统计量较大仅仅表明模型具有一定的拟合不佳，却无法告诉我们拟合不佳的本质。而比较工作模型与更复杂模型却可以让我们知道，如果拟合不佳的话，是什么地方出了问题。

当预测变量为属性变量的时候，我们可以利用残差来比较观测频数与拟合频数。这个想法我们在之前来也提到过——当总体拟合不佳的时候，通过残差来看看哪里出了问题。我们令yiy_i 表示在解释变量的第 ii 组的 nin_i 次试验当中成功的次数， π^i\hat{\pi}_i 则表示由模型拟合得到的成功概率的估计。二项分布的均值的估计 niπ^in_i\hat{\pi}_i 就是成功次数的拟合值。

残差的使用最好是对分组数据，不然残差的意义有限。

若GLM的随机部分服从二项分布，比较yiy_i 与 niπ^in_i\hat{\pi}_i 的：

皮尔逊残差为 ei=yi−niπ^i[niπ^i(1−π^i)]e_i=\frac{y_i-n_i\hat{\pi}_i}{\sqrt{[n_i\hat{\pi}_i(1-\hat{\pi}_i)]}} ，当 nin_i 较大的时候， eie_i 近似服从正态分布。标准化残差为 yi−niπ^iSE=yi−niπ^i[niπ^i(1−π^i)(1−hi)]\frac{y_i-n_i\hat{\pi}_i}{SE}=\frac{y_i-n_i\hat{\pi}_i}{\sqrt{[n_i\hat{\pi}_i(1-\hat{\pi}_i)(1-h_i)]}} , hih_i 来自帽子矩阵对角元。

我们更愿意使用标准化残差，当模型成立的时候，标准化残差也近似服从正态分布。绝对值大于 22 或 33 我们就可以认为拟合不佳了。

栗子：Florida 大学研究生入学情况

> data("UFAdmissions") > UFAdmissions_df <- spread(as.data.frame(UFAdmissions), Decision, Freq) > head(UFAdmissions_df) Dept Gender Admitted Rejected 1 anth Female 32 81 2 anth Male 21 41 3 astr Female 6 0 4 astr Male 3 8 5 chem Female 12 43 6 chem Male 34 110

我们的数据大致长成这样：Florida 大学97-98年文科与理科学院的23个系的研究生院申请情况。一共有4个变量：Dept(系别)，Gender (性别)，Admitted (接收人数)，Rejected ( 拒绝人数)。对于系kk 中性别为 ii 的 nikn_{ik} 个申请者，令 yiky_{ik} 表示接受人数， πik\pi_{ik} 表示接受概率。

在其它条件相同的前提下，我们希望接受的决定是独立于性别的。给定系别下，没有性别效应的模型为：logit(πik)=α+βkDlogit(\pi_{ik})=\alpha+\beta_k^D。但由于某些系列里或许存在性别效应，这个模型也许并不充分。于是我们可以看一下这个超级简单模型的拟合效果：

> m=glm(cbind(Admitted,Rejected)~Dept,data=UFAdmissions_df,family=binomial()) > df=nrow(UFAdmissions_df)-length(coef(m)) > #X2 > X2=sum(resid(m,type="pearson")^2) > X2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=X2_pvalue) X2 pvalue 40.85235935 0.01231099 > #G2 > G2=sum(resid(m,type="deviance")^2) > G2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=G2_pvalue) G2 pvalue 44.735164687 0.004282328

显然模型是糟糕的。这其实侧面说明了性别效应存在的可能。于是我们计算这个模型中女性被接收者的标准化残差。奇数项是女性的残差，男女标准化残差仅仅相差了一个负号。

> #standardized residuals > rstandard(m,type="pearson") 1 2 3 4 5 6 -0.76456728 0.76456728 2.87096225 -2.87096225 -0.26830432 0.26830432 7 8 9 10 11 12 -1.06904497 1.06904497 -0.63260081 0.63260081 1.15751874 -1.15751874 13 14 15 16 17 18 0.94208925 -0.94208925 2.16641024 -2.16641024 -0.26082027 0.26082027 19 20 21 22 23 24 1.88730092 -1.88730092 -0.17627090 0.17627090 1.64564062 -1.64564062 25 26 27 28 29 30 1.37297695 -1.37297695 1.28843796 -1.28843796 1.34164079 -1.34164079 31 32 33 34 35 36 1.32457581 -1.32457581 -0.23317590 0.23317590 -2.27221543 2.27221543 37 38 39 40 41 42 1.26491106 -1.26491106 0.13969719 -0.13969719 0.30122617 -0.30122617 43 44 45 46 -0.01229099 0.01229099 -1.75872555 1.75872555

有大的标准化残差的系可能是造成拟合不佳的罪魁祸首。对女性来说，残差绝对值超过 22 的有天文学系（ 2.872.87 ）、地理学系（ 2.172.17 ）以及心理学系 (−2.27)(-2.27) 。倘若我们把这三个系删了：

> UFAdmission_df_e=UFAdmissions_df[-c(3,4,15,16,35,36),] > m_e=glm(cbind(Admitted,Rejected)~Dept,data=UFAdmission_df_e,family=binomial()) > df=nrow(UFAdmission_df_e)-length(coef(m_e)) > #X2 > X2=sum(resid(m_e,type="pearson")^2) > X2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=X2_pvalue) X2 pvalue 22.7536388 0.3010527 > #G2 > G2=sum(resid(m_e,type="deviance")^2) > G2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=G2_pvalue) G2 pvalue 24.3687511 0.2266554

模型的拟合效果好了很多。

接下来我们可以考虑把性别效应加进去，可是模型与我们最初的相比似乎并没有得到提升。

> m_g=glm(cbind(Admitted,Rejected)~Dept+Gender,data=UFAdmissions_df,family=binomial()) > df=nrow(UFAdmissions_df)-length(coef(m_g)) > #X2 > X2=sum(resid(m_g,type="pearson")^2) > X2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=X2_pvalue) X2 pvalue 38.99080225 0.01414625 > #G2 > G2=sum(resid(m_g,type="deviance")^2) > G2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=G2_pvalue) G2 pvalue 42.360051041 0.005652222

没有提升的原因我们可以来看看模型性别与是否被接受的条件优势比与边缘优势比：

> m_g$coefficients[24] GenderMale -0.1729698 > exp(-m_g$coefficients[24]) 1.18883

也就是说，给定系别，女性接收的优势比男性接收的优势高 18%18\%。可是，如果我们计算边缘优势比：

> library(cdabookfunc) > oddsratio(margin.table(UFAdmissions,margin=c(2,3))) oddsratio 1 0.935911

可以发现女性接受的优势比男性接受的优势低 6%6\% ，辛普森悖论再现了！条件关联与边缘关联方向相反使得性别效应得不到体现。

logistic 回归的影响诊断

与普通的回归一样，某些观测可能会对参数估计产生过大的影响。如果删除这些观测，拟合结果将会非常不同。当残差表明模型对于某些观测拟合得不好时，我们应当将这些观测删去，再对剩下的数据进行拟合。

每个观测的影响诊断有以下几个方面：

当删除观测时分析模型中每个参数估计的变化。这个变化量除以它的标准误得到的统计量为Dfbeta。删除观测后得到的参数联合置信区间的变化的度量。这个置信区间偏移诊断记为 cc。删除观测后的拟合优度统计量X2X^2 或者 G2G^2 的变化。

每一种度量值越大观测的影响就越大。

栗子：心脏病与血压的关系

> data("blood_pressure") > blood_pressure BloodPressure SampleSize ObservedDisease 1 111.5 156 3 2 121.5 252 17 3 131.5 284 12 4 141.5 271 16 5 151.5 139 12 6 161.5 85 8 7 176.5 99 16 8 191.5 43 8

这个数据将 40−5940-59岁的男性按照x=BloodPressurex=BloodPressure (血压水平的得分)，以及 y=HeartDiseasey=HeartDisease 进行了交叉分类。令 πi\pi_i 表示血压类别 ii患心脏病的概率。我们拟合的模型是：logit[πi]=α+βxilogit[\pi_i]=\alpha+\beta x_i 。

> result=influence_logit_sas(m,"data.frame") > result hat resid.pearson resid.deviance dfbetas..Intercept. dfbetas.BloodPressure cvalues cbar difchisq difdev 1 0.2154758 -0.9794311 -1.0616803 -0.5328501641 0.49239757 0.335840815 0.263475258 1.22276052 1.39064034 2 0.2865667 2.0057103 1.8501114 1.2777007175 -1.14316418 2.264933187 1.615878738 5.63875247 5.03879076 3 0.2596674 -0.8133348 -0.8419625 -0.3930528813 0.32844242 0.313402501 0.232022083 0.89353560 0.94092295 4 0.2172430 -0.5067270 -0.5162271 -0.1245126431 0.08164186 0.091041596 0.071263451 0.32803572 0.33775388 5 0.1303602 0.1175833 0.1170033 -0.0004670062 0.00782318 0.002383184 0.002072512 0.01589834 0.01576228 6 0.1293979 -0.3042459 -0.3087740 0.0481729344 -0.06513218 0.015802915 0.013758051 0.10632361 0.10909944 7 0.3797042 0.5134721 0.5049655 -0.3487145849 0.40086540 0.260184479 0.161391345 0.42504495 0.41638154 8 0.3815848 -0.1394648 -0.1402441 0.1127163676 -0.12377298 0.019407074 0.012001629 0.03145206 0.03167005

我们计算出了每一个观测值对应的影响诊断。hat表示杠杆值，resid.pearson是皮尔逊残差，resid.deviance是LR残差，dfbeta.Intercept是截距项的Dfbeta， dfbetas.BloodPressure则是变量的Dfbeta，cvalue是置信区间的偏移，cbar是预测值的偏移，difchisq是除去这个观测之后皮尔逊统计量的差，difdev则是除去这个统计量之后LR统计量的差。

可以看出，删除第二个观测的影响是最大的。但是这并不会让人惊奇——当有很多观测时，一小部分诊断结果较大可能仅仅是因为偶然。

模型整体的拟合优度检验还是不错的。

> #X2 and G2s df > df=nrow(blood_pressure)-length(coef(m)) > #X2 > X2=sum(resid(m,type="pearson")^2) > X2_pvalue=1-pchisq(X2,df) > c(X2=X2,pvalue=X2_pvalue) X2 pvalue 6.2899402 0.3915068 > #G2 > G2=sum(resid(m,type="deviance")^2) > G2_pvalue=1-pchisq(G2,df) > c(G2=G2,pvalue=G2_pvalue) G2 pvalue 5.9091582 0.4334429

我们在诊断性分析中把出现的模式归因于来自模型的偶然变化是需要谨慎的。我们如果要删除，通常都是删除某一类别的整个二项样本，仅仅删去一个单独的个体作用微乎其微。

另外我们还可以画出观测比例与拟合比例的散点图，由图形来看拟合的效果。