网站创建软件,织梦如何做二级网站,华为seo诊断及优化分析,传奇霸业官网下载这里写目录标题 加权标准化生存分析总结个体层面的生存曲线 加权标准化生存分析
我们还可以将加权与标准化结合起来#xff0c;使用 WeightedStandardizedSurvival 模块。在这里#xff0c;我们将逆倾向得分加权模型#xff08;根据基线协变量重新加权人群#xff09;与加… 这里写目录标题 加权标准化生存分析总结个体层面的生存曲线 加权标准化生存分析
我们还可以将加权与标准化结合起来使用 WeightedStandardizedSurvival 模块。在这里我们将逆倾向得分加权模型根据基线协变量重新加权人群与加权回归以及标准化模型相结合
from causallib.survival.weighted_standardized_survival import WeightedStandardizedSurvivalipw IPW(learnerLogisticRegression(max_iter2000))
poly_transform_pipeline Pipeline([(transform, PolynomialFeatures(degree2)), (LR, LogisticRegression(max_iter8000, C1.5))]
)
weighted_standardized_survival WeightedStandardizedSurvival(survival_modelpoly_transform_pipeline, weight_modelipw
)
weighted_standardized_survival.fit(X, a, t, y)population_averaged_survival_curves weighted_standardized_survival.estimate_population_outcome(X, a, t
)plot_survival_curves(population_averaged_survival_curves,labels[non-quitters, quitters],titleWeighted standardized survival of smoke quitters vs. non-quitters in a 10 years observation period,
)或者我们也可以使用 lifelines 包中的 RegressionFitter 类例如 Cox 比例风险拟合器。这是一种加权的 Cox 分析。
ipw IPW(learnerLogisticRegression(max_iter1000))
weighted_standardized_survival WeightedStandardizedSurvival(survival_modellifelines.CoxPHFitter(), weight_modelipw)# Note the fit_kwargs (passed to CoxPHFitter.fit() method)
weighted_standardized_survival.fit(X, a, t, y, fit_kwargs{robust: True})# Without setting robustTrue, well get the following warning:
StatisticalWarning: It appears your weights are not integers, possibly propensity or sampling scores then?
Its important to know that the naive variance estimates of the coefficients are biased. Instead a) set robustTrue in the call to fit, or b) use Monte Carlo to
estimate the variances.population_averaged_survival_curves weighted_standardized_survival.estimate_population_outcome(X, a, t)plot_survival_curves(population_averaged_survival_curves, labels[non-quitters, quitters], titleWeighted standardized survival of smoke quitters vs. non-quitters in a 10 years observation period)总结
不同模型的并列比较。
import itertoolsdef plot_multiple_models(models_dict):grid_dims (int(np.round(np.sqrt(len(models_dict)))), int(np.ceil(np.sqrt(len(models_dict)))))grid_indices itertools.product(range(grid_dims[0]), range(grid_dims[1]))fig, ax plt.subplots(*grid_dims)models_names list(models_dict.keys())for model_name, plot_idx in zip(models_names, grid_indices):model models_dict[model_name]model.fit(X, a, t, y)curves model.estimate_population_outcome(X, a, t, y)ax[plot_idx].plot(curves[0])ax[plot_idx].plot(curves[1])ax[plot_idx].set_title(model_name)ax[plot_idx].set_ylim(0.7, 1.02)ax[plot_idx].grid()plt.tight_layout()plt.show()
MODELS_DICT {MarginalSurvival Kaplan-Meier: MarginalSurvival(survival_modelNone),MarginalSurvival LogisticRegression: MarginalSurvival(survival_modelLogisticRegression(max_iter2000)),MarginalSurvival PiecewiseExponential: MarginalSurvival(survival_modellifelines.PiecewiseExponentialFitter(breakpointsrange(1, 120, 10))),WeightedSurvival Kaplan-Meier: WeightedSurvival(weight_modelIPW(LogisticRegression(max_iter2000)), survival_modelNone),WeightedSurvival LogisticRegression: WeightedSurvival(weight_modelIPW(LogisticRegression(max_iter2000)),survival_modelLogisticRegression(max_iter2000),),WeightedSurvival WeibullFitter: WeightedSurvival(weight_modelIPW(LogisticRegression(max_iter2000)),survival_modellifelines.WeibullFitter(),),StandardizedSurvival LogisticRegression: StandardizedSurvival(survival_modelLogisticRegression(max_iter2000)),StandardizedSurvival Cox: StandardizedSurvival(survival_modellifelines.CoxPHFitter()),WeightedStandardizedSurvival: WeightedStandardizedSurvival(weight_modelIPW(LogisticRegression(max_iter2000)),survival_modelLogisticRegression(max_iter2000),),
}plot_multiple_models(MODELS_DICT)个体层面的生存曲线
在使用直接结果模型StandardizedSurvival 和 WeightedStandardizedSurvival时可以在 causallib 中生成个体层面的效果估计和生存曲线。
%matplotlib inline
import matplotlib as mpl
import seaborn.objects as so
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from causallib.survival import StandardizedSurvival
from causallib.datasets import load_nhefs_survival
data load_nhefs_survival(augmentFalse, onehotFalse)
data.t data.t.rename(longevity)
data.X.join(data.a).join(data.t).join(data.y)现在让我们创建一个基于公式的数据转换器以便轻松指定以下两点
使用样条灵活地建模连续变量创建与所有变量的治疗交互项以允许效应修正。
from formulaic import Formula
from sklearn.base import BaseEstimator, TransformerMixinclass FormulaTransformer(BaseEstimator, TransformerMixin):def __init__(self, formula):super().__init__()self.formula formuladef fit(self, X, yNone):return selfdef transform(self, X, yNone):X_ Formula(self.formula).get_model_matrix(X)return X_
formula f~ 1 {data.a.name}*(C(exercise) C(active) C(education) sex race bs(age, degree5) bs(smokeintensity) bs(smokeyrs) bs(wt71) bs({data.t.name}, degree5) )estimator make_pipeline(FormulaTransformer(formula),LogisticRegression(penaltynone, max_iter1000)
)model StandardizedSurvival(estimator,stratifyFalse,
)
model.fit(data.X, data.a, data.t, data.y)
po model.estimate_individual_outcome(data.X, data.a, data.t)
po遵循 lifelines 的惯例结果的维度将不同的时间点作为行个体作为列。 列进一步按照治疗分配索引因为这些值是潜在结果。 这种结构使我们能够像在非生存分析中那样获得个体层面的效果生存差异
effect po[1] - po[0]
# effect我们现在将结果转置使其变为长格式以便后续绘图
effect effect.reset_index(namestime).melt(id_varstime, var_nameid, value_nameeffect)
effectf mpl.figure.Figure()# Plot inidividual lines:
p so.Plot(effect,xtime,yeffect,groupid,
).add(so.Lines(linewidth.5, alpha0.1, color#919090)
).label(titleSpaghetti plot of the effect difference,
).on(f).plot()# Plot average effect:
avg_effect effect.groupby(time)[effect].mean().reset_index()
ax f.axes[0]
ax.plot(avg_effect[time], avg_effect[effect], color#062f80)
ax.text(0, 0, ATE,verticalalignmentbottom,color#062f80
)
f一旦我们得到了个体级别的生存曲线我们可以任意聚合它们来观察效应在不同的协变量分层中是如何变化的。
f mpl.figure.Figure()
effectX effect.merge(data.X, left_onid, right_indexTrue)
strata racep_eff_strat so.Plot(effectX,xtime,yeffect,colorstrata, # Stratify the effect curves bygroupid,
).add(so.Lines(linewidth.5, alpha0.1)
).scale(colorso.Nominal([#1f77b4, #ff7f0e]),
).label(titleSpaghetti plot for stratified effects,
).on(f).plot()
p_eff_stratavg_effect effectX.groupby([time, strata])[effect].mean().reset_index()
ax f.axes[0]
for s, stratum_data in avg_effect.groupby(strata):ax.plot(stratum_data[time], stratum_data[effect], colorblack, linestyle--,)ax.text(stratum_data[time].iloc[-1], stratum_data[effect].iloc[-1],f{strata}:{s},verticalalignmentcenter,)f
