JMbayes2包学习笔记——Day1

学习途径：JMbayes2 • JMbayes2 (drizopoulos.github.io)     动态预测模型最新R包-JMbayes2包教程(1)-get started

在JMbayes2中拟合联合模型的过程可以概括为以下几个步骤：

拟合纵向结果的混合效应模型：

混合效应模型是一种考虑随机效应和固定效应的统计模型。在JMbayes2中，可以使用函数nlme()或lme()（来自GLMM自适应包）来拟合混合效应模型。

拟合事件过程的Cox或加速失效时间（AFT）模型：

Cox模型和AFT模型都是用于分析生存数据的方法。在JMbayes2中，可以使用函数coxph()或survreg()（来自生存包）来拟合Cox模型或AFT模型。

 拟合联合模型：

在完成上述两个步骤后，可以将生成的模型对象作为参数传递给jm()函数，以拟合联合模型。

案例（pbc数据集）

pbc数据集，描述了关于受肝原发性胆汁性肝硬化影响的人的生存数据点 

id：病人的唯一标识。

years：可能是病人与某种疾病或状况共存的年数。

status：病人的健康状况，如“alive”（存活），“dead”（已故）“transplant”（做了移植手术）。

drug：病人是否使用了特定的药物，这里是“D-penicil”或“placebo”。

age：病人的年龄。

sex：病人的性别。

year：可能是病人的确诊年份或其他相关年份。

ascites：腹水，一种由于腹内体液过多而引起的病症。

hepatomegaly：肝肿大，肝脏体积增大。

spiders：可能是蜘蛛痣，一种由于肝脏疾病引起的皮肤上的小红点。

edema：水肿，身体的某部分异常肿胀。

serBilir：血清胆红素，一种衡量血液中胆红素水平的指标。

serChol：血清胆固醇，一种衡量血液中胆固醇水平的指标。

albumin：白蛋白，一种重要的血液蛋白质。

alkaline：可能是指碱性磷酸酶，一种在骨骼和肝脏中产生的酶。

SGOT：天冬氨酸氨基转移酶（也称为AST），一种肝脏酶。

platelets：血小板，血液中的一种成分，有助于血液凝固。

prothrombin：凝血酶原，一种帮助血液凝固的蛋白质。

histologic status2：可能是关于组织或细胞样本的病理学状态的描述。

从数据中可以看出，这些病人大多数是女性，并且有一些与肝病相关的症状和实验室异常。其中一些病人接受了药物（如“D-penicil”）治疗，而另一些病人则接受了移植手术。

pbc2.id是数据框,status是原始因变量,表示病情("alive","died","transplant")，status2用作新的二值因变量,表示生存情况，status != 'alive' 判断原值不是"alive"，只要不是“alive”，该病人在status2里面的状态均为1。

首先，我们拟合一个考虑性别作为基线协变量的复合事件移植或死亡的Cox模型：

代码：

CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)

coxph()：这是R中用于拟合Cox比例风险模型的函数。

Surv(years, status2)：这是生存对象，其中years是观察时间，status2是事件状态。通常，事件状态是一个二元或多元变量，表示观察是否以事件（如死亡）结束。

~ sex：这表示我们想要研究sex变量对事件发生时间的影响。

data = pbc2.id：这指定了包含模型所需数据的数据框pbc2.id。

运行结果：
Call:
coxph(formula = Surv(years, status2) ~ sex, data = pbc2.id)

  n= 312, number of events= 140 

             coef exp(coef) se(coef)      z Pr(>|z|)   
sexfemale -0.6502    0.5219   0.2183 -2.979   0.0029 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
sexfemale    0.5219      1.916    0.3403    0.8006

Concordance= 0.537  (se = 0.016 )
Likelihood ratio test= 7.73  on 1 df,   p=0.005
Wald test            = 8.87  on 1 df,   p=0.003
Score (logrank) test = 9.18  on 1 df,   p=0.002

Call: 这是运行coxph()函数的命令，它定义了模型。

n= 312, number of events= 140: 这表示你有312个观察值和140个事件（通常是死亡）。

coef: 这是模型估计的系数。在这个例子中，sexfemale的系数是-0.6502。

exp(coef): 这是系数的指数，即exp(-0.6502)，它给出了一个变量变化对风险比的影响。在这个例子中，女性性别的风险比是0.5219，意味着女性的风险较低。

se(coef): 这是系数的标准误差。

z: 这是系数与0之间的z分数，用于测试系数的显著性。

Pr(>|z|): 这是z分数的p值，用于测试系数的显著性。在这个例子中，p值是0.0029，表示性别对事件发生时间有显著影响。

lower .95, upper .95: 这是系数的95%置信区间。

Concordance= 0.537: 这是模型的Cox-Snell决定系数，表示模型的整体拟合程度。这个值介于0和1之间，值越高表示模型拟合越好。

