《属性数据分析引论》部分课后习题R语言实践（第三章、第四章）

最新推荐文章于 2022-11-02 11:32:52 发布

摸鱼的Kaz哥

最新推荐文章于 2022-11-02 11:32:52 发布

阅读量4.6k

点赞数 9

文章标签： r语言数据分析概率论

本文链接：https://blog.csdn.net/LuKenkazia/article/details/121886067

版权

前言

习题选自高等教育出版社译制，Alan Agresti著的《属性数据分析引论(第二版)》中，第三章广义线性模型、第四章Logistic回归中的课后习题。具体题目在文中给出。

本人目前是一位在读的应用统计学专业本科生，这些题目是在课前进行的练习，所给出的思路和答案可能有所错误，欢迎大家批评指正。

第三章广义线性模型

习题3.18

表3.8列出了英乙联赛一个赛季每支球队观赛总人数（千人）和被捕总人数.

表3.8 习题3.18关于足球队逮捕人数数据
球队	观众数	逮捕数	球队	观众数	逮捕数
阿斯顿维拉	404	308	什鲁斯伯里	108	68
布拉德福城	286	197	史云顿	210	67
利兹联	443	184	谢菲尔德联	224	60
伯恩茅斯	169	149	斯托克城	211	57
西布朗维奇	222	132	巴恩斯利	168	55
汉德斯菲德	150	126	米尔沃尔	185	44
米德尔斯堡	321	110	侯城	158	38
伯明翰	189	101	曼彻斯特城	429	35
伊普斯维奇	258	99	普利茅斯	226	29
莱切斯特城	223	81	雷丁	150	20
布莱克本	211	79	奥威	148	19
水晶宫	215	78

a小题

令Y表示观赛总人数为t的球队被捕球迷人数。说明为什么模型E(Y)=μt是可行的。它有等价形式log[E(Y)/t]=α，其中α=log(μ)，给出带位移项的模型表达式。

题解：（本题的解释我并不确定正确）

题目指出了计数响应Y(被捕球迷数)有指标t(总观众数)，那么我们关心的是样本的比率Y/t

若设样本比率的期望值为μ，即，那么两边同乘t便有模型

样本比率的对数模型应为，x是效应因子，在本题中将每个球队的数据视作一次观测，并无效应因子，于是模型表示为，与比较就可得出样本比率的期望值为常数的结论

可以给出带位移的模型表达式为：

b小题

假设样本为泊松样本，拟合模型。给出并解释。

题解：

本题用R进行模型拟合，先将表3.8的数据输入进Excel保存为csv文件，以下是实现的代码

> data3.8=read.csv('table3.8.csv')    #读取数据，并对数据框进行一些处理
> rownames(data3.8)=data3.8[,1]
> data3.8=data3.8[,-1]
> colnames(data3.8)=c('t','Y')
> head(data3.8)                        #展示数据框前6行
###
             t   Y
阿斯顿维拉 404 308
布拉德福城 286 197
利兹联     443 184
伯恩茅斯   169 149
西布朗维奇 222 132
汉德斯菲德 150 126
###
#接着用glm()函数进行拟合，offset表示位移项
> model3.8=glm(Y~NULL,data=data3.8,family=poisson(link='log'),offset=log(t))
> summary(model3.8)
###
Call:
glm(formula = Y ~ NULL, family = poisson(link = "log"), data = data3.8, 
    offset = log(t))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-12.789   -3.426   -0.938    3.079   10.137  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.91028    0.02164  -42.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 669.45  on 22  degrees of freedom
Residual deviance: 669.45  on 22  degrees of freedom
AIC: 812.62

Number of Fisher Scoring iterations: 5
###

从summary(model3.8)返回的模型摘要中可以获得，根据计算得到

这表明被捕球迷数与总观众数的比率期望在0.4左右，即被捕球迷数预计为总观众数的40%

c小题

画出被捕人数与观众人数的散点图以及预测方程。利用残差区分比期望被捕人数更大和更小的球队。

题解：

预测方程即为，模型图像为直线，R中可以通过abline()添加直线

> attach(data3.8)
> mu=exp(model3.8$coe)
> plot(Y~t)
> abline(0,mu)

在直线之上的是被捕人数大于期望值的球队，在直线之下的是被捕人数小于期望值的球队，可以通过残差的正负来判断，也可以用以下命令可以返回被捕人数大于期望值和小于期望值的球队的队名，这种方式和利用残差正负进行判断的方式是等价的

> rownames(data3.8)[Y<mu*t]    #比期望值小，在直线之下，残差小于0的球队
###
 [1] "米德尔斯堡" "伊普斯维奇" "莱切斯特城" "布莱克本"   "水晶宫"     "史云顿"    
 [7] "谢菲尔德联" "斯托克城"   "巴恩斯利"   "米尔沃尔"   "侯城"       "曼彻斯特城"
