この文書は、Allison (2009) （以下「テキスト」と呼ぶ）の中で示されているサンプルデータの分析結果を、R (R Core Team 2021) で再現したものである。テキストには Stata と Mplus のプログラムが掲載されているが、これらの商用プログラムを持っていないユーザもいるだろうし、R の普及のため、という目的もある。なお、この文書は R のプログラミングの経験がある程度ある読者を想定しており、R についてはほとんど解説しない。また、読者はテキストまたはテキストの翻訳を手元に持っていて、すでに読んでいるという想定で書かれているので、この文書単体で読むと意味のわからないところもあるかもしれない。
テキスト中の分析に用いられているデータは、Stata や SAS のデータ形式で公開されているので、これらを使って結果を再現している。ただし、一部データが公開されていないものもあるし、R では簡単に計算できないものもあるので、それらについては割愛している1。なお、この文書の pdf 版が京都大学学術情報リポジトリ KURENAI にも公開される予定である。
まずデータを読み込む。以下では最初に dplyr パッケージ (Wickham et al. 2021) をロードしているが、dplyr を使うのは次の節である（データの整形に使う）。ここでロードしているのは、dplyr を plm パッケージ (Croissant and Millo 2018) （あとで使う）よりあとにロードすると、どちらのパッケージも lag() という関数を持っているため、plm パッケージの lag() が使えなくなってしまうからである（plm パッケージの lag() もあとで使う）。
以下で用いる foreign パッケージ (R Core Team 2020) は、Stata や SAS, SPSS のような主要な統計パッケージで用いられている形式のデータを読み込むための関数を多く含んでいる。read.dta() もその一つ。
library(dplyr) # 3章あたりで使う plm パッケージの lag() と干渉するので先に読み込む
library(foreign)
nlsy <- read.dta("nlsy.dta") # Stata のデータを読み込む関数
head(nlsy) # 最初の6行だけ表示。ワイド形式のパネルデータ## momage anti90 anti92 anti94 gender childage hispanic black momwork married self90 self92 self94 pov90 pov92
## 1 21 1 1 1 1 8.000000 0 0 0 0 21 24 23 1 1
## 2 22 0 0 0 1 8.416667 0 0 1 0 20 24 24 0 0
## 3 18 5 5 5 0 8.083333 0 0 0 1 21 24 24 0 0
## 4 24 2 3 1 0 8.250000 0 0 0 0 23 21 21 0 0
## 5 22 1 0 0 1 9.333333 0 0 0 0 22 23 24 0 0
## 6 24 1 1 1 0 8.583333 0 0 0 0 19 21 24 0 0
## pov94
## 1 1
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0dim(nlsy) # 行数と列数を表示## [1] 581 16以下は 1990年と 1994年のデータを使った場合の回帰分析と 1990年と 1994年のスコアの差を用いた回帰分析。スコアの差を用いることで、時間的に不変な変数を統制した分析が可能。個々の年のデータで回帰分析を行う場合には、独立変数と時間的に不変な変数の間に相関がある場合にはバイアスが生じる。
# 1990年と1994年の間のスコアの差を計算
nlsy$ antidiff <- nlsy$ anti94 - nlsy$ anti90
nlsy$ selfdiff <- nlsy$ self94 - nlsy$ self90
nlsy$ povdiff <- nlsy$ pov94 - nlsy$ pov90
# 反社会的行動指標を従属変数、自尊心と貧困を独立変数とした回帰分析を 3種類
lm90 <- lm(data = nlsy, anti90 ~ self90 + pov90)
lm94 <- lm(data = nlsy, anti94 ~ self94 + pov94)
lmdiff <- lm(data = nlsy, antidiff ~ selfdiff + povdiff)
# 以下のようにすれば詳しい結果が見られる（割愛）
# summary(lm90)
# summary(lm94)
# summary(lmdiff)
# 3つの回帰分析の結果をまとめて表示
library(texreg)
knitreg(
list(lm90, lm94, lmdiff), # 引数は回帰分析の結果のリスト
digits = 3, # 表示桁の指定。省略可
custom.model.names = c("1990年", "1994年", "スコアの差"), # 表頭のラベル。省略可
caption.above = TRUE, # キャプションを上につける、という指定、省略可
caption = "OLSによる回帰分析の結果（表2.1の再現）", # キャプション、省略可
custom.coef.names = c(NA, "self", "pov", "self", "pov", "self", "pov") # 表2.1 と同じように係数を整理するおまじない。
) # 詳しくは help を見よ。省略可| 1990年 | 1994年 | スコアの差 | |
|---|---|---|---|
| (Intercept) | 2.375*** | 2.888*** | 0.209*** |
| (0.384) | (0.447) | (0.063) | |
| self | -0.050** | -0.064** | -0.056*** |
| (0.019) | (0.021) | (0.015) | |
| pov | 0.595*** | 0.547*** | -0.036 |
| (0.126) | (0.148) | (0.128) | |
| R2 | 0.050 | 0.041 | 0.023 |
| Adj. R2 | 0.047 | 0.037 | 0.019 |
| Num. obs. | 581 | 581 | 581 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
上の texreg パッケージ (Leifeld 2013) の knitreg() は複数の回帰分析の結果をきれいにまとめて HTML や Word 等の形式で表示するための関数。この表の最初の 3つの行が回帰係数で、カッコ内は標準誤差を示す。R^2 は決定係数、Adj. R^2 は調整済み決定係数、Num. obs. が計算に用いたサンプルサイズである。
上のスコアの差を使ったモデル (lmdiff) に、時間的に不変の変数を加えたモデルの推定。
lmdiff2 <- update(
lmdiff,
~ . + self90 + pov90 + black + hispanic + childage + married + gender + momage + momwork
)
# 表2.2 の再現
print( # print() は省略可。こうすると、表示桁を減らしてすっきりした出力になる
summary(lmdiff2),
digits = 2 # 表示する有効桁の指定
)##
## Call:
## lm(formula = antidiff ~ selfdiff + povdiff + self90 + pov90 +
## black + hispanic + childage + married + gender + momage +
## momwork, data = nlsy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.32 -0.91 -0.10 0.84 5.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.550 1.360 -0.4 0.686
## selfdiff -0.060 0.020 -3.0 0.002 **
## povdiff 0.031 0.156 0.2 0.845
## self90 -0.018 0.025 -0.7 0.483
## pov90 0.121 0.178 0.7 0.499
## black -0.100 0.155 -0.6 0.518
## hispanic 0.084 0.164 0.5 0.611
## childage 0.220 0.107 2.0 0.041 *
## married -0.206 0.154 -1.3 0.181
## gender 0.101 0.126 0.8 0.426
## momage -0.040 0.030 -1.3 0.184
## momwork -0.153 0.140 -1.1 0.274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.5 on 569 degrees of freedom
## Multiple R-squared: 0.043, Adjusted R-squared: 0.025
## F-statistic: 2.3 on 11 and 569 DF, p-value: 0.0081まず、1990-1992年、1992-1994年の一階差分を使って回帰分析する。
# 1990年と1992年の差を計算
nlsy$antidiff1 <- nlsy$anti92 - nlsy$anti90
nlsy$selfdiff1 <- nlsy$self92 - nlsy$self90
nlsy$povdiff1 <- nlsy$pov92 - nlsy$pov90
# 1992年と1994年の差を計算
nlsy$antidiff2 <- nlsy$anti94 - nlsy$anti92
nlsy$selfdiff2<- nlsy$self94 - nlsy$self92
nlsy$povdiff2 <- nlsy$pov94 - nlsy$pov92
# 回帰分析の結果はあとでまとめて表示するが、summary(lmdif1) などとすれば簡単に結果は見られる
lmdif1 <- lm(data = nlsy, antidiff1 ~ selfdiff1 + povdiff1) # 1990-1992
lmdif2 <- lm(data = nlsy, antidiff2 ~ selfdiff2 + povdiff2) # 1992-1994次に 1990-1992 と 1992-1994 で変数の傾きなどに変化がなかったと仮定して両者を同時に推定する。これは上記の 1990-1992 の差分データと 1992-1994 の差分データを縦に統合し、一人あたり2行の記録のあるデータをまず作り、そのデータを使って回帰分析する2。まずデータをロング形式に変換する。
nlsy$id <- 1 : nrow(nlsy) # 個人 id をデータに加える
dat90_92 <- subset(nlsy, select = c(id, antidiff1, selfdiff1, povdiff1) ) # 1990-1992 の変化のデータ
dat92_94 <- subset(nlsy, select = c(id, antidiff2, selfdiff2, povdiff2) ) # 1992-1994 の変化のデータ
dat90_92$year <- "1990-1992"
dat92_94$year <- "1992-1994"
head(dat90_92)## id antidiff1 selfdiff1 povdiff1 year
## 1 1 0 3 0 1990-1992
## 2 2 0 4 0 1990-1992
## 3 3 0 3 0 1990-1992
## 4 4 1 -2 0 1990-1992
## 5 5 -1 1 0 1990-1992
## 6 6 0 2 0 1990-1992head(dat92_94)## id antidiff2 selfdiff2 povdiff2 year
## 1 1 0 -1 0 1992-1994
## 2 2 0 0 0 1992-1994
## 3 3 0 0 0 1992-1994
## 4 4 -2 0 0 1992-1994
## 5 5 0 1 0 1992-1994
## 6 6 0 3 0 1992-1994# 二つのデータフレームを縦に統合するためには、変数名がすべて一致していなければならないので、同じにする
names(dat90_92) <- c("id", "anti", "self", "pov", "year")
names(dat92_94) <- c("id", "anti", "self", "pov", "year")
# 二つのデータフレームを縦に結合
dat.combined <- rbind(dat90_92, dat92_94)
dat.combined <- dat.combined [order(dat.combined$id), ] # id 順にソート。必要ないがデータがうまく結合されたか確認しやすい
head(dat.combined, 12)## id anti self pov year
## 1 1 0 3 0 1990-1992
## 1100 1 0 -1 0 1992-1994
## 2 2 0 4 0 1990-1992
## 2100 2 0 0 0 1992-1994
## 3 3 0 3 0 1990-1992
## 3100 3 0 0 0 1992-1994
## 4 4 1 -2 0 1990-1992
## 4100 4 -2 0 0 1992-1994
## 5 5 -1 1 0 1990-1992
## 582 5 0 1 0 1992-1994
## 6 6 0 2 0 1990-1992
## 610 6 0 3 0 1992-1994上で作ったデータを使って回帰分析。GLS には gee パッケージ (Carey 2019) の gee() を使う。
# OLS で推定
lm.combined <- lm(anti~ self + pov + year, data = dat.combined)
# GEE で上と同じモデルを推定
# Stata 内の xtreg 関数における pa オプションは、gee 関数における constr = "exchangeable"と同じである
library(gee)
gls <- gee(data = dat.combined, anti ~ self + pov + year, id = id, family = "gaussian", corstr = "exchangeable")## (Intercept) self pov year1992-1994
## 0.04506966 -0.05509885 0.21293142 0.12195171knitreg(list(lmdif1, lmdif2, lm.combined, gls),
digits = 3,
include.scale = FALSE, # GEE ではスケールパラメータもデフォルトで表示されるがその省略を指定
custom.model.names = c("1990-92年", "1992-94年", "統合 OLS", "統合 GLS"),
caption.above = TRUE,
caption = "反社会的行動の一階差分回帰（表2.3の再現）",
custom.coef.names= c("切片", "自尊心", "貧困", "自尊心", "貧困", "自尊心", "貧困", "1992-1994ダミー")
) | 1990-92年 | 1992-94年 | 統合 OLS | 統合 GLS | |
|---|---|---|---|---|
| 切片 | 0.040 | 0.171** | 0.045 | 0.045 |
| (0.053) | (0.059) | (0.056) | (0.054) | |
| 自尊心 | -0.039** | -0.072*** | -0.055*** | -0.055*** |
| (0.014) | (0.016) | (0.010) | (0.011) | |
| 貧困 | 0.197 | 0.216 | 0.213* | 0.139 |
| (0.133) | (0.136) | (0.095) | (0.097) | |
| 1992-1994ダミー | 0.122 | 0.122 | ||
| (0.080) | (0.094) | |||
| R2 | 0.017 | 0.040 | 0.029 | |
| Adj. R2 | 0.013 | 0.036 | 0.027 | |
| Num. obs. | 581 | 581 | 1162 | 1162 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||||
“1990-92年,” “1992-94年,” “統合 OLS” はテキストと同じ推定値が得られている。統合 GLSの結果は、自尊心や切片の標準誤差などが少し異なる以外は概ね同じである。
まず、データをワイド形式からロング形式に変換した表（テキストの表2.4）を作る。ロング形式とワイド形式のあいだの変換には、reshape() が使える。
nlsy.long <- reshape(
nlsy,
direction = "long", # ロング形式への変換
varying = list(c( "anti90", "anti92", "anti94"), # 時間とともに変化する変数の指定
c("self90", "self92","self94"),
c("pov90", "pov92", "pov94") ),
v.names = c("anti", "self", "pov") # ロング形式での時変変数の名前の指定。
) # 省略すると最初の時点の変数名になる。
nlsy.long <- nlsy.long [order(nlsy.long$id), ] # id 順にデータを並べ替え（表2.4の再現のため）
# 表2.4 の再現
head(nlsy.long [, c("id", "time", "anti", "self", "pov", "gender")], 15)## id time anti self pov gender
## 1.1 1 1 1 21 1 1
## 1.2 1 2 1 24 1 1
## 1.3 1 3 1 23 1 1
## 2.1 2 1 0 20 0 1
## 2.2 2 2 0 24 0 1
## 2.3 2 3 0 24 0 1
## 3.1 3 1 5 21 0 0
## 3.2 3 2 5 24 0 0
## 3.3 3 3 5 24 0 0
## 4.1 4 1 2 23 0 0
## 4.2 4 2 3 21 0 0
## 4.3 4 3 1 21 0 0
## 5.1 5 1 1 22 0 1
## 5.2 5 2 0 23 0 1
## 5.3 5 3 0 24 0 1次にロング形式のデータを使って分析（表2.5）。
# 固定効果モデルの推定。今回は、各人のIDとTIMEをファクターとして回帰分析に投入。
lmfix <- lm(data = nlsy.long, anti ~ self + pov + factor(time) + factor(id))
lm.conventional <- update(lmfix, ~. - factor(id)) # 上のモデルから factor(id) を除去して推定
# summary(lmfix) などとするだけで結果は見られるが、表2.5と同じにするために長いスクリプトを書いた
knitreg(list(lmfix, lm.conventional),
digits = 3,
single.row = TRUE, # 標準誤差を係数の横に置く指定
custom.model.names = c("固定効果", "通常の OLS"),
caption.above = TRUE, # キャプションを上に
caption = "反社会的行動を自尊心と貧困に回帰させた結果（ダミー変数法、表2.5）",
custom.coef.map = list( # 表2.5 のように係数を配置するためのまじない。
"self" = NA, "pov" = NA, # 詳しくは help を見よ。
"factor(time)2" = "Time 2", "factor(time)3" = "Time 3",
"factor(id)2" = "id 2", "factor(id)3" = "id 3", "factor(id)4" = "id 4",
"factor(id)5" = "id 5", "factor(id)6" = "id 6", "factor(id)7" = "id 7",
"factor(id)8" = "id 8", "factor(id)9" = "id 9")
) | 固定効果 | 通常の OLS | |
|---|---|---|
| self | -0.055 (0.011)*** | -0.067 (0.011)*** |
| pov | 0.112 (0.093) | 0.518 (0.079)*** |
| Time 2 | 0.044 (0.059) | 0.051 (0.090) |
| Time 3 | 0.211 (0.059)*** | 0.223 (0.091)* |
| id 2 | -0.888 (0.819) | |
| id 3 | 4.131 (0.819)*** | |
| id 4 | 1.057 (0.820) | |
| id 5 | -0.536 (0.819) | |
| id 6 | 0.039 (0.820) | |
| id 7 | 2.170 (0.821)** | |
| id 8 | 0.910 (0.820) | |
| id 9 | -0.276 (0.820) | |
| R2 | 0.734 | 0.048 |
| Adj. R2 | 0.600 | 0.046 |
| Num. obs. | 1743 | 1743 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
次に、個人内平均からの偏差を使って固定効果モデルを推定する。この場合、plm パッケージの plm() を使うと簡単だろう。
# 個人内偏差を用いた固定効果モデルの推定。今回は plm 関数を使用。
library(plm) # データを pdata.frame クラスに変換。plm() はこのクラスのデータでないとエラーになる
nlsy.long.p <- pdata.frame(nlsy.long, index = c("id", "time"))
fix.within <- plm(
anti ~ self + pov,
data = nlsy.long.p,
effect = "twoways" # 時間もダミー変数として統制する指定。この引数を削除してモデル式に + time を加えても同じ
)
print(
summary(fix.within),
digits = 2
)## Twoways effects Within Model
##
## Call:
## plm(formula = anti ~ self + pov, data = nlsy.long.p, effect = "twoways")
##
## Balanced Panel: n = 581, T = 3, N = 1743
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.921 -0.477 -0.025 0.484 3.134
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## self -0.055 0.011 -5.2 2e-07 ***
## pov 0.112 0.093 1.2 0.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1200
## Residual Sum of Squares: 1200
## R-Squared: 0.024
## Adj. R-Squared: -0.47
## F-statistic: 14.4038 on 2 and 1158 DF, p-value: 6.6e-07表2.5 と同じ係数、標準誤差になる。
テキスト（翻訳）の pp.26～27 に解説されている統計量を得たければ、以下のようにすればよい。F検定ではなく、ラグランジュ乗数検定なので、多少の違いはあるがほぼ同じ結果が得られると思う (Long 1997)。
plmtest(fix.within) # H0: 個人間に異質性はない、という帰無仮説の検定##
## Lagrange Multiplier Test - (Honda) for balanced panels
##
## data: anti ~ self + pov
## normal = 23.807, p-value < 2.2e-16
## alternative hypothesis: significant effectsplmtest(fix.within, effect = "time") # H0: 時点ダミーの効果はすべてゼロという帰無仮説の検定##
## Lagrange Multiplier Test - time effects (Honda) for balanced panels
##
## data: anti ~ self + pov
## normal = 1.4772, p-value = 0.06982
## alternative hypothesis: significant effectsplmtest(fix.within, effect = "twoways") # H0: 時点ダミーも個人ダミーもすべてゼロという帰無仮説の検定##
## Lagrange Multiplier Test - two-ways effects (Honda) for balanced panels
##
## data: anti ~ self + pov
## normal = 17.879, p-value < 2.2e-16
## alternative hypothesis: significant effectsF検定したい場合は、以下のようにするとよい。pFtest() は二つの階層的な関係にあるモデルの適合度を比較する F検定をする関数である。
# plm パッケージを使う場合
pFtest( # H0: 個人間の異質性 = 0 という帰無仮説の検定
update(fix.within, ~. + time, effect = "individual"), # いわゆる固定効果モデル
update(fix.within, ~. + time, model = "pooling") # ロング形式のデータでただの OLS 推定
) ##
## F test for individual effects
##
## data: anti ~ self + pov + time
## F = 5.157, df1 = 580, df2 = 1158, p-value < 2.2e-16
## alternative hypothesis: significant effects# lm() で推定した場合
anova(lm.conventional, lmfix) # 表2.5 で推定した個人ダミーを投入したモデルと投入しなかったモデル## Analysis of Variance Table
##
## Model 1: anti ~ self + pov + factor(time)
## Model 2: anti ~ self + pov + factor(time) + factor(id)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1738 4124.8
## 2 1158 1151.2 580 2973.6 5.157 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1テキストには3種類の決定係数が示されているが、いずれも plm() の決定係数とは一致しない。テキストの誤植が疑われるが、詳細は不明。
交互作用を用いることによって、時定変数の効果をモデルに組み込むことが可能になる。
fix.time.interaction <- plm(
anti ~ time * (self + pov + black +
hispanic + childage + married + gender +
momage + momwork),
data = nlsy.long.p,
model = "within")
# 表2.6 の再現
print( # 単に summary(fix1) としても同じ結果が得られるが、こうすると四捨五入されて見やすい
summary(fix.time.interaction),
digits = 2)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = anti ~ time * (self + pov + black + hispanic +
## childage + married + gender + momage + momwork), data = nlsy.long.p,
## model = "within")
##
## Balanced Panel: n = 581, T = 3, N = 1743
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.80 -0.48 -0.02 0.48 3.13
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## time2 0.291 1.245 0.2 0.81
## time3 -0.444 1.258 -0.4 0.72
## self -0.034 0.017 -2.1 0.04 *
## pov 0.097 0.130 0.7 0.46
## time2:self -0.026 0.020 -1.3 0.20
## time3:self -0.023 0.021 -1.1 0.28
## time2:pov -0.112 0.152 -0.7 0.46
## time3:pov 0.099 0.155 0.6 0.52
## time2:black 0.250 0.144 1.7 0.08 .
## time3:black -0.110 0.144 -0.8 0.44
## time2:hispanic 0.190 0.154 1.2 0.22
## time3:hispanic 0.075 0.153 0.5 0.62
## time2:childage 0.076 0.100 0.8 0.45
## time3:childage 0.227 0.101 2.3 0.02 *
## time2:married -0.095 0.143 -0.7 0.51
## time3:married -0.176 0.143 -1.2 0.22
## time2:gender 0.041 0.118 0.3 0.73
## time3:gender 0.107 0.118 0.9 0.37
## time2:momage -0.027 0.028 -1.0 0.34
## time3:momage -0.043 0.028 -1.5 0.13
## time2:momwork 0.137 0.131 1.0 0.29
## time3:momwork -0.144 0.130 -1.1 0.27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1200
## Residual Sum of Squares: 1100
## R-Squared: 0.053
## Adj. R-Squared: -0.45
## F-statistic: 2.92116 on 22 and 1140 DF, p-value: 7.8e-06テキストの本文中（翻訳 p.28）で \(H_0:\) 交互作用効果はすべてゼロ、という帰無仮説を検定してあるので、ここでもやってみる。どういう検定なのか明示されていないが、Wald 検定する (Long 1997)。計算には aod パッケージ (Lesnoff and Lancelot 2012) の wald.test() を使った。
library(aod)
wald.test(
vcov(fix.time.interaction), # 検定したいモデルの係数の分散共分散行列
coef(fix.time.interaction), # 検定したいモデルの係数、全部
Terms = 5 : 22 # すべてゼロという帰無仮説を検定したい係数が何番目かを示すベクトル
)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 24.4, df = 18, P(> X2) = 0.14もちろん、主効果のみのモデルと交互作用付きのモデルの適合度の差を F検定しても上と同じことができる。
fix.main <- update(fix.time.interaction, ~ time + self + pov) # 主効果のみのモデル
pFtest(fix.time.interaction, fix.main)##
## F test for individual effects
##
## data: anti ~ time * (self + pov + black + hispanic + childage + married + ...
## F = 1.3545, df1 = 18, df2 = 1140, p-value = 0.1457
## alternative hypothesis: significant effectsplm() で model = “random” を指定すればランダム効果モデルも推定できる。lme4 パッケージ (Bates et al. 2015) の lmer() のほうが R では有名だろうが、固定効果モデルを plm() で推定するなら、ランダム効果モデルも plm() で推定したほうが簡単だろう。
# plm関数を利用してランダム効果を推定。
rdm <- plm(
anti ~ self + pov + time + black + hispanic + childage +
married + gender + momage + momwork,
data = nlsy.long.p,
model = "random"
)
rdm1 <- update(rdm, ~ self + pov + time) # 上で作った fix.main のランダム効果モデル版
knitreg(list(rdm, rdm1),
digits = 3,
custom.model.names = c( "時定変数含む", "時変変数のみ"),
caption.above = TRUE,
caption = "ランダム効果による推定（表2.7）"
) | 時定変数含む | 時変変数のみ | |
|---|---|---|
| (Intercept) | 2.531* | 2.668*** |
| (1.095) | (0.204) | |
| self | -0.062*** | -0.060*** |
| (0.010) | (0.010) | |
| pov | 0.247** | 0.296*** |
| (0.080) | (0.077) | |
| time2 | 0.047 | 0.047 |
| (0.059) | (0.059) | |
| time3 | 0.216*** | 0.216*** |
| (0.059) | (0.059) | |
| black | 0.227 | |
| (0.126) | ||
| hispanic | -0.218 | |
| (0.138) | ||
| childage | 0.088 | |
| (0.091) | ||
| married | -0.049 | |
| (0.126) | ||
| gender | -0.483*** | |
| (0.106) | ||
| momage | -0.022 | |
| (0.025) | ||
| momwork | 0.261* | |
| (0.115) | ||
| s_idios | 0.997 | 0.997 |
| s_id | 1.136 | 1.169 |
| R2 | 0.057 | 0.036 |
| Adj. R2 | 0.051 | 0.034 |
| Num. obs. | 1743 | 1743 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
s_idios は standard error of idiosyncratic でいわゆる残差の標準偏差である。s_id はランダム効果の標準偏差である。
以下のようにすればハウスマン検定もできる。テキストより p値がやや小さいが、検定しているモデルが若干違うのかもしれない。
phtest(rdm1, fix.main) # 変数などの指定はすべて同じで固定効果か、ランダム効果か、だけが違う二つのモデルを指定##
## Hausman Test
##
## data: anti ~ self + pov + time
## chisq = 12.959, df = 4, p-value = 0.01148
## alternative hypothesis: one model is inconsistentchisq がカイ二乗値であり、df が自由度、p-value がいわゆる p値。
# 個人内平均の計算
nlsy.long.p$ mself <- ave(nlsy.long.p$ self, nlsy.long.p$ id)
nlsy.long.p$ mpov <- ave(nlsy.long.p$ pov, nlsy.long.p$ id)
# 個人内偏差の計算
nlsy.long.p$ dself <- nlsy.long.p$ self - nlsy.long.p$ mself
nlsy.long.p$ dpov <- nlsy.long.p$ pov - nlsy.long.p$ mpov
head( # うまく計算できているか確認
nlsy.long.p [, c("id", "self", "mself", "dself", "pov", "mpov", "dpov")],
15)## id self mself dself pov mpov dpov
## 1-1 1 21 22.66667 -1.6666667 1 1 0
## 1-2 1 24 22.66667 1.3333333 1 1 0
## 1-3 1 23 22.66667 0.3333333 1 1 0
## 2-1 2 20 22.66667 -2.6666667 0 0 0
## 2-2 2 24 22.66667 1.3333333 0 0 0
## 2-3 2 24 22.66667 1.3333333 0 0 0
## 3-1 3 21 23.00000 -2.0000000 0 0 0
## 3-2 3 24 23.00000 1.0000000 0 0 0
## 3-3 3 24 23.00000 1.0000000 0 0 0
## 4-1 4 23 21.66667 1.3333333 0 0 0
## 4-2 4 21 21.66667 -0.6666667 0 0 0
## 4-3 4 21 21.66667 -0.6666667 0 0 0
## 5-1 5 22 23.00000 -1.0000000 0 0 0
## 5-2 5 23 23.00000 0.0000000 0 0 0
## 5-3 5 24 23.00000 1.0000000 0 0 0hyb <- update(rdm,
~ dself + dpov + mself + mpov + black + hispanic +
childage + married + gender + momage + momwork + time)
print(
summary(hyb),
digits = 2
)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = anti ~ dself + dpov + mself + mpov + black + hispanic +
## childage + married + gender + momage + momwork + time, data = nlsy.long.p,
## model = "random")
##
## Balanced Panel: n = 581, T = 3, N = 1743
##
## Effects:
## var std.dev share
## idiosyncratic 0.9942 0.9971 0.435
## individual 1.2896 1.1356 0.565
## theta: 0.5479
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.19 -0.74 -0.17 0.63 3.51
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 2.908 1.160 2.5 0.01 *
## dself -0.055 0.011 -5.2 2e-07 ***
## dpov 0.112 0.093 1.2 0.23
## mself -0.090 0.022 -4.1 4e-05 ***
## mpov 0.616 0.157 3.9 8e-05 ***
## black 0.111 0.132 0.8 0.40
## hispanic -0.280 0.139 -2.0 0.04 *
## childage 0.086 0.091 0.9 0.34
## married -0.128 0.129 -1.0 0.32
## gender -0.508 0.107 -4.8 2e-06 ***
## momage -0.011 0.025 -0.4 0.66
## momwork 0.164 0.119 1.4 0.17
## time2 0.044 0.059 0.8 0.45
## time3 0.211 0.059 3.6 3e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 1800
## Residual Sum of Squares: 1700
## R-Squared: 0.062
## Adj. R-Squared: 0.055
## Chisq: 114.878 on 13 DF, p-value: <2e-16上の出力中に theta という数値が出てくるが、これは時間的に変化しない観察されない異質性も独立変数だとみなしたときの決定係数である。
テキスト本文中の同時検定は、以下のようにする。
# 検定したい帰無仮説を示す行列を作る。
(coef.hyb <- coef(hyb)) # 係数もいるので、まず名前をつける## (Intercept) dself dpov mself mpov black hispanic childage married
## 2.90808043 -0.05515140 0.11247489 -0.09002608 0.61642756 0.11092723 -0.27990345 0.08574845 -0.12841272
## gender momage momwork time2 time3
## -0.50815795 -0.01133842 0.16411912 0.04439337 0.21073657# 帰無仮説を示すデザイン行列
design.matrix <- matrix(0, 2, length(coef.hyb))
colnames(design.matrix) <- names(coef.hyb)
design.matrix [1, "dself"] <- 1
design.matrix [1, "mself"] <- -1
design.matrix [2, "dpov"] <- 1
design.matrix [2, "mpov"] <- -1
design.matrix # こうすると、dself - mself = 0 & dpov - mpov = 0 という帰無仮説を検定してくれる## (Intercept) dself dpov mself mpov black hispanic childage married gender momage momwork time2 time3
## [1,] 0 1 0 -1 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 1 0 -1 0 0 0 0 0 0 0 0 0wald.test(
vcov(hyb),
b = coef.hyb,
L = design.matrix
)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 9.9, df = 2, P(> X2) = 0.0072X2 がカイ二乗値、P(> X2) が p値である。
上のようなデザイン行列を使った帰無仮説の指定には、行列の掛け算の知識が必要である。wald.test() は L の引数であるデザイン行列 (\(L\)) と b の引数（係数のベクトル）を要素とする列ベクトル（1列の行列, \(b\)）の掛け算の値がゼロである、という帰無仮説を検定する。上の例の場合、dself, mself, dpov, mpov の係数を \(d_s, \ m_s, \ d_p, \ m_p\)とすると、 \[ L \times b = \begin{pmatrix} d_s - m_s\\ d_p - m_p \end{pmatrix} = \begin{pmatrix} 0\\ 0 \end{pmatrix} \] が帰無仮説とされている。この例では二つの等式が帰無仮説となっているが、どちらもゼロである、というのが帰無仮説であり、対立仮説は、いずれかはゼロではない、ということになる。
データを読み込む。
library(foreign)
teenpov <- read.dta("teenpov.dta")
head(teenpov)## id pov1 mother1 spouse1 inschool1 hours1 pov2 mother2 spouse2 inschool2 hours2 pov3 mother3 spouse3
## 1 22 1 0 0 1 21 0 0 0 1 15 0 0 0
## 2 75 0 0 0 1 8 0 0 0 1 0 0 0 0
## 3 92 0 0 0 1 30 0 0 0 1 27 0 0 0
## 4 96 0 0 0 0 19 1 1 0 0 54 1 1 0
## 5 141 0 0 0 1 0 0 0 0 1 6 0 0 0
## 6 161 0 0 0 1 0 0 0 0 1 15 0 0 0
## inschool3 hours3 pov4 mother4 spouse4 inschool4 hours4 age black pov5 mother5 spouse5 inschool5 hours5
## 1 1 3 0 0 0 1 0 16 0 0 0 0 1 0
## 2 1 0 0 0 0 1 4 17 0 1 0 0 1 0
## 3 1 24 1 1 0 0 31 16 0 1 1 0 0 0
## 4 0 0 1 1 0 0 26 17 0 1 1 0 0 36
## 5 1 0 0 0 0 1 0 16 0 0 0 0 1 0
## 6 0 37 0 0 0 0 0 17 0 0 0 0 0 0dim(teenpov)## [1] 1151 28addmargins(
table(teenpov$pov1, teenpov$pov5)
)##
## 0 1 Sum
## 0 516 234 750
## 1 211 190 401
## Sum 727 424 1151ここでは tidyverse パッケージ (Wickham et al. 2019) を使ってデータを整形している。パイプ %>% などデータ整形時に便利な関数が含まれている。最後に使われている modelsummary パッケージ (Arel-Bundock 2021) の msummary() も複数の回帰分析の結果をまとめて表示する関数。
# データフレームを整える
library(tidyverse)
for_regression <- teenpov %>%
filter(pov1 != pov5) %>% # pov1 と pov5 の値が異なる行に限定
mutate( # 独立変数の差分を示す変数を計算
dmother = mother5 - mother1,
dspouse = spouse5 - spouse1,
dinschool = inschool5 - inschool1,
dhours = hours5 - hours1
)
# うまく計算できたか確認
head(for_regression [, c("mother5", "mother1", "dmother",
"spouse5", "spouse1", "dspouse",
"inschool5", "inschool1", "dinschool",
"hours5", "hours1", "dhours")])## mother5 mother1 dmother spouse5 spouse1 dspouse inschool5 inschool1 dinschool hours5 hours1 dhours
## 1 0 0 0 0 0 0 1 1 0 0 21 -21
## 2 0 0 0 0 0 0 1 1 0 0 8 -8
## 3 1 0 1 0 0 0 0 1 -1 0 30 -30
## 4 1 0 1 0 0 0 0 0 0 36 19 17
## 9 0 0 0 0 0 0 0 1 -1 40 0 40
## 12 1 0 1 1 0 1 0 1 -1 0 0 0ここから二項ロジットモデルの推定。
m1 <- glm(formula = pov5 ~ dmother + dspouse + dinschool + dhours,
data = for_regression,
family = binomial)
m2 <- update(m1,
formula. = ~ . + black + age,
data = for_regression)
m3 <- update(m2,
formula. = ~ . + mother1 + spouse1 + inschool1 + hours1,
data = for_regression)
library(texreg)
knitreg(
list(m1, m2, m3),
single.row = TRUE,
caption.above = TRUE,
caption = "表3.2: 2 時点間の差分値によるロジスティック回帰",
digits = 3
)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| (Intercept) | 0.539 (0.162)*** | 4.899 (1.644)** | 3.052 (1.826) |
| dmother | 0.730 (0.250)** | 0.744 (0.254)** | 0.909 (0.270)*** |
| dspouse | -1.002 (0.283)*** | -1.032 (0.292)*** | -1.022 (0.301)*** |
| dinschool | 0.343 (0.212) | 0.339 (0.218) | 0.639 (0.251)* |
| dhours | -0.034 (0.006)*** | -0.034 (0.006)*** | -0.034 (0.007)*** |
| black | -0.526 (0.216)* | -0.662 (0.226)** | |
| age | -0.258 (0.103)* | -0.196 (0.111) | |
| mother1 | 0.457 (0.460) | ||
| spouse1 | 0.442 (0.726) | ||
| inschool1 | 1.184 (0.471)* | ||
| hours1 | -0.002 (0.013) | ||
| AIC | 562.839 | 554.636 | 555.621 |
| BIC | 583.330 | 583.322 | 600.700 |
| Log Likelihood | -276.420 | -270.318 | -266.811 |
| Deviance | 552.839 | 540.636 | 533.621 |
| Num. obs. | 445 | 445 | 445 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
teenpov_long <- reshape( # 表2.4 を再現するのに reshape() を使っているので
teenpov, # そちらにも少し解説あり
direction = "long",
varying = list(
paste("pov", 1 : 5, sep = ""), # c("pov1", "pov2", "pov3", "pov4", "pov5") と書くのと同じ
paste("mother", 1 : 5, sep = ""),
paste("spouse", 1 : 5, sep = ""),
paste("inschool", 1 : 5, sep = ""),
paste("hours", 1 : 5, sep = "")
),
v.names = c("pov", "mother", "spouse", "inschool", "hours")
)
teenpov_long <- teenpov_long [order(teenpov_long$ id), ] # id 順にソート
head(teenpov_long, 15) # 中身の確認## id age black time pov mother spouse inschool hours
## 22.1 22 16 0 1 1 0 0 1 21
## 22.2 22 16 0 2 0 0 0 1 15
## 22.3 22 16 0 3 0 0 0 1 3
## 22.4 22 16 0 4 0 0 0 1 0
## 22.5 22 16 0 5 0 0 0 1 0
## 75.1 75 17 0 1 0 0 0 1 8
## 75.2 75 17 0 2 0 0 0 1 0
## 75.3 75 17 0 3 0 0 0 1 0
## 75.4 75 17 0 4 0 0 0 1 4
## 75.5 75 17 0 5 1 0 0 1 0
## 92.1 92 16 0 1 0 0 0 1 30
## 92.2 92 16 0 2 0 0 0 1 27
## 92.3 92 16 0 3 0 0 0 1 24
## 92.4 92 16 0 4 1 1 0 0 31
## 92.5 92 16 0 5 1 1 0 0 0条件付きロジットモデルは、survival パッケージ (Therneau 2021) の clogit() で計算できる。ランダム効果モデルの推定には、lme4 パッケージの glmer() 関数を使った。pglm パッケージ (Croissant 2021) の pglm() でも同様の計算ができる。
# 条件付き最尤法で推定
fixed.logit <- list() # 回帰分析の結果を代入するための空のリストを作成
library(survival)
fixed.logit [[1]] <- clogit(
pov ~ mother + spouse + inschool + hours + factor(time) + strata(id),
data = teenpov_long
) # summary(fixed.logit[[1]]) で詳しい結果が見られる
# GEE で推定
library(gee)
fixed.logit [[2]] <- gee(pov ~ mother + spouse + inschool + hours + factor(time),
id = id,
data = teenpov_long,
family = binomial,
corstr = "unstructured") # summary(fixed.logit [[2]]) で詳しい結果が見られる。## (Intercept) mother spouse inschool hours factor(time)2 factor(time)3
## -0.34610094 0.98621642 -1.25518917 -0.24016286 -0.02657698 0.21601628 0.14757282
## factor(time)4 factor(time)5
## 0.15275281 0.05967264# ランダム効果モデルの推定
library(lme4)
fixed.logit [[3]] <- glmer(
pov ~ mother + spouse + inschool + I(hours/10) + factor(time) + (1 | id),
data = teenpov_long,
family = binomial,
control = glmerControl(optimizer = "bobyqa") # アルゴリズムの指定。普通は必要ないが、デフォルトでは収束しなかったので。
)
library(texreg)
knitreg(
fixed.logit,
digits = 3,
single.row = TRUE,
include.bic = FALSE, # 以下8行はデフォルトで表示される
include.aic = FALSE, # 適合度指標などの表示を抑止する指定
include.scale =FALSE, # 詳しく見たければ summary(fixed.logit [[1]]) などとせよ
include.rsq = FALSE,
include.maxrs = FALSE,
include.missing = FALSE,
include.loglik = FALSE,
caption.above = TRUE,
caption = "条件付き最尤法とその他の方法によるロジットモデルの推定（表3.4）",
custom.model.names = c("固定効果（条件付き最尤法)", "GEE", "ランダム効果")
)| 固定効果（条件付き最尤法) | GEE | ランダム効果 | |
|---|---|---|---|
| mother | 0.582 (0.160)*** | 0.850 (0.093)*** | 1.083 (0.116)*** |
| spouse | -0.748 (0.175)*** | -0.930 (0.127)*** | -1.250 (0.150)*** |
| inschool | 0.272 (0.113)* | -0.045 (0.080) | -0.079 (0.097) |
| hours | -0.020 (0.003)*** | -0.021 (0.002)*** | |
| factor(time)2 | 0.332 (0.102)** | 0.223 (0.073)** | 0.284 (0.100)** |
| factor(time)3 | 0.335 (0.108)** | 0.172 (0.079)* | 0.222 (0.103)* |
| factor(time)4 | 0.433 (0.117)*** | 0.196 (0.081)* | 0.254 (0.108)* |
| factor(time)5 | 0.403 (0.128)** | 0.122 (0.091) | 0.165 (0.114) |
| (Intercept) | -0.543 (0.099)*** | -0.661 (0.124)*** | |
| hours/10 | -0.269 (0.028)*** | ||
| Num. events | 2169 | ||
| Num. obs. | 5755 | 5755 | 5755 |
| Num. groups: id | 1151 | ||
| Var: id (Intercept) | 1.288 | ||
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
Num. groups はグループの数（この例では、有効回答者数）、Var: id (Intercept) は、ランダム効果の分散の推定値である。
ランダム効果モデルの結果が若干テキストとは異なっている。アルゴリズムが異なるからかもしれない。
glmer() では、hours は 10 で割った値をモデルに投入している。独立変数のバラツキが変数によって顕著に異なる場合、固定効果の標準誤差も変数によって顕著に異なることが多いが、このような条件下では、lmer() や glmer() の推定が収束しにくくなることがある。このような場合、以下のような警告が出るので、他の独立変数よりも顕著にバラツキが大きい独立変数（そのような変数は係数も標準誤差も小さくなりがち）を上のように適当な定数で割ってバラツキを小さくするか、逆に顕著にバラツキの小さな独立変数に適当な定数をかけてバラツキを大きくするか（いずれも定数は1よりも大きな値）、するとうまく推定できることもある。
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?library(survival)
clogit.interaction <- list()
clogit.interaction [[1]] <- update(
fixed.logit [[1]],
~ . + mother:black
)
clogit.interaction [[2]] <- update(
fixed.logit [[1]],
~ . + time : (inschool + hours + black + age)
)
# 表3.5
library(modelsummary)
msummary(clogit.interaction, stars = TRUE)| Model 1 | Model 2 | |
|---|---|---|
| mother | 0.982*** | 0.684*** |
| (0.253) | (0.164) | |
| spouse | －0.783*** | －0.737*** |
| (0.178) | (0.178) | |
| inschool | 0.267* | －0.573* |
| (0.113) | (0.261) | |
| hours | －0.019*** | －0.001 |
| (0.003) | (0.008) | |
| factor(time)2 | 0.332** | 1.081** |
| (0.102) | (0.403) | |
| factor(time)3 | 0.334** | 1.875* |
| (0.108) | (0.783) | |
| factor(time)4 | 0.430*** | 2.873* |
| (0.117) | (1.157) | |
| factor(time)5 | 0.400** | 3.828* |
| (0.128) | (1.527) | |
| mother × black | －0.599* | |
| (0.290) | ||
| inschool × time | 0.251*** | |
| (0.069) | ||
| hours × time | －0.005* | |
| (0.002) | ||
| time × black | －0.182*** | |
| (0.048) | ||
| time × age | －0.056* | |
| (0.024) | ||
| Num.Obs. | 5755 | 5755 |
| R2 | 0.018 | 0.025 |
| AIC | 3053.9 | 3015.9 |
| BIC | 3105.0 | 3084.0 |
| Log.Lik. | －1517.948 | －1495.929 |
| n | 5755.000 | 5755.000 |
| nevent | 2169.000 | 2169.000 |
| statistic.log | 101.613 | 145.652 |
| p.value.log | 0.000 | 0.000 |
| statistic.sc | 98.417 | 139.225 |
| p.value.sc | 0.000 | 0.000 |
| statistic.wald | 93.640 | 131.170 |
| p.value.wald | 0.000 | 0.000 |
| statistic.robust | ||
| p.value.robust | ||
| r.squared.max | 0.420 | 0.420 |
| concordance | 0.597 | 0.619 |
| std.error.concordance | 0.014 | 0.014 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
teenpov_for_hybrid <- mutate(teenpov_long,
mmother = ave(teenpov_long$mother, teenpov_long$id),
mspouse = ave(teenpov_long$spouse, teenpov_long$id),
minschool = ave(teenpov_long$inschool, teenpov_long$id),
mhours = ave(teenpov_long$hours, teenpov_long$id),
dmother = mother - mmother,
dspouse = spouse - mspouse,
dinschool = inschool - minschool,
dhours = hours - mhours
)
library(pglm)
teenpov_for_hybrid <- pdata.frame(teenpov_for_hybrid, index = c("id", "time"))
hybrid <- pglm(
pov ~ dmother + dspouse + dinschool + I(dhours/10) +
mmother + mspouse + minschool + I(mhours/10) + black + age + time,
data = teenpov_for_hybrid,
model = "random",
family = binomial
)
library(texreg)
knitreg(
hybrid,
single.row = TRUE,
digits = 3,
caption.above = TRUE,
caption = "ハイブリッド法による推定（M で始まる変数は、個人内平均である。他方で、D で始まる変数は、個人内平均からの偏差を示す。表3.6）"
)| maximum | |
|---|---|
| (Intercept) | 1.893 (0.819)* |
| dmother | 0.594 (0.158)*** |
| dspouse | -0.807 (0.179)*** |
| dinschool | 0.275 (0.113)* |
| dhours/10 | -0.210 (0.032)*** |
| mmother | 1.080 (0.181)*** |
| mspouse | -2.147 (0.255)*** |
| minschool | -1.362 (0.202)*** |
| mhours/10 | -0.468 (0.058)*** |
| black | 0.572 (0.097)*** |
| age | -0.123 (0.050)* |
| time2 | 0.333 (0.101)** |
| time3 | 0.330 (0.107)** |
| time4 | 0.431 (0.115)*** |
| time5 | 0.391 (0.125)** |
| sigma | 1.118 (0.057)*** |
| DF | 16 |
| Log Likelihood | -3363.536 |
| AIC | 6759.072 |
| nobs | 5755 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |
R ではランダム効果モデルは、lme4 パッケージの glmer() 関数が有名だと思うが、計算が収束しなかったので、pglm パッケージの pglm() 関数を用いた。係数の最後にある sigma はおそらくランダム効果の標準偏差の推定値と思われるが、マニュアルにもヘルプにも特に記載がない。
上の分析結果の係数が本当は等しいといえる余地があるかどうかを Wald 検定する。ここでは car パッケージ (Fox and Weisberg 2019) の linearHypothesis() を使って計算した。aod パッケージの wald.test() でも同じことができる。linearHypothesis() のほうが直感的に帰無仮説を指定できるという利点があるが、扱えるモデルが限られるという難点がある。逆に wald.test() のほうは自分でデザイン行列を作らなければならないという難点があるが、係数のベクトルと係数の分散共分散行列さえあればどんなモデルで推定された係数でも扱えるという利点がある。
# Wald test
library(car)
wald.hybrid <- list()
(
wald.hybrid [[1]] <- linearHypothesis(hybrid, "dmother = mmother")
)## Linear hypothesis test
##
## Hypothesis:
## dmother - mmother = 0
##
## Model 1: restricted model
## Model 2: pov ~ dmother + dspouse + dinschool + I(dhours/10) + mmother +
## mspouse + minschool + I(mhours/10) + black + age + time
##
## Df Chisq Pr(>Chisq)
## 1
## 2 1 4.1606 0.04138 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1wald.hybrid [[2]] <- linearHypothesis(hybrid, "dspouse = mspouse")
wald.hybrid [[3]] <- linearHypothesis(hybrid, "dinschool = minschool")
wald.hybrid [[4]] <- linearHypothesis(hybrid, "I(dhours/10) = I(mhours/10)" )
# Combined test
wald.hybrid [[5]] <- linearHypothesis(
hybrid,
c("dmother = mmother",
"dspouse = mspouse",
"dinschool = minschool",
"I(dhours/10) = I(mhours/10)")
)
# あとは wald.hybrid とタイプすれば十分だが、表3.7を作ってみる
names(wald.hybrid) <- c("MOTHER", "SPOUSE", "SCHOOL", "HOURS", "同時検定（自由度 4）")
# 検定結果 (wald.hybrid) の中からカイ二乗値とp値だけとってくる
tab.wald <- sapply(wald.hybrid, function(x) c(x$ Chisq [2], x $ "Pr(>Chisq)" [2]))
library(knitr)
kable(
t(tab.wald),
digits = c(2, 3),
col.names = c("カイ2乗値", "p値"),
caption = "個人内偏差と個人内平均の係数の等値性の検定（表3.7）"
)| カイ2乗値 | p値 | |
|---|---|---|
| MOTHER | 4.16 | 0.041 |
| SPOUSE | 19.30 | 0.000 |
| SCHOOL | 49.88 | 0.000 |
| HOURS | 15.70 | 0.000 |
| 同時検定（自由度 4） | 79.06 | 0.000 |
上の結果を見るに、係数が等しいとはいえない。なお knitr パッケージ (Xie 2021) の kable() 関数は、表をきれいにして latex や html, word などの形式で出力するための関数である。
テキストでは、当時ランダム効果を仮定した順序ロジットモデルを推定するソフトウェアがなかったために、ランダム効果を仮定しない通常の順序ロジットモデルで係数を推定し、誤差相関を仮定したロバスト標準誤差を推定している。しかし、現在の R にはランダム効果を仮定した順序ロジットモデルを推定する関数がすでに実装されているので、これを使う。ordinal パッケージ (Christensen 2019) の clmm2() や clmm() がそうである。
# 個人内平均と個人内偏差を計算
nlsy.long <- nlsy.long %>%
mutate(mself = ave(nlsy.long$self, nlsy.long$id),
mpov = ave(nlsy.long$pov, nlsy.long$id),
dself = self - mself,
dpov = pov - mpov)
# clmm2() の仕様上、従属変数を因子型に変換する必要がある
nlsy.long$anti.c <- as.factor(nlsy.long$anti)library(ordinal)
ordered.logit <- clmm2(
anti.c ~ dself + dpov + mself + mpov +
black + hispanic + childage + married + gender + momage + momwork +
time + (1 | id),
data = nlsy.long
)
summary(ordered.logit, digits = 3)## Call:
## clm2(location = anti.c ~ dself + dpov + mself + mpov + black +
## hispanic + childage + married + gender + momage + momwork +
## time + (1 | id), data = nlsy.long)
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## dself -0.064 0.019 -3.343 0.00082885
## dpov 0.116 0.168 0.694 0.48786767
## mself -0.108 0.018 -5.936 2.9161e-09
## mpov 0.696 0.127 5.461 4.7309e-08
## black 0.153 0.107 1.424 0.15436839
## hispanic -0.309 0.114 -2.699 0.00694699
## childage 0.084 0.074 1.133 0.25739533
## married -0.188 0.106 -1.774 0.07607765
## gender -0.599 0.088 -6.820 9.1281e-12
## momage -0.017 0.021 -0.823 0.41043779
## momwork 0.189 0.097 1.945 0.05174536
## time 0.083 0.053 1.572 0.11605410
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 0|1 -2.567 0.953 -2.695
## 1|2 -1.516 0.951 -1.594
## 2|3 -0.679 0.951 -0.714
## 3|4 0.185 0.952 0.194
## 4|5 1.184 0.955 1.240
## 5|6 2.334 0.965 2.419
##
## log-likelihood: -2867.459
## AIC: 5770.918
## Condition number of Hessian: 2736146.97上の出力中、mself の p値が 2.9161e-09 となっているが、これは、\(2.9161 \times 10^{-9} = 0.0000000029161\) を短縮表記したもの。Threshold coefficients が切片に対応するものである。
ordinal パッケージには clmm() という、やはりランダム効果を仮定した順序ロジットモデルを推定する関数もある。clmm() のほうが新しく開発されており、複数のランダム効果を仮定できるが、clmm2() のほうがテキストに近い推定値をえられたので、そちらの結果を示した。clmm2() と clmm() のどちらのほうが正確かは不明。標準誤差はテキストよりも clmm2() のほうが大きくなっている。
テキストで利用されていたteenpov2データが見当たらないため、表3.9 は再現不可能であった。
まず、2期の分析に必要なデータを読み込む。patents.dta を d0 と名付ける。d0 は 346 行 24 列のデータである。時定変数として、科学分野ダミーを表す science と、企業の純資産を表す logsize がある。時変変数として、logr70 から logr79 までの研究開発費と、pat70 から pat79 までの特許数がある。cusip は企業ID として使われる。ardssic はテキストでは使われていないし、欠損値を含むので、リストワイズ削除の前に分析から除外する必要がある。
library(foreign)
d0 <- read.dta("patents.dta")
dim(d0)## [1] 346 24head(d0)## cusip ardssic science logsize logr70 logr71 logr72 logr73 logr74 logr75 logr76 logr77
## 1 800 15 0 6.08360 0.99684 0.88311 0.94196 1.06678 1.02901 0.92327 1.02309 0.97240
## 2 1030 14 1 1.97492 -0.45815 -0.21637 0.08434 -0.15087 -0.68464 -1.48519 -1.19495 -0.60968
## 3 2824 4 1 5.64899 3.39054 3.40697 3.44199 3.52962 3.58542 3.67343 3.77871 3.82205
## 4 4644 13 0 0.68377 0.54340 0.48454 0.58779 0.48840 0.53714 0.43436 0.33836 0.36561
## 5 7842 16 0 2.06433 0.39511 -0.88895 -0.31471 -0.21752 -0.36157 -1.29773 -1.67471 -2.14971
## 6 10202 3 0 5.98276 2.05881 2.37646 2.42241 2.47475 2.59761 2.56063 2.56625 2.60261
## logr78 logr79 pat70 pat71 pat72 pat73 pat74 pat75 pat76 pat77 pat78 pat79
## 1 1.09500 1.07624 30 28 22 34 31 32 41 60 57 77
## 2 -0.58082 -0.60915 3 1 2 1 2 3 2 1 1 1
## 3 3.88021 3.90665 48 43 61 63 58 49 42 63 77 80
## 4 0.43860 0.42459 1 0 0 1 0 0 0 1 0 0
## 5 -1.32526 -2.83380 2 3 2 1 8 0 0 0 1 0
## 6 2.67044 2.68509 32 32 24 38 36 41 52 68 74 602期のデータでポアソンモデルを推定する場合には、glm() 関数で、family = binomial という引数を指定して推定する。従属変数は cbind(2期目のカウント, 1期目のカウント) とする。まずテキスト本文中（翻訳 p.68）の定数項だけのモデルの結果を再現する。
m0 <- glm(formula = cbind(pat79, pat75) ~ 1,
family = binomial,
data = d0)
summary(m0)##
## Call:
## glm(formula = cbind(pat79, pat75) ~ 1, family = binomial, data = d0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.7456 -1.1192 0.0000 0.8994 5.0897
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.13858 0.01298 -10.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1001.4 on 299 degrees of freedom
## Residual deviance: 1001.4 on 299 degrees of freedom
## AIC: 1953
##
## Number of Fisher Scoring iterations: 3100 * (exp(coef(m0)) - 1) # 1975-1979 の減少率 (%)## (Intercept)
## -12.9409ブートストラップ標準誤差の推定には、car パッケージの Boot() を使う。実行するたびに表示される値が変化するが、テキストとよく似た結果になる。
library(car)
m0b <- Boot(m0)
sd(m0b$t) # ブートストラップ標準誤差## [1] 0.03608946ブートストラップ標準誤差は単に m0b とタイプするだけでも表示される（出力中の BootSE がブートストラップ標準誤差である）が、なぜか original の値が上の切片の値と一致しないので（本来なら一致するはず）3、単純に標準誤差だけ出力した。
研究開発費の差分として、rd0からrd5までの変数を作成する。
d0$rd0 <- d0$logr79 - d0$logr75
d0$rd1 <- d0$logr78 - d0$logr74
d0$rd2 <- d0$logr77 - d0$logr73
d0$rd3 <- d0$logr76 - d0$logr72
d0$rd4 <- d0$logr75 - d0$logr71
d0$rd5 <- d0$logr74 - d0$logr70rd0 から rd5 までの変数を独立変数とするモデルを m1.1 と名付ける。加えて、science と logsize も含むモデルを m1.2 と名付ける。ブートストラップ標準誤差は、実行するたびに表示される値が変化するが、テキストとよく似た結果になる。
m1.1 <- glm(formula = cbind(pat79, pat75) ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5,
family = binomial,
data = d0)
m1.2 <- glm(formula = cbind(pat79, pat75) ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5 + science + logsize,
family = binomial,
data = d0)
library(texreg)
knitreg(
list(m1.1, m1.2),
single.row = TRUE,
digits = 3,
caption = "2 期の特許データについての条件付きポアソン推定値（表4.1）",
caption.above = TRUE)| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | -0.222 (0.018)*** | -0.347 (0.062)*** |
| rd0 | 0.521 (0.084)*** | 0.533 (0.085)*** |
| rd1 | -0.207 (0.113) | -0.192 (0.113) |
| rd2 | -0.118 (0.111) | -0.137 (0.112) |
| rd3 | 0.060 (0.096) | 0.062 (0.096) |
| rd4 | 0.181 (0.090)* | 0.183 (0.091)* |
| rd5 | -0.093 (0.069) | -0.100 (0.069) |
| science | 0.023 (0.028) | |
| logsize | 0.017 (0.008)* | |
| AIC | 1912.916 | 1912.510 |
| BIC | 1939.841 | 1947.128 |
| Log Likelihood | -949.458 | -947.255 |
| Deviance | 949.303 | 944.898 |
| Num. obs. | 300 | 300 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
m1.1b <- Boot(m1.1)
# 上の Model 1 のブートストラップ標準誤差
apply(m1.1b$t, 2, sd) # 単に summary(m1.1b) とタイプするだけでもブートストラップ標準誤差 (BootSE) は得られるが、やはり original の値が上の係数と一致しない。## (Intercept) rd0 rd1 rd2 rd3 rd4 rd5
## 0.05178397 0.17594388 0.20391905 0.26734192 0.25270359 0.20015806 0.15275587m1.2b <- Boot(m1.2)
apply(m1.2b$t, 2, sd)## (Intercept) rd0 rd1 rd2 rd3 rd4 rd5 science logsize
## 0.13784630 0.18182436 0.22451408 0.27621738 0.25324927 0.20514817 0.15144522 0.06950166 0.01759517patents.dta を d1 と名付ける。d1 をロング形式に変換したものを d1r と名付ける。d1r は 3,460行 6列のデータである。時変変数として、研究開発費を表す rd0 と、特許数を表す pat がある。time は時点を表す変数である。
d1 <- read.dta("patents.dta")
d1r <- reshape(
d1,
drop = "ardssic", # 削除する変数名
direction = "long",
idvar = "cusip",
varying = list(5:14, 15:24),
v.names = c("rd0", "pat")
)
dim(d1r)## [1] 3460 6head(d1r)## cusip science logsize time rd0 pat
## 800.1 800 0 6.08360 1 0.99684 30
## 1030.1 1030 1 1.97492 1 -0.45815 3
## 2824.1 2824 1 5.64899 1 3.39054 48
## 4644.1 4644 0 0.68377 1 0.54340 1
## 7842.1 7842 0 2.06433 1 0.39511 2
## 10202.1 10202 0 5.98276 1 2.05881 32最初に、d1r を pdata.frame 形式に変換したものを d1p と名付ける。d1p は 3,460行 6列のデータで、index に指定された cusip と time でソートされている。次に、研究開発費のラグ変数として、rd1 から rd5 までの変数を作成する。最後に、d1p にリストワイズ削除を適用したものを d2 と名付ける。d2 は 1,730行 11列のデータである。
library(plm)
d1p <- pdata.frame(d1r, index = c("cusip", "time"))
d1p$rd1 <- lag(d1p$rd0, k = 1)
d1p$rd2 <- lag(d1p$rd0, k = 2)
d1p$rd3 <- lag(d1p$rd0, k = 3)
d1p$rd4 <- lag(d1p$rd0, k = 4)
d1p$rd5 <- lag(d1p$rd0, k = 5)
d2 <- na.omit(d1p)
dim(d2)## [1] 1730 11kable(
head(d2 [, c(1, 4, 6:5, 7:11)], n = 20L),
caption = "再構築されたデータセットの最初の4 つの企業についての観測値（表4.2）"
)| cusip | time | pat | rd0 | rd1 | rd2 | rd3 | rd4 | rd5 |
|---|---|---|---|---|---|---|---|---|
| 800 | 6 | 32 | 0.92327 | 1.02901 | 1.06678 | 0.94196 | 0.88311 | 0.99684 |
| 800 | 7 | 41 | 1.02309 | 0.92327 | 1.02901 | 1.06678 | 0.94196 | 0.88311 |
| 800 | 8 | 60 | 0.97240 | 1.02309 | 0.92327 | 1.02901 | 1.06678 | 0.94196 |
| 800 | 9 | 57 | 1.09500 | 0.97240 | 1.02309 | 0.92327 | 1.02901 | 1.06678 |
| 800 | 10 | 77 | 1.07624 | 1.09500 | 0.97240 | 1.02309 | 0.92327 | 1.02901 |
| 1030 | 6 | 3 | -1.48519 | -0.68464 | -0.15087 | 0.08434 | -0.21637 | -0.45815 |
| 1030 | 7 | 2 | -1.19495 | -1.48519 | -0.68464 | -0.15087 | 0.08434 | -0.21637 |
| 1030 | 8 | 1 | -0.60968 | -1.19495 | -1.48519 | -0.68464 | -0.15087 | 0.08434 |
| 1030 | 9 | 1 | -0.58082 | -0.60968 | -1.19495 | -1.48519 | -0.68464 | -0.15087 |
| 1030 | 10 | 1 | -0.60915 | -0.58082 | -0.60968 | -1.19495 | -1.48519 | -0.68464 |
| 2824 | 6 | 49 | 3.67343 | 3.58542 | 3.52962 | 3.44199 | 3.40697 | 3.39054 |
| 2824 | 7 | 42 | 3.77871 | 3.67343 | 3.58542 | 3.52962 | 3.44199 | 3.40697 |
| 2824 | 8 | 63 | 3.82205 | 3.77871 | 3.67343 | 3.58542 | 3.52962 | 3.44199 |
| 2824 | 9 | 77 | 3.88021 | 3.82205 | 3.77871 | 3.67343 | 3.58542 | 3.52962 |
| 2824 | 10 | 80 | 3.90665 | 3.88021 | 3.82205 | 3.77871 | 3.67343 | 3.58542 |
| 4644 | 6 | 0 | 0.43436 | 0.53714 | 0.48840 | 0.58779 | 0.48454 | 0.54340 |
| 4644 | 7 | 0 | 0.33836 | 0.43436 | 0.53714 | 0.48840 | 0.58779 | 0.48454 |
| 4644 | 8 | 1 | 0.36561 | 0.33836 | 0.43436 | 0.53714 | 0.48840 | 0.58779 |
| 4644 | 9 | 0 | 0.43860 | 0.36561 | 0.33836 | 0.43436 | 0.53714 | 0.48840 |
| 4644 | 10 | 0 | 0.42459 | 0.43860 | 0.36561 | 0.33836 | 0.43436 | 0.53714 |
3期以上のポアソンモデルを推定するには、glm() で family = poisson と指定して推定する。表4.3 の固定効果モデルの結果が再現される。ブートストラップ標準誤差は、計算に時間がかかるが、テキストとよく似た結果になる。
m3.1 <- glm(formula = pat ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5 + time + cusip,
family = poisson,
data = d2)
# summary(m3.1) で詳しい結果は見られる
# m3.1b <- Boot(m3.1) # ブートストラップ標準誤差を得る
# summary(m3.1b) 結果は割愛ランダム効果モデルを推定するには、familyをpoissonとするglmer関数を使う。デフォルトの推定法では計算が収束しなかったので、control = glmerControl(optimizer = “bobyqa”) という引数を指定して、別の推定法で推定してある。
library(lme4)
m3.2 <- glmer(
formula = pat ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5 + time + (1 | cusip),
family = poisson,
data = d2,
control = glmerControl(optimizer = "bobyqa")
) # summary(m3.2) で詳しい結果は見られるGEEモデルを推定するには、gee() で family = poisson という引数を指定する。表4.3 の GEEモデルの結果が再現される。gee() 関数では、ロバスト標準誤差がデフォルトで表示される。
library(gee)
d2$ time <- factor(d2$ time)
m3.3 <- gee(
formula = pat ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5 + time,
id = cusip,
data = d2,
family = poisson,
corstr = "unstructured")## (Intercept) rd0 rd1 rd2 rd3 rd4 rd5 time7 time8
## 1.82830787 0.15251014 0.02198137 0.04368987 0.08271559 0.10397102 0.30109113 -0.04399151 -0.06043235
## time9 time10
## -0.18924632 -0.22978995knitreg(
list(m3.1, m3.2, m3.3) ,
custom.model.names = c("固定効果", "ランダム効果", "GEE"),
single.row = TRUE,
digits = 3,
caption = "5 期の特許データについてのポアソン推定値（表4.3）",
caption.above = TRUE,
include.dev = FALSE,
include.aic = FALSE,
omit.coef = "(Intercept)|(cusip)"
)| 固定効果 | ランダム効果 | GEE | |
|---|---|---|---|
| rd0 | 0.322 (0.046)*** | 0.488 (0.042)*** | 0.303 (0.053)*** |
| rd1 | -0.087 (0.049) | -0.003 (0.048) | 0.049 (0.056) |
| rd2 | 0.079 (0.045) | 0.140 (0.045)** | 0.167 (0.051)** |
| rd3 | 0.001 (0.041) | 0.063 (0.041) | 0.085 (0.062) |
| rd4 | -0.005 (0.038) | 0.028 (0.037) | 0.050 (0.042) |
| rd5 | 0.003 (0.032) | 0.086 (0.031)** | 0.038 (0.043) |
| time7 | -0.043 (0.013)** | -0.047 (0.013)*** | -0.048 (0.017)** |
| time8 | -0.040 (0.013)** | -0.057 (0.013)*** | -0.052 (0.026) |
| time9 | -0.157 (0.014)*** | -0.192 (0.014)*** | -0.178 (0.043)*** |
| time10 | -0.198 (0.015)*** | -0.256 (0.014)*** | -0.234 (0.041)*** |
| BIC | 11502.717 | 10634.672 | |
| Log Likelihood | -4424.213 | -5272.601 | |
| Num. obs. | 1730 | 1730 | 1730 |
| Num. groups: cusip | 346 | ||
| Var: cusip (Intercept) | 0.969 | ||
| Scale | 22.469 | ||
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
研究開発費のラグ変数の代わりに rd0 と science の交互作用を含むモデルを m4 と名付ける。表4.4 のモデルの結果が再現される。ブートストラップ標準誤差は、計算に時間がかかるが、テキストとよく似た結果になる。
m4 <- glm(formula = pat ~ rd0 + rd0 : science + time + cusip,
family = poisson,
data = d2)
knitreg(
m4,
single.row = TRUE,
digits = 3,
caption = "時定共変量を含む条件付きポアソン推定値（表4.4）",
caption.above = TRUE,
include.aic = FALSE, # 以下4行は、要らない統計量の出力を省略
include.bic = FALSE,
include.loglik = FALSE,
include.dev = FALSE,
omit.coef = "(Intercept)|(cusip)" # 固定効果の推定量の出力を省略
)| Model 1 | |
|---|---|
| rd0 | 0.375 (0.048)*** |
| time7 | -0.034 (0.013)** |
| time8 | -0.034 (0.013)** |
| time9 | -0.151 (0.014)*** |
| time10 | -0.189 (0.015)*** |
| rd0:science | -0.204 (0.067)** |
| Num. obs. | 1730 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |
# m4b <- Boot(m4)
# summary(m4b)連続変数としての time と science の交互作用を含むモデルを m5 と名付ける。表4.5 のモデルの結果が再現される。ブートストラップ標準誤差は、計算に時間がかかるが、テキストとよく似た結果になる。
m5 <- glm(
formula = pat ~ rd0 + as.numeric(time) + as.numeric(time) : science + cusip,
family = poisson,
data = d2)
knitreg(
m5,
single.row = TRUE,
digits = 3,
caption = "時点との交互作用を含む条件付きポアソン推定値（表4.5）",
caption.above = TRUE,
include.aic = FALSE, # 以下4行は、要らない統計量の出力を省略
include.bic = FALSE,
include.loglik = FALSE,
include.dev = FALSE,
omit.coef = "(Intercept)|(cusip)")| Model 1 | |
|---|---|
| rd0 | 0.276 (0.039)*** |
| as.numeric(time) | -0.049 (0.005)*** |
| as.numeric(time):science | -0.001 (0.006) |
| Num. obs. | 1730 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |
# m5b <- Boot(m5)
# summary(m5b) # 結果は割愛固定効果負の二項モデルを推定するには、MASS パッケージ (Venables and Ripley 2002) の glm.nb() 関数を使う。表4.3 の固定効果ポアソンモデルと同じ変数を含むモデルを m6 と名付ける。加えて、勾配の外積（OPG）標準誤差を用いた検定の結果を m6opg と名付ける。表4.6 のモデルの結果が再現される。lmtest パッケージ (Zeileis and Hothorn 2002) の coeftest() 関数では、sandwich パッケージ (Zeileis, Köll, and Graham 2020) を読み込むと、vcov = vcovOPG という引数を指定することができる。
library(MASS)
m6 <- glm.nb(formula = pat ~ rd0 + rd1 + rd2 + rd3 + rd4 + rd5 + time + cusip, data = d2)
library(lmtest)
library(sandwich)
m6opg <- coeftest(m6, vcov = vcovOPG)
knitreg(
list(m6, m6opg),
single.row = TRUE,
digits = 3,
caption.above = TRUE,
caption = "固定効果負の二項モデル（表4.6）",
include.aic = FALSE, # 以下4行は、要らない統計量の出力を省略
include.bic = FALSE,
include.loglik = FALSE,
include.dev = FALSE,
omit.coef = "cusip")| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 3.677 (0.118)*** | 3.677 (0.101)*** |
| rd0 | 0.371 (0.063)*** | 0.371 (0.072)*** |
| rd1 | -0.083 (0.068) | -0.083 (0.073) |
| rd2 | 0.064 (0.064) | 0.064 (0.075) |
| rd3 | 0.014 (0.059) | 0.014 (0.071) |
| rd4 | 0.034 (0.056) | 0.034 (0.060) |
| rd5 | 0.002 (0.046) | 0.002 (0.052) |
| time7 | -0.049 (0.023)* | -0.049 (0.027) |
| time8 | -0.051 (0.023)* | -0.051 (0.029) |
| time9 | -0.159 (0.024)*** | -0.159 (0.027)*** |
| time10 | -0.224 (0.025)*** | -0.224 (0.028)*** |
| Num. obs. | 1730 | |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
表4.7 は条件付き負の二項推定の結果を示しているが、R には対応する関数がない。しかし、テキストでは条件付き負の二項推定は真の固定効果推定でないことが示されており、表4.7 のモデルの結果を再現してもあまり意味はない。
研究開発費の平均変数として、mrd0 から mrd5 までの変数を作成する。加えて、研究開発費の偏差変数として、drd0 から drd5 までの変数を作成する。
# 企業内平均の計算
d2$mrd0 <- ave(d2$rd0, d2$cusip)
d2$mrd1 <- ave(d2$rd1, d2$cusip)
d2$mrd2 <- ave(d2$rd2, d2$cusip)
d2$mrd3 <- ave(d2$rd3, d2$cusip)
d2$mrd4 <- ave(d2$rd4, d2$cusip)
d2$mrd5 <- ave(d2$rd5, d2$cusip)
# 企業内偏差の計算
d2$drd0 <- d2$rd0 - d2$mrd0
d2$drd1 <- d2$rd1 - d2$mrd1
d2$drd2 <- d2$rd2 - d2$mrd2
d2$drd3 <- d2$rd3 - d2$mrd3
d2$drd4 <- d2$rd4 - d2$mrd4
d2$drd5 <- d2$rd5 - d2$mrd5ランダム効果モデルを推定するには、lme4 パッケージの glmer.nb() が使えるが、計算が収束しなかったので、pglm パッケージの pglm() を使って推定した。
library(pglm)
m8.1 <- pglm(
pat ~ drd0 + drd1 + drd2 + drd3 + drd4 + drd5 +
mrd0 + mrd1 + mrd2 + mrd3 + mrd4 + mrd5 +
science + logsize + time,
family = negbin,
model = "random",
data = d2
) # 結果はあとでまとめて表示。summary(m8.1) で詳しい結果は見られる前述のとおり、gee関数を使うことによって、GEEモデルを推定することができる。しかし、gee関数では、負の二項モデルが利用可能でない。結果として、表4.8のGEEモデルの結果は再現されない。他方で、ポアソンモデルは利用可能である。負の二項モデルの代理としてのポアソンモデルをm8.2と名付ける。ロバスト標準誤差がデフォルトで表示されており、偏差変数の係数は、テキストとよく似た結果になる。
library(gee)
m8.2 <- gee(
pat ~ drd0 + drd1 + drd2 + drd3 + drd4 + drd5 + mrd0 + mrd1 + mrd2 + mrd3 + mrd4 + mrd5 +
science + logsize + factor(time),
id = cusip,
data = d2,
family = poisson,
corstr = "exchangeable"
)## (Intercept) drd0 drd1 drd2 drd3 drd4 drd5
## 0.794809948 0.252326326 -0.066331272 0.077330616 -0.013367326 0.001166999 0.003558630
## mrd0 mrd1 mrd2 mrd3 mrd4 mrd5 science
## -0.137634734 2.241298048 -3.742662355 1.627737243 0.102901539 0.407035311 0.464056447
## logsize factor(time)7 factor(time)8 factor(time)9 factor(time)10
## 0.245856609 -0.041090800 -0.035646767 -0.151105229 -0.187039425knitreg(
list(m8.1, m8.2),
custom.model.names = c("ランダム効果", "GEE"),
single.row = TRUE,
caption = "負の二項モデルのハイブリッド推定値（表4.8）",
caption.above = TRUE,
custom.coef.names = c(rep(NA, 21), "time7", "time8", "time9", "time10"),
omit.coef= "a|b", # pglm() の係数に a, b という値があり、おそらくスケールパラメータか
digits = 3 # ランダム効果の分散あたりと思われるが、不明なので、出力から削除
)| ランダム効果 | GEE | |
|---|---|---|
| (Intercept) | 1.104 (0.179)*** | 0.870 (0.247)*** |
| drd0 | 0.320 (0.071)*** | 0.281 (0.062)*** |
| drd1 | -0.054 (0.075) | -0.075 (0.059) |
| drd2 | 0.083 (0.068) | 0.078 (0.051) |
| drd3 | -0.015 (0.064) | -0.008 (0.066) |
| drd4 | 0.011 (0.058) | -0.000 (0.058) |
| drd5 | 0.025 (0.050) | 0.003 (0.063) |
| mrd0 | -0.059 (0.731) | -0.463 (1.051) |
| mrd1 | 1.286 (1.492) | 2.969 (2.438) |
| mrd2 | -1.281 (1.665) | -4.146 (3.470) |
| mrd3 | -0.220 (1.458) | 1.388 (3.007) |
| mrd4 | 0.740 (1.101) | 0.534 (2.286) |
| mrd5 | 0.228 (0.521) | 0.235 (0.998) |
| science | 0.044 (0.109) | 0.441 (0.173)* |
| logsize | 0.077 (0.047) | 0.227 (0.061)*** |
| time7 | -0.042 (0.021)* | |
| time8 | -0.049 (0.022)* | |
| time9 | -0.170 (0.023)*** | |
| time10 | -0.207 (0.024)*** | |
| DF | 21 | |
| Log Likelihood | -4941.246 | |
| AIC | 9924.491 | |
| nobs | 1730 | |
| Scale | 18.999 | |
| Num. obs. | 1730 | |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
library(foreign)
nsfg <- read.dta("nsfg.dta")
head(nsfg)## caseid pregordr age lbw married passt dur birth nobreast caesar multiple college _st _d _t _t0
## 1 10000002 4 28.50 1 1 1 27 0 0 1 0 1 1 0 27 0
## 2 10000002 3 22.17 1 1 1 76 1 0 1 0 1 1 1 76 0
## 3 10000002 2 21.25 1 1 1 11 1 1 1 0 1 1 1 11 0
## 4 10000002 1 16.42 1 0 1 58 1 1 0 0 1 1 1 58 0
## 5 10000004 1 29.08 2 0 1 30 0 1 0 0 1 1 0 30 0
## 6 10000007 4 32.92 2 1 0 41 0 0 1 0 1 1 0 41 0dim(nsfg)## [1] 14932 16nsfg <- nsfg [, 1 : 12] # 最後の4つの列は使わないし、意味がわからないので削除R で連続時間イベントヒストリー分析をする場合、ほとんど lm() とスクリプトの書き方は同じだが、以下のように従属変数の指定の仕方が異なる。
coxph(Surv(イベント発生までの時間, イベント発生ダミー) ~ 独立変数群, data = データフレーム)
イベント発生ダミーとは、イベントが発生したケースに関しては 1、発生する前に観察を打ち切られた（センサーされた）場合には 0 をとる変数である。また、クラスター・ロバスト標準誤差も計算したい場合は、
cluster = クラスター変数
という引数を追加すればよい。
まったく同じ時点にイベントが発生すると、コックス回帰分析の推定を歪ませることが知られている。このように同時にイベントが起きることをタイと呼ぶことがある。連続時間モデルなので、十分な精度でイベントの発生時間が記録されていれば、タイが生じることはほとんどない。しかし、社会科学では年単位のように粗く時間についてたずねることもよくあるため、タイは頻繁に生じる。そのような場合、コックス回帰の推定の歪みを補正する必要がある (Allison 2014)。この補正の方法として、Breslow の近似が比較的よく知られているように思うが、Efron の近似のほうが有効である（誤差が小さい）と言われている (Therneau and Grambsch 2000)。R の survival パッケージ所収の coxph() 関数のデフォルト補正法は、この Efron の近似である。ほとんどの場合デフォルトのままで良いと思うが、Breslow の近似を用いたい場合は、
ties = “breslow”
という引数を追加すればよい。
テキストの表5.1 のコックス回帰分析では Stata の stcox のデフォルトの Breslow の近似で計算されているようなので、まず Efron の近似で計算し、次に Breslow の近似でも計算した。
### 表 5.1 ###
library(survival)
# Efron 近似でコックス回帰
cox1 <- coxph(
Surv(dur, birth) ~ pregordr + age + married + passt + nobreast + lbw + caesar + multiple + college,
data = nsfg,
cluster = caseid
)
# Breslow 近似でコックス回帰
cox1.breslow <- update(cox1, ties = "breslow") # 近似法だけが上と異なるので、update() で修正
library(texreg)
knitreg(list(cox1, cox1.breslow),
digits = 3,
custom.model.names = c("Efron 近似", "Breslow 近似"),
caption.above = TRUE,
caption = "タイ・データ補正の方法によるコックス回帰の推定結果の違い") | Efron 近似 | Breslow 近似 | |
|---|---|---|
| pregordr | -0.164*** | -0.163*** |
| (0.016) | (0.016) | |
| age | -0.066*** | -0.065*** |
| (0.003) | (0.003) | |
| married | 0.223*** | 0.221*** |
| (0.030) | (0.029) | |
| passt | 0.138*** | 0.137*** |
| (0.030) | (0.030) | |
| nobreast | -0.272*** | -0.270*** |
| (0.023) | (0.023) | |
| lbw | -0.002 | -0.003 |
| (0.043) | (0.043) | |
| caesar | -0.117*** | -0.116*** |
| (0.028) | (0.028) | |
| multiple | -0.707*** | -0.702*** |
| (0.144) | (0.144) | |
| college | -0.208*** | -0.207*** |
| (0.026) | (0.026) | |
| AIC | 142893.229 | 143013.610 |
| R2 | 0.108 | 0.107 |
| Max. R2 | 1.000 | 1.000 |
| Num. events | 8021 | 8021 |
| Num. obs. | 14932 | 14932 |
| Missings | 0 | 0 |
| PH test | 0.000 | 0.000 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
Efron と Breslow のタイ・データ処理法の差は、標準誤差よりもずっと小さく、無視しうる程度のものである。ただし、同時発生イベントがもっと多くなれば差はもっと大きくなる場合もあろう。
なお、上の表の PH test とは Proportional Hazard、いわゆる平行性の仮定に関する検定の p値である。帰無仮説は、「すべての独立変数が平行性の仮定を満たしている」である。この場合は p値が 0.000 なので、独立変数のうちの一つ以上が平行性の仮定を満たしていないことになる。どの変数が問題か知りたい場合は、以下のように cox.zph() を使うと詳しい検定結果が見られる。
cox.zph(cox1)## chisq df p
## pregordr 184.3242 1 < 2e-16
## age 46.3431 1 9.9e-12
## married 4.0585 1 0.0439
## passt 0.9545 1 0.3286
## nobreast 0.1835 1 0.6684
## lbw 8.3499 1 0.0039
## caesar 0.0169 1 0.8967
## multiple 0.3947 1 0.5299
## college 15.3763 1 8.8e-05
## GLOBAL 210.9296 9 < 2e-16以下はテキストの表5.1の再現である。
### 表 5.1 の再現 ###
tab5.1 <- summary(cox1.breslow) $ coefficients [, c("coef", "se(coef)", "robust se", "exp(coef)")]
library(knitr)
kable(tab5.1,
digits = 3,
col.names = c("係数", "通常の標準誤差", "ロバスト標準誤差", "ハザード比"),
caption = "通常モデルのコックス回帰推定結果（表5.1） ")| 係数 | 通常の標準誤差 | ロバスト標準誤差 | ハザード比 | |
|---|---|---|---|---|
| pregordr | -0.163 | 0.011 | 0.016 | 0.849 |
| age | -0.065 | 0.003 | 0.003 | 0.937 |
| married | 0.221 | 0.029 | 0.029 | 1.247 |
| passt | 0.137 | 0.029 | 0.030 | 1.147 |
| nobreast | -0.270 | 0.023 | 0.023 | 0.763 |
| lbw | -0.003 | 0.042 | 0.043 | 0.997 |
| caesar | -0.116 | 0.031 | 0.028 | 0.890 |
| multiple | -0.702 | 0.143 | 0.144 | 0.495 |
| college | -0.207 | 0.026 | 0.026 | 0.813 |
固定効果モデルの場合は、モデル式の末尾に
- strata(クラスター変数)
という項を追加すればよい。なお “cluster = クラスター変数” という引数は指定しない。以下は表5.2の再現（ただし完全に再現するのは手間がかかるので、ハザード比だけ別表に）。
### 表5.2 ###
cox2.1 <- coxph(
Surv(dur, birth) ~ pregordr + age + married + passt + nobreast + lbw + caesar + multiple + college + strata(caseid),
data = nsfg
)
cox2.2 <- update(cox2.1, ~. + college : nobreast)
# ハザード比以外の結果をまとめる（ハザード比も含めて表5.2を完全にR上で作るのはちょっとめんどうだし、こちらのほうが一般的なまとめ方なので）
library(texreg)
knitreg(list(cox2.1, cox2.2),
digits = 3,
single.row = TRUE,
caption.above = TRUE,
include.zph = FALSE,
include.missing = FALSE,
include.maxrs = FALSE,
caption = "表5.2 固定効果コックス回帰モデル（ハザード比除く、おまけの統計量が下につく）")| Model 1 | Model 2 | |
|---|---|---|
| pregordr | -0.717 (0.034)*** | -0.717 (0.034)*** |
| age | 0.008 (0.011) | 0.008 (0.011) |
| married | 0.183 (0.070)** | 0.183 (0.070)** |
| passt | 0.076 (0.069) | 0.075 (0.069) |
| nobreast | -0.128 (0.060)* | 0.042 (0.100) |
| lbw | -0.236 (0.081)** | -0.242 (0.081)** |
| caesar | -0.078 (0.093) | -0.079 (0.093) |
| multiple | -0.607 (0.219)** | -0.589 (0.219)** |
| college | ||
| nobreast:college | -0.266 (0.125)* | |
| AIC | 9397.550 | 9395.005 |
| R2 | 0.162 | 0.162 |
| Num. events | 8021 | 8021 |
| Num. obs. | 14932 | 14932 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
# ハザード比も別表にまとめる
kable(exp(coef(cox2.1)), digits = 3)| x | |
|---|---|
| pregordr | 0.488 |
| age | 1.008 |
| married | 1.201 |
| passt | 1.079 |
| nobreast | 0.880 |
| lbw | 0.789 |
| caesar | 0.925 |
| multiple | 0.545 |
| college | NA |
表5.3 では Gompertz 分布を仮定した連続時間比例ハザード・モデル（ランダム効果付き）が推定されているが、これを再現する R のパッケージを見つけることができなかった。そこで、ランダム効果コックス回帰モデルを推定するスクリプトと、推定結果を以下に示す。そもそも Allison はコックス回帰の shared frailty model を推定したかったが、計算が収束しなかったので、Gompertz 分布を仮定して shared frailty model を推定したわけだが、coxph() 関数で簡単に shared frailty model は推定できる。単にモデル式の最後に、
- frailty(クラスター変数)
という項を追加すればよい。計算例とその結果を示す前にもう少しだけ関連するパッケージについて概観しておく。
Gompertz 分布を仮定した連続時間比例ハザード・モデルは parfm パッケージ (Munda, Rotolo, and Legrand 2012) に実装されていることになっているが、計算結果がテキストの表5.3とはまったく異なっていた。また eha パッケージ (Brostrom 2020) の phreg() 関数は Gompertz 分布を仮定した連続時間比例ハザード・モデルを推定できるが、ランダム効果 (shared frailty) についてはマニュアルに記載がない。coxph() と同じようにモデル式に + frailty(クラスター変数) を加えてみたが、ただの共変量として扱われてしまっているようで、ランダム効果は推定されなかった。
# 時間はかかったが、コックス回帰でランダム効果モデル (shared frailty model と呼ばれる) も推定できる
cox3.1 <- update(cox1.breslow, ~ . + frailty(caseid), cluster = NULL)
cox3.1## Call:
## coxph(formula = Surv(dur, birth) ~ pregordr + age + married +
## passt + nobreast + lbw + caesar + multiple + college + frailty(caseid),
## data = nsfg, ties = "breslow")
##
## coef se(coef) se2 Chisq DF p
## pregordr -0.16318 0.01150 0.01149 201.50834 1.0 < 2e-16
## age -0.06531 0.00306 0.00306 456.41566 1.0 < 2e-16
## married 0.22103 0.02867 0.02867 59.42495 1.0 1.3e-14
## passt 0.13726 0.02868 0.02868 22.89891 1.0 1.7e-06
## nobreast -0.27000 0.02332 0.02332 133.99992 1.0 < 2e-16
## lbw -0.00258 0.04204 0.04204 0.00377 1.0 0.95101
## caesar -0.11636 0.03054 0.03054 14.51660 1.0 0.00014
## multiple -0.70242 0.14258 0.14258 24.27012 1.0 8.4e-07
## college -0.20666 0.02599 0.02598 63.25062 1.0 1.8e-15
## frailty(caseid) 0.39326 0.4 0.29272
##
## Iterations: 5 outer, 23 Newton-Raphson
## Variance of random effect= 5e-05 I-likelihood = -71497.8
## Degrees of freedom for terms= 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.4
## Likelihood ratio test=1683 on 9.4 df, p=<2e-16
## n= 14932, number of events= 8021# 上で十分だが、表5.3 と同じ形に分析結果を整形
se.cox3.1 <- sqrt(diag(vcov(cox3.1))) # 標準誤差
tab5.3 <- cbind(
"係数" = coef(cox3.1),
"標準誤差" = se.cox3.1,
"ハザード比" = exp(coef(cox3.1))
)
# 表5.3をゴンペルツ・モデルではなく、コックス回帰で再現（推定値は若干異なる）
kable(tab5.3, digits = 3, caption = "表5.3改 ランダム効果コックス回帰モデルの推定結果")| 係数 | 標準誤差 | ハザード比 | |
|---|---|---|---|
| pregordr | -0.163 | 0.011 | 0.849 |
| age | -0.065 | 0.003 | 0.937 |
| married | 0.221 | 0.029 | 1.247 |
| passt | 0.137 | 0.029 | 1.147 |
| nobreast | -0.270 | 0.023 | 0.763 |
| lbw | -0.003 | 0.042 | 0.997 |
| caesar | -0.116 | 0.031 | 0.890 |
| multiple | -0.702 | 0.143 | 0.495 |
| college | -0.207 | 0.026 | 0.813 |
表5.4 と表5.5 は用いられているデータを入手できなかったので、再現不可能。以下では Allison (2014) で用いられていたサンプルデータを用いて、繰り返しのないイベントの場合の固定効果イベントヒストリーモデルの推定を R でやってみた結果を示す。生化学の助教がいつ准教授に昇進したかを調べたデータで、助教がリスク・セットの構成員で、イベントは准教授への昇進、以下では所属研究機関の威信と累積論文出版数という時変独立変数の効果を検討する。データの詳細は、Long, Allison, and McGinnis (1979) を見よ。
離散時間モデルが用いられるが、これについては、Allison (2014) を参照。離散時間モデルの場合、データを準備するのに手間がかかる。Rでの準備例として 太郎丸 (2018) を参照。計算そのものは、survival パッケージの clogit() 関数を使えばよい。これについては 3章のスクリプトとして紹介されているので、そちらを参照。
library(foreign) # Stata のデータを読み込むためのパッケージ
rank.dta <- read.dta("rank.dta") # read.dta() で dta 形式のファイルを読み込める
head(rank.dta) # dur は観察が何年目に打ち切られたかを示す## dur promo undgrad phdmed phdprest art1 art2 art3 art4 art5 art6 art7 art8 art9 art10 cit1 cit2 cit3 cit4
## 1 10 0 7.00 0 2.21 0 0 2 2 2 2 2 2 2 2 0 0 1 1
## 2 4 1 6.00 0 2.21 8 10 14 18 NA NA NA NA NA NA 27 44 57 63
## 3 4 0 4.95 0 2.21 0 0 0 2 NA NA NA NA NA NA 0 0 0 2
## 4 10 0 2.00 1 4.54 2 3 3 3 3 4 6 6 6 6 11 17 17 17
## 5 7 1 5.00 1 2.15 1 1 1 2 2 3 5 NA NA NA 14 14 14 19
## 6 6 1 4.95 0 4.54 0 2 3 5 5 6 NA NA NA NA 0 0 2 6
## cit5 cit6 cit7 cit8 cit9 cit10 prest1 prest2 jobtime
## 1 1 1 1 1 1 1 2.36 2.36 NA
## 2 NA NA NA NA NA NA 1.84 2.88 3
## 3 NA NA NA NA NA NA 2.40 2.40 NA
## 4 17 20 24 24 24 24 2.56 2.56 NA
## 5 19 20 23 NA NA NA 2.40 2.40 NA
## 6 6 11 NA NA NA NA 3.36 3.36 NA# jobtime は勤め先の異動があった場合、何年目から2つ目の職場に異動したかを示し、NA は異動なし
# promo が 1 は昇進した人、0 はしなかった人
rank.dta$jobtime [is.na(rank.dta$jobtime)] <- 0 # NA を 0 に置き換える。あとでデータの整形をするさいに役立つ
# データをロングフォーマットに変換
rank.long <- reshape(
rank.dta,
direction = "long",
varying = list(
paste("art", 1:10, sep = ""),
paste("cit", 1:10, sep = "")
),
v.names = c("art", "cit") # 時変変数名の定義
)
rank.long <- rank.long [order(rank.long$id), ] # id 順に並べ替える（こっちのほうが分かりやすい）
head(rank.long) # time という新しい変数ができ、何年目のパーソンイヤーかを示す## dur promo undgrad phdmed phdprest prest1 prest2 jobtime time art cit id
## 1.1 10 0 7 0 2.21 2.36 2.36 0 1 0 0 1
## 1.2 10 0 7 0 2.21 2.36 2.36 0 2 0 0 1
## 1.3 10 0 7 0 2.21 2.36 2.36 0 3 2 1 1
## 1.4 10 0 7 0 2.21 2.36 2.36 0 4 2 1 1
## 1.5 10 0 7 0 2.21 2.36 2.36 0 5 2 1 1
## 1.6 10 0 7 0 2.21 2.36 2.36 0 6 2 1 1nrow(rank.long)## [1] 3010離散時間イベントヒストリーモデルの推定のために、下記のようにデータを整形する。
rank.long <- subset(rank.long, dur >= time) # 観察が打ち切られた後のパーソンイヤーは削除
rank.long $ promotion <- rank.long $ promo # promo のバックアップ。あとで使う
rank.long$promo [rank.long$time < rank.long$dur] <- 0 # イベントが起きる前は従属変数をゼロに
rank.long$jobpres <- rank.long$prest1 # 勤め先の威信を示す jobpres という時変変数を作る
rank.long$jobpres [rank.long$time >= rank.long$jobtime] <- # 勤め先が変わった後は
rank.long$prest2 [rank.long$time >= rank.long$jobtime] # 2番目の勤め先の威信を代入まずは、通常の離散時間イベントヒストリー分析をしてみる。
# 繰り返しのない場合の固定効果モデルでは、独立変数はダミー変数でなければならないので、ダミー変数を作る。
rank.long $ job.prestige <- rank.long $ jobpres > quantile(rank.long $ jobpres, 0.2) # 勤め先威信が 20パーセンタイルより大きいかどうかのダミー変数
rank.long $ n.articles <- rank.long $ art > median(rank.long $ art) # 論文数が中央値より大きいかどうかのダミー変数
# 時間以外の主効果のみのモデル
logit0 <- glm(promo ~ undgrad + phdmed + phdprest + job.prestige + n.articles + cit + time + I(time ^ 2),
data = rank.long,
family = binomial)
summary(logit0)##
## Call:
## glm(formula = promo ~ undgrad + phdmed + phdprest + job.prestige +
## n.articles + cit + time + I(time^2), family = binomial, data = rank.long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6069 -0.5466 -0.2388 -0.0732 3.5317
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.7462923 0.7713668 -11.339 < 2e-16 ***
## undgrad 0.1546917 0.0616757 2.508 0.01214 *
## phdmed -0.2184992 0.1717745 -1.272 0.20337
## phdprest 0.0144057 0.0922596 0.156 0.87592
## job.prestigeTRUE -0.5173064 0.2001169 -2.585 0.00974 **
## n.articlesTRUE 0.8795313 0.2024845 4.344 1.40e-05 ***
## cit 0.0022643 0.0008833 2.564 0.01036 *
## time 2.0948823 0.2340815 8.949 < 2e-16 ***
## I(time^2) -0.1567487 0.0202249 -7.750 9.17e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1309.5 on 1740 degrees of freedom
## Residual deviance: 1006.4 on 1732 degrees of freedom
## AIC: 1024.4
##
## Number of Fisher Scoring iterations: 7勤め先の威信 job.prestige も論文の数 n.articles も有意である。次が固定効果モデルの推定結果である。
d0 <- subset(rank.long, subset = promotion == 1, select = c(id, time, promo, n.articles, job.prestige)) # 昇進した人だけに限定
library(survival)
# 勤め先威信の係数の推定
clogit1 <- clogit(job.prestige ~ promo + time + I(time ^ 2) + n.articles + strata(id), data = d0)
# 論文数の係数の推定
clogit2 <- clogit(n.articles ~ promo + time + I(time ^ 2) + job.prestige + strata(id),
data = d0) # update() を使うとエラーになった
knitreg(
list(clogit1, clogit2),
custom.model.names = c("勤め先威信", "論文数"),
caption.above = TRUE,
caption = "テキストとは別のデータで繰り返しのないイベントの固定効果モデルの推定をした結果")| 勤め先威信 | 論文数 | |
|---|---|---|
| promo | -1.92* | 1.35 |
| (0.84) | (3999.55) | |
| time | 0.91 | 19.98 |
| (0.49) | (3046.95) | |
| time^2 | -0.08 | -0.20 |
| (0.05) | (389.79) | |
| n.articlesTRUE | 0.96 | |
| (1.01) | ||
| job.prestigeTRUE | -16.43 | |
| (20049.29) | ||
| AIC | 61.82 | 8.00 |
| R2 | 0.01 | 0.33 |
| Max. R2 | 0.05 | 0.33 |
| Num. events | 947 | 679 |
| Num. obs. | 1203 | 1203 |
| Missings | 0 | 0 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||
論文数の係数の推定は、Allisonが指摘するように、計算が収束しなかった。この変数は、その時点までに出版した論文数なので、必ず時間とともに上昇するからであろう。いちおう参考のためにこの結果も示してある。上の表の promo の行が、それぞれ勤め先威信と論文数の昇進に及ぼす効果である。この行以外は、実質的には無意味な数値ばかりだが、やはり参考のために表記した。勤め先威信は通常のイベントヒストリー分析と同様に有意であった。テキストと同じように、1, 2, 3 年以内に論文数が全体の中央値を超えたかどうかを示すダミー変数も作ってそれらの効果を推定してみたが、計算は収束するもののやはり有意な結果はえられなかった。
データを読み込む。
library(foreign)
nlsy <- read.dta("nlsy.dta")
head(nlsy)## momage anti90 anti92 anti94 gender childage hispanic black momwork married self90 self92 self94 pov90 pov92
## 1 21 1 1 1 1 8.000000 0 0 0 0 21 24 23 1 1
## 2 22 0 0 0 1 8.416667 0 0 1 0 20 24 24 0 0
## 3 18 5 5 5 0 8.083333 0 0 0 1 21 24 24 0 0
## 4 24 2 3 1 0 8.250000 0 0 0 0 23 21 21 0 0
## 5 22 1 0 0 1 9.333333 0 0 0 0 22 23 24 0 0
## 6 24 1 1 1 0 8.583333 0 0 0 0 19 21 24 0 0
## pov94
## 1 1
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0dim(nlsy)## [1] 581 16共分散構造分析あるいは構造方程式モデリング（structural equation modeling, SEM）は sem (Fox, Nie, and Byrnes 2021) などいくつかのパッケージに実装されているが、本稿では lavaan パッケージ (Rosseel 2012) を用いる。lavaan パッケージ で SEM を実行する場合、sem() という関数もあるものの、以下では、
lavaan(model=モデル,data=データフレーム)
と指定している。lavaan() の引数にあるモデルは次のように ’ ’ （バッククォート）の間に指定する。
model <- ‘ ’
また、SEM でのモデルの指定は、大きく、構造方程式の部分と測定方程式の部分に分けられる。構造方程式の部分は、一般的な回帰モデルと同じように指定する。たとえば、被説明変数が y で、説明変数が x1、x2、x3 の回帰モデルは、
y ~ x1 + x2 + x3
と指定する。
SEM では、複数の回帰式や、ある回帰式の説明変数が別の回帰式の被説明変数になっているモデルを指定することができる。たとえば、
y1 ~ x1 + x2 + x3
y2 ~ x4 + x5
y3 ~ y1 + y2
というような場合である。
測定方程式は =~ を用いて記述する。観測変数 x1, x2, x3 に、潜在因子 f1 を仮定した場合、
f1 =~ x1 + x2 + x3
と指定する。固定効果モデルやランダム効果モデルにおける時間的に変化しない観察されない異質性も、SEM ではこの潜在変数として位置づけられる。
また、観測変数や誤差の分散や変数間の共分散（相関）は、~~ を用いて記述する。たとえば、
x1 ~~ x3
と指定する。
以下では、
auto.var = TRUE
という引数を使用している。これは外生潜在共変量（潜在変数だが、構造方程式で独立変数と指定されているが、従属変数とはなっていない変数）の分散を推定させるための引数である。この場合だと、切片＝時不変の観察されない異質性（以下のスクリプト中では f1 と名付けられている）の分散を推定させるためのものである。
summary() の出力には多くの情報が含まれているが、表6.1の内容はregressionの部分を確認する。
library(lavaan)
m_random <- '
f1 =~ 1 * anti90 + 1 * anti92 + 1 * anti94 # 測定方程式（負荷量を1で固定）。これがランダム切片（時不変の観測されない異質性）に対応
anti90 ~ a * self90 + b * pov90 + c * black + d * hispanic + e * childage + # ここから構造方程式
f * married + g * gender + h * momage + i * momwork # 3時点ですべての変数に等値制約
anti92 ~ a * self92 + b * pov92 + c * black + d * hispanic + e * childage +
f * married + g * gender + h * momage + i * momwork
anti94 ~ a * self94 + b * pov94 + c * black + d * hispanic + e * childage +
f * married + g * gender + h * momage + i * momwork
anti90 ~~ j * anti90 # 誤差分散を一定にする（j）
anti92 ~~ j * anti92
anti94 ~~ j * anti94
'
fit1 <- lavaan(
model = m_random,
data = nlsy,
auto.var = TRUE
)
# summary(fit1) 固定効果モデルについては、固定効果と時変変数のあいだに共変関係を仮定する。さらに、外生測定変数間の相関をすべて明示的に指定する必要がある。そのため、モデルの指定は以下に示すようにとても長くなる。表6.1の内容はregressionの部分で確認する。
m_fixed <- paste( # lavaan() のモデルは、lm() や glm() と違い、formula クラス・オブジェクトで指定するのではなく、
m_random, # ' ' でくくった単なる文字列で指定する。そのため update() でモデルを修正することはできない
' # しかし、単にモデルの記述を書き加えるだけならば、paste() で追加の指定が可能
#### 固定効果f1と時変変数とのあいだに相関を仮定する ####
f1 ~~ pov90 + self90 + pov92 + self92 + pov94 + self94
# 時変変数と時不変変数を含むすべての外生変数間に共変関係があることを指定する
pov90 ~~ self90 + black + hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
self90 ~~ black + hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
black ~~ hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
hispanic ~~ childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
childage ~~ married + gender + momage + momwork + pov92 + self92 + pov94 + self94
married ~~ gender + momage + momwork + pov92 + self92 + pov94 + self94
gender ~~ momage + momwork + pov92 + self92 + pov94 + self94
momage ~~ momwork + pov92 + self92 + pov94 + self94
momwork ~~ pov92 + self92 + pov94 + self94
pov92 ~~ self92 + pov94 + self94
self92 ~~ pov94 + self94
pov94 ~~ self94
'
)
fit2 <- update(fit1, m_fixed)
# summary(fit2)固定効果モデル同様、外生変数間の共変関係を指定する。
m_comp <- paste( # lavaan() のモデルは、lm() や glm() と違い、formula クラス・オブジェクトで指定するのではなく、
m_random, # ' ' でくくった単なる文字列で指定する。そのため update() でモデルを修正することはできない
' # しかし、単にモデルの記述を書き加えるだけならば、paste() で追加の指定が可能
# 固定効果f1と時変変数とのあいだに相関を仮定する(ただし、self変数はここに含めない)
f1 ~~ pov90 + pov92 + pov94
# 時変変数と時不変変数を含むすべての外生変数間に共変関係があることを指定する
pov90 ~~ self90 + black + hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
self90 ~~ black + hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
black ~~ hispanic + childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
hispanic ~~ childage + married + gender + momage + momwork + pov92 + self92 + pov94 + self94
childage ~~ married + gender + momage + momwork + pov92 + self92 + pov94 + self94
married ~~ gender + momage + momwork + pov92 + self92 + pov94 + self94
gender ~~ momage + momwork + pov92 + self92 + pov94 + self94
momage ~~ momwork + pov92 + self92 + pov94 + self94
momwork ~~ pov92 + self92 + pov94 + self94
pov92 ~~ self92 + pov94 + self94
self92 ~~ pov94 + self94
pov94 ~~ self94
'
)
fit3 <- update(fit1, m_comp)
# summary(fit3)library(texreg)
knitreg( # 表6.1 の再現
list(fit1, fit2, fit3),
custom.model.names = c("ランダム効果", "固定効果", "折衷"),
single.row = TRUE,
custom.coef.map = list( # 係数がたくさん表示されすぎるので、
"anti90 ~ self90" = "SELF", # 表6.1 に記載の係数だけ選んでラベルを指定
"anti90 ~ pov90" = "POV",
"anti90 ~ black" = "BLACK",
"anti90 ~ hispanic" = "HISPANIC",
"anti90 ~ childage" = "CHILDAGE",
"anti90 ~ married" = "MARRIED",
"anti90 ~ gender" = "GENDER",
"anti90 ~ momage" = "MOMAGE",
"anti90 ~ momwork" = "MOMWORK"
),
digits = 3,
caption.above = TRUE,
caption = "表6.1 の再現: 構造方程式モデルで NLSY データを分析した結果"
)| ランダム効果 | 固定効果 | 折衷 | |
|---|---|---|---|
| SELF | -0.062 (0.009)*** | -0.055 (0.011)*** | -0.062 (0.009)*** |
| POV | 0.247 (0.080)** | 0.112 (0.093) | 0.111 (0.093) |
| BLACK | 0.227 (0.125) | 0.269 (0.126)* | 0.269 (0.126)* |
| HISPANIC | -0.218 (0.137) | -0.198 (0.138) | -0.201 (0.138) |
| CHILDAGE | 0.088 (0.091) | 0.089 (0.091) | 0.090 (0.091) |
| MARRIED | -0.050 (0.126) | -0.022 (0.127) | -0.025 (0.127) |
| GENDER | -0.483 (0.106)*** | -0.476 (0.106)*** | -0.479 (0.106)*** |
| MOMAGE | -0.022 (0.025) | -0.026 (0.025) | -0.025 (0.025) |
| MOMWORK | 0.261 (0.114)* | 0.296 (0.115)* | 0.295 (0.115)* |
| agfi | 0.644 | 0.932 | 0.936 |
| AIC | 5876.398 | 23695.634 | 23692.563 |
| BIC | 5924.411 | 24167.027 | 24150.862 |
| cfi | 0.928 | 0.978 | 0.978 |
| chisq | 84.564 | 66.565 | 69.494 |
| npar | 11.000 | 108.000 | 105.000 |
| rmsea | 0.051 | 0.049 | 0.046 |
| rmsea.conf.high | 0.064 | 0.064 | 0.061 |
| srmr | 0.029 | 0.022 | 0.024 |
| tli | 0.912 | 0.907 | 0.916 |
| converged | 1 | 1 | 1 |
| ngroups | 1 | 1 | 1 |
| nobs | 581 | 581 | 581 |
| norig | 581 | 581 | 581 |
| nexcluded | 0 | 0 | 0 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | |||
agfi～tli までは構造方程式モデリングの適合度指標である。converged は計算が収束したかどうかを示す（1: 収束、0: 収束せず）、ngroups はグループ数（多母集団用のモデルだとこれが2以上になる）、norig は読み込んだデータの行数。デフォルトだとリストワイズで欠損値処理するので、欠損値があれば、nobs （計算に使用した事例の数）は、そのぶんだけ norig より少なくなる。このリストワイズによりモデルから除去された事例の数が nexcluded である。
尤度比検定をする場合、まず、ランダム効果モデルと固定効果モデルの場合、
anova(fit1, fit2)## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit2 28 23695.6 24167.0 66.565
## fit1 34 5876.4 5924.4 84.564 17.998 6 0.006236 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1と指定すればよい。Chisq diff が尤度比統計量、Df diff がその自由度、Pr(>Chisq) が p値である。
続いて、固定効果モデルと折衷モデルの比較の場合、
anova(fit2, fit3)## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit2 28 23696 24167 66.565
## fit3 31 23693 24151 69.494 2.9293 3 0.4027と指定すればよい。
standardizedSolution (lavaanクラスオブジェクト) で標準化係数がえられる。外生変数間の相関係数もこれでえられる。
std.coefs <- standardizedSolution(
fit2,
pvalue = FALSE, # p値の出力を抑止
ci = FALSE # 信頼区間の出力を抑止
)
nrow(std.coefs) # SEM ではたくさんの係数が推定されることも多いので何行あるか見てみる## [1] 131head(std.coefs) # 中身を少しだけ見てみる## lhs op rhs label est.std se z
## 1 f1 =~ anti90 0.724 0.017 43.794
## 2 f1 =~ anti92 0.722 0.016 44.384
## 3 f1 =~ anti94 0.725 0.016 44.469
## 4 anti90 ~ self90 a -0.112 0.022 -5.184
## 5 anti90 ~ pov90 b 0.034 0.028 1.207
## 6 anti90 ~ black c 0.082 0.039 2.131# 必要な行（時間的に変化しない観察されない異質性と時変変数の相関）だけとってくる
std.coefs <- std.coefs [std.coefs$ lhs == "f1" & # 左側の変数が f1 で、
std.coefs$ op == "~~" & # オペレーターが分散・共分散で（回帰係数ではない）
std.coefs$ rhs != "f1", ] # 右側の変数が f1 ではない行だけ
std.coefs## lhs op rhs label est.std se z
## 34 f1 ~~ pov90 0.151 0.045 3.372
## 35 f1 ~~ self90 -0.036 0.048 -0.757
## 36 f1 ~~ pov92 0.060 0.045 1.320
## 37 f1 ~~ self92 -0.082 0.049 -1.696
## 38 f1 ~~ pov94 0.116 0.046 2.494
## 39 f1 ~~ self94 -0.048 0.048 -0.999# 上で十分だが、表6.2 どおりに整形
kable(
std.coefs [c(2, 4, 6, 1, 3, 5),
c("rhs", "est.std", "z")],
col.names = c("", "相関係数", "z統計量"),
row.names = FALSE,
caption = "アルファと時変予測変数との相関（表6.2）",
digits = 3
)| 相関係数 | z統計量 | |
|---|---|---|
| self90 | -0.036 | -0.757 |
| self92 | -0.082 | -1.696 |
| self94 | -0.048 | -0.999 |
| pov90 | 0.151 | 3.372 |
| pov92 | 0.060 | 1.320 |
| pov94 | 0.116 | 2.494 |
相関係数の値はテキストとは結構ずれている印象だが、z値の大きさは大差ない感じである。原因は不明。
データを読み込む。
occ <- read.dta("occ.dta")
head(occ)## pf1 pf2 pf3 pf4 mdwgf1 mdwgf2 mdwgf3 mdwgf4
## 1 0.3658537 0.3971743 0.4810659 0.5161744 13.96397 17.05506 16.75154 17.54808
## 2 0.4054054 0.4704830 0.4977444 0.5803981 14.59221 16.79543 17.86831 17.78846
## 3 0.1570248 0.2868217 0.3629630 0.4234234 12.26452 14.34610 20.66023 20.19231
## 4 0.1926978 0.3228963 0.3511111 0.4223512 16.08589 15.09781 16.33275 19.23077
## 5 0.3402923 0.4380165 0.5682159 0.6046852 16.31181 17.66186 16.77636 17.93785
## 6 0.6470588 0.6178344 0.7621440 0.7454545 18.17193 18.60539 15.96980 16.82692dim(occ)## [1] 178 8まず賃金中央値を1期前の賃金中央値と女性割合で予測する固定効果モデルを推定する。
m_recip1 <- ' # 測定方程式：固定効果に対応
f1 =~ 1 * mdwgf4 + 1 * mdwgf3 + 1 * mdwgf2 # 負荷量を1で固定
mdwgf4 ~ a * mdwgf3 + b * pf3 # 構造方程式
mdwgf3 ~ a * mdwgf2 + b * pf2 # 3時点ですべての変数に等値制約 a, b
mdwgf2 ~ a * mdwgf1 + b * pf1
f1 ~~ mdwgf1 + pf1 + pf2 + pf3 # 固定効果と外生変数間に相関を仮定
mdwgf1 ~~ pf1 + pf2 + pf3 # 外生変数間に相関を仮定
pf1 ~~ pf2 + pf3 # 〃
pf2 ~~ pf3 # 〃
mdwgf2 ~~ pf3 # テキスト付録の Mplus のプログラムで指定してあったので同じように指定
'
fit4 <- lavaan(m_recip1,
data = occ,
auto.var = TRUE)
# summary(fit4)次に女性割合を1期前の女性割合と賃金中央値で予測する。
m_recip2 <- ' # 測定方程式（負荷量を1で固定している）
f1 =~ 1 * pf4 + 1 * pf3 + 1 * pf2
# 構造方程式（3時点ですべての変数に等値制約 c, d）
pf4 ~ c * mdwgf3 + d * pf3
pf3 ~ c * mdwgf2 + d * pf2
pf2 ~ c * mdwgf1 + d * pf1
f1 ~~ pf1 + mdwgf1 + mdwgf2 + mdwgf3 # 固定効果と外生変数の相関を仮定
pf1 ~~ mdwgf1 + mdwgf2 + mdwgf3 # 外生変数間に共変関係を指定
mdwgf1 ~~ mdwgf2 + mdwgf3
mdwgf2 ~~ mdwgf3
pf2 ~~ mdwgf3 # これも理由はよくわからないが Mplus の指定通りに共変関係を指定
'
fit5 <- lavaan(m_recip2,
data = occ,
auto.var=TRUE)
# summary(fit5)
knitreg(
list(fit4, fit5),
single.row = TRUE,
custom.coef.map = list( # 必要な係数だけとって表6.3の通りに並べるまじない
"mdwgf4 ~ mdwgf3" = "賃金中央値",
"mdwgf4 ~ pf3" = "女性割合",
"pf4 ~ mdwgf3" = "賃金中央値",
"pf4 ~ pf3" = "女性割合"
),
custom.model.names = c("賃金中央値", "女性割合"),
caption.above = TRUE,
caption = "双方向効果モデルの推定値（表6.3）",
digits = 3
)| 賃金中央値 | 女性割合 | |
|---|---|---|
| 賃金中央値 | 0.344 (0.064)*** | -0.001 (0.002) |
| 女性割合 | -0.159 (2.436) | 0.299 (0.082)*** |
| agfi | 0.918 | 0.880 |
| AIC | 2639.983 | 1389.364 |
| BIC | 2706.801 | 1456.181 |
| cfi | 0.996 | 0.994 |
| chisq | 13.815 | 20.398 |
| npar | 21.000 | 21.000 |
| rmsea | 0.074 | 0.104 |
| rmsea.conf.high | 0.131 | 0.157 |
| srmr | 0.034 | 0.011 |
| tli | 0.989 | 0.981 |
| converged | 1 | 1 |
| ngroups | 1 | 1 |
| nobs | 178 | 178 |
| norig | 178 | 178 |
| nexcluded | 0 | 0 |
| ***p < 0.001; **p < 0.01; *p < 0.05 | ||