We use data from the previous example and also convert the concentrations into a factor:
We can fit parametric survival time models using the survreg() function, where the dist argument specifies the assumed distribution. The Gamma distribution is not supported by survreg (but see below).
From the models we can access the Log-Likelihood via object$loglig. Note that survreg gives 2 likelihoods, one without covariables (intercept only) and one with covariables. We are interested in the second one.
and extract the AIC-values (extractAIC) for model comparisons:
Like with SAS the weibull-model has the lowest AIC. Not that the Log-Likelihood values differ by a constant of ~920 from SAS. Therefore also the AIC deviates by a constant of ~1840 (=2*920).
From the summary we get the estimates and their standard errors, as well as a Wald test for the coefficients.
Note that the estimate of $0.0602 \pm 0.2566$ for WETWT in the book must be a typo: The standard error is way too big to be statistically significant, so I think the R result is right here.
We see that the concentration has a negative effect on time to death and that bigger fish survive longer.
From our model we can predict for any salt concentration or fish wet weight the median time to death.
To produce Figure 4.12 we first generate some new data for which we want to predict the median time to death. I use expand.grid here, which gives every combination of wet weight from 0 to 1g and four salt concentrations.
Then I use the predict function to predict for all these combinations the median-time-to-death from our model (type='quantile', p = 0.5 specifies the median time).
Afterwards this is plotted using the ggplot2 plotting system.
We can also plot the data and the model into one plot. Here I fit a model only with concentration, for demonstration.
We can plot the Kaplan-Meier curve as shown in previous example. Then add the model as curve using the predict function:
To produce a plot similar to Figure 4.11 we can use the following code:
To fit a gamma distribution we have to use another package: flexsurv.
Usage is (intended) similar to survival. However the loglogistic distribution is not build in. But this can be fixed easily, see ?flexsurvreg:
Comparing the results from flexsurvreg with survreg, we see that the estimates are identical for all models. Therefore I conclude that flexsurv is an alternative when fitting with gamma distribution.
Code and data are available on my github-repo under file name ‘p180’.