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
虛無設計矩陣嵌套在對立設計矩陣中。它的秩比對立設計矩陣少 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，而在對立假設下平均值低得多。類似地，我們預期在虛無假設下，約 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