[R] survit function and cox model with frailty

[Thomas Lumley]

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