前言
纵向数据分析有多种分析方法，之前讲过了混合线性效应模型，这次讨论Landmark analysis（LM），后续讨论joint model方法。学无止境，但作为应用派，我们要知道什么时候用什么方法。
LM方法的应用
LM方法有多种应用场景
- 癌症研究：在肿瘤研究中，研究人员可能会使用landmark analysis来评估特定治疗方案或药物对患者生存时间的影响。通过分析患者在特定时间点的存活率，研究人员可以更好地了解治疗效果和预后。
- 临床试验设计：在临床试验设计中，landmark analysis可以用来评估特定时间点上治疗干预的效果。研究人员可以使用这种方法来确定治疗在不同时间点的效果，并为临床试验的结果提供更深入的理解。
- 慢性病进展监测：在慢性病研究中，如心血管疾病或糖尿病，研究人员可以使用landmark analysis来监测患者在特定时间点上的疾病进展情况。这有助于确定治疗干预是否能够延缓病情恶化或降低并发症风险。
我们来看具体的例子：两个消化顶刊的例子
文中就运用了LM法来分析肝纤维化指标的动态变化与脂肪肝中肝硬化及肝癌发生风险的关系。
文中采用了LM法来分析转氨酶动态变化对自身免疫性肝炎病人的生存影响。
在更多的肿瘤学研究中，LM法经常作为比较两种或多种治疗手段差异的重要办法，此处不再举例。
LM和混合线性效应模型的区别
landmark analysis和混合线性效应模型在纵向数据分析中有着不同的应用和作用。
-
Landmark Analysis（里程碑分析）：Landmark analysis是一种统计分析方法，用于在研究中对特定的时间点进行分析。它通常用于处理生存分析或时间至事件数据，以评估特定事件发生前后的效应。在医学研究中，例如在临床试验或观察性研究中，研究者可能会在患者达到特定时间点时进行里程碑分析，以了解治疗效果或预后的变化。
-
混合线性效应模型：混合线性效应模型是一种用于处理纵向数据的统计模型，它可以考虑数据中个体间或观测间的相关性，并允许引入随机效应以捕捉个体差异。这种模型通常用于探究数据中的随时间变化的趋势、处理重复测量数据，并进行预测和推断。
总的来说，landmark analysis侧重于在特定时间点上进行分析，而混合线性效应模型更侧重于整体纵向数据的模式和变化趋势。在某些情况下，这两种方法可能会结合使用，以全面理解纵向数据的动态变化和效应。
实战LM
此处我们借用一个新的R包，dynamicLM包来进行实战，这个包简单易用，容易上手，作者配备了齐全的例子分析。 并且作者发表了这个包的正式论文，大家可以引用。
示例数据
基本信息
包中自带了示例数据，该数据集包含了3,332名曾吸烟的患者，他们是多种族队列研究（Multiethnic Cohort Study，MEC）的一部分，这些患者在1993年至2007年之间被诊断为初发原发性肺癌（IPLC），并在此后至2017年进行了随访以观察是否发生了次发性肺癌（SPLC）。
通过与SEER登记处的关联，确定了SPLC病例。未在最后的随访日期（2017年12月31日）发生任何事件的患者将被删除。根据Martini和Melamed的临床标准来确定了SPLC的发生：如果随后的肿瘤是在IPLC诊断后两年内诊断的，或者具有不同的组织学类型。就在队列入组时（1993-1996年）收集了全面的基线信息，通过自填问卷获得；更新后的与吸烟有关的变量在10年随访调查中得到（2003-2008年），占研究队列的21.1%（n=704）。如果10年随访信息是在IPLC诊断之前或同时获得的，我们将使用更新后的数据作为基线信息；如果随访数据在IPLC诊断后收集的，我们将使用基线和随访数值，以捕捉IPLC诊断后随时间变化的协变量变化。
研究终点是从IPLC诊断到SPLC发生的时间（cause 1）；在发生SPLC之前可以出现的两种竞争事件包括肺癌死亡（cause 2）和其他原因死亡（cause 3）
具体指标
dynamicLM包需要设置一个super dataset
super dataset需要：（1）the outcome columns（终点事件）（2）variable types (fixed vs. time-varying) （变量类型（固定 vs. 随时间变化））。
固定和随时间变化的变量都被视为模型开发的候选预测因子。
静态（固定）变量包括：性别（“male”）、IPLC诊断前的癌症病史（“ph”）、IPLC分期（“stage.ix”）、IPLC组织学（“hist_AD”、“hist_LC”、“hist_NSCLC_NOS”、“hist_SC”、“hist_OTH”）、首次治疗（“surgery.ix”、“radiation.ix”、“chemo.ix”）。
随时间变化的协变量——即在基线（即IPLC诊断）后随时间可能发生变化的状态或数值——包括以下内容：吸烟状态（“smkstatus”）、吸烟包年数（“packyears”）、每日吸烟量（“cigday”），以及在10年随访调查中收集的戒烟年份，并可按年度自动或手动更新。
表格具体的含义如下：
代码实战
1.加载R包，设置变量
devtools::install_github("thehanlab/dynamicLM")
library(dynamicLM)
data(splc) ##加载示例数据
outcome <- list(time = "Time", status = "event") ##结局变量
fix_covs <- c("age.ix", "male", "fh", "ph", "bmi", "stage.ix",
"hist_AD", "hist_LC", "hist_NSCLC_NOS", "hist_SC",
"hist_OTH", "surgery.ix", "radiation.ix",
"chemo.ix", "quityears") ##固定变量
vary_covs <- c("smkstatus", "cigday", "packyears") ##时依协变量
covs <- list(fixed = fix_covs, varying = vary_covs) 2.设置landmark，创建dataset
要在IPLC诊断后0至3年的任何评估时间点上预测SPLC的5年风险，我们必须设置里程碑时间点、预测窗口和数据格式。基于数据可用性以及基线暴露可能影响感兴趣结果的合理时间窗口，可以确定里程碑和预测时间窗口（“w”）。我们设置了从IPLC诊断后的0至3年（0、1、2和3年）每年一个里程碑（“lms”）。
数据格式在“format”参数中进行了指定。原始数据集（splc）采用长格式，如果患者有10年随访调查数据，则可能包括患者的多次观测。列“T.fup”表示10年随访时间，即基线和10年随访测量之间的持续时间，包括随时间变化的协变量（“smkstatus”、“cigday”和“packyears”）。这个10年随访时间在“rtime”参数中进行了指定。
从在给定里程碑创建的里程碑数据集生成一个堆叠数据集(stack data)（“lmdata”）。时间变化的协变量的值在数据集中进行了更新（如果可用）。堆叠数据的列“LM”表示每个LM。
> w <- 5 # 5-year prediction window
> lms <- seq(0, 3, by = 1)# landmarks
> lmdata <- stack_data(data = splc, outcome = outcome, lms = lms, w = w,
covs = covs, format = "long", rtime = "T.fup", id = "ID")
> table(lmdata$data$LM) # number of patients per landmark dataset (Corresponding to Supp. Fig 1A, but please keep in mind that the example dataset (splc) in this R package is synthetic because the MEC data is only available under a data use agreement. Readers can apply the code but cannot replicate the same estimates.)
在一些时变协变量中，随着时间的推移可以自动更新。在我们的数据中，曾吸烟者的戒烟年限呈线性增长。 也即随着LM的增加而增加，LM加1，戒烟相当于多1年，所以需要手动设置一下。
former_smokers <- lmdata$data$smkstatus == 2 ##选出曾吸烟者
##如果是曾吸烟者则在quityears基础上加上LM设置的年限，生成新的quityears
lmdata$data[former_smokers, "quityears"] <-
lmdata$data[former_smokers, "quityears"] +
lmdata$data[former_smokers, "LM"]
3.创建交互项以捕捉时效作用，并整合平滑效果
啥叫时效作用，最简单的做法，变量*时间，就是变量和时间的交互作用项。
在这一步中有两个步骤：1）创建选定协变量（“covars”）与LM time（LM time可以是线性、二次方或其他形式，使用func_covars参数设置）之间的交互项，以检查协变量的时效影响；2）线性和二次标志时间变量（“func_LMs”）在模型拟合中整合平滑效果。
与时变协变量不同，时效作用是指某个变量对特定时间内IPLC风险的影响在时间上发生变化。例如，据报道，用于治疗IPLC的放疗的风险随时间增加。这一步在捕捉协变量对结果的真实影响以及检验比例风险假设方面至关重要。
如果协变量与LM time之间存在显著的交互项，则违反了比例风险假设（即风险比随标志时间变化）。因此，有必要同时包括主要项和交互项（也即前面提到的选择变量及其与时间的交互项）。
在当前SPLC模型开发中，基于文献，我们提前选择了一组变量（IPLC放疗、IPLC分期、IPLC组织学、先前的癌症病史和吸烟包年数）以检查时效影响（但读者可能希望测试所有可用的协变量）；基于专业知识我们选定了LM time的两种转换形式（线性和二次）。通常情况下，一个变量的时效影响可以以线性或二次方形式表示，但如果必要还可以添加其他转换形式（如对数）。
新创建的交互项由变量名称后跟下划线和数字来识别（如：stage.ix_1 = stage.ix x LM，以及stage.ix_2 = stage.ix x LM2）。
> # A priori variables to be examined for time-dependent effect
> interaction_covs <- c("packyears", "radiation.ix", "ph", "stage.ix",
"hist_AD", "hist_LC", "hist_NSCLC_NOS", "hist_SC", "hist_OTH")
> lmdata <- add_interactions(lmdata = lmdata, lm_covs = interaction_covs,
func_covars = c("linear", "quadratic"),
func_lms = c("linear", "quadratic"))
> colnames(lmdata$data) # Check the newly created variables
4.识别对SPLC风险具有时效影响的变量
在创建了交互项之后，我们要确定具有显著时效影响的协变量。
为了检测特定变量（例如“吸烟包年数”）的时效影响，我们拟合一个包含主协变量（“吸烟包年数”）和交互项（“吸烟包年数_1”，“吸烟包年数_2”对应LM和LM的2次方）的模型。
如下面的示例代码所示。
如果添加交互项使主协变量失去显著性，那么我们认为该协变量在标志时间上没有时效影响，符合比例风险假设。
如果在添加交互项后主协变量仍然显著，我们会检验交互项的显著性，仅保留显著的交互项。
如果没有显著的交互效应，那么我们也不认为该协变量有时效影响。
这种情况下，我们在模型选择过程中只使用主协变量（例如“吸烟包年数”）。如果主协变量和交互项都显著，我们认为该协变量在标志时间上有时效影响，并在模型选择中同时包括主协变量和交互项（例如“放疗.ix”和“放疗.ix_1”）。
为了更全面地评估，读者还可以考虑添加交互项如何改变模型的性能（例如校准）。
> formula <- "Hist(Time, event, LM) ~ packyears + packyears_1 + packyears_2 + cluster(ID)"
> check_td <- dynamic_lm(lmdata, as.formula(formula), "CSC")
> print(check_td, cause = 1)
> # packyears showed a significant effect, but packyears_1 and packyear_2 did not.
> formula <- "Hist(Time, event, LM) ~ packyears + packyears_1 + cluster(ID)"
> check_td <- dynamic_lm(lmdata, as.formula(formula), "CSC")
> print(check_td, cause = 1)
> # packyears showed a significant effect, but packyears_1 did not.
> # packyears does not have a time-dependent effect.
> formula <- "Hist(Time, event, LM) ~ radiation.ix + radiation.ix_1 + radiation.ix_2 + cluster(ID)"
> check_td <- dynamic_lm(lmdata, as.formula(formula), "CSC")
> print(check_td, cause = 1)
> # radiation.ix showed a significant effect, but radiation.ix_1 and radiation.ix_2 did not.
> formula <- "Hist(Time, event, LM) ~ radiation.ix + radiation.ix_1 + cluster(ID)"
> check_td <- dynamic_lm(lmdata, as.formula(formula), "CSC")
> print(check_td, cause = 1)
> # Both radiation.ix and radiation.ix_1 showed a significant effect.
> # radiation.ix has a time-dependent effect.
> # include both radiation.ix and radiation.ix_1 together in model development.
5.变量选择及建模
批量单因素LM分析
选择p<0.1的变量
> predictors <- c(fix_covs, vary_covs)
> for (i in predictors) {
formula <- paste("Hist(Time, event, LM) ~ ", i[[1]], "+ cluster(ID)")
supermodel_uni <- dynamic_lm(lmdata, as.formula(formula), "CSC")
print(supermodel_uni, cause = 1)
}
多因素LM分析
在检验了单变量模型之后，我们拟合一个全模型，包括所有表现出P值小于0.1的特征，以及文献中已知与SPLC风险相关的（或者是潜在混杂因素）附加特征。对于确定具有时效影响的协变量，我们添加了交互项。此外，我们还包括线性和二次的标志时间（LM_1和LM_2）变量，以整合在stack data上的平滑效果；这使得拟合的超模型的基线风险能依赖于LM time。
> formula <- "Hist(Time, event, LM) ~ age.ix + bmi + male + ph + ph_1 +
surgery.ix + cigday + packyears + quityears +
hist_AD + hist_LC + hist_NSCLC_NOS + hist_OTH + hist_SC +
hist_AD_1 + hist_LC_1 + hist_NSCLC_NOS_1 + hist_OTH_1 +
hist_SC_1 + radiation.ix + radiation.ix_1 + stage.ix +
stage.ix_1 + stage.ix_2 + LM_1 + LM_2 + cluster(ID)"
> supermodel_full <- dynamic_lm(lmdata, as.formula(formula), "CSC")
> print(supermodel_full)
在进行变量选择时，我们从全模型开始，逐个剔除一个变量，通过评估估计的风险比的解释性、预测性能指标（区分度、校准和布里尔分数）以及统计显著性来进行评估。下面是所选的最终模型。在附表S2中显示了所选模型的动态特异原因危险比。具有时效影响的协变量包括IPLC放疗，其效应随标志时间增加（附图S2）（见原文附件）。
> formula <- "Hist(Time, event, LM)~ male + ph + ph_1 +
packyears + quityears + hist_AD + hist_LC + hist_NSCLC_NOS +
hist_OTH + hist_SC + hist_AD_1 + hist_LC_1 + hist_NSCLC_NOS_1 +
hist_OTH_1 + hist_SC_1 + radiation.ix + radiation.ix_1 +
stage.ix + stage.ix_1 + stage.ix_2 + LM_1 + LM_2 + cluster(ID)"
> supermodel <- dynamic_lm(lmdata, as.formula(formula), "CSC", x = TRUE)
> print(supermodel)
作图
可以绘制所选模型中包含的变量的动态风险比。
par(mfrow = c(2, 3))
plot(supermodel,
logHR = FALSE,
covars = c("ph", "stage.ix", "hist_LC", "hist_SC", "radiation.ix"),
main = c("Prior history of cancer", "IPLC stage (Advanced)",
"IPLC histology (LC)", "IPLC histology (SC)",
"IPLC radiotherapy"),
xlab = "Risk Assessment Time \n (Landmark Prediction Time)",
ylab = "Adjusted csHR")
基于所选模型，可以使用以下方法来估计个体肺癌患者发展SPLC的5年预测风险。
p <- predict(supermodel)
predicted_risk <- p$preds
summary(predicted_risk$risk)
下图显示了两位患者的动态5年SPLC风险的示例。
par(mfrow = c(1, 1))
inds <- lmdata$data[lmdata$data$ID %in% c("ID2", "ID7"), ]
plotrisk(supermodel, inds, format = "long", ylim = c(0, 0.2))
在不同的风险评估时间点（标志时间）上，显示了两位示例个体预测的动态5年SPLC风险轨迹。在一个预定义的5年预测窗口（w=5）内，在每个标志时间点（即在水平坐标上的风险评估时间或预测时间点），都估计了在5年内患SPLC的风险。在本例中，患者P1在IPLC诊断时（标志=0）的5年SPLC发展风险为0.97％，在IPLC诊断后1年（标志=1）预测的风险更新为1.43％，在IPLC诊断后2年（标志=2）预测的风险为3.06％，在IPLC诊断后3年（标志=3）预测的风险为4.38％。这些都是基于截至每个标志时间的最新患者信息来预测的。
6.模型评价
为了评估预测模型的性能，我们创建了校准图并计算了AUC和Brier分数。总体而言，所提出的模型在不同的标志时间上表现良好（校准斜率范围为[0.9-1.2]）
par(mfrow = c(2, 2))
cal <- calplot(list("CSC" = p), method = "quantile", q = 10,
ylim = c(0, 0.25), xlim = c(0, 0.25))
scores <- score(list("CSC" = p), cause = 1)
print(scores)
par(mfrow = c(1, 2))
plot(scores)
基准时刻的AUC大于85％，随着标志时间的推移而下降，因为用于模型拟合的大多数可用数据来自MEC的基线。如果随着时间的推移，患者特征得到更新和补充，预测精度可能会保持一致。
7.模型的验证和泛化性
Bootstrapped内部验证
par(mfrow = c(2, 2))
cal <- calplot(list("CSC" = supermodel), method = "quantile", q = 10,
ylim = c(0, 0.25), xlim = c(0, 0.25),
split.method = "bootcv", B = 10)
score(list("CSC" = supermodel), cause = 1, split.method = "bootcv", B = 10)
外部验证
# New stacked data with a LM column
data(splc_test)
table(splc_test$LM)
# lms argument gives the landmark column
par(mfrow = c(2,2))
cal <- calplot(list("CSC" = supermodel), data = splc_test, lms = "LM",
method = "quantile", q = 10,
ylim = c(0, 0.25), xlim = c(0, 0.25))
scores_external <- score(list("CSC" = supermodel), cause = 1,
data = splc_test, lms = "LM")
print(scores_external)
258
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
