拓端tecdat|R語言泊松Poisson回歸模型預測人口死亡率和期望壽命

原文鏈接：http://tecdat.cn/?p=18782

 

本文我們討論了期望壽命的計算。人口統計模型的起點是死亡率表。但是，這種假設有偏差，因為它假設生活條件不會得到改善。為了正確處理問題，我們使用了更完整的數據，其中死亡人數根據x歲而定，還包括日期t。

1.  
    
2.  
   DE=read.table("DE.txt",skip = 3,header=TRUE)
3.  
   EXPS=read.table("EXPS.txt",skip = 3,header=TRUE)

 我們用 Dx，t表示死亡人數，Ex，t表示暴露人數。因此，對於在日期t上x歲的某人，在該年死亡的概率為 qx，t = Dx，t / Ex，t。這些數據存儲在矩陣中進行可視化，存儲在數據庫中進行回歸。

1.  
    
2.  
   QF[QF==0]=NA
3.  
   QH[QH==0]=NA

 必須進行一些修改以避免出現零值的問題，因為（i）我們求出比率（ii）然后我們對數化）。我們可以可視化為x和t的函數。

persp(log(QF))

 或

1.  
    
2.  
   persp3d(ages,annees,log(QH),col="light blue")

 

為了模擬qx，t的演化，我們可以從Lee＆Carter（1992）的模型中獲得啟發，該模型  假設log （qx，t）= Ax + Bx⋅Kt。A =（A0，A1，⋯，A110）在某種程度上是log（qx，t）。K =（K1816，K1817，⋯，K2015）使我們了解生活條件的改善，一年內死亡的可能性降低。這些改善不是均勻的，因此我們使用B =（B0，B1，⋯，B110）來使改善取決於l '年齡。

為了估計參數A，B和K，我們嘗試使用二項式模型。B（Ex，t，qx，t），這是人壽保險的基本模型。這里Dx，t〜B（Ex，t，exp [ Ax + Bx⋅Kt]）。

另一個線索是使用小數定律，即如果概率低（一年中的死亡概率就是這種情況），則二項式定律可以近似由泊松分布。我們在這里用到了Poisson回歸，其解釋變量為年齡x，年t和暴露量為偏移變量。唯一的問題是它不是線性回歸。我們這里有非線性模型，因為E [ Dx，t] =（exp[log（Ex，t）+ Ax + Bx⋅Kt]）。

1.  
    
2.  
   gnm( DH ~ offset(log(EH) + as.factor(age) +
3.  
   Multas.factor(age,as.factor(annee),
4.  
   family = poisson(link="log")

 我們有估計系數A ^，B ^和K ^。

1.  
    
2.  
   Ax=reg$coefficients[2:111]
3.  
   Bx=reg$coefficients[112:222]
4.  
   Kt=reg$coefficients[223:length(reg$coefficients)]

我們可以表示三組系數。首先 A ^表示平均變化，

plot(ages[-1],Ax)

 

我們還可以用 K ^來繪制時間。

 

 

同樣，該模型不可被識別。簡而言之，改善沒有任何意義。我們可以表示-K ^，它的優點是描述了生活條件的改善。最后，讓我們作圖-B ^

 

 

困難在於，為了預測期望壽命，我們需要針對t的大值（尚未觀察到）計算qt，x。例如，某人可能想知道q50,2020（對於1970年出生的人）。我們要使用q50,2020 = exp（A ^ 50 + B ^ 50 K ^ 2020）。問題是K ^ 2020不屬於估計數量K ^。

這個想法是Lee＆Carter（1992）的初衷，我們可以嘗試指數模型或線性模型（在1950年以后的原始K ^序列上）

1.  
    
2.  
   lm(log(Kt[idx])~ann[idx])
3.  
   futur=2016:2125
4.  
    
5.  
   lm(Kt[idx]~ann[idx])
6.  
    
7.  
   points(futur,pr,col="blue")

 

然后，我們可以根據過去的數據建立一系列預測，q ^ x，t = exp [A ^ x + B ^ x K ^ t]，以及未來數據q〜x，t = exp [A ^ x + B ^ x K〜t]。

我們保留過去的數據，這里是1880年死亡的概率

1.  
    
2.  
   plot(BASE$x[BASE$t==1880],BASE$pred[BASE$t==1880],
3.  
   log="y")

 

 

同樣，我們在未來（此處為2050年）使用這兩種模型

1.  
    
2.  
   BASE2$Qpred1=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt1)
3.  
    
4.  
    
5.  
   plot(BASE2$x[BASE2$t==2050],BASE2$Qpred1[BASE2$t==
6.  
   2050],log="y")

 

 

用於指數預測

 對於線性預測，對1968年出生的人，我們有第二年死亡的概率

1.  
    
2.  
   if(sbase$t[i]<= 2015)
3.  
   {vq[i]=BASE[ BASE$x==sbase$x[i]) & BASE$t==sbase$t[i]),"Qpred"]
4.  
   if(sbase$t[i] <2015)
5.  
   {vq[i]=BASE2[(BASE2$x==sbase$x[i]) & (BASE2$t==sbase$t[i]),"Qpred2"]
6.  
    

 

 

左邊是我們模型估算值，右邊是預測值。

要計算出生時的期望壽命，我們使用以下代碼

1.  
   sum(cumprod(exp(-vq[1:110])))
2.  
   [1] 77.62047

 然后，我們可以做函數可視化這種期望壽命的演變

1.  
    
2.  
   vP = cumprod(exp(-(sbase$vq[1:110])))
3.  
   sum(vP)}

 

1.  
    
2.  
   ANN =1930:2010
3.  
   plot(ANN ,E2)

 

 

如果我們看一下變化，我們發現每年（大約）有0.25的變化

 

 

 

另一方面，如果我們采用保留Kt指數變化的預測，則可以得出

 

 

結果不符合實際，它更少地考慮曲線的變化。

 

 

---

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應用案例

2.面板平滑轉移回歸(PSTR)分析案例實現

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優度檢驗

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現

7.在R語言中實現Logistic邏輯回歸

8.python用線性回歸預測股票價格

9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標