StatsmodelsはPythonというプログラミング言語上で動く統計解析ソフトです。statsmodelsのサンプルを動かすにはPCにPythonがインストールされている必要があります。まだインストールされていない方はJupyter notebookのインストールを参照してください。Jupyter notebookはstatsmodelsを動かすのに大変便利です。

多くの部分はstatsmodelsのリファレンスなどを参考にしています。

一般化推定方程式 (Generalized Estimating Equations)

一般化推定方程式は、パネル、クラスタ、反復測定データという観測値がクラスター内で相関を持ち、クラスター外で相関をもたないデータの集団における平均的な効果を推定することを目的としたモデルで、一般化線形モデルの1つである。クラスター内では相関をもつとしているが、相関が生じる原因を特定していない、または特定できないともいえる。平均構造と相関構造を別々に扱える。GLMは

従属変数$y_i$が指数分布族にしたがう。
従属変数がそれぞれ独立である。

という仮定のもとに

確率変数$y_i$の確率分布を仮定する。
$y_i$の期待値$\mu$の線形回帰モデルを設定する。
$y_i$の仮定した確率分布の尤度関数を作りそれを最大化する$\beta(\hat{\beta})$を求める。
最尤推定量$\hat{\beta}$は一致性と漸近正規性をもつ。

として、モデルを推定する。

GEEでは、

従属変数$y_i$が指数分布族にしたがうという仮定を緩める。
従属変数が独立であるという仮定も緩める。
疑似尤度を用いて未知な分布に対応する。
作業用相関行列(working correlation matrix) $R_i(\alpha)$を相関構造を表現するために導入する。
一般化推定方程式を2段階に分けて数値計算で解く。

という順序でモデルを推定する。したがって、データの多様性に対処できるが、情報量も必要になる。

例: てんかん患者の発作の回数

Thall と Vail (1990) は59人のてんかん患者の8週間のベースライン期間と4回の連続した2週間ごとの発作回数について分析した。ベースライン期間においては患者への処置は施さない。ベースライン期間後に患者は無作為にプラセボまたはプロガビドを処方される。処方後の観測は４ｘ２＝８週間行われる。

データの詳細：
使用:てんかん
フォマット：236行9列

y：2週間の発作の数
trt: 処置、プラセボ、またはプロガビド
base：8週間のベースライン期間の発作の数
age：被験者の年齢、歳
V4：4期間の0/1指示変数
subject：被験者の番号, 1 to 59.
period : 期間, 1 to 4.
lbase : ベースライン期間での発作数の対数、平均をゼロに標準化
lage : 年齢の対数、平均をゼロに標準化

てんかんの発作は生存時間解析の再発事象に相当する。再発事象データは，同一被験者の発作の発生を繰り返し測定することで得られ、経時データ (longitudinal data) や反復測定データ (repeated measures data) の一種としてとらえられる。被験者が異なればてんかんの発作の回数に相関がないと考えられるが，同一被験者についての測定であれば結果には相関があると考えるのが自然である。

薬剤投与前にベースライン期間として8週間，薬剤投与後は2週間間隔で計 4 回の観測を行い、区間ごとに発作回数を調べている。被験者の年齢だけが共変量である。

まず、初期化してデータを得る。

import statsmodels.formula.api as smf
data = sm.datasets.get_rdataset('epil', package='MASS').data

記述統計

データの様子を調べる。

data.head(10)

y	trt	base	age	V4	subject	period	lbase	lage
0	5	placebo	11	31	0	1	1	-0.756354	0.114204
1	3	placebo	11	31	0	1	2	-0.756354	0.114204
2	3	placebo	11	31	0	1	3	-0.756354	0.114204
3	3	placebo	11	31	1	1	4	-0.756354	0.114204
4	3	placebo	11	30	0	2	1	-0.756354	0.081414
5	5	placebo	11	30	0	2	2	-0.756354	0.081414
6	3	placebo	11	30	0	2	3	-0.756354	0.081414
7	3	placebo	11	30	1	2	4	-0.756354	0.081414
8	2	placebo	6	25	0	3	1	-1.362490	-0.100908
9	4	placebo	6	25	0	3	2	-1.362490	-0.100908

プラセボという有効成分を含まない薬とプロガビドというGABA類似体系抗てんかん薬を処方された後のてんかんの発作回数の分布を比べてみる。データの数が前者が112で後者が124である。

