survit function and cox model with frailty
On Tue, 20 May 2003 gc4 at duke.edu wrote:
In this regard, I offer a brief clarification on what I'm attempting to do. The example I submitted is a simplified working example I am using to make sure I am able to make the code work. My assumption is that if it does not work on something simple, it will not work on something more complicated. I am estimating a model measuring the survival of political leaders in office in a sample of 1992 leaders from 166 countries from 1919 to 1999. My interest in on the effect of victory and defeat in war on leaders' survival in office. I include variables measuring economic conditions, domestic political institutions, domestic unrest, leaders' age and previous times in office, and war participation and war outcomes. I also include a country-level frailty term on the assumption that my covariates only partially capture the range of country-specific conditions affecting leaders' political survival. If I interpret a frailty term as a latent effect that enters multiplicatively into the specification of the hazard function, I can consider it as an additional covariate associated with a hidden coefficient of 1. Then, for example, I would ask what the survival probabilities are given a covariate path for a political leader in a high frailty country or in a low frailty country. And if this is statistically legitimate, can the survfit function accommodate a frailty term?
Yes, it is statistically legitimate. No, survfit can't do it. You could do
it yourself by extracting the baseline cumulative hazard and multiplying
it by the coefficients
An example, using the data in example(frailty)
data(kidney)
kfit <- coxph(Surv(time, status)~ age + sex + disease + frailty(id),
kidney
)
H<-basehaz(kfit,centered=FALSE)
# three new data points
temp <- kidney[1:3,c("age","sex","disease")]
# model.matrix without intercept
Mtemp <- model.matrix(~age+sex+disease,temp)[,-1]
# fitted log hazard ratio
logHR <- Mtemp%*%coef(kfit)[,1]
# turn it into a vector, not matrix
logHR <- drop(logHR)
# Hazard ratio
HR <- exp(logHR)
# survival curves with frailty=1
frail1 <- exp( -outer(H$hazard,HR))
# survival curves with frailty=2
frail2 <- exp( -outer(H$hazard,HR)*2)
# survival curves with frailty=0.5
frail.5 <- exp( -outer(H$hazard,HR)*0.5)
-thomas