[R] AFT model
tsn4867
tsn4867 at rice.edu
Wed Mar 18 21:58:28 CET 2009
Thank you, Dr. Lumley. I have implemented the following code for
the pareto distribution (see below). However, the estimates obtained
from survreg are very small & inaccurate. What I need help with is the
function for the deviance (the code below is wrong). I just don't
understand how to obtain this since the manual is not clear.
Thank you,
T.N.
library(survival)
#pareto: f(x) = a/(x+1)^(a+1), x>0, a>2#
# S(t) = 1/(t+1)^a, a>2, t>0#
survreg.distributions$pareto$name<-"pareto"
survreg.distributions$pareto$variance<-function(a)
a/((a-1)^2*(a-2))
survreg.distributions$pareto$parms$a<-3
survreg.distributions$pareto$init<-function(x,weights,a){
if(a<=2)
stop("a needs to be greater than 2")
mean <- sum(x*weights)/sum(weights)
var <- sum(weights*(x-mean)^2)/sum(weights)
c(mean+1/(a-1),var*(a-1)^2*(a-2)/a)
}
survreg.distributions$pareto$deviance<-function(y,scale,parms){
a<-parms; status <- y[,ncol(y)]
width <- ifelse(status==3, (y[,2]-y[,1])/scale,1)
temp <- width/(exp(width)-1)
center <- ifelse(status==3, y[,1]-log(temp),y[,1])
temp2 <- (-temp) + log(1-1/(width+1)^a)
best <- ifelse(status==1, -log(a*scale),
ifelse(status==3,temp2,0))
list(center=center,loglik=best)
}
survreg.distributions$pareto$density<-function(x,a){
w<-1/(x+1)^a; ww<-a/(x+1)
cbind(1-w, w, w * ww, -(a+1)/(x+1), ((a+1)*(a+2))/(x+1)^2)
}
survreg.distributions$pareto$quantile<-function(p,a){
1/(1-p)^(1/a) - 1
}
#example
#pareto (a=3) dist.#
my.dist<-survreg.distributions$pareto
my.dist$name<-"pareto"
X<-matrix(rnorm(20*5,0,.1),20,5); beta1<-matrix(runif(5,.5,.5),5,1)
e<-1/(1-runif(20,0,1))^(1/3)-1
y<-X%*%beta1+e
status<-c(rep(0,4),rep(0,16))
survreg(Surv(y,status),dist="extreme")
survreg(Surv(y,status)~X,dist=my.dist,parms=3)
Thomas Lumley wrote:
> On Tue, 17 Mar 2009, tsn4867 wrote:
>
>> Hi,
>> In the package survival, using the function survreg for AFT model, I
>> only see 4 distributions for the response y: weibull, gaussian,
>> logistic, lognormal and log-logistic, which correspond to certain
>> distributions for the error terms. I'm wondering if there is a
>> package or how to obtain the parameter estimates (the beta's are of
>> great interest) from the AFT model (maximizing log-likelihood with
>> censored data as in Klein and Moeschberger book) if we have
>> 1) error ~ Cauchy(0,1)
>> 2) error comes from a contaminated gaussian dist., for example,
>> density of error, f(x) = 0.9*phi(x, mean=0, sd=1) + 0.1*phi(x,
>> mean=0, sd=3), where phi(.) is the standard gaussian pdf.
>
> If you look at
> ?survreg.distributions
> it has an example showing how to specify your own distributions.
>
> -thomas
>
> Thomas Lumley Assoc. Professor, Biostatistics
> tlumley at u.washington.edu University of Washington, Seattle
>
>
>
More information about the R-help
mailing list