GEE スコア検定¶
このノートブックは、シミュレーションを使用して頑健な GEE スコア検定を実証しています。これらの検定は、GEE 分析において、平均構造に関する入れ子になった仮説を比較するために使用できます。これらの検定は、ワーキング相関モデルの誤指定、および分散構造のある種の誤指定（例えば、準ポアソン分析における尺度パラメータで捉えられるもの）に対して頑健です。
データは、クラスタ内で依存性があり、クラスタ間では依存性がないクラスタとしてシミュレートされます。クラスタごとの依存性は、コプラアプローチを使用して誘導されます。データは周辺的に負の二項（ガンマ/ポアソン）混合に従います。
検定の性能を評価するために、検定のレベルと検出力について以下で検討します。
[1]:
import pandas as pd
import numpy as np
from scipy.stats.distributions import norm, poisson
import statsmodels.api as sm
import matplotlib.pyplot as plt
次のセルで定義されている関数は、コプラアプローチを使用して、周辺的に負の二項分布に従う相関のあるランダム値をシミュレートします。入力パラメータuは(0, 1)の値の配列です。uの要素は、(0, 1)上で周辺的に一様に分布している必要があります。uの相関は、返される負の二項値に相関を誘導します。配列パラメータmuは周辺平均を与え、スカラーパラメータscaleは平均/分散の関係を定義します（分散はscaleと平均の積です）。uとmuの長さは同じでなければなりません。
[2]:
def negbinom(u, mu, scale):
p = (scale - 1) / scale
r = mu * (1 - p) / p
x = np.random.gamma(r, p / (1 - p), len(u))
return poisson.ppf(u, mu=x)
以下は、シミュレーションで使用されるデータを制御するいくつかのパラメータです。
[3]:
# Sample size
n = 1000
# Number of covariates (including intercept) in the alternative hypothesis model
p = 5
# Cluster size
m = 10
# Intraclass correlation (controls strength of clustering)
r = 0.5
# Group indicators
grp = np.kron(np.arange(n/m), np.ones(m))
シミュレーションでは、固定された設計行列を使用します。
[4]:
# Build a design matrix for the alternative (more complex) model
x = np.random.normal(size=(n, p))
x[:, 0] = 1
帰無仮説の設計行列は、対立仮説の設計行列に入れ子になっています。それは、対立仮説の設計行列よりも2ランク小さいです。
[5]:
x0 = x[:, 0:3]
GEE スコア検定は、依存性と過分散に対して頑健です。ここでは、過分散パラメータを設定します。各観測値の負の二項分布の分散は、その平均値のscale倍に等しくなります。
[6]:
# Scale parameter for negative binomial distribution
scale = 10
次のセルでは、帰無仮説と対立仮説の平均構造を設定します。
[7]:
# The coefficients used to define the linear predictors
coeff = [[4, 0.4, -0.2], [4, 0.4, -0.2, 0, -0.04]]
# The linear predictors
lp = [np.dot(x0, coeff[0]), np.dot(x, coeff[1])]
# The mean values
mu = [np.exp(lp[0]), np.exp(lp[1])]
以下は、シミュレーションを実行する関数です。
[8]:
# hyp = 0 is the null hypothesis, hyp = 1 is the alternative hypothesis.
# cov_struct is a statsmodels covariance structure
def dosim(hyp, cov_struct=None, mcrep=500):
# Storage for the simulation results
scales = [[], []]
# P-values from the score test
pv = []
# Monte Carlo loop
for k in range(mcrep):
# Generate random "probability points" u that are uniformly
# distributed, and correlated within clusters
z = np.random.normal(size=n)
u = np.random.normal(size=n//m)
u = np.kron(u, np.ones(m))
z = r*z +np.sqrt(1-r**2)*u
u = norm.cdf(z)
# Generate the observed responses
y = negbinom(u, mu=mu[hyp], scale=scale)
# Fit the null model
m0 = sm.GEE(y, x0, groups=grp, cov_struct=cov_struct, family=sm.families.Poisson())
r0 = m0.fit(scale='X2')
scales[0].append(r0.scale)
# Fit the alternative model
m1 = sm.GEE(y, x, groups=grp, cov_struct=cov_struct, family=sm.families.Poisson())
r1 = m1.fit(scale='X2')
scales[1].append(r1.scale)
# Carry out the score test
st = m1.compare_score_test(r0)
pv.append(st["p-value"])
pv = np.asarray(pv)
rslt = [np.mean(pv), np.mean(pv < 0.1)]
return rslt, scales
独立ワーキング共分散構造を使用してシミュレーションを実行します。帰無仮説の下では平均は約0になり、対立仮説の下でははるかに低くなると予想されます。同様に、帰無仮説の下では、p値の約10％が0.1未満であり、対立仮説の下では0.1未満のp値の割合がはるかに大きくなると予想されます。
[9]:
rslt, scales = [], []
for hyp in 0, 1:
s, t = dosim(hyp, sm.cov_struct.Independence())
rslt.append(s)
scales.append(t)
rslt = pd.DataFrame(rslt, index=["H0", "H1"], columns=["Mean", "Prop(p<0.1)"])
print(rslt)
Mean Prop(p<0.1)
H0 0.482817 0.106
H1 0.051356 0.858
次に、尺度パラメータの推定値が妥当であることを確認します。GEE スコア検定の依存性と過分散に対する頑健性を評価しているので、ここでは、予想通りに過分散が存在することを確認しています。
[10]:
_ = plt.boxplot([scales[0][0], scales[0][1], scales[1][0], scales[1][1]])
plt.ylabel("Estimated scale")
[10]:
Text(0, 0.5, 'Estimated scale')
次に、交換可能なワーキング相関モデルを使用して同じ分析を行います。これは、独立したワーキング相関を使用した上記の例よりも遅くなるため、モンテカルロ反復回数を少なくします。
[11]:
rslt, scales = [], []
for hyp in 0, 1:
s, t = dosim(hyp, sm.cov_struct.Exchangeable(), mcrep=100)
rslt.append(s)
scales.append(t)
rslt = pd.DataFrame(rslt, index=["H0", "H1"], columns=["Mean", "Prop(p<0.1)"])
print(rslt)
Mean Prop(p<0.1)
H0 0.503855 0.10
H1 0.055889 0.81