Desejo criar dados de sobrevivência de brinquedos (tempo até o evento) que sejam corretamente censurados e sigam alguma distribuição com riscos proporcionais e riscos constantes da linha de base.
Criei os dados da seguinte maneira, mas não consigo obter taxas de risco estimadas próximas dos valores reais depois de ajustar um modelo de riscos proporcionais de Cox aos dados simulados.
O que eu fiz errado?
Códigos R:
library(survival)
#set parameters
set.seed(1234)
n = 40000 #sample size
#functional relationship
lambda=0.000020 #constant baseline hazard 2 per 100000 per 1 unit time
b_haz <-function(t) #baseline hazard
{
lambda #constant hazard wrt time
}
x = cbind(hba1c=rnorm(n,2,.5)-2,age=rnorm(n,40,5)-40,duration=rnorm(n,10,2)-10)
B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)
hist(x %*% B) #distribution of scores
haz <-function(t) #hazard function
{
b_haz(t) * exp(x %*% B)
}
c_hf <-function(t) #cumulative hazards function
{
exp(x %*% B) * lambda * t
}
S <- function(t) #survival function
{
exp(-c_hf(t))
}
S(.005)
S(1)
S(5)
#simulate censoring
time = rnorm(n,10,2)
S_prob = S(time)
#simulate events
event = ifelse(runif(1)>S_prob,1,0)
#model fit
km = survfit(Surv(time,event)~1,data=data.frame(x))
plot(km) #kaplan-meier plot
#Cox PH model
fit = coxph(Surv(time,event)~ hba1c+age+duration, data=data.frame(x))
summary(fit)
cox.zph(fit)
Resultados:
Call:
coxph(formula = Surv(time, event) ~ hba1c + age + duration, data = data.frame(x))
n= 40000, number of events= 3043
coef exp(coef) se(coef) z Pr(>|z|)
hba1c 0.236479 1.266780 0.035612 6.64 3.13e-11 ***
age 0.351304 1.420919 0.003792 92.63 < 2e-16 ***
duration 0.356629 1.428506 0.008952 39.84 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
hba1c 1.267 0.7894 1.181 1.358
age 1.421 0.7038 1.410 1.432
duration 1.429 0.7000 1.404 1.454
Concordance= 0.964 (se = 0.006 )
Rsquare= 0.239 (max possible= 0.767 )
Likelihood ratio test= 10926 on 3 df, p=0
Wald test = 10568 on 3 df, p=0
Score (logrank) test = 11041 on 3 df, p=0
mas os valores verdadeiros são definidos como
B = c(1.1,1.2,1.3) # hazard ratios (model coefficients)
survival
cox-model
monte-carlo
stats_newb
fonte
fonte
Respostas:
Não está claro para mim como você gera seus horários de evento (que, no seu caso, podem ser ) e indicadores de evento:< 0
Então, aqui está um método genérico, seguido por algum código R.
Gerando tempos de sobrevivência para simular modelos de riscos proporcionais de Cox
Para gerar tempos de eventos a partir do modelo de riscos proporcionais, podemos usar o método de probabilidade inversa (Bender et al., 2005) : se for uniforme em e se é a função de sobrevivência condicional derivada do modelo de riscos proporcionais, ou seja, então é fato que a variável aleatória possui a função de sobrevivênciaV ( 0 , 1 ) S( ⋅|x )
Exemplo [risco de linha de base Weibull]
Deixe com a forma e a escala . Então e . Seguindo o método de probabilidade inversa, uma realização de é obtida calculando com uma variável uniforme em . Usando resultados em transformações de variáveis aleatórias, pode-se notar que tem uma distribuição Weibull condicional (dadoh0(t)=λρtρ−1 ρ>0 λ>0 H0(t)=λtρ T∼S(⋅H−10(t)=(tλ)1ρ t = ( - log ( v )T∼S(⋅|x) v(0,1)Txρλexp(x′β)
Código R
A função R a seguir gera um conjunto de dados com uma única covariável binária (por exemplo, um indicador de tratamento). O risco da linha de base tem uma forma Weibull. Os tempos de censura são sorteados aleatoriamente a partir de uma distribuição exponencial.x
Teste
fonte
flexsurvreg(Surv(time, status) ~ x, data=dat, dist = "weibull")
os mesmos dados simulados, o coeficiente aparece como0.6212
. Por que é isso?então eu modifiquei assim
se rho = 1, o resultado será o mesmo.
fonte