以下はWindows10、R version 4.2.1で実行しています。
データを読み込むにあたって、作業場所の指定とそこにデータが置いてあることが前提になります。 R本体を操作している場合は、「ファイル」→「ディレクトリの変更」でデータの置いてある場所を指定するのが簡単です。 RStudioを操作している場合は、例えばデスクトップの「R」という名前のフォルダにデータがあるとすると
setwd("C:/Users/ユーザー名/Desktop/R")を最初に実行するのがよいでしょう。ユーザー名は各自異なるので注意です。
使用するデータを読み込み、表示します。なお、R version 4.1以前で読み込む際には
data1 <- read.csv("ファイル名.csv", fileEncoding = "UTF-8-BOM")のように、エンコードのオプションを指定する必要があります。version 4.2以降では必要ありませんので、以下ではオプションを指定していません。
data1 <- read.csv("ch10-12ex.csv")使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、RやRStudioからクリックでインストールすることもできます。
install.packages(“stargazer”, dependencies = TRUE)パッケージの機能が使えるように読み込みます。
library(dplyr)##
## 次のパッケージを付け加えます: 'dplyr'## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, unionlibrary(car)## 要求されたパッケージ carData をロード中です##
## 次のパッケージを付け加えます: 'car'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## recodelibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(mfx)## 要求されたパッケージ sandwich をロード中です## 要求されたパッケージ lmtest をロード中です## 要求されたパッケージ zoo をロード中です##
## 次のパッケージを付け加えます: 'zoo'## 以下のオブジェクトは 'package:base' からマスクされています:
##
## as.Date, as.Date.numeric## 要求されたパッケージ MASS をロード中です##
## 次のパッケージを付け加えます: 'MASS'## 以下のオブジェクトは 'package:dplyr' からマスクされています:
##
## select## 要求されたパッケージ betareg をロード中ですlibrary(erer)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(MASS)使用する変数の準備をします。
data1 <- data1 %>% mutate(
age=q10_5,
agecl=as.factor(recode(age,'20:29=1;30:39=2;40:49=3;50:59=4;60:69=5')),
fem=recode(gender,'2=1; 1=0'),
mar=recode(q10_7,'1=0; 2=1; 3:9=NA'),
edu=as.factor(recode(q10_12,'1:2=1; 3:4=2; 5:6=3; 7:9=NA')),
wrk=as.factor(recode(q10_11,'1:3=1; 4:7=2; 9:10=3; 8=NA; 11:99=NA')),
smk=recode(q10_1,'1=1; 2:3=0; 9=NA'),
hap=as.factor(recode(q04_1,'1=5;2=4;3=3;4=2;5=1;9=NA')))リコードの確認をします。
table(data1$gender, data1$fem, useNA = "ifany")##
## 0 1
## 1 1556 0
## 2 0 1191table(data1$q10_7, data1$mar, useNA = "ifany")##
## 0 1 <NA>
## 1 892 0 0
## 2 0 1734 0
## 3 0 0 93
## 9 0 0 28table(data1$q10_12, data1$edu, useNA = "ifany")##
## 1 2 3 <NA>
## 1 25 0 0 0
## 2 597 0 0 0
## 3 0 177 0 0
## 4 0 327 0 0
## 5 0 0 1391 0
## 6 0 0 184 0
## 7 0 0 0 12
## 9 0 0 0 34table(data1$q10_11, data1$wrk, useNA = "ifany")##
## 1 2 3 <NA>
## 1 77 0 0 0
## 2 329 0 0 0
## 3 773 0 0 0
## 4 0 303 0 0
## 5 0 65 0 0
## 6 0 207 0 0
## 7 0 20 0 0
## 8 0 0 0 113
## 9 0 0 435 0
## 10 0 0 335 0
## 11 0 0 0 49
## 99 0 0 0 41table(data1$age, data1$agecl)##
## 1 2 3 4 5
## 20 20 0 0 0 0
## 21 24 0 0 0 0
## 22 35 0 0 0 0
## 23 33 0 0 0 0
## 24 35 0 0 0 0
## 25 48 0 0 0 0
## 26 61 0 0 0 0
## 27 62 0 0 0 0
## 28 69 0 0 0 0
## 29 69 0 0 0 0
## 30 0 42 0 0 0
## 31 0 42 0 0 0
## 32 0 48 0 0 0
## 33 0 53 0 0 0
## 34 0 40 0 0 0
## 35 0 56 0 0 0
## 36 0 60 0 0 0
## 37 0 90 0 0 0
## 38 0 66 0 0 0
## 39 0 88 0 0 0
## 40 0 0 55 0 0
## 41 0 0 48 0 0
## 42 0 0 46 0 0
## 43 0 0 44 0 0
## 44 0 0 41 0 0
## 45 0 0 40 0 0
## 46 0 0 52 0 0
## 47 0 0 46 0 0
## 48 0 0 47 0 0
## 49 0 0 38 0 0
## 50 0 0 0 55 0
## 51 0 0 0 68 0
## 52 0 0 0 67 0
## 53 0 0 0 52 0
## 54 0 0 0 58 0
## 55 0 0 0 49 0
## 56 0 0 0 46 0
## 57 0 0 0 38 0
## 58 0 0 0 54 0
## 59 0 0 0 36 0
## 60 0 0 0 0 83
## 61 0 0 0 0 88
## 62 0 0 0 0 80
## 63 0 0 0 0 78
## 64 0 0 0 0 94
## 65 0 0 0 0 85
## 66 0 0 0 0 79
## 67 0 0 0 0 62
## 68 0 0 0 0 49
## 69 0 0 0 0 28table(data1$q10_1, data1$smk, useNA = "ifany")##
## 0 1 <NA>
## 1 0 419 0
## 2 721 0 0
## 3 1581 0 0
## 9 0 0 26table(data1$q04_1, data1$hap, useNA = "ifany")##
## 1 2 3 4 5 <NA>
## 1 0 0 0 0 310 0
## 2 0 0 0 1289 0 0
## 3 0 0 577 0 0 0
## 4 0 352 0 0 0 0
## 5 187 0 0 0 0 0
## 9 0 0 0 0 0 32ロジット・モデルを推定します。
logit <- glm(smk~agecl+fem+mar+edu+wrk, family=binomial(link=logit), data=data1)
stargazer(logit, type="text",star.cutoffs=c(0.1,0.05,0.01),keep.stat = c("n","ll","rsq"))##
## ==========================================
## Dependent variable:
## ---------------------------
## smk
## ------------------------------------------
## agecl2 0.124
## (0.227)
##
## agecl3 0.454*
## (0.233)
##
## agecl4 0.409*
## (0.235)
##
## agecl5 0.081
## (0.249)
##
## fem -1.087***
## (0.148)
##
## mar 0.166
## (0.143)
##
## edu2 -0.052
## (0.177)
##
## edu3 -0.593***
## (0.140)
##
## wrk2 -0.129
## (0.152)
##
## wrk3 -0.490***
## (0.173)
##
## Constant -1.156***
## (0.223)
##
## ------------------------------------------
## Observations 2,416
## Log Likelihood -986.252
## ==========================================
## Note: *p<0.1; **p<0.05; ***p<0.0120代に比べて40代、50代でいずれも10％水準で有意に喫煙確率が高いことがわかります。
また、女性、大学・大学院、無職の場合、いずれもベースカテゴリーに比べて1％水準で有意に喫煙確率が低いことがわかります。
既婚であることは有意ではなく、影響があるかどうかわからないという結果となっています。
マクファデンの疑似決定係数を計算します。
1-(logit$deviance/logit$null.deviance)## [1] 0.05886078平均的な限界効果を計算します。
logitmfx(smk~agecl+fem+mar+edu+wrk, data=data1, atmean=FALSE) ## Call:
## logitmfx(formula = smk ~ agecl + fem + mar + edu + wrk, data = data1,
## atmean = FALSE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## agecl2 0.0159058 0.0296909 0.5357 0.592158
## agecl3 0.0621881 0.0346040 1.7971 0.072314 .
## agecl4 0.0550559 0.0339508 1.6216 0.104881
## agecl5 0.0102393 0.0319074 0.3209 0.748282
## fem -0.1263771 0.0158269 -7.9849 1.406e-15 ***
## mar 0.0203696 0.0172088 1.1837 0.236541
## edu2 -0.0064067 0.0217286 -0.2948 0.768109
## edu3 -0.0769470 0.0186583 -4.1240 3.723e-05 ***
## wrk2 -0.0158316 0.0181964 -0.8700 0.384277
## wrk3 -0.0573417 0.0187257 -3.0622 0.002197 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "agecl2" "agecl3" "agecl4" "agecl5" "fem" "mar" "edu2" "edu3"
## [9] "wrk2" "wrk3"平均的な限界効果で見ると、女性は男性に比べて12.6％ポイント喫煙確率が低いことがわかります。また、大学・大学院卒は高校以下卒に比べて7.7％ポイント喫煙確率が低くなっています。無職の場合は、正規職と比べて5.7％ポイント確率が低くなります。
平均での限界効果を計算します。
logitmfx(smk~agecl+fem+mar+edu+wrk, data=data1, atmean=TRUE) ## Call:
## logitmfx(formula = smk ~ agecl + fem + mar + edu + wrk, data = data1,
## atmean = TRUE)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## agecl2 0.0151617 0.0283837 0.5342 0.593224
## agecl3 0.0601689 0.0340133 1.7690 0.076897 .
## agecl4 0.0531418 0.0332130 1.6000 0.109591
## agecl5 0.0097479 0.0304242 0.3204 0.748664
## fem -0.1241191 0.0156490 -7.9314 2.166e-15 ***
## mar 0.0192825 0.0162425 1.1872 0.235165
## edu2 -0.0060777 0.0205860 -0.2952 0.767815
## edu3 -0.0733143 0.0178385 -4.1099 3.958e-05 ***
## wrk2 -0.0149767 0.0171492 -0.8733 0.382488
## wrk3 -0.0546154 0.0178874 -3.0533 0.002264 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "agecl2" "agecl3" "agecl4" "agecl5" "fem" "mar" "edu2" "edu3"
## [9] "wrk2" "wrk3"平均での限界効果もほぼ同じ結果となっています。
オッズ比を計算します。
logitor(smk~agecl+fem+mar+edu+wrk, data=data1) ## Call:
## logitor(formula = smk ~ agecl + fem + mar + edu + wrk, data = data1)
##
## Odds Ratio:
## OddsRatio Std. Err. z P>|z|
## agecl2 1.132344 0.257227 0.5471 0.584283
## agecl3 1.574769 0.366790 1.9497 0.051216 .
## agecl4 1.504734 0.354294 1.7354 0.082662 .
## agecl5 1.084191 0.270044 0.3245 0.745530
## fem 0.337309 0.050016 -7.3291 2.318e-13 ***
## mar 1.180416 0.168837 1.1596 0.246192
## edu2 0.949642 0.168130 -0.2918 0.770403
## edu3 0.552493 0.077344 -4.2382 2.253e-05 ***
## wrk2 0.878884 0.133290 -0.8513 0.394620
## wrk3 0.612367 0.105668 -2.8421 0.004482 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1オッズ比も推定された係数の正負に従って、１より大きくなったり小さくなったりしていることがわかります。
具体的には、女性は男性に比べてオッズが0.34倍、大学・大学院卒は高校以下卒に比べてオッズが0.55倍、無職の場合は、正規職と比べてオッズが0.61倍となっていることがわかります。
順序ロジット・モデルを推定します。
ologit <- polr(hap~agecl+fem+mar+edu+wrk, data=data1, method = "logistic")
summary(ologit)##
## Re-fitting to get Hessian## Call:
## polr(formula = hap ~ agecl + fem + mar + edu + wrk, data = data1,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## agecl2 -0.231305 0.13514 -1.71160
## agecl3 -0.244836 0.14532 -1.68480
## agecl4 -0.544231 0.14589 -3.73046
## agecl5 -0.005085 0.14991 -0.03392
## fem 0.696633 0.08917 7.81199
## mar 1.419292 0.09717 14.60633
## edu2 0.348271 0.12129 2.87150
## edu3 0.388102 0.09692 4.00455
## wrk2 -0.222659 0.10290 -2.16393
## wrk3 -0.386496 0.10716 -3.60688
##
## Intercepts:
## Value Std. Error t value
## 1|2 -1.6598 0.1535 -10.8146
## 2|3 -0.3817 0.1419 -2.6903
## 3|4 0.7898 0.1425 5.5412
## 4|5 3.3983 0.1585 21.4392
##
## Residual Deviance: 6262.628
## AIC: 6290.628
## ( 344 個の観測値が欠損のため削除されました )疑似決定係数の計算で使用するために、推定に使用された対象をマークします。
data1$used <- TRUE;
data1$used[na.action(ologit)] <- FALSE;
table(data1$used)##
## FALSE TRUE
## 344 2403独立変数のないモデルを推定します。
ologitn <- polr(hap~1, data=data1, subset=used==1, method = "logistic")
summary(ologitn)##
## Re-fitting to get Hessian## Call:
## polr(formula = hap ~ 1, data = data1, subset = used == 1, method = "logistic")
##
## No coefficients
##
## Intercepts:
## Value Std. Error t value
## 1|2 -2.6205 0.0811 -32.3015
## 2|3 -1.4274 0.0516 -27.6417
## 3|4 -0.3885 0.0416 -9.3454
## 4|5 2.0018 0.0630 31.7745
##
## Residual Deviance: 6621.112
## AIC: 6629.112マクファデンの疑似決定係数を計算します。
1-(logLik(ologit)/logLik(ologitn))## 'log Lik.' 0.05414252 (df=14)stargazer(ologit, type="text",star.cutoffs=c(0.1,0.05,0.01),keep.stat = c("n","ll","rsq"))##
## ========================================
## Dependent variable:
## ---------------------------
## hap
## ----------------------------------------
## agecl2 -0.231*
## (0.135)
##
## agecl3 -0.245*
## (0.145)
##
## agecl4 -0.544***
## (0.146)
##
## agecl5 -0.005
## (0.150)
##
## fem 0.697***
## (0.089)
##
## mar 1.419***
## (0.097)
##
## edu2 0.348***
## (0.121)
##
## edu3 0.388***
## (0.097)
##
## wrk2 -0.223**
## (0.103)
##
## wrk3 -0.386***
## (0.107)
##
## ----------------------------------------
## Observations 2,403
## ========================================
## Note: *p<0.1; **p<0.05; ***p<0.0120代に比べて30代、40代は10％水準で有意に幸福度が低く、50代は1％水準で有意に幸福度が低くなっていることがわかります。 また、女性、既婚者、高学歴の人はいずれもベースカテゴリーに比べて1％水準で有意に幸福度が高く、正規以外の職や無職の人は、正規の人に比べてそれぞれ5％水準、1％水準で有意に幸福度が低くなっていることもわかります。
念のため、閾値の有意性の確認です。
ctable <- coef(summary(ologit))##
## Re-fitting to get Hessianp <- pnorm(abs(ctable[, "t value"]), lower.tail=FALSE)*2
(ctable <- cbind(ctable, "p value"=p))## Value Std. Error t value p value
## agecl2 -0.231305362 0.13514011 -1.71159672 8.697101e-02
## agecl3 -0.244836175 0.14532082 -1.68479770 9.202761e-02
## agecl4 -0.544231084 0.14588836 -3.73046269 1.911284e-04
## agecl5 -0.005084938 0.14991363 -0.03391911 9.729417e-01
## fem 0.696632790 0.08917485 7.81198735 5.629321e-15
## mar 1.419291939 0.09716964 14.60633062 2.559360e-48
## edu2 0.348270898 0.12128530 2.87150142 4.085269e-03
## edu3 0.388102368 0.09691525 4.00455402 6.213459e-05
## wrk2 -0.222658980 0.10289572 -2.16392851 3.046984e-02
## wrk3 -0.386495793 0.10715528 -3.60687570 3.099060e-04
## 1|2 -1.659773268 0.15347475 -10.81463427 2.934608e-27
## 2|3 -0.381736606 0.14189456 -2.69028355 7.139133e-03
## 3|4 0.789814558 0.14253417 5.54122967 3.003550e-08
## 4|5 3.398273374 0.15850729 21.43922508 5.756422e-102オッズ比
係数の推定値が正の場合、オッズ比は１より大きく、負の場合は１より小さくなっています。 女性は、男性に比べてより高い幸福度を選ぶオッズが2.01倍となっていることがわかります。 また、既婚の人は未婚の人に比べてより高い幸福度を選ぶオッズが4.13倍となっています。
exp(coef(ologit))## agecl2 agecl3 agecl4 agecl5 fem mar edu2 edu3
## 0.7934971 0.7828328 0.5802878 0.9949280 2.0069834 4.1341921 1.4166160 1.4741807
## wrk2 wrk3
## 0.8003878 0.6794336限界効果
ここでは平均での限界効果が示されています。 女性は、最も低い幸福度（1）を選ぶ確率が男性に比べて3.4％ポイント低く、最も高い幸福度（5）を選ぶ確率が6.6％ポイント高くなっていることがわかります。 また、既婚の人は、最も低い幸福度（1）を選ぶ確率が未婚の人に比べて9.7％ポイント低く、最も高い幸福度（5）を選ぶ確率が10.7％ポイント高くなっていることがわかります。
me_ol <- ocME(w = ologit, digits=4)##
## Re-fitting to get Hessianme_ol## effect.1 effect.2 effect.3 effect.4 effect.5
## agecl2 0.0124 0.0214 0.0221 -0.0360 -0.0200
## agecl3 0.0134 0.0229 0.0232 -0.0385 -0.0208
## agecl4 0.0322 0.0527 0.0481 -0.0895 -0.0435
## agecl5 0.0003 0.0005 0.0005 -0.0008 -0.0005
## fem -0.0343 -0.0607 -0.0685 0.0973 0.0662
## mar -0.0969 -0.1402 -0.1021 0.2324 0.1068
## edu2 -0.0161 -0.0294 -0.0356 0.0466 0.0345
## edu3 -0.0203 -0.0353 -0.0375 0.0587 0.0345
## wrk2 0.0119 0.0206 0.0213 -0.0346 -0.0193
## wrk3 0.0210 0.0360 0.0365 -0.0604 -0.0331