統計解析R演習6

注意事項

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

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

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

install.packages("ggplot2")
install.packages("foreign")
install.packages("car")
install.packages("sandwich")

library(ggplot2)
library(foreign)
library(car)
library(sandwich)

データのダウンロード

• Stata pressで公開されているデータを使用します。
• プログラムはGermán Rodríguez氏やここのWebサイトを参考にしています。
• 以下の要領でデータを読み込みます。

npub <- read.dta("http://www.stata-press.com/data/lf2/couart2.dta")

• データの要約は以下の通り。

summary(npub)

      art            fem           mar           kid5             phd       
 Min.   : 0.000   Men  :494   Single :309   Min.   :0.0000   Min.   :0.755  
 1st Qu.: 0.000   Women:421   Married:606   1st Qu.:0.0000   1st Qu.:2.260  
 Median : 1.000                             Median :0.0000   Median :3.150  
 Mean   : 1.693                             Mean   :0.4951   Mean   :3.103  
 3rd Qu.: 2.000                             3rd Qu.:1.0000   3rd Qu.:3.920  
 Max.   :19.000                             Max.   :3.0000   Max.   :4.620  
      ment       
 Min.   : 0.000  
 1st Qu.: 3.000  
 Median : 6.000  
 Mean   : 8.767  
 3rd Qu.:12.000  
 Max.   :77.000  

• 平均と分散を比較すると、過分散な状態であることがわかります。

c(mean=mean(npub$art), var=var(npub$art), ratio=var(npub$art)/mean(npub$art))

    mean      var    ratio 
1.692896 3.709742 2.191358 

• ヒストグラムは以下のようになります。

ggplot(npub, aes(x=art)) + geom_histogram()

• 変数mentとartの散布図は以下のようになります。

ggplot(npub) + geom_point(aes(x = ment, y = art, col = 2)) + 
  guides(colour = "none")

ポアソン回帰モデルの推定

• まず、全変数を投入したモデルnpub.poisson1を推定します。

npub.poisson1 <- glm(art~fem+mar+kid5+phd+ment, data=npub, 
                    family = poisson(link = "log"))
summary(npub.poisson1)

Call:
glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson(link = "log"), 
    data = npub)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.5672  -1.5398  -0.3660   0.5722   5.4467  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.304617   0.102981   2.958   0.0031 ** 
femWomen    -0.224594   0.054613  -4.112 3.92e-05 ***
marMarried   0.155243   0.061374   2.529   0.0114 *  
kid5        -0.184883   0.040127  -4.607 4.08e-06 ***
phd          0.012823   0.026397   0.486   0.6271    
ment         0.025543   0.002006  12.733  < 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: 1817.4  on 914  degrees of freedom
Residual deviance: 1634.4  on 909  degrees of freedom
AIC: 3314.1

Number of Fisher Scoring iterations: 5

• 変数mentとartの散布図にモデル推定結果を重ねます。

ggplot(npub, aes(x = ment, y = art, col = 2)) + geom_point() + 
  geom_smooth(method = "glm", method.args = list(family = "poisson")) +
  guides(colour = "none")

推定結果の評価

• 最大対数尤度

logLik(npub.poisson1)

'log Lik.' -1651.056 (df=6)

• 自由度909での5%水準の\(\chi^2\)値は、以下のように得られます。

qchisq(0.95, df.residual(npub.poisson1))

[1] 980.2518

• モデルがデータに当てはまるかどうかをみるために、devianceとPearson’s \(\chi^2\)値を計算します。
• Deviance

npub.poisson1$deviance

[1] 1634.371

• Pearson’s \(\chi^2\)値

sum((residuals(npub.poisson1,"pearson"))^2)

[1] 1662.547

• これらの結果から、モデルは当てはまりが良くないことがわかります。
• 過分散とdispersion parameterを以下のように計算します。
• 過分散

npub.poisson1$deviance/npub.poisson1$df.res

[1] 1.797988

• dispersion parameter

sum(residuals(npub.poisson1, type="pearson")^2)/npub.poisson1$df.res

[1] 1.828984

モデル選択

