Eleni Rapsomaniki
2010-Dec-10 16:21 UTC
[R] survreg vs. aftreg (eha) - the relationship between fitted coefficients?
Dear R-users,
I need to use the aftreg function in package 'eha' to estimate failure
times for left truncated survival data. Apparently, survreg still cannot fit
such models. Both functions should be fitting the accelerated failure time
(Weibull) model. However, as G?ran Brostr?m points out in the help file for
aftreg, the parameterisation is different giving rise to different coefficients.
The betas for adjusted covariates are opposite in sign but otherwise identical,
whereas the intercept is quite different in a non-obvious way. The
log-likelihoods are similar also, but not identical. I would like to find out
how I can convert one set of coefficients to the other so as to obtain the same
linear predictors using either model. Any ideas???
#the example below uses right-censored data for simplicity (the principle should
be the same with left truncation I hope)
library(survival)
library(eha)
# COMPARE coefs between survreg ('survival' pkg) and aftreg
('eha' pkg)
#Fitting NULL models (no covariates) results in (approximately) the same coefs
(which is good!)
m1_NULL=survreg(Surv(futime/365, status==1) ~ 1, data=pbcseq)
m2_NULL=aftreg(Surv(futime/365, status==1) ~ 1, data=pbcseq)
c(m1_NULL$coef, 1/m1_NULL$scale) #--> intercept= 3.878656 , shape =
1.478177
c(m2_NULL$coef[1], exp(m2_NULL$coef[2])) #--> intercept= 3.878859 ,
shape=1.478150
# NOW I adjust for covariates
m1=survreg(Surv(futime/365, status==1) ~ chol+stage, data=pbcseq)
m2= aftreg(Surv(futime/365, status==1) ~ chol+stage, data=pbcseq)
### m2 #######
#Coefficients:
# (Intercept) chol stage
# 5.944641913 -0.001692574 -0.470861324
#Scale= 0.6416744
#Loglik(model)= -483.9 Loglik(intercept only)= -506.8
# Chisq= 45.91 on 2 degrees of freedom, p= 1.1e-10
#n=1124 (821 observations deleted due to missingness)
### m2 #######
#Covariate W.mean Coef Exp(Coef) se(Coef) Wald p
#chol 303.777 0.002 1.002 0.000 0.000
#stage 3.298 0.460 1.584 0.119 0.000
#
#log(scale) 5.029 152.807 0.477 0.000
#log(shape) 0.467 1.595 0.095 0.000
#
#Events 92
#Total time at risk 9017
#Max. log. likelihood -484.31
#LR test statistic 45.0
#Degrees of freedom 2
#Overall p-value 1.64669e-10
Many thanks for any help you may be able to provide.
Eleni Rapsomaniki
Research Associate
University of Cambridge
Institute of Primary and Public Health
Göran Broström
2010-Dec-14 17:00 UTC
[R] survreg vs. aftreg (eha) - the relationship between fitted coefficients?
On Fri, Dec 10, 2010 at 5:21 PM, Eleni Rapsomaniki <er339 at medschl.cam.ac.uk> wrote:> Dear R-users, > > I need to use the aftreg function in package 'eha' to estimate failure times for left truncated survival data.Be careful! This is only possible under strong assumptions about the (unobserved) covariate vector. This is not (yet) reflected on the help page for aftreg, but has been discussed here earlier. Apparently, survreg still cannot fit such models. Both functions should be fitting the accelerated failure time (Weibull) model. However, as G?ran Brostr?m points out in the help file for aftreg, the parameterisation is different giving rise to different coefficients. The betas for adjusted covariates are opposite in sign but otherwise identical, whereas the intercept is quite different in a non-obvious way. The log-likelihoods are similar also, but not identical. I would like to find out how I can convert one set of coefficients to the other so as to obtain the same linear predictors using either model. Any ideas??? Yes. Since you work with Weibull data, try 'phreg' instead. You will have to change sign of estimated coefficients, and divide them by estimated shape, see the help page for 'weibreg', to compare to parameters from survreg, but they are directly comparable with those from coxph and coxreg. G?ran> > #the example below uses right-censored data for simplicity (the principle should be the same with left truncation I hope) > library(survival) > library(eha) > > # ?COMPARE coefs between survreg ('survival' pkg) and aftreg ('eha' pkg) > #Fitting NULL models (no covariates) results in (approximately) the same coefs (which is good!) > > m1_NULL=survreg(Surv(futime/365, status==1) ~ 1, data=pbcseq) > m2_NULL=aftreg(Surv(futime/365, status==1) ~ 1, data=pbcseq) > > c(m1_NULL$coef, 1/m1_NULL$scale) #--> intercept= 3.878656 ?, ?shape = 1.478177 > c(m2_NULL$coef[1], exp(m2_NULL$coef[2])) #--> intercept= 3.878859 , ?shape=1.478150 > > > # NOW I adjust for covariates > > m1=survreg(Surv(futime/365, status==1) ~ chol+stage, data=pbcseq) > m2= aftreg(Surv(futime/365, status==1) ~ chol+stage, data=pbcseq) > > ### ? ? ?m2 ?####### > #Coefficients: > # (Intercept) ? ? ? ? chol ? ? ? ?stage > # 5.944641913 -0.001692574 -0.470861324 > > #Scale= 0.6416744 > > #Loglik(model)= -483.9 ? Loglik(intercept only)= -506.8 > # ? ? ? ?Chisq= 45.91 on 2 degrees of freedom, p= 1.1e-10 > #n=1124 (821 observations deleted due to missingness) > > ### ? ? ?m2 ?####### > > #Covariate ? ? ? ? ?W.mean ? ? ?Coef Exp(Coef) ?se(Coef) ? ?Wald p > #chol ? ? ? ? ? ? ?303.777 ? ? 0.002 ? ? 1.002 ? ? 0.000 ? ? 0.000 > #stage ? ? ? ? ? ? ? 3.298 ? ? 0.460 ? ? 1.584 ? ? 0.119 ? ? 0.000 > # > #log(scale) ? ? ? ? ? ? ? ? ? ?5.029 ? 152.807 ? ? 0.477 ? ? 0.000 > #log(shape) ? ? ? ? ? ? ? ? ? ?0.467 ? ? 1.595 ? ? 0.095 ? ? 0.000 > # > #Events ? ? ? ? ? ? ? ? ? ?92 > #Total time at risk ? ? ? ? ?9017 > #Max. log. likelihood ? ? ?-484.31 > #LR test statistic ? ? ? ? 45.0 > #Degrees of freedom ? ? ? ?2 > #Overall p-value ? ? ? ? ? 1.64669e-10 > > Many thanks for any help you may be able to provide. > > Eleni Rapsomaniki > Research Associate > University of Cambridge > Institute of Primary and Public Health > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >-- G?ran Brostr?m