Я пытаюсь построить сплайн-график, который показывает отношение рисков на оси Y и переменную риска на оси X. Я думаю, что я, вероятно, нашел решение, но я не совсем уверен, действительно ли ось Y показывает отношение рисков или скорее показывает что-то другое. В моем примере переменная риска - это возраст. Коэффициент опасности рассчитывается с помощью регрессии Кокса, скорректированной с учетом переменных риска «курение» и «пол». Вот код:
library(dplyr)
library(pspline)
library(Greg)
library(survival)
library(rms)
############ Creating dataset ####################
n <- 1000
set.seed(731)
age <- c(rnorm(50, mean=20, sd=1), round(50 + 12*rnorm(900), 1),rnorm(50, mean=80, sd=1))
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n,
rep=TRUE, prob=c(.6, .4)))
cens <- c(17*runif(n-60),1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
smoking <- factor(sample(c('Yes','No'), n,
rep=TRUE, prob=c(.2, .75)))
h <- .02*exp(.02*(age-50)+.1*((age-50)/10)^3+.8*(sex=='Female')+2*(smoking=='Yes'))
dt <- -log(runif(n))/h
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
# Add missing data to smoking
smoking[sample(1:n, round(n*0.05))] <- NA
sex<- as.factor(sex)
smoking<-as.factor(smoking)
age<- as.numeric(age)
ds <- data.frame(
dt = dt,e = e,
age=age,
smoking=smoking,
sex=sex)
ds <- na.omit(ds)
###################################################
library(splines)
Srv <- Surv(ds$dt,ds$e)
fit.death <- coxph(Srv ~ pspline(age,4) + sex + smoking, data=ds, x=TRUE, y=TRUE)
## get predicted values for fitted spline
# what do the values of predicted( type="terms") mean? Are those the Hazard ratios? or the log(Hazard ratio)?
predicted <- predict(fit.death , type = "terms" , se.fit = TRUE , terms = 1)
#Termplot delivers the same values sorted by age?
termplot(fit.death, se=FALSE, plot=FALSE,term=1)
## Set the median age as center in order to get hazard ratio= 1 for the median age
ptemp <- termplot(fit.death, se=TRUE, plot=FALSE)
ageterm <- ptemp$age # this will be a data frame
center <- ageterm$y[which.min(abs(ageterm$x-median(age)))]
## Plot 1
plot( ds$age, exp(predicted$fit-center) , type="n" , xlim = c(20,83),ylim = c(0.1, 20), log="y")
#
lines( sm.spline(ds$age , exp(predicted$fit-center)) , col = "black" , lty = 1 )
lines( sm.spline(ds$age , exp(predicted$fit + 1.96 * predicted$se-center)) , col = "red" , lty = 2 )
lines( sm.spline(ds$age , exp(predicted$fit - 1.96 * predicted$se-center)) , col = "red" , lty = 2 )
##Plot 2
# This delivers the same plot by using matplot instead
ytemp <- ageterm$y + outer(ageterm$se, c(0, -1.96, 1.96), "*" ) #
matplot(ageterm$x, exp(ytemp - center), log= "y" ,
type= "l" , lty=c(1,2,2), col=1,
xlab="Age", ylab="Hazard ratio", ylim = c(0.1,20), xlim = c(20,83)) #
##Plot 3:
plotHR(fit.death,se=TRUE, term="age", plot.bty="o", xlim=c(20,83),ylim=c(0.1,20), ylog=TRUE, rug="ticks", xlab="age") # alternativ für term=1 kann man wie hier zu sehen auch den Variablennahmen hier Peguero_mV einfügen
##Plot 4
#what do the y values of the plot created by termplot mean?
termplot(fit.death,term=1, se=FALSE, plot=TRUE)
# What does the y-axis of the plot without the exponential of the predict function mean?
##Plot 5
noexp<-predicted$fit-center
noexpSE1<-predicted$fit + 1.96 * predicted$se-center
noexpSE2<-predicted$fit - 1.96 * predicted$se-center
plot( ds$age, noexp , type="n" , xlim = c(20,83),ylim = c(0.1, 20))
#
lines( sm.spline(ds$age , noexp) , col = "black" , lty = 1 )
lines( sm.spline(ds$age , noexpSE1) , col = "red" , lty = 2 )
lines( sm.spline(ds$age , noexpSE2) , col = "red" , lty = 2 )
что означают значения прогнозируемых (тип = "условия")? Это те крысы Hazard ios (+ SE-значения)? или журнал (коэффициент опасности)?
predict(fit.death , type = "terms" , se.fit = TRUE , terms = 1)
Termplot выдает те же «коэффициенты опасности?» - значения, отсортированные по возрасту с соответствующим x (= возраст ) и y (= коэффициент опасности?) значение?
termplot(fit.death, se=FALSE, plot=FALSE, term=1)
Графики 1 и 2, похоже, показывают один и тот же результирующий график с коэффициентом опасности при y- Ось и возраст на оси х?
График 1 слева, График 2 справа, ось X: возраст, ось Y: Коэффициент опасности?
Почему строится exp(predicted$fit-center)), а не просто predicted$fit-center, как показано на графике 5?
На графике 3 создается один и тот же график с двумя проблемами:
- средний возраст не имеет отношения риска 1
- средний рассчитывается для значений возраста в пределах xlim, а не для всех значений возраста
График 3, созданный с помощью команды Greg package plotHR, на правом графике добавлена вертикальная линия для среднего возраста и горизонтальная линия для коэффициента опасности = 1
Что осевые значения o f сюжет, созданный с помощью термплота? (График 4) termplot(fit.death,term=1, se=FALSE, plot=TRUE)
График 4
Что означает ось Y графика без exp () прогнозирования -значение функции значит? (Участок 5)
участок 5