Standard errors from flexsurvreg

I'm using the flexsurv package (in R) to fit an exponential distribution to the veteran dataset in survival. I'd like to estimate the rate parameter of the exponential distribution, as well as an estimate of the standard error of the MLE.

Here's some code (and output):

library(flexsurv)

testPatients <- subset(veteran,trt==2)
testPatSurvObj <- with(data=testPatients,expr={Surv(time,status,type="right")})

### exponential fit
expFit <- flexsurvreg(testPatSurvObj ~ 1, dist="exp")

expFit$res ### est L95% U95% ### rate 0.007341177 0.005746003 0.009379195 expFit$res.t
###           est      L95%      U95%
### rate -4.914256 -5.159251 -4.669261


Now, I see that expFit$cov is the covariance matrix of the parameter estimates, with positive parameters on the log scale, but I don't see how I can use this fact to transform expFit$cov into something I can use to make a (say) normal-based 95% confidence interval.

Getting the confidence interval directly

Firstly, printing an flexsurvreg object (or its res element) already shows the 95% confidence interval:

> expFit
Estimates:
est       L95%      U95%      se
rate  0.007341  0.005746  0.009379  0.000918


Reproducing the the confidence interval manually

So I guess you’re asking how to reproduce the above CI manually using the estimated covariance matrix for the parameter estimator. The help page (?flexsurvreg) says that ‘Parameters defined to be positive are estimated on the log scale.’ So we get log-transformed values when we extract the estimates:

> lrate = coef(expFit)
> lrate
[1] -4.914262


To get the actual estimated rate parameter, we need to exponentiate this number:

> rate = exp(lrate)
> rate
[1] 0.007341133


Now to generate a confidence interval, we just need the standard error, and then we can use the normal approximation (Wald confidence interval). The standard error is:

> se = sqrt(vcov(expFit))
> se
rate
rate 0.125


So the confidence interval for the log rate is:

> z_alpha = qnorm(1-.05/2)
> ci = lrate + c(-1,1)*z_alpha*se
> ci
[1] -5.159258 -4.669267


To get a 95% CI for the actual rate parameter, we need only exponentiate the two CI limits:

> exp(ci)
[1] 0.005745964 0.009379146


The effect of the number of observations in practice

In the example data we had 68 observations (of which 4 were censored), which is quite a lot. The distributions of 1) the mean of 68 independent exponential observations (let’s just ignore the censored ones for now) with rate 0.0073, 2) the rate estimate (68 divided by this mean) and 3) the logarithm of this rate estimate look approximately like this:

The R code for this image is:

library(MASS)
sim = function(n, rate) {
x = replicate(10^5, mean(rexp(n, rate)))
par(mfrow = c(1,3), mai = c(0.9,0.4,0.2,0))
truehist(x, col = "brown", border = NA)
truehist(n / x, col = "brown", border = NA)
truehist(log(n / x), col = "brown", border = NA)
}
sim(n = 68, rate = rate)


We see that there is some skewness in the distribution for the rate estimator (the middle panel): there’s a tail to the right. But the distribution is approximately symmetric (and normal), so the estimate based on the normal approximation without a log-transform should be quite good. Still, we see from the right-most panel that it is better to use the log-transform even here.

If we had only 15 observations, the picture would change:

sim(n=15, rate=rate)


Here we have very much skewness in the middle panel, and we need to use the log-transform (right panel) to get a good CI.

Summary and closing words

We calculate the (approximate) confidence interval for the rate parameter on the log scale because the sample distribution of the estimator is more symmetrical and normal on this scale.

Note that this is common practice for parameters that must be positive, since the sample distribution of the estimators are then typically right-skewed. (When the standard error is large, using the normal approximation without log-transforming first can in some cases even result in the lower confidence limit for the parameter being less than 0.)

Also note there are also other transformations besides the log transformation that can be useful for transforming estimators to be more symmetrical/normal, and that can be more suitable for some distributions/parameters. But the log transform is the most commonly used.