Likelihood ratio test= 7.73: 对数似然比检验的结果，用于测试模型的整体显著性。

Wald test= 8.87: Wald检验的结果，用于测试系数的显著性。

Score (logrank) test= 9.18: Logrank检验的结果，用于测试模型的整体显著性。

总的来说，这个模型表明性别对事件发生时间有显著影响（p<0.01），女性性别的风险较低（风险比为0.5219）。
 

 代码：

fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)

给出的代码是一个使用lme函数（来自nlme包）来拟合线性混合效应模型的R代码片段。
解释这个代码的每一部分：
lme: 这是函数的名字，它用于拟合线性混合效应模型。
log(serBilir) ~ year * sex: 这是模型的公式。
log(serBilir)：响应变量是serBilir的的对数。
year * sex：这是两个预测变量的交互项。也就是说，模型将考虑year和sex的交互效应。
data = pbc2: 这告诉函数数据集的名字是pbc2。
random = ~ year | id: 这定义了模型的随机效应部分。
year | id: 表示year是作为随机截距，而id是随机截距的类别变量。换句话说，对于每个不同的id值，我们都预期有一个不同的截距。
总结：这段代码的目的是使用线性混合效应模型来分析pbc2数据集中的数据，其中响应变量是serBilir的对数，主要解释变量是year和sex，并且考虑了这两个变量的交互效应。模型还考虑了随机效应，其中每个不同的id值都有一个不同的截距。

运行结果：

 

Linear mixed-effects model fit by REML
  Data: pbc2 
       AIC      BIC    logLik
  3080.439 3125.006 -1532.219

Random effects:
 Formula: ~year | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 0.9973978 (Intr)
year        0.1706167 0.395 
Residual    0.3487097       

Fixed effects:  log(serBilir) ~ year * sex 
                    Value  Std.Error   DF   t-value p-value
(Intercept)     0.7383347 0.17073919 1631  4.324342  0.0000
year            0.2398386 0.03618717 1631  6.627726  0.0000
sexfemale      -0.2737765 0.18152767  310 -1.508180  0.1325
year:sexfemale -0.0717076 0.03851019 1631 -1.862042  0.0628
 Correlation: 
               (Intr) year   sexfml
year            0.235              
sexfemale      -0.941 -0.221       
year:sexfemale -0.221 -0.940  0.235

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.31522853 -0.49860164 -0.01329973  0.45353616  5.29370408 

Number of Observations: 1945
Number of Groups: 312 

 

这是一个线性混合效应模型的输出，通过REML（限制极大似然）方法进行拟合。详细解释输出中的各个部分和数字所对应的统计学意义。

1. AIC, BIC, logLik:
   AIC和BIC是用于模型选择的准则，值越小表示模型拟合越好。logLik是对数似然值，描述了模型对数据的拟合好坏。
2. Random effects:
   这部分描述了模型的随机效应部分。
   Formula: ~year | id 表示随机效应由年份和id两个因素决定。
   Structure: General positive-definite, Log-Cholesky parametrization 描述了随机效应的参数化方式。
   StdDev & Corr: 分别给出了标准差和相关性。
3. Fixed effects:
   这部分描述了模型的固定效应部分。
   log(serBilir) ~ year * sex 表示响应变量log(serBilir)与年份和性别两个因素之间的关系。
   Value & Std.Error: 分别给出了固定效应的估计值和标准误差。
   DF: 自由度，表示用于估计固定效应的独立数据点数量。
   t-value & p-value: t值是t检验的统计量，p值是对应的p值，用于检验固定效应的显著性。
4. Correlation:
   这部分描述了固定效应部分的各个系数之间的相关性。
5. Standardized Within-Group Residuals:
   这部分给出了标准化后的残差分布情况。
   Min, Q1, Med, Q3, Max: 分别表示残差分布的最小值、第一四分位数、中位数、第三四分位数和最大值。
6. Number of Observations & Groups:
   Number of Observations: 观察值的数量，即数据点的数量。
   Number of Groups: 组的数量，即随机效应的类别数量。
   综上所述，这个线性混合效应模型通过REML方法进行了拟合，并给出了模型的详细输出结果。从输出结果来看，模型拟合良好，固定效应显著，残差分布合理。

多变量联合模型的搭建

代码：

 fForms <- list(
  "log(serBilir)" = ~ slope(log(serBilir)) + slope(log(serBilir)):sex,
  "prothrombin" = ~ area(prothrombin)
)
jointFit3 <- update(jointFit2, functional_forms = fForms)
summary(jointFit3) 

 jm()函数的默认设置是将混合模型个体专属线性预测变量纳入生存比例风险模型。jm()函数的functional_forms参数提供了额外的选项。

• 随时间变化的斜率:混合模型个体线性预测随时间变量的一阶导数；

• 随时间变化的正则化累积效应:混合模型个体线性预测从0到当前时间t的积分,除以t。积分是主体特定纵向轨迹下的面积；将积分除以t，我们得到相应时期（0，t）内主体特定纵向轨迹的平均值。