[13] "普利茅斯"   "雷丁"       "奥威"  
###
> rownames(data3.8)[Y>=mu*t]    #比期望值大，在直线之上，残差大于0的球队
###
[1] "阿斯顿维拉" "布拉德福城" "利兹联"     "伯恩茅斯"   "西布朗维奇" "汉德斯菲德"
[7] "伯明翰"     "什鲁斯伯里"
###

d小题

用负二项分布拟合模型 $log[E(Y)/t]=\alpha$ . 将 $\hat{\alpha }$ 及其SE与(b)中结果比较。基于这个信息和散布参数及其SE的估计值，泊松假设合适吗？

题解：

负二项对数模型用MASS包内的glm.nb()函数，不采用glm()进行负二项对数拟合的原因是我们暂不知晓样本的散布参数，虽然可以用logtrans()函数确定散布参数倒数θ的取值，但是用glm.nb()可以一步到位，比较方便。不过glm.nb()没有offset参数(位移)，但是我们可以调整formula参数的表达进行带位移项的拟合，这个调整也适用于glm()函数

> library(MASS)
> model3.8_nb=glm.nb(Y~offset(log(t)),data=data3.8,init.theta=1,link='log')
> summary(model3.8_nb)
###
Call:
glm.nb(formula = Y ~ offset(log(t)), data = data3.8, init.theta = 3.135631071, 
    link = "log")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2049  -0.7464  -0.1857   0.6129   1.5568  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9052     0.1200  -7.546 4.49e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(3.1356) family taken to be 1)

    Null deviance: 24.15  on 22  degrees of freedom
Residual deviance: 24.15  on 22  degrees of freedom
AIC: 244.24

Number of Fisher Scoring iterations: 1


              Theta:  3.136 
          Std. Err.:  0.920 

 2 x log-likelihood:  -240.236
###
#比较两个模型截距的估计值和标准误
> summary(model3.8)$coe
###
  	            Estimate Std. Error   z value Pr(>|z|)
(Intercept) -0.9102802 0.02163712 -42.07031        0
###
> summary(model3.8_nb)$coe
###
 	             Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -0.9051888  0.1199579 -7.545888 4.492147e-14
###

从模型摘要可见负二项对数模型的θ的估计值为3.136，标准误为0.920，则散布参数的估计值为，说明样本具有一定的超散布性

两个模型对α的估计值相似，但是负二项对数模型的α的标准误相对较高，用模型的偏差进行比较，也可以得出负二项对数模型的拟合效果更好。综上，泊松假设并不适合。

习题3.19

表3.4给出了火车事故数据

表3.4 英国涉及火车的碰撞事故
年份	火车里程	火车碰撞	火车-道路碰撞	年份	火车里程	火车碰撞	火车-道路碰撞
2003	518	0	3	1988	443	2	4
2002	516	1	3	1987	397	1	6
2001	508	0	4	1986	414	2	13
2000	503	1	3	1985	418	0	5
1999	505	1	2	1984	389	5	3
1998	487	0	4	1983	401	2	7
1997	463	1	1	1982	372	2	3
1996	437	2	2	1981	417	2	2
1995	423	1	2	1980	430	2	2
1994	415	2	4	1979	426	3	3
1993	425	0	4	1978	430	2	4
1992	430	1	4	1977	425	1	8
1991	439	2	6	1976	426	2	12
1990	431	1	2	1975	436	5	2
1989	436	4	4

a小题

比较只有截距项的撞车比率的泊松GLM和具有时间趋势项的GLM，这两个模型的偏差分别是35.1和23.5。通过上述结果，能将这29年里每年的撞车事件数看作具有相同参数的独立泊松变量吗？

题解：

不带时间效应的模型偏差为35.1，加入时间效应的模型偏差为23.5，这其实已经说明带时间效应的模型拟合效果更好

另外我们可以通过对两个模型的偏差做差，得到的值近似服从卡方分布，自由度是两个模型的参数数量差，对该题来说，这正是β=0的似然比检验，自由度df=1

用R辅助计算P值

> Dev1=35.1;Dev2=23.5
> p.value=1-pchisq(Dev1-Dev2,df=1);p.value 
###
[1] 0.0006595182
###

β显著性的似然比检验P值很小，说明时间对撞击次数的影响还是存在的，即使模型的偏差并没有减少很多。这样来看这29年的撞车事件数并不能看作具有相同参数的独立泊松变量。

b小题

