統計解析R演習7

注意事項

• 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
• 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
• Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

パッケージのインストール

• 今回の演習では、以下のパッケージを使います。

install.packages("ggplot2")
install.packages("AER")
install.packages("MASS")
install.packages("pscl")

library(ggplot2)
library(AER)
library(MASS)
library(pscl)

データの読み込み

• 応用計量経済に関するAERパッケージに含まれる高齢者の医療機関診察回数のデータNMES1988を用います。
• プログラムはこちらのWebサイトを参考にしています。

data("NMES1988")
# help(NMES1988)

• データの要約

summary(NMES1988)

     visits          nvisits           ovisits            novisits       
 Min.   : 0.000   Min.   :  0.000   Min.   :  0.0000   Min.   :  0.0000  
 1st Qu.: 1.000   1st Qu.:  0.000   1st Qu.:  0.0000   1st Qu.:  0.0000  
 Median : 4.000   Median :  0.000   Median :  0.0000   Median :  0.0000  
 Mean   : 5.774   Mean   :  1.618   Mean   :  0.7508   Mean   :  0.5361  
 3rd Qu.: 8.000   3rd Qu.:  1.000   3rd Qu.:  0.0000   3rd Qu.:  0.0000  
 Max.   :89.000   Max.   :104.000   Max.   :141.0000   Max.   :155.0000  
   emergency          hospital           health        chronic           adl      
 Min.   : 0.0000   Min.   :0.000   poor     : 554   Min.   :0.000   normal :3507  
 1st Qu.: 0.0000   1st Qu.:0.000   average  :3509   1st Qu.:1.000   limited: 899  
 Median : 0.0000   Median :0.000   excellent: 343   Median :1.000                 
 Mean   : 0.2635   Mean   :0.296                    Mean   :1.542                 
 3rd Qu.: 0.0000   3rd Qu.:0.000                    3rd Qu.:2.000                 
 Max.   :12.0000   Max.   :8.000                    Max.   :8.000                 
       region          age          afam         gender     married        school     
 northeast: 837   Min.   : 6.600   no :3890   female:2628   no :2000   Min.   : 0.00  
 midwest  :1157   1st Qu.: 6.900   yes: 516   male  :1778   yes:2406   1st Qu.: 8.00  
 west     : 798   Median : 7.300                                       Median :11.00  
 other    :1614   Mean   : 7.402                                       Mean   :10.29  
                  3rd Qu.: 7.800                                       3rd Qu.:12.00  
                  Max.   :10.900                                       Max.   :18.00  
     income        employed   insurance  medicaid  
 Min.   :-1.0125   no :3951   no : 985   no :4004  
 1st Qu.: 0.9122   yes: 455   yes:3421   yes: 402  
 Median : 1.6982                                   
 Mean   : 2.5271                                   
 3rd Qu.: 3.1728                                   
 Max.   :54.8351                                   

• 平均と分散を計算すると、過分散なデータであることがわかります。

c(mean=mean(NMES1988$visits), var=var(NMES1988$visits), ratio=var(NMES1988$visits)/mean(NMES1988$visits))

     mean       var     ratio 
 5.774399 45.687117  7.912013 

• ヒストグラムを表示すると、0の多いデータであることがわかります。

ggplot(NMES1988, aes(x=visits)) + geom_histogram()

ポアソン回帰モデル

• まず、ポアソン回帰モデルを推定します。

nmes.poi <- glm(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                data=NMES1988, 
                family = poisson(link = "log"))
summary(nmes.poi)

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、ポアソン回帰モデルでは0を過小推計していることがわかります。

c(obs=sum(NMES1988$visits==0), poi=sum(dpois(0,exp(predict(nmes.poi)))))

      obs       poi 
683.00000  48.98773 

負の二項分布モデル

• 負の二項分布モデルを推定します。

nmes.nb <- glm.nb(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                  data=NMES1988)
summary(nmes.nb)

