逆概率加权（R和Python案例）

news2024/11/20 0:27:00

逆概率加权（Inverse Probability Weighting, IPW）是一种统计技术，用于观察性研究中调整混杂变量的影响，以便更准确地估计因果关系。这种方法特别有用于在无法进行随机化实验的情况下，通过给予不同个体不同的权重，从而使得样本在处理组和对照组之间更加平衡，达到类似于随机对照试验的效果。

基本原理

逆概率加权的核心思想是通过给予每个个体一个权重，这个权重是其接受治疗或暴露的概率的倒数。具体来说：

对于接受治疗的个体，权重是 $\frac{1}{\text{Pr}(T=1 | X)}$ ；
对于未接受治疗的个体，权重是 $\frac{1}{\text{Pr}(T=0 | X)}$ 。

其中， $T$ 表示治疗或暴露的指派， $X$ 是可能影响 $T$ 分配的协变量。

详细公式

逆概率加权的公式可以表示为：
$\text{权重} = \left( \frac{T}{\text{Pr}(T=1 | X)} + \frac{1-T}{\text{Pr}(T=0 | X)} \right)$
这里， $T$ 是治疗指示变量（接受治疗为1，未接受为0），而 $\text{Pr}(T | X)$ 是给定协变量 $X$ 下接受治疗 $T$ 的概率。

案例

假设我们研究一种新药对疾病恢复的影响。我们有以下数据：

$T$ ：治疗指示（1=接受新药，0=接受常规治疗）
$X$ ：患者的年龄、性别等协变量
$Y$ ：恢复情况（1=完全恢复，0=未完全恢复）

首先，我们需要使用逻辑回归或其他适当的统计方法来估计 $\text{Pr}(T | X)$ 。然后，我们计算每个患者的权重，并使用这些权重来估计新药的因果效应。

实现

在实际应用中，逆概率加权可以通过以下步骤实现：

使用逻辑回归或其他模型估计 $\text{Pr}(T | X)$ 。
计算每个个体的权重。
使用加权的观测数据进行因果效应分析。

例如，在R语言中，可以使用如下代码进行逆概率加权的计算和分析：

# 模拟数据
n <- 1000
simdat <- data.frame(l = rnorm(n, 10, 5))
a.lin <- simdat$l - 10
pa <- exp(a.lin)/(1 + exp(a.lin))
simdat$a <- rbinom(n, 1, prob = pa)
simdat$y <- 10*simdat$a + 0.5*simdat$l + rnorm(n, -10, 5)

# 计算逆概率权重
library(ipw)
temp <- ipwpoint(exposure = a, family = "binomial", link = "logit", numerator = ~ 1, denominator = ~ l, data = simdat)

# 应用权重进行分析
simdat$sw <- temp$ipw.weights
msm <- svyglm(y ~ a, design = svydesign(~ 1, weights = ~ sw, data = simdat))
coef(msm)
confint(msm)

在Python中，可以使用statsmodels库来实现逆概率加权（IPW）。以下是一个简单的示例，展示了如何使用Python进行IPW的计算和分析：

import numpy as np
import pandas as pd
from statsmodels.miscmodels.ordinal_model import OrderedModel
from statsmodels.genmod.generalized_linear_model import GLM
from statsmodels.genmod.families import Binomial

# 模拟数据
np.random.seed(10)
n = 1000
data = pd.DataFrame({
    'age': np.random.normal(10, 5, n),
    'treatment': np.random.binomial(1, 0.5, n),
    'recovery': np.random.binomial(1, 0.5, n)
})

# 估计治疗分配概率模型
glm = GLM(data['treatment'], data[['age']], family=Binomial())
result = glm.fit()
predictions = result.predict(data[['age']])

# 计算逆概率权重
data['weight'] = data['treatment'] / predictions + (1 - data['treatment']) / (1 - predictions)

# 使用权重估计治疗效果
weighted_glm = GLM(data['recovery'] * data['treatment'], data[['age']], family=Binomial(), weights=data['weight'])
weighted_result = weighted_glm.fit()
print(weighted_result.summary())

# 使用权重估计控制组效果
weighted_glm_control = GLM(data['recovery'] * (1 - data['treatment']), data[['age']], family=Binomial(), weights=data['weight'])
weighted_result_control = weighted_glm_control.fit()
print(weighted_result_control.summary())

# 计算平均治疗效果 (ATE)
ate = (weighted_result.params[0] - weighted_result_control.params[0])
print(f'Estimated Average Treatment Effect (ATE): {ate}')