Standard errors from flexsurvreg

Question

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.

Answer 1

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

Which is what expFit$res gave us.

Why use the log scale?

To get this answer back to statistics instead of just R, why do we calculate the confidence interval on the log scale? Can’t we just use the standard error returned by expFit?

Well, we can (in this case):

> expFit$res[1] + c(-1,1)*z_alpha*expFit$res[4]
[1] 0.005542589 0.009139678

This CI is very similar to the log-based CI. The reason is the large sample size. The actual formula for estimating the rate parameter in an exponential distribution is actually very simple; it’s the number of events divided by the total amount of time:

> with(testPatients, sum(status)/sum(time))
[1] 0.007341133

If none of the events were censored, this is 1 divided by the average time to event, which makes sense intuitively. For example, if the time to event is on average 15 minutes (1/4 of an hour), you would expect the events to happen with a rate of 4 per hour.

The times were assumed to be i.i.d. exponential with a fixed ‘rate’ parameter. These time variables are thus all right-skewed. The sum of many such independent time variables has a gamma distribution – or approximately a normal distribution for large n (per the central limit theorem).

But we’re looking at a constant divided by this sum, i.e. divided by something with a distribution that is a (right-skewed) gamma (for small $n$) or approximately symmetric (a normal distribution in the limit). It’s right-skewed in both cases. For large $n$, the amount of skewness is small (the value of the rate parameter is actually not important for the skewness).

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.