[R] AFT model

[tsn4867]

     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
>
>
>