这段代码是用于更新之前创建的联合模型 `jointFit2`，主要包括添加了函数形式（functional forms）的定义。具体解释如下：

1. `fForms` 是一个列表，其中包含了两个函数形式的定义：
   - `"log(serBilir)" = ~ slope(log(serBilir)) + slope(log(serBilir)):sex`：表示 `log(serBilir)` 这个纵向变量的函数形式。`slope(log(serBilir))` 表示该变量的斜率，而 `slope(log(serBilir)):sex` 表示该斜率在不同性别之间的交互效应。
   - `"prothrombin" = ~ area(prothrombin)`：表示 `prothrombin` 这个纵向变量的函数形式，其中 `area(prothrombin)` 表示 `prothrombin` 的面积效应。

2. `update(jointFit2, functional_forms = fForms)`：使用 `update` 函数更新之前拟合的联合模型 `jointFit2`，通过指定 `functional_forms = fForms` 来添加新的函数形式。

3. `summary(jointFit3)`：对更新后的联合模型 `jointFit3` 进行总结和展示。

总体来说，这段代码的作用是更新联合模型，加入了对 `log(serBilir)` 和 `prothrombin` 的函数形式定义，以更全面地描述这两个变量在模型中的效应。函数形式的定义允许模型更加灵活地捕捉不同变量之间的关系和交互效应。

使用收缩先验惩罚系数

代码：

jointFit4 <- update(jointFit3, priors = list("penalty_alphas" = "horseshoe"))
cbind("un-penalized" = unlist(coef(jointFit3)), 
      "penalized" = unlist(coef(jointFit4))) 

这段代码执行了模型的更新，使用了 "horseshoe" 先验对模型进行惩罚。接着，通过 `cbind` 函数，将未经惩罚的模型 (`jointFit3`) 和经过惩罚的模型 (`jointFit4`) 的系数进行合并，以便比较两者的效果。

解释：

- **un-penalized:** 这列包含了 `jointFit3` 模型中的未惩罚的系数。这是在没有额外先验（惩罚）的情况下，由模型学到的系数。

- **penalized:** 这列包含了 `jointFit4` 模型中的经过 "horseshoe" 先验惩罚的系数。 "horseshoe" 先验是一种用于稀疏模型的先验分布，它可以推动一些系数变得非常接近于零，从而实现变量选择（模型中只有一部分变量对响应有贡献）。

通过比较这两列，你可以看到在使用 "horseshoe" 先验后，一些系数可能会被推向零，从而实现模型的稀疏性。这可以在变量选择问题中有用，因为它可以使模型更简单，更容易解释，同时提高泛化性能。

 运行结果：

                                           un-penalized   penalized
gammas.Mean                                  -0.2100947 -0.71280039
association.slope(log(serBilir))              4.1105892  2.67481801
association.slope(log(serBilir)):sexfemale   -2.2855265 -0.49996888
association.area(prothrombin)                 0.1354759  0.09110531
association.value(ascites)                    0.7797344  0.71585830

在给定的代码中，运行结果中的 "un-penalized" 和 "penalized" 对应了两种不同的模型配置，分别是未使用惩罚先验（`jointFit3`）和使用了惩罚先验（`jointFit4`）的情况。以下是对结果的解释：

1. **"gammas.Mean"：**
   - 这个参数可能是 Cox 模型中的性别（sex）的效应。在未惩罚的模型中，其估计为 -0.2100947，而在惩罚的模型中为 -0.71280039。这表明在引入惩罚先验的情况下，性别的效应变得更加负向。

2. **"association.slope(log(serBilir))" 和 "association.slope(log(serBilir)):sexfemale"：**
   - 这两个参数可能是纵向数据（log(serBilir)）的斜率（slope）和交互效应（interaction term）。在未惩罚的模型中，斜率的估计为 4.1105892，而在惩罚的模型中为 2.67481801。这表示引入惩罚先验后，模型对纵向数据的某些效应进行了压缩。
   - 交互效应 "association.slope(log(serBilir)):sexfemale" 在未惩罚的模型中为 -2.2855265，在惩罚的模型中为 -0.49996888。同样，引入惩罚先验后，模型对交互效应进行了压缩。

3. **"association.area(prothrombin)" 和 "association.value(ascites)"：**
   - 这两个参数可能涉及纵向数据中 prothrombin 和二元数据中 ascites 的效应。在未惩罚的模型中，这两者的估计为 0.1354759 和 0.7797344，而在惩罚的模型中为 0.09110531 和 0.71585830。这表明引入惩罚先验后，模型对这两个效应进行了压缩。

总体来说，"un-penalized" 和 "penalized" 结果的比较显示了在引入惩罚先验的情况下，模型参数的估计值相对于未惩罚模型更加稀疏。这有助于防止模型过拟合，提高模型的泛化能力。