In the past few months I've posted a few questions about survival
analysis... all of the replies have come from T. Therneau... a summary of
the questions and responses are below. Might be helpful to some just
getting started....
Cheers
andrew
-----------------------------------
#Question: Can one actually code interactions in Parametric models and Are
there any example of how to interpret interactions in survival models....
VR has a little on this...
#Response
Yes, interactions work fine in the survival focused regression routines.
Why no examples?
I can answer for the Splus manuals, since I'm responsible for most of their
examples. They all happen to use fairly small data sets, not surprizing if
you knew that they all started as a part of my test suite (bugs are easier to
narrow down in a small data set). Because they are small, there just isn't
a big enough N to entertain interactions from a statistical/analysis point of
view.
------------------------------------------------
#Question
1) How does one "really" interpret the predict( ,type="terms") predictions
- this plus the SPLUS guide examples p336-8 were helpful
2) How to fit a model to each "unit" in an experimental design... i.e. if
your treatments are applied plot or hospital wide, how could you get an
estimate per plot.... The answer is useful, but I found very slow for large
data sets... use a for loop indexed over the unit, with the model being
intercept only produced a much faster result (3 minutes vs. >2days for
95000 observations in 192 units).
#Response
Here are some thoughts on your questions:
1. survReg fits models of the type t(y) ~ D(mean= Xb, std=sigma), where
X are the covariates, b=the coefficients, D is some named distribution, and
t is a transform.
By default (Wiebull model) t =log and D= extreme value.
predict(fit, type='terms') will give you the components of Xb, but when
looking at the terms you need to always be aware of t(), at least to understand
what "partial prediction" means.
2. fit <- survReg(Surv(time, status) ~ factor(inst) + strata(inst), data=lung)
is a way to fit a separate mean and variance to each institution, i.e.,
it is the same as fitting a separate model to each institution.
--------------------------------------------------------------
#Question:Following a strategy of beginning with a full model (4-way
interaction), and getting a signficant p-value in coxph or survReg for that
4 way, I am curious as how to:
1) be confident of the 4-way interaction
2) be justified in stratifying or factoring
3) interpret the interactions
4) If i stratify under censorReg, the assumption is that a seperate model
is being fit under each strata. Can you modify the model under each strata
independently-i.e.if a 3 way interaction is insignificant in one strata
level, i'd like to remove it and refit that portion of the model.
I have examined the proportional hazard assumptions via the log and cloglog
plots and beta(t) plots and they are not too good (the cox.zph results are
almost all significant; the log plots are parallel but not linear) - plus,
there are convergence errors with an unstratified model in coxph.
#Response:
1. The non interaction model:
There is a separate baseline death rate for each density level, but
"fed", for instance, decreases the death rate by 12%, the same "treatment"
effect within each density level. The interaction model says that the
effect of "fed" might differ for different densities. The interaction model
is essentially equivalent to fitting each density as a separate run.
The classic interaction: time to maturity for corn, stratified by amount
water (irrigated/not), the effect of fertilizer. Fertilizer could help the
irrigated and hurt the non-irrigated.
2. With lots of interactions I find that Splus doesn't always count
correctly, and puts in more dummy variables than it needs. This appears to
be the cause of most of your NA coefficients.
The issue isn't unique to coxph.
3. Non-ph: Like all regression models, the success of the Cox model depends
on how accurately it can capture the effect of a variable in a single
coefficient. For analysis of covariance models, a key question is common
slopes (drug effect controlling for age, for instance). If common slopes
holds, then a possibly complex situation is well desribed by a small set
of intercepts.
In the Cox model, assume the coef of "fed" is 0.5. This says that the
fed subjects have an exp(0.5) = 1.6 fold higher death rate than non-fed.
The PH assumption is that it is 1.6 fold higher on day 1, 1.6 fold higher on
day 2, 1.6 fold higher on day 3, .... Very similar in interpretation to
the parallel line assumption in ancova.
If the assumption is not true -- perhaps it's time to try another model.
The survReg/censorReg pair are a good alternate choice, and much closer
to glm's in their interpretation.
(The programs are nearly identical, just a biomedical slant vs an industrial
reliability slant).
---------------------------------------------------
#Question
Basically:
1)how does one use weights?
2) what really is the unit of observation in a SA; if there are two
replicates of a treatment in a exp design, does replicate need to be coded....
#Response
To collapse the data, you can use case weights within a tube: if
27 mites die on day 1, 4 are lost to fu on day 1 (you actually lost them
when counting), 83 die on day 2, 35 on day 3, ... this would be only 4
observations with weights of 27, 4, 83 and 35.
In a Cox model, there is no "residual" variance (it's essentially a
Poisson model), so the tube vs individual distinction is completely lost.
You can't summarize a tube in a single observation, due to the different
death days, but that is just a technical issue of how data is represented
to the fitting function.
The model you wrote has density as a strata -- which means something
quite different from "Error" strata in a linear model. The two statistical
areas used the word for different things. Your model says: there is a
baseline death rate for each of the 3 densities (much like multiple
intercepts in an analysis of covariance model), and the other factors
multiply this baseline rate. No assumption is made about how the baseline
rates relate to each other, or their shape over time. You can examine
them, after a Cox fit, with the survfit() function.
-------------------------------------------------------------------------------
#Question
How to get Fitted Curves for a cox model with strata and interactions
#Response
You need to use the 'survfit' function to get baseline hazard estimates.
fit <- coxph(Surv(event,cens) ~ strata(cjdens) + jft*jfa, data=juv.work)
temp <- data.frame(jft=c(0,0,1,1), jfa=c(0,1,0,1)) #four combinations
curves <- survfit(fit, newdata=temp)
Now, "curves" will contain 8 survival curves. Print and plot them as you
would a Kaplan-Meier created by survfit(). They are the predictions for
the four covariate combinations, within strata 1 and then within strata 2.
------------------------------------------------------------------
#Question
A coding question about using survReg to accomplish what censorReg does...
which was all along in the SPLUS manual...
#Response
fit <- survReg(Surv(time, status) ~ (x1 + x2 +x3) * strata(group) +
factor(group),...
This is a model with separate coefficients for x1,x2 and x3 in each group,
and separate scale factors for each group, and a separate intercept for
each group (the last term); that is, it is the same as
fitting a separate model to each group.
Thus, the survreg strata() is more general, because it allows models, e.g.
separate scale but not separate intercept, that censorReg does not.
------------------------------------------------------------------------------
Dr. Andrew Beckerman
Institute of Biological Science
University of Stirling, Stirling FK9 4LA, Scotland, UK
phone:+44 (0)1786 467808 fax: +44 (0)1786 464994
|