3.3.6节拟合了负二项模型。1975年之后第x年撞车比率的估计值为 $e^{-4.20}\cdot (e^{-0.0337})^{x}=(0.015)(0.967)^{x}$ . ML估计 $\hat{\beta }=-0.0337$ 的SE=0.0130。建立 $H_{0}:\beta =0$ 对 $H_{a}:\beta \neq 0$ 的Wald检验.

题解：

题目要求的检验便是参数β的显著性Wald检验

β的估计值除以其标准误便是显著性检验的Wald统计量，其近似服从标准正态分布

我们可以在R中进行相同的拟合，3.3.6节中给出散布参数D=0.099

> data3.19=read.csv('table3.4.csv')    #数据只取了年份、火车里程、火车-道路碰撞次数
> data3.19[,1]=data3.19[,1]-1975
> head(data3.19)
###
  年份 里程 碰撞
1   28  518    3
2   27  516    3
3   26  508    4
4   25  503    3
5   24  505    2
6   23  487    4
###
> model3.19_3=glm(碰撞~年份+offset(log(里程)),data=data3.19,family=negative.binomial(theta=1/0.099,link='log'))
> summary(model3.19_3)$coe        #获取模型的参数估计和检验
###
               Estimate Std. Error    t value     Pr(>|t|)
(Intercept) -4.19997478 0.20170528 -20.822334 3.658918e-18
年份        -0.03366993 0.01326265  -2.538703 1.720186e-02
###

R中得到的参数估计和标准误SE与书中一致，这里进行的显著性检验就是Wald检验，可见P值约为0.0172，并没有比0.05小很多，但是依然能够作为拒绝原假设的依据。

c小题

β的似然比95%置信区间为(-0.060,-0.008).求出事故率的年乘积效应的区间，解释结果。

题解：

变换模型的表达形式，有

于是事故率的年乘积效应就是，β的似然比95%置信区间题目给出为(-0.060,-0.008)，通过指数变换可以得到的95%置信区间

用R辅助计算

> c(exp(-0.06),exp(-0.008))
###
[1] 0.9417645 0.9920319
###

计算出的95%置信区间约为(0.942,0.992)，说明每到下一年，有95%的把握估计该年的事故率相比上一年事故率减少0.8%至5.8%

第四章 Logistic回归

习题4.1

一项研究利用logistic回归确定与Y=癌症是否缓解（1=是）相关联的特征量。最重要的解释变量是通过对病人注射氚，标记胸苷后，测量细胞繁殖的标记指数（LI）。该研究给出被“标记”细胞的百分比。表4.8给出了分组数据，表4.9是以LI 预测 $\pi =P(Y=1)$ 的logistic回归模型的结果。

表4.8 习题4.1关于癌症缓解的数据
LI	案例数	缓解数	LI	案例数	缓解数
8	2	0	22	2	1
10	2	0	24	1	0
12	3	0	26	1	1
14	3	0	28	1	1
16	3	0	32	1	0
18	1	1	34	1	1
20	3	2	38	3	2

表4.9 习题4.1的电脑输出结果

Standard Likelihood Ratio 95% Chi-

Parameter Estimate Error Confidence Limits Square

Intercept -3.7771 1.3786 -6.9946 -1.4097 7.51

li 0.1449 0.0593 0.0425 0.2846 5.96

Scale 1.0000 0.0000 1.0000 1.0000

LR Statistics For Type 3 Analysis

Chi-

Source DF Square Pr > Chisq

li 1 8.30 0.0040

Obs li nc nr pi_hat lower upper

1 8 2 0 0.06797 0.01121 0.31925

2 10 2 0 0.08879 0.01809 0.34010

...

a小题

说明当LI=8时，软件如何得到 $\hat{\pi }=0.068$ .

题解：

由于表4.9已经给出了模型的拟合结果，接下来R只用于辅助计算，不再次拟合模型。

根据表4.9得到的结果，模型可表示为

将LI=8代入模型，可以得出的logit值，通过反解公式

就可以得到当LI=8时的值，在R中可以进行如下计算得到

> T.logit=function(x){exp(x)/(1+exp(x))}
> T.logit(-3.7771+0.1449*8)
###
[1] 0.06799525
###

于是得到

b小题

证明当LI=26.0时， $\hat{\pi}=0.50$ .

题解：

当时，logit值为0

令模型的线性部分为0，反解出LI的值即可

即解方程

解得

c小题

证明当LI=8时 $\hat{\pi}$ 的变化率为0.009，当LI=26时为0.036

题解：

将模型表示为与LI的函数，有

求LI在各个取值时的变化率，可以对上式求导，有

将LI=8和LI=26分别代入上式便可得到变化率，在R中可以进行如下计算得到

