5
$\begingroup$

I am currently estimating a survival model (specifically, accelerated failure time model) with a log-logistic distribution using the survreg function in the survival R apckage. I want to simulate expected survival times in line with King et al. (2001), but I am unsure of the link function needed to calculate the expected survival time for the log-logistic distribution from the survreg regression output. I have included a minimal working example below:

library(survival)
data(kidney)
survreg(formula = Surv(time, status) ~ 
                  age + 
                  cluster(id), 
                  data = kidney, 
                  dist = "loglogistic", 
                  robust = TRUE)

##                Value Std. Err (Naive SE)     z        p
## (Intercept)  4.38127   0.6783     0.5338  6.46 1.05e-10
## age         -0.00298   0.0135     0.0114 -0.22 8.26e-01
## Log(scale)  -0.23009   0.0732     0.1038 -3.14 1.67e-03

## Scale= 0.794 

## Log logistic distribution
## Loglik(model)= -342   Loglik(intercept only)= -342
##  Chisq= 0.07 on 1 degrees of freedom, p= 0.79 
## (Loglikelihood assumes independent observations)
## Number of Newton-Raphson Iterations: 3 
## n= 76 

I simply want to know how I can calculate expected survival time from the estimated parameters from the survreg output.

Reference: King, G., Tomz, M., & Wittenberg, J. (2000). Making the most of statistical analyses: Improving interpretation and presentation. American journal of political science, 347-361.

$\endgroup$
3
$\begingroup$

Seems that this was more of a coding question and might have gotten a more prompt coding response on StackOverflow, but since no close votes have been offered I put in a belated CV response. Most R regression functions have an associated predict method and survival::survreg is no exception. You need to assign the output of the model call to a named object and then run predict:

sreg.model <- .Last.value
predict(sreg.model)
#_________
 [1] 73.54492 73.54492 69.29323 69.29323 72.67421 72.67421
 [7] 72.89091 72.67421 77.59404 77.59404 76.22015 75.99355
snipped

So those are the expected values for each individual combinations of covariates in the original dataset. If you wanted now predictions that might be more suitable for constructing plots, you would supply a newdata argument in the form of dataframe. See ?predict.survival for a worked example. It also shows how to plot expected survival curves and if you picked out the expected 50% survival you would have the predicted median values.

pct <- 1:98/100   
ptime <- predict(sreg.model, newdata=data.frame(age=10:69, id=1) , type='quantile',   p=pct, se=TRUE)
ptime$fit[ ,50]
 [1] 77.59404 77.36335 77.13335 76.90403 76.67539 76.44743
 [7] 76.22015 75.99355 75.76762 75.54236 75.31777 75.09384
[13] 74.87059 74.64800 74.42606 74.20479 73.98418 73.76422
#snipped
$\endgroup$
1
  • $\begingroup$ Thanks! Sorry for the much belated reply. This solved it. $\endgroup$ – chrstnsn Jun 1 '15 at 9:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.