data[data.trt!="placebo"].y.hist()
data[data.trt=="placebo"].y.hist()

data[data.trt!="placebo"].y.count(),data[data.trt=="placebo"].y.count()

(124, 112)

標準化したベースライン期間の発作回数の対数と年齢の対数を頻度図にしてみる。

data.lbase.hist()

data.lage.hist()

つぎに記述統計を得る。ベースライン期間の発作の回数の平均と標準偏差

data.base.mean(),data.base.std()

(31.220338983050848, 26.705051109822254)

プロガビドを施した後の発作の回数の平均と標準偏差

data[data.trt!="placebo"].y.mean(),data[data.trt!="placebo"].y.std()

(7.959677419354839, 13.92978863002629)

プロガビドを施した患者のベースライン期間の発作回数の平均と標準偏差

data[data.trt!="placebo"].base.mean(),data[data.trt!="placebo"].base.std()

(31.612903225806452, 27.638405492528324)

プラセボを施した後の発作の回数の平均と標準偏差

data[data.trt=="placebo"].y.mean(),data[data.trt=="placebo"].y.std()

(8.580357142857142, 10.369426989352696)

プラセボを施した患者のベースライン期間の発作回数の平均と標準偏差

data[data.trt=="placebo"].base.mean(),data[data.trt=="placebo"].base.std()

(30.785714285714285, 25.749111266541444)

被験者の年齢とベースライン期間の発作の回数の対数をプロットしている。

plt.scatter(data.lage,data.lbase)
rg=sm.OLS(data.lbase,data.lage).fit()
plt.plot(data.lage,rg.fittedvalues)

推測統計

モデルを推定する。発作発生回数をポアソン分布にしたがうとする。

作業用相関行列の違い

independence()

作業用相関行列について、それぞれの患者の反復測定の対となるデータは独立であると仮定する。
"subject"はgroupsを示している。

fam = sm.families.Poisson()
ind = sm.cov_struct.Independence()
mod = smf.gee("y ~ age + trt + base", "subject", data,
              cov_struct=ind, family=fam)

res = mod.fit()
print(res.summary())

GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  236
Model:                                 GEE   No. clusters:                       59
Method:                        Generalized   Min. cluster size:                   4
                      Estimating Equations   Max. cluster size:                   4
Family:                            Poisson   Mean cluster size:                 4.0
Dependence structure:         Independence   Num. iterations:                     2
Date:                     Sat, 02 Nov 2019   Scale:                           1.000
Covariance type:                    robust   Time:                         13:06:54
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept            0.5730      0.361      1.589      0.112      -0.134       1.280
trt[T.progabide]    -0.1519      0.171     -0.888      0.375      -0.487       0.183
age                  0.0223      0.011      1.960      0.050    2.11e-06       0.045
base                 0.0226      0.001     18.451      0.000       0.020       0.025
==============================================================================
Skew:                          3.7823   Kurtosis:                      28.6672
Centered skew:                 2.7597   Centered kurtosis:             21.9865
==============================================================================

つぎに被験者ごとの発作の回数の平均値と分散の対数をとってブロットしてみる。

yg = mod.cluster_list(np.asarray(data["y"]))
ymn = [x.mean() for x in yg]
yva = [x.var() for x in yg]
plt.grid(True)
plt.plot(np.log(ymn), np.log(yva), 'o')
plt.plot([-2, 6], [-2, 6], '-', color='grey')
plt.xlabel("Log variance", size=17)
plt.ylabel("Log mean", size=17)

exchangeable()

この場合には、作業用相関行列についてそれぞれの患者の反復測定の対となるデータは同じ相関をもつと仮定する。

fam = sm.families.Poisson()
ex = sm.cov_struct.Exchangeable()
mod = smf.gee("y ~ age + trt + base", "subject", data,
              cov_struct=ex, family=fam)

res = mod.fit()
print(res.summary())

 GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  236
Model:                                 GEE   No. clusters:                       59
Method:                        Generalized   Min. cluster size:                   4
                      Estimating Equations   Max. cluster size:                   4
