GEE 评分检验¶
此笔记本使用模拟来演示稳健的 GEE 评分检验。这些检验可以在 GEE 分析中使用，用于比较关于均值结构的嵌套假设。这些检验对工作相关模型的错误指定以及对方差结构的某些形式的错误指定（例如，在准泊松分析中由尺度参数捕获）具有鲁棒性。
数据模拟为聚类，其中聚类内部存在依赖性，但聚类之间不存在依赖性。聚类依赖性是使用 copula 方法引起的。数据边缘服从负二项式（伽马/泊松）混合。
以下是考虑的检验水平和功效，以评估检验的性能。
[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
以下单元格中定义的函数使用 copula 方法来模拟边缘服从负二项式分布的相关随机值。输入参数 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
零设计矩阵嵌套在备择设计矩阵中。它的秩比备择设计矩阵少两个。
[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，在备择假设下均值要低得多。同样，我们预计在零假设下，约 10% 的 p 值小于 0.1，而在备择假设下， p 值小于 0.1 的比例要大得多。
[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