前言

第三章? 廣義線性模型

習題3.18

a小題

?b小題

?c小題

d小題

習題3.19

a小題

b小題

c小題

第四章? Logistic回歸

習題4.1

a小題

b小題

c小題

d小題

e小題

習題4.2

a小題

b小題

c小題

?d小題

?小結

前言

????????習題選自高等教育出版社譯制，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.62Number 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.24Number of Fisher Scoring iterations: 1Theta:  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回歸模型的習題還是剛上手的，具體的知識還沒仔細思考過，用的還都是在第三章廣義線性模型中所了解的知識。

? ? ? ? 再次聲明本人只是一名小小的本科生，題目可能做錯，歡迎批評指正和交流。希望能幫到大家。