Я выполняю анализ выживания, и я не уверен, что делаю это правильно.Мой набор данных является результатом эксперимента по прорастанию семян.Основная переменная интереса - «лечить» (категориальная с 3 уровнями).В своем сценарии я пытаюсь выяснить, есть ли разница между обработками, какая из них является лучшей, и в какой степени, сравнивая проценты коэффициента PH.Может ли кто-нибудь помочь мне с некоторыми проблемами, с которыми я сталкиваюсь?
1) Нужно ли объявлять мои переменные как .factor (), чтобы использовать их?Или целое число интерпретируется одинаково?
2) Если допущена пропорциональность предположения об опасностях (PH), что мне делать с моими данными, чтобы перейти к построению модели Кокса?Я интенсивно исследовал, но не смог понять программирование, чтобы добавить ковариацию * временное взаимодействие или стратификацию в мою модель.
3) Как включить хрупкие термины в модель Кокса и обнаружить случайный эффект (например, пластину)в котором семена прорастали, категориальная переменная с 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 - тот же, но другой набор данных, контрольная обработка1 без событий!
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, кажется, работает нормально, и интерпретация коэффициентов опасности в порядке, но в сценарии 2 нет понятия, как это интерпретировать или решить (в контрольной обработке1 не было события).