Family:                            Poisson   Mean cluster size:                 4.0
Dependence structure:         Exchangeable   Num. iterations:                     2
Date:                     Thu, 31 Oct 2019   Scale:                           1.000
Covariance type:                    robust   Time:                         07:23:48
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept            0.5730      0.361      1.589      0.112      -0.134       1.280
trt[T.progabide]    -0.1519      0.171     -0.888      0.375      -0.487       0.183
age                  0.0223      0.011      1.960      0.050    2.11e-06       0.045
base                 0.0226      0.001     18.451      0.000       0.020       0.025
==============================================================================
Skew:                          3.7823   Kurtosis:                      28.6672
Centered skew:                 2.7597   Centered kurtosis:             21.9865
==============================================================================

fig = res.plot_isotropic_dependence()
fig.get_axes()[0].set_ylim(0, 1)

rslts = res.params_sensitivity(0., 0.9, 5)
dp = [x.cov_struct.dep_params for x in rslts]
age_params = np.asarray([x.params[2] for x in rslts])
age_se = np.asarray([x.bse[2] for x in rslts])
plt.grid(True)
plt.plot(dp, age_params, '-', color='orange', lw=4)
plt.plot(dp, age_params + age_se, '-', color='grey')
plt.plot(dp, age_params - age_se, '-', color='grey')
plt.axvline(res.cov_struct.dep_params, ls='--', color='black')
plt.xlabel("Autocorrelation parameter", size=14)
plt.ylabel("Age effect", size=14)

plt.plot(res.fittedvalues, res.resid, 'o', alpha=0.5)
plt.xlabel("Fitted values", size=17)
plt.ylabel("Residuals", size=17)

Autoregressive()

この場合には、作業用相関行列についてそれぞれの患者の反復測定の対となるデータは自己相関をもつと仮定する。

fam = sm.families.Poisson()
au = sm.cov_struct.Autoregressive()
mod = smf.gee("y ~ age + trt + base", "subject", data,
              cov_struct=au, family=fam)

res = mod.fit()
print(res.summary())

GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  236
Model:                                 GEE   No. clusters:                       59
Method:                        Generalized   Min. cluster size:                   4
                      Estimating Equations   Max. cluster size:                   4
Family:                            Poisson   Mean cluster size:                 4.0
Dependence structure:       Autoregressive   Num. iterations:                    11
Date:                     Sat, 02 Nov 2019   Scale:                           1.000
Covariance type:                    robust   Time:                         13:08:52
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
Intercept            0.4123      0.422      0.977      0.329      -0.415       1.240
trt[T.progabide]    -0.0704      0.164     -0.430      0.667      -0.391       0.250
age                  0.0274      0.013      2.108      0.035       0.002       0.053
base                 0.0223      0.001     18.252      0.000       0.020       0.025
==============================================================================
Skew:                          3.9903   Kurtosis:                      30.1753
Centered skew:                 2.7597   Centered kurtosis:             21.9865
==============================================================================

fig = res.plot_isotropic_dependence()
fig.get_axes()[0].set_ylim(0, 1)

rslts = res.params_sensitivity(0., 0.9, 5)
dp = [x.cov_struct.dep_params for x in rslts]
age_params = np.asarray([x.params[2] for x in rslts])
age_se = np.asarray([x.bse[2] for x in rslts])
plt.grid(True)
plt.plot(dp, age_params, '-', color='orange', lw=4)
plt.plot(dp, age_params + age_se, '-', color='grey')
plt.plot(dp, age_params - age_se, '-', color='grey')
plt.axvline(res.cov_struct.dep_params, ls='--', color='black')
plt.xlabel("Autocorrelation parameter", size=14)
plt.ylabel("Age effect", size=14)

試験薬の投与がプラセボと比較して発作の回数を減らしているかどうかを、作業用相関行列を変えながら検定したが、この臨床試験から一般化推定方程式を用いた場合には有意だという結果は得られない。

例：てんかん発作（Leppik et.al.1987）

同じ59人のてんかん発作のデータを用いて、一般推定方程式を用いて臨床試験の結果を評価した場合を見てみる。

まずはデータの構造を変えて、ベースライン期間の発作回数を$y_{ij}$の中に入れる。
$$log E(Y_{ij})=log(T_{ij})+\beta_0+\beta_1x_{ij1}+\beta_2x_{ij2}+\beta_3x_{ij1}x_{ij2}$$
ここで$x_{ij1}$はランダム化の前(0)と後(1)の時期を二値で表現している。$x_{ij2}$はプラセボ(0)、試験薬(1)の治療群を二値で示している。$T_{ij}$はベースライン期間の8週間と治療後の各測定期間の2週間の区間の長さである。