Call:
glm.nb(formula = visits ~ hospital + health + chronic + adl + 
    region + age + afam + gender + married + school + income + 
    insurance, data = NMES1988, init.theta = 1.217628984, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0269  -0.9883  -0.2981   0.3160   5.4784  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.246772   0.206612   6.034 1.60e-09 ***
hospital         0.215474   0.020175  10.680  < 2e-16 ***
healthpoor       0.279125   0.050012   5.581 2.39e-08 ***
healthexcellent -0.346856   0.060976  -5.688 1.28e-08 ***
chronic          0.173100   0.012184  14.207  < 2e-16 ***
adllimited       0.094163   0.042546   2.213  0.02688 *  
regionnortheast  0.115242   0.043328   2.660  0.00782 ** 
regionmidwest   -0.014302   0.039829  -0.359  0.71953    
regionwest       0.141049   0.044639   3.160  0.00158 ** 
age             -0.045946   0.026043  -1.764  0.07770 .  
afamyes         -0.042551   0.051885  -0.820  0.41216    
gendermale      -0.099978   0.033922  -2.947  0.00321 ** 
marriedyes      -0.047230   0.035684  -1.324  0.18565    
school           0.024527   0.004606   5.325 1.01e-07 ***
income          -0.001069   0.005534  -0.193  0.84686    
insuranceyes     0.237042   0.041343   5.734 9.83e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 5778.8  on 4405  degrees of freedom
Residual deviance: 5045.5  on 4390  degrees of freedom
AIC: 24346

Number of Fisher Scoring iterations: 1

              Theta:  1.2176 
          Std. Err.:  0.0340 

 2 x log-likelihood:  -24312.2180 

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、負の二項分布モデルでは0を過小推計していることがわかります。

c(obs=sum(NMES1988$visits==0), nb=sum(dpois(0,exp(predict(nmes.nb)))))

      obs        nb 
683.00000  60.14426 

ハードルポアソン回帰モデル

• psclパッケージのhurdle()関数を用いて、ハードルポアソン回帰モデルを推定します。

nmes.hpoi <- hurdle(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                    data=NMES1988,
                    dist="poisson", zero.dist="binomial")
summary(nmes.hpoi)

Call:
hurdle(formula = visits ~ hospital + health + chronic + adl + region + age + 
    afam + gender + married + school + income + insurance, data = NMES1988, 
    dist = "poisson", zero.dist = "binomial")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-5.1651 -1.1572 -0.4635  0.5874 21.5303 

Count model coefficients (truncated poisson with log link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.999723   0.087833  22.767  < 2e-16 ***
hospital         0.157436   0.006105  25.789  < 2e-16 ***
healthpoor       0.231037   0.018427  12.538  < 2e-16 ***
healthexcellent -0.301275   0.031256  -9.639  < 2e-16 ***
chronic          0.098874   0.004770  20.729  < 2e-16 ***
adllimited       0.098458   0.016668   5.907 3.48e-09 ***
regionnortheast  0.097790   0.017845   5.480 4.26e-08 ***
regionmidwest   -0.030219   0.016943  -1.784   0.0745 .  
regionwest       0.102099   0.018131   5.631 1.79e-08 ***
age             -0.082740   0.010959  -7.550 4.37e-14 ***
afamyes          0.005235   0.022646   0.231   0.8172    
gendermale      -0.024091   0.014419  -1.671   0.0948 .  
marriedyes      -0.077677   0.014910  -5.210 1.89e-07 ***
school           0.018090   0.001971   9.179  < 2e-16 ***
income          -0.003003   0.002336  -1.286   0.1986    
insuranceyes     0.105460   0.017990   5.862 4.57e-09 ***
Zero hurdle model coefficients (binomial with logit link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -1.198795   0.590182  -2.031 0.042232 *  
hospital         0.310636   0.092133   3.372 0.000747 ***
healthpoor       0.035582   0.167337   0.213 0.831610    
healthexcellent -0.315586   0.143892  -2.193 0.028291 *  
chronic          0.540807   0.046108  11.729  < 2e-16 ***
adllimited      -0.168423   0.129917  -1.296 0.194842    
regionnortheast  0.127225   0.124238   1.024 0.305818    
regionmidwest    0.077052   0.113200   0.681 0.496081    
regionwest       0.228673   0.133732   1.710 0.087276 .  
age              0.176535   0.075318   2.344 0.019085 *  
afamyes         -0.298554   0.129067  -2.313 0.020713 *  
gendermale      -0.507012   0.095884  -5.288 1.24e-07 ***
marriedyes       0.225581   0.102415   2.203 0.027621 *  
school           0.049193   0.012682   3.879 0.000105 ***
income           0.004608   0.018165   0.254 0.799742    
insuranceyes     0.642775   0.106945   6.010 1.85e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 22 
Log-likelihood: -1.604e+04 on 32 Df

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、ハードルポアソン回帰モデルでは0をうまく推計できている（0のカウント数が同じ）ことがわかります。

c(obs=sum(NMES1988$visits==0), hpoi=sum(predict(nmes.hpoi,type="prob")[,1]))

 obs hpoi 
 683  683 

ハードル負の二項分布モデル

• ハードル負の二項分布モデルを推定します。

nmes.hnb <- hurdle(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                    data=NMES1988,
                    dist="negbin", zero.dist="binomial")
summary(nmes.hnb)

Call:
hurdle(formula = visits ~ hospital + health + chronic + adl + region + age + 
    afam + gender + married + school + income + insurance, data = NMES1988, 
    dist = "negbin", zero.dist = "binomial")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.1824 -0.7139 -0.2688  0.3457 15.4518 

Count model coefficients (truncated negbin with log link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.787308   0.217027   8.235  < 2e-16 ***
hospital         0.209671   0.021218   9.882  < 2e-16 ***
healthpoor       0.280878   0.049363   5.690 1.27e-08 ***
healthexcellent -0.335310   0.065881  -5.090 3.59e-07 ***
chronic          0.124951   0.012479  10.013  < 2e-16 ***
adllimited       0.118884   0.042107   2.823  0.00475 ** 
regionnortheast  0.108557   0.044530   2.438  0.01478 *  
regionmidwest   -0.028365   0.041130  -0.690  0.49042    
regionwest       0.121093   0.045518   2.660  0.00781 ** 
age             -0.084018   0.027194  -3.090  0.00200 ** 
afamyes          0.017383   0.055452   0.313  0.75391    
gendermale      -0.017356   0.035800  -0.485  0.62782    
marriedyes      -0.099650   0.037308  -2.671  0.00756 ** 
school           0.019428   0.004720   4.116 3.85e-05 ***
income          -0.001309   0.005523  -0.237  0.81261    
insuranceyes     0.136054   0.044666   3.046  0.00232 ** 
Log(theta)       0.353544   0.042644   8.291  < 2e-16 ***
Zero hurdle model coefficients (binomial with logit link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -1.198795   0.590182  -2.031 0.042232 *  
hospital         0.310636   0.092133   3.372 0.000747 ***
healthpoor       0.035582   0.167337   0.213 0.831610    
healthexcellent -0.315586   0.143892  -2.193 0.028291 *  
chronic          0.540807   0.046108  11.729  < 2e-16 ***
adllimited      -0.168423   0.129917  -1.296 0.194842    
regionnortheast  0.127225   0.124238   1.024 0.305818    
regionmidwest    0.077052   0.113200   0.681 0.496081    
regionwest       0.228673   0.133732   1.710 0.087276 .  
age              0.176535   0.075318   2.344 0.019085 *  
afamyes         -0.298554   0.129067  -2.313 0.020713 *  
gendermale      -0.507012   0.095884  -5.288 1.24e-07 ***
marriedyes       0.225581   0.102415   2.203 0.027621 *  
school           0.049193   0.012682   3.879 0.000105 ***
income           0.004608   0.018165   0.254 0.799742    
insuranceyes     0.642775   0.106945   6.010 1.85e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.4241
Number of iterations in BFGS optimization: 24 
Log-likelihood: -1.206e+04 on 33 Df

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、ハードル負の二項分布モデルでは0をうまく推計できている（0のカウント数が同じ）ことがわかります。

c(obs=sum(NMES1988$visits==0), hnb=sum(predict(nmes.hnb,type="prob")[,1]))

obs hnb 
683 683 

ゼロ過剰ポアソン回帰モデル

• psclパッケージのzeroinfl()関数を用いてゼロ過剰ポアソン回帰モデルを推定します。

nmes.zip <- zeroinfl(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                    data=NMES1988,
                    dist="poisson")
summary(nmes.zip)

Call:
zeroinfl(formula = visits ~ hospital + health + chronic + adl + region + age + 
    afam + gender + married + school + income + insurance, data = NMES1988, 
    dist = "poisson")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-5.1537 -1.1578 -0.4631  0.5875 21.5527 

Count model coefficients (poisson with log link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.998750   0.088015  22.709  < 2e-16 ***
hospital         0.157479   0.006103  25.802  < 2e-16 ***
healthpoor       0.230867   0.018428  12.528  < 2e-16 ***
healthexcellent -0.300908   0.031218  -9.639  < 2e-16 ***
chronic          0.098987   0.004773  20.739  < 2e-16 ***
adllimited       0.098643   0.016667   5.919 3.25e-09 ***
regionnortheast  0.098031   0.017845   5.493 3.94e-08 ***
regionmidwest   -0.029868   0.016946  -1.763   0.0780 .  
regionwest       0.102042   0.018136   5.627 1.84e-08 ***
age             -0.082717   0.010978  -7.535 4.88e-14 ***
afamyes          0.005249   0.022665   0.232   0.8169    
gendermale      -0.024015   0.014424  -1.665   0.0959 .  
marriedyes      -0.077535   0.014914  -5.199 2.01e-07 ***
school           0.018137   0.001973   9.194  < 2e-16 ***
income          -0.002955   0.002330  -1.268   0.2047    
insuranceyes     0.105049   0.017999   5.836 5.34e-09 ***

Zero-inflation model coefficients (binomial with logit link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.257736   0.605357   2.078 0.037739 *  
hospital        -0.300904   0.092342  -3.259 0.001120 ** 
healthpoor      -0.025139   0.168428  -0.149 0.881352    
healthexcellent  0.268029   0.150898   1.776 0.075696 .  
chronic         -0.537627   0.046890 -11.466  < 2e-16 ***
adllimited       0.183301   0.131575   1.393 0.163581    
regionnortheast -0.116960   0.126137  -0.927 0.353800    
regionmidwest   -0.081213   0.115962  -0.700 0.483714    
regionwest      -0.221981   0.136222  -1.630 0.103196    
age             -0.190539   0.077341  -2.464 0.013754 *  
afamyes          0.303607   0.131008   2.317 0.020479 *  
gendermale       0.514769   0.097886   5.259 1.45e-07 ***
marriedyes      -0.237879   0.104473  -2.277 0.022790 *  
school          -0.047477   0.012984  -3.657 0.000256 ***
income          -0.004535   0.018505  -0.245 0.806417    
insuranceyes    -0.644546   0.108957  -5.916 3.31e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 39 
Log-likelihood: -1.604e+04 on 32 Df

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、ゼロ過剰ポアソン回帰モデルでは0をうまく推計できていることがわかります。

c(obs=sum(NMES1988$visits==0), zip=sum(predict(nmes.zip,type="prob")[,1]))

     obs      zip 
683.0000 682.9868 

ゼロ過剰負の二項分布モデル

• ゼロ過剰負の二項分布モデルを推定します。

nmes.zinb <- zeroinfl(visits~hospital+health+chronic+adl+region+age+afam+gender+married+school+income+insurance, 
                      data=NMES1988,
                      dist="negbin")
summary(nmes.zinb)

Call:
zeroinfl(formula = visits ~ hospital + health + chronic + adl + region + age + 
    afam + gender + married + school + income + insurance, data = NMES1988, 
    dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.2023 -0.7094 -0.2692  0.3482 15.0422 

Count model coefficients (negbin with log link):
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.8120470  0.2072655   8.743  < 2e-16 ***
hospital         0.1995827  0.0202149   9.873  < 2e-16 ***
healthpoor       0.2544212  0.0470370   5.409 6.34e-08 ***
healthexcellent -0.3016651  0.0624819  -4.828 1.38e-06 ***
chronic          0.1297708  0.0120666  10.755  < 2e-16 ***
adllimited       0.0842867  0.0403009   2.091 0.036489 *  
regionnortheast  0.1044950  0.0425585   2.455 0.014076 *  
regionmidwest   -0.0335426  0.0394323  -0.851 0.394970    
regionwest       0.1198416  0.0434925   2.755 0.005861 ** 
age             -0.0849993  0.0258513  -3.288 0.001009 ** 
afamyes         -0.0062532  0.0534015  -0.117 0.906782    
gendermale      -0.0299478  0.0344446  -0.869 0.384603    
marriedyes      -0.0919262  0.0357454  -2.572 0.010120 *  
school           0.0187213  0.0045494   4.115 3.87e-05 ***
income          -0.0005312  0.0052921  -0.100 0.920049    
insuranceyes     0.1534746  0.0436530   3.516 0.000438 ***
Log(theta)       0.4093154  0.0353695  11.573  < 2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      4.492282   1.504956   2.985  0.00284 ** 
hospital        -0.854501   0.474445  -1.801  0.07169 .  
healthpoor      -0.169905   0.504974  -0.336  0.73652    
healthexcellent  0.298400   0.280089   1.065  0.28671    
chronic         -1.175042   0.170809  -6.879 6.02e-12 ***
adllimited      -0.036422   0.371401  -0.098  0.92188    
regionnortheast -0.172214   0.263138  -0.654  0.51281    
regionmidwest   -0.378864   0.284730  -1.331  0.18332    
regionwest      -0.332403   0.304282  -1.092  0.27465    
age             -0.604502   0.196083  -3.083  0.00205 ** 
afamyes          0.444891   0.256862   1.732  0.08327 .  
gendermale       0.965100   0.232445   4.152 3.30e-05 ***
marriedyes      -0.617430   0.242156  -2.550  0.01078 *  
school          -0.087157   0.028580  -3.050  0.00229 ** 
income           0.001024   0.044678   0.023  0.98172    
insuranceyes    -0.977300   0.237684  -4.112 3.93e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 1.5058 
Number of iterations in BFGS optimization: 44 
Log-likelihood: -1.206e+04 on 33 Df

• 0をうまく推計できているかどうかを、実際に0の数をカウントして確認しましょう。その結果、ゼロ過剰ポアソン回帰モデルでは0を過大推計していることがわかります。

c(obs=sum(NMES1988$visits==0), zinb=sum(predict(nmes.zinb,type="prob")[,1]))

     obs     zinb 
683.0000 709.3136 

AICによるモデル間の比較

• AICを用いて、推定した６つのモデルを比較します。
• この結果、ハードル負の二項分布モデルのAICが最も小さくなります。このモデルは、0もうまく推計できているモデルとなっています。

c(AIC.model1=AIC(nmes.poi), AIC.model2=AIC(nmes.nb),
  AIC.model3=AIC(nmes.hpoi), AIC.model4=AIC(nmes.hnb),
  AIC.model5=AIC(nmes.zip), AIC.model5=AIC(nmes.zinb))

AIC.model1 AIC.model2 AIC.model3 AIC.model4 AIC.model5 AIC.model5 
  35807.64   24346.22   32134.63   24182.64   32134.23   24189.51