Cox模型，coxph()，无事件对照处理，种子萌发

Question

我正在进行生存分析，我不确定我做得是否正确。我的数据集是种子萌发实验的结果。感兴趣的主要变量是“对待”变量(具有3个级别的分类变量)。在我的脚本中，我试图通过比较PH系数百分比来找出不同处理之间是否存在差异，哪一个是最好的，以及在多大程度上。有人能帮我解决一些我正在处理的问题吗？

1)我需要声明变量as.factor()才能使用它们吗？还是integer被同等解释？

2)如果违反了危险假设(PH)的相称性，我应该如何处理我的数据才能继续构建cox模型？我进行了深入的研究，但还不能理解如何在我的模型中添加协变量*时间交互或分层。

3)如何在cox模型中加入脆弱项并检测随机效应(例如，种子发芽的图版，4个水平的分类变量，代表重复)。

4)我也不能解释打印(.*(cox.fra))摘要

*见下文

请看下面我的两个完整的脚本和注释。

脚本1

    rd01 <- read.table("sa_kb01.txt", header = T) # raw dataset, seed 
    survival
    rd01

    str(rd01) 

    rd01$begin <- as.factor(rd01$begin) # integers to factors
    rd01$spp <- as.factor(rd01$spp)
    rd01$cit <- as.factor(rd01$cit)
    rd01$treat <- as.factor(rd01$treat)
    rd01$plate <- as.factor(rd01$plate)

    str(rd01) 

    summary(rd01)

    names(rd01) # headers

    ### Survival analysis

    # install.packages("survival")

    library(survival)
    library (survminer)

    ?survfit
    ?survfit.formula
    ?survfit.coxph
    ?ggsurvplot

    ## Fit Kaplan-Meier survivor function

    km.fit <- survfit(Surv(day, status) ~ treat, data= rd01, type="kaplan-meier")
    km.fit
    print(summary(km.fit))

    plot(km.fit, conf.int= T, fun = "event", mark.time = c(140), pch = c("S", "W", "A"), col = c("darkred","darkblue","darkgreen"), lty = c("solid","dotted","longdash"),lwd = 1.5, xlab = "time [days]", ylab = "germination probability [%]")

    print(summary(km.fit))

    ## Comparison of Survivor Functions

    # Log-rank tests

    ?survdiff

    # Log-rank or Mantel-Haenszel test in "rho = 0" OR 
    # Peto & Peto modification of the Gehan-Wilcoxon test in "rho = 1"
    # ... Assess all groups for heterogeneity
    lrmh.123 <- survdiff(Surv(day,status) ~ treat, data= rd01, rho= 0) 

    print(lrmh.123) # If p<0.05 there are difference between all groups!

    # ... Comparing groups pairwise

    lrmh.120 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset= {treat!=3}, rho= 0)
    lrmh.103 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset= {treat!=2}, rho= 0)
    lrmh.023 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset= {treat!=1}, rho= 0)

    print(lrmh.120)
    print(lrmh.103)
    print(lrmh.023) # If p<0.05 there are difference pairwised groups!

    ## Checking Proportional Hazard (PH) assumption

    # Define function mlogmlog() to calculate -log(-log(S(t)))
    mlogmlog <- function(y){-log(-log(y))}

    # Use estimated Kaplan-Meier survivor functions
    km.fit

    # ... to plot -log(-log(S(t))) versus log(t)
    plot(km.fit, fun= mlogmlog, log="x", mark.time= c(140), pch = c("S", "W", "A"), col = c("darkred","darkblue","darkgreen"), lty = c("solid","dotted","longdash"), lwd = 1.5, xlab="time [days]", ylab= "-log(-log(S(t)))") # If lines do not cross, PH assumption is plausible!

    # Interpretarion: http://www.sthda.com/english/wiki/cox-model-assumptions#testing-proportional-hazards-assumption

    ## Checking for multicollinearity

    # install.packages("HH")
    library(HH)

    # Fit a generalized linear model predicting days from treatment
    ?glm
    mc.glm <- glm(day ~ treat, data=rd01)
    print(mc.glm) # doesn't need interpretation, only used to create object to VIF function

    # Check for multicollinearity among covariates throught variance inflation factor (VIF)
    ?vif
    mc.vif <- vif(mc.glm)
    print(mc.vif) # VIF can determine what proportion of the variation in each covariate 
    # is explained by the other covariates:
    # VIF > 10, serious multicollinearity; VIF = 5, evidence of multicollinearity;
    # VIF < 1, no evidence of multicollinearity

    ## Adding covariates to the Cox model

    # Create a Cox model
    cox.mod <- coxph(Surv(day, status) ~ treat, data= rd01)
    print(summary(cox.mod)) 

    # Interpretation: http://www.sthda.com/english/wiki/cox-proportional-hazards-model

    # Double check for PH assumption now with Cox model built
    dc.ph <- cox.zph(cox.mod)
    dc.ph  
    ggcoxzph(dc.ph) # if global and individual p-vale > 0.05, PH assumption is plausible! 

    ## Including random effects
    ?frailty

    # Adding plate variable as frailty term 
    cox.fra <- coxph(Surv(day, status) ~ treat + frailty(plate), data= rd01)
    print(summary(cox.fra)) # if global and individual p-vale < 0.05, 
    # maintain frailty term while adding covariates 1 at a time in cox model!`