• まずANOVAを用いて、モデルnpub.poisson1の推定結果を評価します。
• ANOVA(Type I)はanova()関数を用いて計算します。

anova(npub.poisson1, test="Chisq")

Analysis of Deviance Table

Model: poisson, link: log

Response: art

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   914     1817.4              
fem   1   23.029       913     1794.4 1.596e-06 ***
mar   1    0.251       912     1794.1  0.616719    
kid5  1   17.388       911     1776.7 3.047e-05 ***
phd   1   10.499       910     1766.2  0.001194 ** 
ment  1  131.868       909     1634.4 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• Type IのANOVAでは、変数marに効果がないと、誤って検出されてしまいます。
• Type IIとType IIIのANOVAはcarパッケージのAnova()関数を用いて研鑽します。
• ANOVA(TypeII)

car::Anova(npub.poisson1, type=c("II"))

Analysis of Deviance Table (Type II tests)

Response: art
     LR Chisq Df Pr(>Chisq)    
fem    17.075  1  3.593e-05 ***
mar     6.431  1    0.01121 *  
kid5   22.082  1  2.613e-06 ***
phd     0.236  1    0.62696    
ment  131.868  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• ANOVA(Type III)

car::Anova(npub.poisson1, type=c("III"))

Analysis of Deviance Table (Type III tests)

Response: art
     LR Chisq Df Pr(>Chisq)    
fem    17.075  1  3.593e-05 ***
mar     6.431  1    0.01121 *  
kid5   22.082  1  2.613e-06 ***
phd     0.236  1    0.62696    
ment  131.868  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• 次に、変数phdを減らしたモデルnpub.poisson2を推定します。

npub.poisson2 <- glm(art~fem+mar+kid5+ment, data=npub, 
                     family = poisson(link = "log"))
summary(npub.poisson2)

Call:
glm(formula = art ~ fem + mar + kid5 + ment, family = poisson(link = "log"), 
    data = npub)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.6436  -1.5408  -0.3583   0.5623   5.3986  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.34517    0.06012   5.741 9.41e-09 ***
femWomen    -0.22530    0.05461  -4.125 3.70e-05 ***
marMarried   0.15218    0.06107   2.492   0.0127 *  
kid5        -0.18499    0.04014  -4.609 4.05e-06 ***
ment         0.02576    0.00195  13.212  < 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: 1817.4  on 914  degrees of freedom
Residual deviance: 1634.6  on 910  degrees of freedom
AIC: 3312.3

Number of Fisher Scoring iterations: 5

• npub.poisson1とnpub.poisson2のAICを比較します。

c(AIC.model1=AIC(npub.poisson1), AIC.model2=AIC(npub.poisson2))

AIC.model1 AIC.model2 
  3314.113   3312.349 

• ANOVAによりnpub.poisson1とnpub.poisson2のAICを比較します。

anova(npub.poisson1, npub.poisson2, test="Chisq")

Analysis of Deviance Table

Model 1: art ~ fem + mar + kid5 + phd + ment
Model 2: art ~ fem + mar + kid5 + ment
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       909     1634.4                     
2       910     1634.6 -1  -0.2362    0.627

参考：Robust SE

• Robust SEは以下の方法で計算できます。

result <- npub.poisson2
cov.result <- vcovHC(result, type="HC0")
std.err <- sqrt(diag(cov.result))
r.est <- cbind(Estimate= coef(result), "Robust SE" = std.err,
               "Pr(>|z|)" = 2 * pnorm(abs(coef(result)/std.err), lower.tail=FALSE),
               LL = coef(result) - 1.96 * std.err,
               UL = coef(result) + 1.96 * std.err)
r.est

               Estimate   Robust SE     Pr(>|z|)          LL          UL
(Intercept)  0.34517404 0.082177620 2.665130e-05  0.18410591  0.50624218
femWomen    -0.22530260 0.071201599 1.554611e-03 -0.36485773 -0.08574746
marMarried   0.15217527 0.083642215 6.885659e-02 -0.01176348  0.31611401
kid5        -0.18499262 0.055799227 9.153637e-04 -0.29435910 -0.07562613
ment         0.02576117 0.003491732 1.609692e-13  0.01891738  0.03260497