У меня возник вопрос по поводу "центрирования" сплайн-графика. Какое значение plotHR принимает за коэффициент риска = 1? Можно ли, например, присвоить коэффициент опасности 1 среднему возрасту (или среднему возрасту) при использовании функции coxph? Можно ли вручную выбрать значение оси x (например, возраст), которое не связано с величиной оси y Коэффициент опасности = 1 при использовании функции coxph и pspline?
plotHR включен вГрег Пакет. Приведенный ниже пример представляет собой слегка измененную версию этого plotHR-примера: https://rdrr.io/cran/Greg/man/plotHR.html
Приведенный ниже код назначит отношение опасности 1 к значению оси x где-то между срединным возрастом 3-го квартиля. Как я могу назначить средний возраст для соотношения рисков 1?
library(Greg)
library(survival)
library(rms)
# Get data for example
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
# Create a data frame since plotHR will otherwise
# have a hard time getting the names of the variables
ds <- data.frame(
dt = dt,e = e,
age=age,
smoking=smoking,
sex=sex)
library(splines)
Srv <- Surv(dt,e)
org_par <- par(xaxs="i", ask=TRUE)
dd <- datadist(ds)
options(datadist='dd')
fit.cph <- coxph(Srv ~ pspline(age,4) + sex + smoking, data=ds, x=TRUE, y=TRUE)
unadjusted_fit <- coxph(Srv ~ pspline(age,4), data=ds, x=TRUE, y=TRUE)
plotHR(list(fit.cph, unadjusted_fit), term="age", xlab="Age",
polygon_ci=c(TRUE, FALSE),
col.term = c("#08519C", "#77777799"),
col.se = c("#DEEBF7BB", grey(0.6)),
lty.term = c(1, 2),
plot.bty="l", xlim=c(30, 90),rug="ticks")
par(org_par)
Результатом кода выше является следующий график (красное поле было добавлено вручную, чтобы показать, что HR = 1, средний возрастне назначены друг другу): Результирующий сплайн-график