> g=expression(exp(-3.7771+0.1449*x)/(1+exp(-3.7771+0.1449*x)))
> D(g,'x')
###
exp(-3.7771 + 0.1449 * x) * 0.1449/(1 + exp(-3.7771 + 0.1449 * 
    x)) - exp(-3.7771 + 0.1449 * x) * (exp(-3.7771 + 0.1449 * 
    x) * 0.1449)/(1 + exp(-3.7771 + 0.1449 * x))^2
###
> x=8
> eval(D(g,'x'))
###
[1] 0.009182588
###
> x=26
> eval(D(g,'x'))
###
[1] 0.03622415
###

计算得到LI=8时的变化率约为0.009；LI=26时的变化率约为0.036

d小题

LI的下四分位数和上四分位数分别为14和28。证明 $\pi$ 在这两个值之间从0.15增加到0.57，增幅为0.42

题解：

依然通过将LI的取值代入模型函数来计算

在R中可以进行如下运算

> g=expression(exp(-3.7771+0.1449*x)/(1+exp(-3.7771+0.1449*x)))
> x=14;a=eval(g);a        #LI=14时的预测概率
###
[1] 0.1482365
###
> x=28;b=eval(g);b        #LI=28时的预测概率
###
[1] 0.5695707
###
> b-a                     #增幅
###
[1] 0.4213342
###

可得到当LI=14时；当LI=28时，增幅为0.42

e小题

证明当LI增加1，缓解的优势的估计值扩大1.16倍

题解：

在logistic模型中，优势可以表示为

x每增加1，优势便扩大 $e^{\beta}$ 倍，于是该题我们要求的便是

根据表4.9可知，则的计算为

> exp(0.1449)
###
[1] 1.155924
###

得到当LI增加1，缓解的优势的估计值扩大约1.16倍

习题4.2

续上题。利用表4.9的信息：

a小题

建立LI效应的Wald检验，并解释结果

题解：

根据表4.9中的信息，LI的效应估计值，标准误

Wald统计量为，表4.9中已经给出了 $z^{2}$ 的值为5.96

在大样本下z近似服从标准正态分布，则 $z^{2}$ 近似服从自由度df=1的卡方分布

LI效应的Wald检验P值计算

> 1-pchisq(5.96,df=1)
###
[1] 0.01463404
###

P值约等于0.015，小于0.05，可以认为LI的效应是有显著性意义的

b小题

建立相应于LI增加1个单位优势比的Wald置信区间，并解释结果

题解：

由上文可知，缓解的优势可以表示为

那么LI增加1个单位的优势比就是

求的95%Wald置信区间可以从求β的95%Wald置信区间开始，再通过指数变换得到

> beta=0.1449
> SE=0.0593
> a=c(beta-qnorm(1-0.05/2)*SE,beta+qnorm(1-0.05/2)*SE);a    #β的置信区间
###
[1] 0.02867414 0.26112586
###
> exp(a)        #exp(β)的置信区间
###
[1] 1.029089 1.298391
###

得到LI增加1个单位的优势比即eβ的95%Wald置信区间约为(1.029,1.299)

这说明LI每增加1个单位，我们有95%的把握认为优势会变为原来的1.029到1.299倍，总的来说优势是随着LI上升的

c小题

建立LI效应的似然比检验，并解释结果

题解：

本题所给出的样本量并不大，Wald检验的功效和可信度不如似然比检验，表4.9已经给出了似然比检验的结果

Chi-

Source DF Square Pr > Chisq

li 1 8.30 0.0040

似然比统计量值为8.30，自由度为1，P值为0.004

检验结果与Wald检验相同，可以认为LI的效应是有显著性意义的，不过似然比检验的结果给出了比Wald检验更强烈的证据(似然比检验的P值更小)

d小题

建立优势比的似然比置信区间，并解释结果

题解：

本题依然是求的置信区间，依然是从β的置信区间入手，不过本次是用β的似然比置信区间

表4.9已经给出了β的95%似然比置信区间为(0.0425,0.2846)，对其进行指数变换即可得出的95%似然比置信区间

> exp(c(0.0425,0.2846))
###
[1] 1.043416 1.329230
###

得的95%似然比置信区间约为(1.0434,1.3292)，与Wald置信区间的结论相似，LI每增加1个单位，我们有95%的把握认为优势会变为原来的1.0434到1.3292倍

小结

以上是从广义线性模型和Logistic回归两章选的习题的练习结果。Logistic回归模型也算是广义线性模型中的一种，其应用比较广泛，所以书上总共用了两个章节讲解Logistic回归模型。本次关于Logistic回归模型的习题还是刚上手的，具体的知识还没仔细思考过，用的还都是在第三章广义线性模型中所了解的知识。

再次声明本人只是一名小小的本科生，题目可能做错，欢迎批评指正和交流。希望能帮到大家。