6 繰り返しのある事象モデル

研究で扱うほとんどの事象には繰り返しがある。例えば、一個体が繰り返し経験しうる転職・出産・離婚・逮捕などがこれに該当する。

第5章で扱った元服役囚の再犯データにも実は繰り返し逮捕された元服役囚がいた。逮捕された回数(arrstcount)ごとに人数を見てみると、以下のようになる(表6.1)。1年間の間に最大6回まで再犯を犯した個体がいる。再犯を犯した服役囚の内、2回以上逮捕された元服役囚は実に約39.5%にのぼる(132/334)。

tarp %>% 
  group_by(arrstcount) %>% 
  summarise(N = n()) %>% 
  kable(caption = "逮捕回数ごとの人数") %>% 
  kable_styling(font_size = 13, full_width = FALSE)

表6.1: 逮捕回数ごとの人数
arrstcount	N
0	598
1	202
2	88
3	28
4	10
5	4
6	2

本章ではこうした繰り返しの事象を扱う方法を解説する。なお、話を簡単にするため、繰り返された逮捕は一種類の犯罪であると仮定する。

6.1 繰り返しのある事象のカウントデータ・モデル

6.1.1 モデルの概要

もっとも簡単な方法は、事象が発生するタイミングを無視して国個体の事象発生数を従属変数にするものである。これは、以下の2つが仮定できる場合には最善の方法である。なぜならば、以下の条件を満たすのであれば。データに含まれている観察期間における事象の発生タイミングは有益な情報をほとんど持たないからである。

時間依存独立変数がない
独立変数が観測期間全てで同じ効果を持っている

カウントデータの分析で最も適切かつ容易に利用できるのは負の二項分布モデルである。ポアソン分布モデルも使えるが、しばしば過分散の問題が生じる¹³。個体\(i\)に発生した事象の数を\(Y_i\)とし、\(Y_i\)が期待値\(\lambda_i\)の負の二項分布に従うと仮定するとき、回帰モデルは以下のように表現される(式(6.1))。つまり、通常の一般化線形モデルと同様である。

\[ \begin{aligned} Y_i &\sim NegBinomial(\lambda_i)\\ log(\lambda_i) &= b_0 + b_1x_{i1} + b_1x_{i2} + \cdots + b_kx_{ik} \end{aligned} \tag{6.1} \]

6.1.2 Rでの実装

負の二項モデルはMASSパッケージのglm.nb()関数で以下のように実行できる。

mod_nb <- glm.nb(arrstcount ~  fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb, data = tarp)

結果は以下の通り(表6.2)。

model_parameters(mod_nb) %>% 
  data.frame() %>% 
  mutate(`95%CI` = str_c("[",sprintf("%.2f",CI_low),", ",sprintf("%.2f",CI_high),"}")) %>% 
  mutate(`Exp(Coef)` = exp(Coefficient)) %>% 
  dplyr::select(-df_error, -CI_low, -CI_high, -CI) %>% 
  dplyr::select(Parameter, Coefficient, `95%CI`, `Exp(Coef)`, everything()) %>% 
  mutate_if(.predicate = is.numeric,.funs = ~sprintf("%.3f",.)) %>% 
  kable(align = "c", caption = "mod1の結果") %>% 
  kable_styling(font_size = 11, full_width = FALSE)

表6.2: mod1の結果
Parameter	Coefficient	95%CI	Exp(Coef)	SE	z	p
(Intercept)	0.201	[-0.69, 1.09}	1.223	0.458	0.439	0.661
fin	0.148	[-0.06, 0.36}	1.160	0.109	1.359	0.174
age	-0.033	[-0.05, -0.02}	0.968	0.008	-4.132	0.000
white	-0.153	[-0.38, 0.07}	0.858	0.113	-1.354	0.176
male	0.372	[-0.14, 0.92}	1.450	0.270	1.379	0.168
married	-0.052	[-0.29, 0.18}	0.949	0.120	-0.435	0.664
paro	-0.323	[-0.55, -0.10}	0.724	0.116	-2.795	0.005
numprop	0.278	[0.13, 0.43}	1.321	0.073	3.827	0.000
crimprop	0.343	[0.09, 0.60}	1.409	0.129	2.658	0.008
numarst	0.013	[0.00, 0.02}	1.013	0.005	2.904	0.004
edcomb	-0.060	[-0.11, -0.01}	0.942	0.024	-2.478	0.013

