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. ______________________________________________ R-help@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.