脚本2-相同但不同的数据集，控制无事件的treat1！

    rd01 <- read.table("sa_hal01.txt", header = T) # raw dataset, seed         survival
    rd01

    str(rd01) 

    rd01$begin <- as.factor(rd01$begin) # integers to factors
    rd01$spp <- as.factor(rd01$spp)
    rd01$cit <- as.factor(rd01$cit)
    rd01$treat <- as.factor(rd01$treat)
    rd01$plate <- as.factor(rd01$plate)

    str(rd01) 

    summary(rd01)

    names(rd01) # headers

    ### Survival analysis

    # install.packages("survival")

    library(survival)
    library (survminer)

    ?survfit
    ?survfit.formula
    ?survfit.coxph
    ?ggsurvplot

    ## Fit Kaplan-Meier survivor function

    km.fit <- survfit(Surv(day, status) ~ treat, data= rd01, type="kaplan-meier")
    km.fit
    print(summary(km.fit))

    plot(km.fit, conf.int= T, fun = "event", mark.time = c(140), pch = c("S", "W", "A"), col = c("darkred","darkblue","darkgreen"), lty = c("solid","dotted","longdash"),lwd = 1.5, xlab = "time [days]", ylab = "germination probability [%]")

    print(summary(km.fit))

    ## Comparison of Survivor Functions

    # Log-rank tests

    ?survdiff

    # Log-rank or Mantel-Haenszel test in "rho = 0" OR 
    # Peto & Peto modification of the Gehan-Wilcoxon test in "rho = 1"
    # ... Assess all groups for heterogeneity
    lrmh.123 <- survdiff(Surv(day,status) ~ treat, data= rd01, rho= 0) 

    print(lrmh.123) # If p<0.05 there are difference between all groups!

    # ... Comparing groups pairwise

    lrmh.120 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset= {treat!=3}, rho= 0)
    lrmh.103 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset= {treat!=2}, rho= 0)
    lrmh.023 <- survdiff(Surv(day,status) ~ treat, data= rd01, subset=         {treat!=1}, rho= 0)

    print(lrmh.120)
    print(lrmh.103)
    print(lrmh.023) # If p<0.05 there are difference pairwised groups!

    ## Checking Proportional Hazard (PH) assumption

    # Define function mlogmlog() to calculate -log(-log(S(t)))
    mlogmlog <- function(y){-log(-log(y))}

    # Use estimated Kaplan-Meier survivor functions
    km.fit

    # ... to plot -log(-log(S(t))) versus log(t)
    plot(km.fit, fun= mlogmlog, log="x", mark.time= c(140), pch =         c("S", "W", "A"), col = c("darkred","darkblue","darkgreen"), lty =         c("solid","dotted","longdash"), lwd = 1.5, xlab="time [days]", ylab= "-        log(-log(S(t)))") # If lines do not cross, PH assumption is plausible!

    # Interpretarion: http://www.sthda.com/english/wiki/cox-model-        assumptions#testing-proportional-hazards-assumption

    ## Checking for multicollinearity

    # install.packages("HH")
    library(HH)

    # Fit a generalized linear model predicting days from treatment
    ?glm
    mc.glm <- glm(day ~ treat, data=rd01)
    print(mc.glm) # doesn't need interpretation, only used to create object to         VIF function

    # Check for multicollinearity among covariates throught variance inflation         factor (VIF)
    ?vif
    mc.vif <- vif(mc.glm)
    print(mc.vif) # VIF can determine what proportion of the variation in each covariate 
    # is explained by the other covariates:
    # VIF > 10, serious multicollinearity; VIF = 5, evidence of                 multicollinearity;
    # VIF < 1, no evidence of multicollinearity

    ## Adding covariates to the Cox model

    # Create a Cox model
    cox.mod <- coxph(Surv(day, status) ~ treat, data= rd01)
    print(summary(cox.mod)) 

    # Interpretation: http://www.sthda.com/english/wiki/cox-proportional-hazards-model

    # Double check for PH assumption now with Cox model built
    dc.ph <- cox.zph(cox.mod)
    dc.ph  
    ggcoxzph(dc.ph) # if global and individual p-vale > 0.05, PH assumption is                         plausible! 

    ## Including random effects
    ?frailty

    # Adding plate variable as frailty term 
    cox.fra <- coxph(Surv(day, status) ~ treat + frailty(plate), data=                 rd01)
    print(summary(cox.fra)) # if global and individual p-vale < 0.05, 
    # maintain frailty term while adding covariates 1 at a time in cox model!

似乎有统计学上的显着差异，并且treat3在这两个脚本中都与其他组不同。在脚本1中，PH被违反了，我现在不知道该怎么做。除此之外，脚本1中的Cox模型似乎工作得很好，风险比的解释也是可以的，但在脚本2中，不知道如何解释或解决这个问题(控制treat1中没有事件)。

Answer 1

1)我需要声明变量as.factor()才能使用它们吗？还是integer被同等解释？

我认为在你的情况下，as.factor是正确的。如果您有连续的数字变量，则可以使用整数-例如，如果您希望在实验之前存储时间种子，则可以使用as.numeric作为时间变量。

2)如果违反PH，我应该如何处理我的数据才能继续构建cox模型？我一直在深入研究，但还不能理解如何将协变量x时间交互或分层添加到我的模型中。

Cox回归，又称Cox比例风险模型，是建立在比例风险假设基础上的。如果违反了这个假设，你就得不到可靠的结果。您可能可以尝试一些数据转换，看看是否会有所帮助。或者，如果在某个子实验/小组中违反了它，您可以直接将其省略。