import pandas as pd
yy=pd.DataFrame(np.arange(59)+1,columns=['subject'])
yy['period']=np.log(8)
yy['y']=np.NaN
yy['trt']=0#np.NaN
yy['btrt']=0
data['btrt']=1
data['period']=np.log(2)

ii=0
for i in range(len(data)):
    if data.subject.iloc[i]==yy.subject.iloc[ii]:
        yy.iloc[ii,2]=data.base.iloc[i]
        yy.iloc[ii,3]=data.trt.iloc[i]
        yy.iloc[ii,0]=data.subject.iloc[i]
        ii+=1
        if ii>=59:
            break
yy=pd.concat([data,yy],axis=0,join='inner')

for i in range(len(yy)):
    if yy.iloc[i,1]=='progabide':
        yy.iloc[i:i+1,1]=1
    if yy.iloc[i,1]=='placebo':
        yy.iloc[i:i+1,1]=0

mod = smf.gee("y ~ btrt + trt ", "subject", yy,
              cov_struct=ex, family=fam,offset='period')

res = mod.fit()
print(res.summary())

  GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  295
Model:                                 GEE   No. clusters:                       59
Method:                        Generalized   Min. cluster size:                   5
                      Estimating Equations   Max. cluster size:                   5
Family:                            Poisson   Mean cluster size:                 5.0
Dependence structure:         Exchangeable   Num. iterations:                     6
Date:                     Sat, 02 Nov 2019   Scale:                           1.000
Covariance type:                    robust   Time:                         23:30:28
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.3240      0.168      7.883      0.000       0.995       1.653
btrt           0.0560      0.106      0.530      0.596      -0.151       0.263
trt            0.0705      0.213      0.331      0.740      -0.346       0.487
==============================================================================
Skew:                          3.3538   Kurtosis:                      15.0311
Centered skew:                 2.0882   Centered kurtosis:              6.6169
==============================================================================

臨床試験の試験群が有意であるという結果は得られない。そこで少し変えて見える。共変量である被験者の年齢を加えてみる。

import pandas as pd
yy=pd.DataFrame(np.arange(59)+1,columns=['subject'])
yy['period']=0
yy['y']=np.NaN
yy['trt']=0#np.NaN
yy['btrt']=0
yy['age']=0
data = sm.datasets.get_rdataset('epil', package='MASS').data
data['btrt']=1

ii=0
for i in range(len(data)):
    if data.subject.iloc[i]==yy.subject.iloc[ii]:
        yy.iloc[ii,2]=data.base.iloc[i]
        #yy.iloc[ii,3]=data.trt.iloc[i]
        yy.iloc[ii,0]=data.subject.iloc[i]
        yy.iloc[ii,5]=data.age.iloc[i]
        ii+=1
        if ii>=59:
            break
#print(yy.head(200))
yy=pd.concat([data,yy],axis=0,join='inner')

for i in range(len(yy)):
    if yy.iloc[i,1]=='progabide':
        yy.iloc[i,1]=1
    if yy.iloc[i,1]=='placebo':
        yy.iloc[i,1]=0
yy[yy.trt==0].count()

mod = smf.gee("y ~ age + btrt + trt ", "subject", yy,
              cov_struct=ind, family=fam,offset='period')

res = mod.fit()
print(res.summary())

GEE Regression Results                              
===================================================================================
Dep. Variable:                           y   No. Observations:                  295
Model:                                 GEE   No. clusters:                       59
Method:                        Generalized   Min. cluster size:                   5
                      Estimating Equations   Max. cluster size:                   5
Family:                            Poisson   Mean cluster size:                 5.0
Dependence structure:         Independence   Num. iterations:                     2
Date:                     Sat, 02 Nov 2019   Scale:                           1.000
Covariance type:                    robust   Time:                         23:48:30
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      3.9937      0.614      6.503      0.000       2.790       5.197
age           -0.0198      0.021     -0.954      0.340      -0.060       0.021
btrt          -4.3316      0.154    -28.050      0.000      -4.634      -4.029
trt           -0.1014      0.335     -0.303      0.762      -0.758       0.555
==============================================================================
Skew:                          2.6891   Kurtosis:                      13.3017
Centered skew:                 0.3066   Centered kurtosis:              3.7456
==============================================================================

有意にはならないという点は変わらない。

参考文献：

再発事象データの解析
 Wiki notebooks for GEE
repeated measures of epileptic seizure counts