Ниже я запустил некоторый код, который смотрит на выполнение регрессии Кокса для нескольких типов результатов (инсульт, рак, респираторные заболевания), которые отображаются в отдельных столбцах. мурр, кажется, делает это довольно хорошо. Но я также хотел бы
Я знаю, что это довольно сложная задача, но она важна, поскольку мой реальный набор данных имеет почти 20 типов результатов. Любая помощь приветствуется!
library(survival)
library(purrr)
mydata <- read.table(header=T,
text = "age Sex survival stroke cancer respiratory
51 2 1.419178082 2 1 1
60 1 5 1 2 2
49 2 1.082191781 2 2 2
83 1 0.038356164 1 1 2
68 2 0.77260274 2 1 2
44 2 2.336986301 1 2 1
76 1 1.271232877 1 2 2")
outcomes <- names(mydata[4:6])
purrr::map(outcomes, ~coxph(as.formula(paste("Surv(survival,", .x, ") ~ Sex + age")),
mydata))
Я не совсем уверен, что это то, что вы ищете, но если вы запустите следующий код:
result <- purrr::map(outcomes, function(x) {
f <- as.formula(paste("Surv(survival,", x, ") ~ Sex + age"))
model <- coxph(f, mydata)
model$call$formula <- f
s <- summary(model)
cat(x, ':\n', paste0(apply(s$coefficients, 1,
function(x) {
paste0("HR : ", round(exp(x[1]), 2),
' (95% CI ', round(exp(x[1] - 1.96 * x[3]), 2),
' - ', round(exp(x[1] + 1.96 * x[3]), 2), ')')}),
collapse = '\n'), '\n\n', sep = '')
invisible(model)
})
Он распечатает:
#> stroke:
#> HR : 650273590159.06 (95% CI 0 - Inf)
#> HR : 1.36 (95% CI 0.75 - 2.49)
#>
#> cancer:
#> HR : 1121.58 (95% CI 0 - 770170911.09)
#> HR : 1.33 (95% CI 0.78 - 2.28)
#>
#> respiratory:
#> HR : 24.1 (95% CI 0.31 - 1884.85)
#> HR : 1.2 (95% CI 0.99 - 1.45)
И ваш список моделей будет храниться с правильным вызовом над ними:
result
#> [[1]]
#> Call:
#> coxph(formula = Surv(survival, stroke) ~ Sex + age, data = mydata)
#>
#> coef exp(coef) se(coef) z p
#> Sex 2.720e+01 6.503e+11 2.111e+04 0.001 0.999
#> age 3.105e-01 1.364e+00 3.066e-01 1.013 0.311
#>
#> Likelihood ratio test=6.52 on 2 df, p=0.03834
#> n= 7, number of events= 3
#>
#> [[2]]
#> Call:
#> coxph(formula = Surv(survival, cancer) ~ Sex + age, data = mydata)
#>
#> coef exp(coef) se(coef) z p
#> Sex 7.0225 1121.5843 6.8570 1.024 0.306
#> age 0.2870 1.3325 0.2739 1.048 0.295
#>
#> Likelihood ratio test=2.58 on 2 df, p=0.2753
#> n= 7, number of events= 4
#>
#> [[3]]
#> Call:
#> coxph(formula = Surv(survival, respiratory) ~ Sex + age, data = mydata)
#>
#> coef exp(coef) se(coef) z p
#> Sex 3.18232 24.10259 2.22413 1.431 0.1525
#> age 0.18078 1.19815 0.09772 1.850 0.0643
#>
#> Likelihood ratio test=5.78 on 2 df, p=0.05552
#> n= 7, number of events= 5
Привет, Спенсер. Функция cat
принимает аргумент file
, поэтому вы можете поместить file = 'HR.txt')
в конец cat
, и она запишет вывод в файл, а не в консоль.
Привет, Аллан, спасибо за этот ответ, это было действительно полезно. Я хотел спросить, есть ли способ сохранить результаты HR в объекте или хотя бы сохранить их в файл??