mixed_cox.Rmd

---
title: "mixed_cox"
author: "liuc"
date: "10/28/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Mixed effects cox regression

> https://stats.oarc.ucla.edu/r/dae/mixed-effects-cox-regression/

Mixed effects cox regression models are used to model survival data when there
are repeated measures on an individual, individuals nested within some other
hierarchy, or some other reason to have both fixed and random effects.

*shared frailty model: *
Shared frailty model是一种有参生存分析模型，它考虑了群组内个体之间的相关性，常用于分析有群组结构的数据，如医院、家庭、社区等。在R中，可以使用coxme包实现shared frailty model。
R allows fitting a frailty model via coxph by adding a frailty() term to the model formula. There is a new and more general approach in Therneau’s coxme library, which includes the coxme() function to fit mixed Cox survival models with Gaussian random effects using a Laplace approximation. In this example the two approaches give very similar answers, but this is not always the case.

```{r}
library(tidyverse)
library(survival)
library(coxme)

```

```{r}
m <- coxme(Surv(time, status) ~ age + sex + (1|ph.ecog), 
           data = lung)
summary(m)
```




# Competing Risks Regression

竞争风险模型是指event有多种，一种发生则另外一种不发生。在竞争事件存在的情况下各个结局事件的累计发病函数是如何随着时间变化的，以及如何来比较不同的累计发病函数，以及如何探讨影响因素.

不得不说还是finalfit简单好用。

```{r}
library(cmprsk)
# library(ccrstep)
library(condSURV)
library(finalfit)
library(ggsurvfit)
```

```{r}

# The crr() syntax differs from survival::coxph() but finalfit brings these together
# https://finalfit.org/articles/survival.html

data(melanoma, package = 'boot')

melanoma <- melanoma %>%
  mutate(
    # Overall survival
    status_os = ifelse(status == 2, 0, # "still alive"
                                         1), # "died of melanoma" or "died of other causes"
    
    # Diease-specific survival
    status_dss = case_when(
        status == 2 ~ 0, # "still alive"
      status == 1 ~ 1, # "died of melanoma"
      TRUE ~  0),     # "died of other causes is censored"

    # Competing risks regression
    status_crr = case_when(
        status == 2 ~ 0, # "still alive"
      status == 1 ~ 1, # "died of melanoma"
      TRUE ~ 2),       # "died of other causes"
    
    # Label and recode other variables
    age = ff_label(age, "Age (years)"), # ff_label to make table friendly var labels
    thickness = ff_label(thickness, "Tumour thickness (mm)"), # ff_label to make table friendly var labels
    sex = factor(sex) %>% 
        fct_recode("Male" = "1", 
                             "Female" = "0") %>% 
        ff_label("Sex"),
    ulcer = factor(ulcer) %>% 
        fct_recode("No" = "0",
                             "Yes" = "1") %>% 
        ff_label("Ulcerated tumour")
  )

explanatory   <- c("age", "sex", "thickness", "ulcer")
dependent_dss <- "Surv(time, status_dss)"
dependent_crr <- "Surv(time, status_crr)"

melanoma %>%
    # Summary table
  summary_factorlist(dependent_dss, explanatory, 
                     column = TRUE, fit_id = TRUE) %>%
    # CPH univariable
      ff_merge(
    melanoma %>%
      coxphmulti(dependent_dss, explanatory) %>%
      fit2df(estimate_suffix = " (DSS CPH univariable)")
    ) %>%
    # CPH multivariable
  ff_merge(
    melanoma %>%
      coxphmulti(dependent_dss, explanatory) %>%
      fit2df(estimate_suffix = " (DSS CPH multivariable)")
    ) %>%
    # Fine and Gray competing risks regression
  ff_merge(
    melanoma %>%
      crrmulti(dependent_crr, explanatory) %>%
      fit2df(estimate_suffix = " (competing risks multivariable)")
    ) %>%
  select(-fit_id, -index) %>%
  dependent_label(melanoma, "Survival")


# 固定时间点的累计发病置信区间
CumIncidence
```

```{r}
# competing HM by cmprsk

attach(melanoma)
fit=cuminc (ftime = time, 
            fstatus = status_crr, sex, cencode = 0) 
plot(fit) # 各个分组的累计发病曲线
fit

detach()
```

```{r}
# 上面介绍的方法相当于没有竞争风险的时候的KM法，通过上面的方法我们知道有了不同风险结局累计发病曲线的差异，继续我们会继续看影响因素，要做的就是竞争风险回归模型了，需要用到的函数就是crr
cov1 = model.matrix(as.formula(paste0("~", paste0(explanatory, 
        collapse = "+"))), melanoma)[, -1]

mod_crr <- crr(
  ftime = time,
  fstatus = status_crr,
  cov1 = cov1
)

summary(mod_crr)
```

```{r}
# equal to above 
melanoma %>%
      crrmulti(dependent_crr, explanatory) %>%
      fit2df(estimate_suffix = "(competing risks multivariable)")
```

### ROC曲线绘制, 列线图，校准曲线
```{r}
library(riskRegression)
library(pec)
library(timeROC)
library(regplot)

fgr1<-FGR(Hist(time,status_crr)~age+thickness+sex+ulcer,
          data = melanoma,cause = 1)
summary(fgr1) #绘制竞争风险模型
```
绘制某指标预测复发的time-dependent-ROC曲线
```{r}
tROC<-timeROC(T=melanoma$time,
              delta = melanoma$status_crr,
              weighting = "aalen",
              marker = melanoma$age,cause = 1,times =c(185,204,232))
tROC

plotAUCcurve(tROC,col="darkcyan")


## 计算C-index(复发，6、12、18...等月复发的C-index)
cindex(fgr,data = melanoma,cause = 1,eval.times = seq(6,131,6))


calPlot(fgr1,time = 131,xlim = c(0,0.8),ylim = c(0,0.8))

# 绘制交互列线图(查看某位患者的生存率)
library(mstate)
dat.w<-crprep("ftime","Status",
              data = melanoma,
              trans = c(1,2),cens = 0,id="ID",
              keep = c("Age","Gender","Diagnosis","phase_cr","source"))
dat.w$Time<-dat.w$Tstop-dat.w$Tstart

crm<-coxph(Surv(Time,status==1)~Age+Gender+Diagnosis+phase_cr+source,
           data = dat.w[dat.w$failcode==1,],weights = weight.cens,subset = failcode==1)
summary(crm)

regplot(crm,observation = dat.w[dat.w$ID==66&dat.w$failcode==1,],
        failtime = c(12,36,60),prfail = T, droplines=T)
```
```{r}
# 绘制校准曲线
xs=Score(list(Cox=fgr1),
         Surv(time,status)~1,
         data=df_test,
         plots="cal",
         metrics=NULL)
plotCalibration(xs, cens.method="local")

###例如评估两年的ROC及AUC值
f1 <- coxph(Surv(time,status==2)~age+sex+ph.ecog+
              ph.karno+pat.karno,
            data=aa,
            x=T)
model<- Score(list(model1=f1), 
            Hist(time, status==2)~1, 
            data = aa,
            times = 24,
            plots = 'roc',
            metrics ="auc")
plotROC(model, 
        xlab="1-Specificity",
        ylab="Sensitivity",
        lty=1,  #线型，2=虚线
        cex=1.1,#字体大小
        pch=2,  #文字格式
        lwd=2,    #线粗
        col="red",
        legend="模型1")

```