[R] Survreg function for loglogistic hazard estimation

[chenjiakai]

I am trying to use R to do loglogistic hazard estimation.  My plan is to
generate a loglogistic hazard sample data and then use survreg to estimate
it.  If everything is correct, survreg should return the parameters I have
used to generate the sample data.  

I have written the following code to do a time invariant hazard estimation. 
The output of summary(modloglog) shows the factor loading of dfico is -0.002
instead of 0.005.  Also I can not find the survreg's output of alpha and
beta of loglogistic regression.  Any comments?  Thanks a lot, 

Jiakai Chen

Code: 

============BEGIN==============
# A time invariant loglogistic estimation
library(survival)

# Clear the workspace
rm(list = ls(all = TRUE))

# Totally 100K samples which may die during 100 periods 
timeline = 100 
samplesize = 100000


dfico <- rnorm(samplesize, mean = 0, sd = 25)
dficoeff <- 0.005

# Baseline loglogistic hazard function stored in bhaz
a<- 20
b <- 4
time <- 1:timeline
bhaz <- ((b/a) * (time /a) ^ (b-1))/(1+(time/a)^b)

# Event time stored in endtime.  Baseline hazard function is controlled by 
# dfico value with factor loading dficoeff

endtime <- numeric(samplesize)
for ( i in 1:samplesize) {
rnos <- runif(timeline)
endtime[i] <- which(rnos <= bhaz * exp(dfico[i]*dficoeff))[1]
}

# Construct data frame
raw <- data.frame(end = endtime, status = 1, fico_demean = dfico)

# Adding censorship
for( i in 1:samplesize) {
if(is.na(raw$end[i])) {raw$status[i] <- 0} 
}
for( i in 1:samplesize) {
if(is.na(raw$end[i])) {raw$end[i] <- timeline }
}

# Output the factory
factory <- factor(raw$end)
plot(factory)

# Use survreg to estimate the coefficents of loglogistic distribution
modloglog <- survreg(Surv(end, status) ~ fico_demean, data = raw, model =
TRUE, dist = 'loglogistic')

summary(modloglog)
============END==============
-- 
View this message in context: http://www.nabble.com/Survreg-function-for-loglogistic-hazard-estimation-tp23907598p23907598.html
Sent from the R help mailing list archive at Nabble.com.