Hi:
I have a question about the use of the survfit function after the
estimation of a cox proportional hazard model with a frailty term. My goal
is to estimate expected survival probabilities while controlling for the
group-specific frailty term.
First, I estimate a model of the following form:
model1 <- coxph(Surv(t0, t, d) ~ x1 + x2 + frailty(id), na.action=na.exclude,
data=My.data)
Then, I prepare a data frame:
temp <- data.frame(t0=c(0,365,730,1095,1460),
t=c(365,730,1095,1460,3000),
d=c(0,0,0,0,0),
x1=c(0,0,0,0,0),
x2=c(1.5,1.5,1.5,1.5,1.5))
I think I would need to enter a statement with respect to the frailty
term, but I don't know how.
Indeed, when I use the survit function with temp data frame, I get an
error message:
Giacomo -
Here's just a guess, since I have no experience with this.
Try adding a sixth variable to the data frame "temp", a
variable named "id" with value c(1,1,1,1,1). Maybe that
will do it ? This is consistent with your intention in
setting individual=TRUE, and it will give the frailty()
term that's already in the model something to work with.
- tom blackwell - u michigan medical school - ann arbor -
On Mon, 19 May 2003 gc4 at duke.edu wrote:
Hi:
I have a question about the use of the survfit function after the
estimation of a cox proportional hazard model with a frailty term. My goal
is to estimate expected survival probabilities while controlling for the
group-specific frailty term.
First, I estimate a model of the following form:
model1 <- coxph(Surv(t0, t, d) ~ x1 + x2 + frailty(id), na.action=na.exclude,
data=My.data)
Then, I prepare a data frame:
temp <- data.frame(t0=c(0,365,730,1095,1460),
t=c(365,730,1095,1460,3000),
d=c(0,0,0,0,0),
x1=c(0,0,0,0,0),
x2=c(1.5,1.5,1.5,1.5,1.5))
I think I would need to enter a statement with respect to the frailty
term, but I don't know how.
Indeed, when I use the survit function with temp data frame, I get an
error message:
Hi:
I have followed through Thomas' suggestions, but unfortunately, to no
avail. I added a statement identifying the id group variable to the temp
data frame which I feed into the survfit function.
temp <- data.frame(t0=c(0,365,730,1095,1460),
t=c(365,730,1095,1460,3000),
d=c(0,0,0,0,0),
x1=c(0,0,0,0,0),
x2=c(1.5,1.5,1.5,1.5,1.5),
id=c(1,1,1,1,1))
However, I obtain the following error message:
fit2 <- survfit(model1, newdata=temp, individual=TRUE)
Error in x2 %*% coef : non-conformable arguments
The survfit function seems to be looking for a coefficient for the frailty
term, which, of course, is not there given that:
h(t)=h0(t)*a*exp(Xb), where a is the frailty term.
I also tried to add a "frail(id)"=c(1,1,1,1,1) line to the statement
that creates the temp data frame. But, again, I obtain the following error
message:
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
variable lengths differ
Thank you very much for any suggestions. Thanks, in particular, to Thomas.
giacomo
On Mon, 19 May 2003, Thomas W Blackwell wrote:
Giacomo -
Here's just a guess, since I have no experience with this.
Try adding a sixth variable to the data frame "temp", a
variable named "id" with value c(1,1,1,1,1). Maybe that
will do it ? This is consistent with your intention in
setting individual=TRUE, and it will give the frailty()
term that's already in the model something to work with.
- tom blackwell - u michigan medical school - ann arbor -
Giacomo -
I'm stumped. Again, I have never had occasion to use the
survival package myself and I am just looking for syntax
irregularities in the code you show. The new error message
"Error in x2 %*% coef" is suggestive. Do take a look at
coef(model1) to see if there is something irregular there.
You might print temp out to the screen to see that x2 has
the exepcted value. Could the original data's data frame
have changed between the time model1 was created and the
time survfit() is run on it ?
(I've had SO much difficulty in the past with alternating
commas and decimal points that I would have defined x2 as
x2 = rep(1.5, 5).)
I reassure you that you CAN get to the bottom of this.
I think you're making progress. One step going forward
might be to look at help("model.frame"), since that's
what causes the initial error message, and think about
how the arguments get matched in survfit().
Aaaah. Here's one more thing you could try. Give the
new data frame of dimensions 5 x 6 some other name than
"temp" and try again. I see that the variable name "temp"
is already used in survfit() for something else. Again,
this is just a guess. Unfortunately, the package code
is very complex.
- tom blackwell - u michigan medical school - ann arbor -
On Mon, 19 May 2003 gc4 at duke.edu wrote:
... However, I obtain the following error message:
fit2 <- survfit(model1, newdata=temp, individual=TRUE)
Error in x2 %*% coef : non-conformable arguments
I also tried to add a "frail(id)"=c(1,1,1,1,1) line to the statement
that creates the temp data frame. But, again, I obtain the following error
message:
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
variable lengths differ
On Mon, 19 May 2003, Thomas W Blackwell wrote:
Giacomo -
Here's just a guess, since I have no experience with this.
Try adding a sixth variable to the data frame "temp", a
variable named "id" with value c(1,1,1,1,1). Maybe that
will do it ? This is consistent with your intention in
setting individual=TRUE, and it will give the frailty()
term that's already in the model something to work with.
- tom blackwell - u michigan medical school - ann arbor -
I reassure you that you CAN get to the bottom of this.
I think you're making progress. One step going forward
might be to look at help("model.frame"), since that's
what causes the initial error message, and think about
how the arguments get matched in survfit().
Actually I think one CAN'T get to the bottom of this. I don't think
survfit works on models with frailties.
I don't actually understand what the intent is (and why the multiple time
periods with identical covariates), but it isn't going to be at all
straightforward to do a proper prediction for a single new case: the
survival curve should be the survival distribution with the frailty
integrated out, which is hard.
It should be possible to do a prediction setting the frailty to 1, but it
isn't. Given that coxph will fit user-defined penalised likelihoods
survfit would have to be fairly clever to guess what it was supposed to do
in each circumstance.
-thomas
Hi:
Special thanks to Thomas B. and Thomas L.
The question that arises is whether it is statistically legitimate to
estimate survival probabilities after fitting a cox model with a frailty
term.
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?
Thanks again to you all.
giacomo
On Tue, 20 May 2003, Thomas Lumley wrote:
I don't actually understand what the intent is (and why the multiple time
periods with identical covariates), but it isn't going to be at all
straightforward to do a proper prediction for a single new case: the
survival curve should be the survival distribution with the frailty
integrated out, which is hard.
It should be possible to do a prediction setting the frailty to 1, but it
isn't. Given that coxph will fit user-defined penalised likelihoods
survfit would have to be fairly clever to guess what it was supposed to do
in each circumstance.
-thomas
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