負の二項分布モデルではオフセット項に観察期間を含めることで、観察期間の違う複数の個体がいるデータについても簡単に分析が可能である¹⁴。

6.2 時間のギャップに基づく方法

6.2.1 モデルの概要

時間依存変数や独立変数の影響が時間によって変化すると仮定される場合は、より複雑な方法が必要になる。この方法では、観察の開始から最初の事象までを一つの記録、また次の事象までを発生までを別の一つの記録、というように事象の発生ごとに観察期間を区切ったデータを作成する。

データはすでに作成されているので読み込む。lengthが事象間の時間のギャップを表している。

arrest <- read_dta("data/arrests.dta")

datatable(arrest,
          options = list(scrollX = 60))

分析では、これらの観察期間の記録全てを異なった個体の観測値とみなし、「時間のギャップ」(= 観測の開始時点と終了時点の差)をコックス回帰で分析する。

6.2.2 Rでの実装

第3で解説したのと同様にコックス回帰を行う。データにはなぜかlengthが負になるものが含まれているので、それらは除く。

arrest <- filter(arrest, length > 0)

mod_tg <- coxph(Surv(length,arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb, data = arrest)

分析の結果は以下の通り。

mod_tg

## Call:
## coxph(formula = Surv(length, arrind) ~ fin + age + white + male + 
##     married + paro + numprop + crimprop + numarst + edcomb, data = arrest)
## 
##               coef exp(coef)  se(coef)      z        p
## fin       0.135958  1.145633  0.090046  1.510 0.131075
## age      -0.030654  0.969811  0.006762 -4.533 5.82e-06
## white    -0.156933  0.854761  0.094189 -1.666 0.095683
## male      0.339312  1.403981  0.230061  1.475 0.140246
## married  -0.093797  0.910468  0.100009 -0.938 0.348303
## paro     -0.303321  0.738362  0.096908 -3.130 0.001748
## numprop   0.258334  1.294771  0.056161  4.600 4.23e-06
## crimprop  0.304769  1.356312  0.109110  2.793 0.005218
## numarst   0.011737  1.011806  0.003474  3.379 0.000728
## edcomb   -0.053929  0.947500  0.019860 -2.715 0.006620
## 
## Likelihood ratio test=93.22  on 10 df, p=1.227e-15
## n= 1444, number of events= 530

6.2.3 データの非独立性の問題

今回の分析では、同じ個体のデータが複数個入っていることがあるにもかかわらず、全てのデータが独立であると仮定されてしまっている。その結果、先ほどのモデル(mod_tg)の偏回帰係数の標準誤差は実際より小さく推定されてしまっており、結果的に\(z\)値は過大に評価されてしまっている。正しい統計的推論を行うためには、これを修正する必要がある。修正方法は主に二つある。

6.2.3.1 頑強推定

これは、「サンドイッチ」推定法と呼ばれる方法を使って標準誤差を計算することである。Rでは以下のようにcloster(id)を式に加えることで実行できる。

mod_tg_rob <- coxph(Surv(length,arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + cluster(id), data = arrest)

結果は以下の通り。偏回帰係数の推定値自体は先ほどのモデルと全く同じである。一方で、頑強推定を行わない場合の標準誤差(se(coef))よりも、頑強推定での標準誤差(robust se)の方が大きくなっており、先ほどより\(z\)値が小さくなっていることが分かる。

mod_tg_rob

## Call:
## coxph(formula = Surv(length, arrind) ~ fin + age + white + male + 
##     married + paro + numprop + crimprop + numarst + edcomb, data = arrest, 
##     cluster = id)
## 
##               coef exp(coef)  se(coef) robust se      z        p
## fin       0.135958  1.145633  0.090046  0.102972  1.320 0.186721
## age      -0.030654  0.969811  0.006762  0.008399 -3.650 0.000263
## white    -0.156933  0.854761  0.094189  0.113864 -1.378 0.168126
## male      0.339312  1.403981  0.230061  0.290408  1.168 0.242646
## married  -0.093797  0.910468  0.100009  0.122438 -0.766 0.443630
## paro     -0.303321  0.738362  0.096908  0.105280 -2.881 0.003963
## numprop   0.258334  1.294771  0.056161  0.065745  3.929 8.52e-05
## crimprop  0.304769  1.356312  0.109110  0.131040  2.326 0.020030
## numarst   0.011737  1.011806  0.003474  0.004411  2.661 0.007789
## edcomb   -0.053929  0.947500  0.019860  0.025856 -2.086 0.037004
## 
## Likelihood ratio test=93.22  on 10 df, p=1.227e-15
## n= 1444, number of events= 530

6.2.3.2 共用フレイルティ・モデル

一般的なコックス回帰の場合、通常は標準誤差を修正するだけでは理想的ではない。「共用フレイルティ」を回帰モデルの構成要素として含むモデルは、頑強モデルだけでは十分でない点も修正することができる。

基本的なモデルは、標準のコックス回帰モデルにランダム切片¹⁵を加えたランダム効果モデル(混合モデル)である。式(6.2)の\(h_ij(t)\)は個体\(i\)が\(j\)番目の事象を経験するハザード率であり、\(t\)は事象が最後に発生した時点からの経過時間である。\(e_i\)は「共用フレイルティ」で、全ての測定されていない個体ごとの異質性を表す。\(e_i\)は独立変数とは独立であると仮定され、平均0で分散\(\theta\)の正規分布かガンマ分布が仮定される。\(\theta\)が大きいほど、同じ個体のデータ同士の依存度が高くなる。

\[ log(h_ij(t)) = a(t) + b_1x_{i1} + b_2x_{i2} + e_i　\tag{6.2} \]

Rでは、以下のようにして簡単に実行できる。\(e_i\)の分布に正規分布を仮定するなら式にfrailty.gaussian(id)をガンマ分布を仮定するならfrailty.gamma(id)を追加する。

mod_tg_frail <- coxph(Surv(length,arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + frailty.gaussian(id), data = arrest)

結果は以下の通り。\(\theta\)は0.49と推定された。偏回帰係数の標準誤差は頑強推定と似通っている。StataやSASとは異なり、統計検定量は\(\chi^2\)値で出力されるようだ。

summary(mod_tg_frail)

## Call:
## coxph(formula = Surv(length, arrind) ~ fin + age + white + male + 
##     married + paro + numprop + crimprop + numarst + edcomb + 
##     frailty.gaussian(id), data = arrest)
## 
##   n= 1444, number of events= 530 
## 
##                      coef     se(coef) se2      Chisq  DF    p      
## fin                   0.13868 0.104475 0.091060   1.76   1.0 1.8e-01
## age                  -0.03250 0.007734 0.006917  17.66   1.0 2.6e-05
## white                -0.15138 0.108495 0.094794   1.95   1.0 1.6e-01
## male                  0.36425 0.260060 0.232914   1.96   1.0 1.6e-01
## married              -0.10136 0.115930 0.101118   0.76   1.0 3.8e-01
## paro                 -0.30459 0.111573 0.098411   7.45   1.0 6.3e-03
## numprop               0.27501 0.070026 0.058390  15.42   1.0 8.6e-05
## crimprop              0.32415 0.124763 0.110972   6.75   1.0 9.4e-03
## numarst               0.01226 0.004419 0.003650   7.70   1.0 5.5e-03
## edcomb               -0.05746 0.022866 0.019871   6.31   1.0 1.2e-02
## frailty.gaussian(id)                            285.23 188.1 6.1e-06
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## fin         1.1488     0.8705    0.9361    1.4098
## age         0.9680     1.0330    0.9535    0.9828
## white       0.8595     1.1634    0.6949    1.0632
## male        1.4394     0.6947    0.8646    2.3964
## married     0.9036     1.1067    0.7199    1.1341
## paro        0.7374     1.3561    0.5926    0.9177
## numprop     1.3165     0.7596    1.1477    1.5102
## crimprop    1.3829     0.7231    1.0829    1.7659
## numarst     1.0123     0.9878    1.0036    1.0211
## edcomb      0.9442     1.0591    0.9028    0.9874
## 
## Iterations: 6 outer, 29 Newton-Raphson
##      Variance of random effect= 0.4898072 
## Degrees of freedom for terms=   0.8   0.8   0.8   0.8   0.8   0.8   0.7   0.8   0.7   0.8 188.1 
## Concordance= 0.829  (se = 0.008 )
## Likelihood ratio test= 537.6  on 195.7 df,   p=<2e-16

なお、coxmeパッケージのcoxme()関数を用いても共用フレイルティ・モデルの推定は行えるようである。ただし、こちらはランダム効果の分布に正規分布のみが仮定できる。

mod_tg_frail_b <- coxme(Surv(length,arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + (1|id), data = arrest)

結果は以下の通り。\(\theta\)は0.52と推定され、偏回帰係数の推定値も少し違うようだ。この違いが何によるものかは不明だが、おそらく推定方法の違いによるものと思われる。

summary(mod_tg_frail_b)

## Cox mixed-effects model fit by maximum likelihood
##   Data: arrest
##   events, n = 530, 1444
##   Iterations= 13 70 
##                     NULL Integrated    Fitted
## Log-likelihood -3663.536  -3598.781 -3386.021
## 
##                    Chisq     df p    AIC     BIC
## Integrated loglik 129.51  11.00 0 107.51   60.51
##  Penalized loglik 555.03 206.67 0 141.68 -741.41
## 
## Model:  Surv(length, arrind) ~ fin + age + white + male + married + paro +      numprop + crimprop + numarst + edcomb + (1 | id) 
## Fixed coefficients
##                 coef exp(coef)    se(coef)     z       p
## fin       0.13891349 1.1490247 0.105410172  1.32 1.9e-01
## age      -0.03262094 0.9679054 0.007796089 -4.18 2.9e-05
## white    -0.15146381 0.8594490 0.109447779 -1.38 1.7e-01
## male      0.36627568 1.4423528 0.262117579  1.40 1.6e-01
## married  -0.10216178 0.9028835 0.116970737 -0.87 3.8e-01
## paro     -0.30459710 0.7374204 0.112528421 -2.71 6.8e-03
## numprop   0.27621349 1.3181292 0.070848073  3.90 9.7e-05
## crimprop  0.32554843 1.3847899 0.125788213  2.59 9.7e-03
## numarst   0.01227966 1.0123554 0.004473929  2.74 6.1e-03
## edcomb   -0.05770664 0.9439268 0.023065816 -2.50 1.2e-02
## 
## Random effects
##  Group Variable  Std Dev   Variance 
##  id    Intercept 0.7222207 0.5216028

6.2.3.3 モデルの比較

「時間のギャップ」を用いた3つのモデルで推定された偏回帰係数とそれを指数変換したもの、\(z\)値をまとめたのが表6.3である。なお、頑強推定(mod_tg_rob)でも通常のコックス回帰(mod_tg)と偏回帰係数の推定値自体は変わらないので、それらは省略している。共用フレイルティのあるモデルでは、統計的に有意であるような変数の偏回帰係数の推定値(Coef(frail))が通常のコックス回帰(Coef)よりも大きくなっている。\(z\)値(Z(frail))は基本的には頑強推定のもの(Z(robust))と似通っているが、有意な変数では頑強推定のものより高くなっている。しかし、どの変数も通常のコックス回帰の\(z\)値(Z)よりは小さくなっている。

tibble(Covariate = rownames(mod_tg$coefficients %>% data.frame()),
       "Coef" = c(sprintf("%.3f",coef(mod_tg))[1:10]),
       "Exp(Coef)" = c(sprintf("%.3f",exp(coef(mod_tg)))[1:10]),
       "Z" = sprintf("%.3f",coef(mod_tg)/sqrt(diag(vcov(mod_tg))))[1:10],
       "Z(robust)" = sprintf("%.3f",coef(mod_tg)/sqrt(diag(vcov(mod_tg_rob))))[1:10],
       "Coef(frail)" = sprintf("%.3f",coef(mod_tg_frail))[1:10],
       "Exp(Coef_frail)" = sprintf("%.3f",exp(coef(mod_tg_frail))[1:10]),
       "Z(frail)" = sprintf("%.3f",coef(mod_tg_frail)/sqrt(diag(vcov(mod_tg_frail))))[1:10]) %>% 
  flextable() %>% 
  add_header_row(colwidth = c(1,4,3),
                 values = c("","コックス回帰\n(通常＋頑強推定)", "コックス回帰\n(共用フレイルティ)")) %>% 
  flextable::align(align = "center", part = "all") %>% 
  set_caption("各モデルの推定値の比較")

表6.3: 各モデルの推定値の比較
	コックス回帰 (通常＋頑強推定)				コックス回帰 (共用フレイルティ)
Covariate	Coef	Exp(Coef)	Z	Z(robust)	Coef(frail)	Exp(Coef_frail)	Z(frail)
fin	0.136	1.146	1.510	1.320	0.139	1.149	1.327
age	-0.031	0.970	-4.533	-3.650	-0.032	0.968	-4.202
white	-0.157	0.855	-1.666	-1.378	-0.151	0.860	-1.395
male	0.339	1.404	1.475	1.168	0.364	1.439	1.401
married	-0.094	0.910	-0.938	-0.766	-0.101	0.904	-0.874
paro	-0.303	0.738	-3.130	-2.881	-0.305	0.737	-2.730
numprop	0.258	1.295	4.600	3.929	0.275	1.317	3.927
crimprop	0.305	1.356	2.793	2.326	0.324	1.383	2.598
numarst	0.012	1.012	3.379	2.661	0.012	1.012	2.775
edcomb	-0.054	0.947	-2.715	-2.086	-0.057	0.944	-2.513

なお、ここで紹介したモデルは全てパラメトリックなモデル(例えばワイブル回帰モデル)に対しても同様に適用できる。

6.3 観察開始からの時間に基づく方法

6.3.1 モデルの概要

ほとんどの場合では、時間のギャップを用いるのが最善であると考えられる。しかし、場合によってはハザード率を他の時間の長さに依存させた方が妥当かもしれない。例えば、再逮捕のハザード率は最後の逮捕からの時間ではなく、最初に刑務所を出所してからの時間に依存しているかもしれない。

この分析でも先ほどと同じデータを使用するが、一つ一つの観察期間を従属変数にするのではなく、時間依存変数があるときのように各行のデータの開始時点と終了時点を明示的に指定する形になる。データarrestには、各行の観察開始時点を表す変数beginと終了を表す変数endが存在する。

6.3.2 Rでの実装

Rでは、従属変数をSurv(begin, end, arrind)とすることでこのような推定ができる。今回も同じ個体のデータが複数あるので、頑強モデルと共用フレイルティモデルの両方を推定する。

mod_st_rob <- coxph(Surv(begin, end, arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + cluster(id), data = arrest)

mod_st_frail <- coxph(Surv(begin, end, arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + frailty.gaussian(id), data = arrest)

両モデルで推定された偏回帰係数とそれを指数変換したもの、\(z\)値をまとめたのが表6.4である。いずれも時間のギャップを用いたモデル(表6.4)と非常によく似ている。

tibble(Covariate = rownames(mod_st_rob$coefficients %>% data.frame()),
       "Coef(rob)" = c(sprintf("%.3f",coef(mod_st_rob))[1:10]),
       "Exp(Coef_rob)" = c(sprintf("%.3f",exp(coef(mod_st_rob)))[1:10]),
       "Z(rob)" = sprintf("%.3f",coef(mod_st_rob)/sqrt(diag(vcov(mod_st_rob))))[1:10],
       "Coef(frail)" = sprintf("%.3f",coef(mod_st_frail))[1:10],
       "Exp(Coef_frail)" = sprintf("%.3f",exp(coef(mod_st_frail))[1:10]),
       "Z(frail)" = sprintf("%.3f",coef(mod_st_frail)/sqrt(diag(vcov(mod_st_frail))))[1:10]) %>% 
  flextable() %>% 
  add_header_row(colwidth = c(1,3,3),
                 values = c("","コックス回帰\n(頑強推定)", "コックス回帰\n(共用フレイルティ)")) %>% 
  flextable::align(align = "center", part = "all") %>% 
  set_caption("各モデルの推定値の比較")

表6.4: 各モデルの推定値の比較
	コックス回帰 (頑強推定)			コックス回帰 (共用フレイルティ)
Covariate	Coef(rob)	Exp(Coef_rob)	Z(rob)	Coef(frail)	Exp(Coef_frail)	Z(frail)
fin	0.140	1.150	1.320	0.140	1.151	1.343
age	-0.031	0.970	-3.638	-0.032	0.968	-4.165
white	-0.153	0.858	-1.311	-0.144	0.866	-1.331
male	0.344	1.411	1.147	0.361	1.434	1.390
married	-0.144	0.866	-1.140	-0.136	0.873	-1.171
paro	-0.313	0.731	-2.938	-0.309	0.734	-2.773
numprop	0.278	1.320	4.092	0.283	1.327	4.055
crimprop	0.303	1.354	2.257	0.316	1.372	2.537
numarst	0.011	1.011	2.466	0.012	1.012	2.630
edcomb	-0.058	0.944	-2.151	-0.059	0.943	-2.583

6.3.3 観察開始時点からの時間を使う利点

この方法の魅力の一つは、独立変数の影響が観察期間全体で変化するかを調べることができる点である。これは、独立変数と開始時点からの時間(\(t\))の交互作用項を含めることでモデリングできる。

\[ log(h_ij(t)) = a(t) + b_1x_{i1} + b_2x_{i2} + b_3x_{i2}t \]

実際に、先ほどの頑強推定モデル(mod_st_rob)で比例ハザード性が満たされているか、シェーンフィールド残差を用いて調べたところ、人種(white)、前科の数(numarst)、教育レベル(edcomb)は満たされていないことが示唆された。

cox.zph(mod_st_rob)

##           chisq df      p
## fin       1.549  1 0.2133
## age       0.606  1 0.4361
## white     6.917  1 0.0085
## male      1.945  1 0.1631
## married   0.713  1 0.3986
## paro      3.570  1 0.0588
## numprop   1.062  1 0.3029
## crimprop  0.192  1 0.6615
## numarst   9.802  1 0.0017
## edcomb    8.450  1 0.0037
## GLOBAL   37.849 10  4e-05

そこで、これらと時間の交互作用を含むモデルを考える。

mod_st_rob_int <- coxph(Surv(begin, end, arrind) ~ fin + age + white + male + married + paro + numprop + 
                       crimprop + numarst + edcomb + white:end + numarst:end + edcomb:end + cluster(id), data = arrest)

結果が以下の通り。交互作用はいずれも5%水準で有意になっていることが分かる。

mod_st_rob_int

## Call:
## coxph(formula = Surv(begin, end, arrind) ~ fin + age + white + 
##     male + married + paro + numprop + crimprop + numarst + edcomb + 
##     white:end + numarst:end + edcomb:end, data = arrest, cluster = id)
## 
##                   coef  exp(coef)   se(coef)  robust se       z        p
## fin          1.484e-01  1.160e+00  9.381e-02  1.098e-01   1.352 0.176359
## age         -5.265e-02  9.487e-01  7.309e-03  9.226e-03  -5.706 1.15e-08
## white        1.072e+00  2.921e+00  2.072e-01  2.813e-01   3.810 0.000139
## male         5.813e-01  1.788e+00  2.368e-01  2.913e-01   1.995 0.045996
## married     -2.620e-01  7.695e-01  1.071e-01  1.217e-01  -2.153 0.031297
## paro        -3.028e-01  7.388e-01  1.005e-01  1.158e-01  -2.614 0.008953
## numprop      2.184e-01  1.244e+00  6.102e-02  7.316e-02   2.985 0.002840
## crimprop     7.244e-02  1.075e+00  1.140e-01  1.206e-01   0.601 0.548112
## numarst      3.972e-02  1.041e+00  6.570e-03  1.177e-02   3.375 0.000739
## edcomb       1.173e+00  3.231e+00  4.338e-02  7.149e-02  16.403  < 2e-16
## white:end   -3.368e-03  9.966e-01  8.888e-04  1.095e-03  -3.075 0.002104
## numarst:end -2.565e-04  9.997e-01  4.673e-05  8.041e-05  -3.189 0.001426
## edcomb:end  -4.789e-03  9.952e-01  1.548e-04  2.492e-04 -19.215  < 2e-16
## 
## Likelihood ratio test=2604  on 13 df, p=< 2.2e-16
## n= 1444, number of events= 530

例えば前科の数と教育レベルに着目すると、それぞれの独立変数の偏回帰係数(\(b+ct\))は時間と共に減少していき、途中で負に転じることが分かる(表6.5)。

tibble(t = seq(0,365,50),
       `偏回帰係数\n(b + ct)` = sprintf("%.3f",mod_st_rob_int$coefficients[[9]] + mod_st_rob_int$coefficients[[12]]*t),
         `Exp(b+ct)` = sprintf("%.3f", exp(mod_st_rob_int$coefficients[[9]] + mod_st_rob_int$coefficients[[12]]*t)),
       `偏回帰係数\n(b + ct) ` = sprintf("%.3f", mod_st_rob_int$coefficients[[10]] + mod_st_rob_int$coefficients[[13]]*t),
         `Exp(b+ct) ` = sprintf("%.3f", exp(mod_st_rob_int$coefficients[[10]] + mod_st_rob_int$coefficients[[13]]*t))) %>% 
  rename(`出所して\nからの日数(t)` = t) %>% 
  flextable() %>% 
  width(width = 2/2) %>% 
  add_header_row(colwidth = c(1,2,2),
                 values = c("","前科の数","教育レベル")) %>% 
  flextable::align(align = "center", part = "all") %>% 
  set_caption("モデルから推定された交互作用項を含む偏回帰係数")

表6.5: モデルから推定された交互作用項を含む偏回帰係数
	前科の数		教育レベル
出所してからの日数(t)	偏回帰係数 (b + ct)	Exp(b+ct)	偏回帰係数 (b + ct)	Exp(b+ct)
0	0.040	1.041	1.173	3.231
50	0.027	1.027	0.933	2.543
100	0.014	1.014	0.694	2.001
150	0.001	1.001	0.454	1.575
200	-0.012	0.988	0.215	1.240
250	-0.024	0.976	-0.025	0.976
300	-0.037	0.963	-0.264	0.768
350	-0.050	0.951	-0.504	0.604

6.4 分析モデルの拡張

繰り返し発生する事象について、これまで検討した方法はさまざまに拡張できる。例えば競合リスク・モデルや離散時間モデルについても同様の方法が適用できる。一つの個体が複数の観察記録を持つ問題については、本章で解説したように頑強推定を用いたり、ランダム効果モデルを用いたりすることで対処できる。

References

久保拓弥. (2012). データ解析のための統計モデリング入門. 岩波書店.

大東健太郎. (2010). 線形モデルから一般化線形モデル（GLM）へ. 雑草研究, 55(4), 268–274.

粕谷英一. (2012). 一般化線形モデル. 共立出版.

過分散については大東 (2010) や粕谷 (2012) を参照。負の二項分布については、こちらも参照。簡単に言えば負の二項分布はポワソン分布の期待値\(\lambda\)が毎回ガンマ分布から得られると仮定する混合分布である。よって、ポワソン分布は負の二項分布の特殊例といえる。↩︎
オフセット項については久保 (2012) を参照。↩︎
ランダム効果については、久保 (2012) や粕谷 (2012) を参照。↩︎