Calculating confidence intervals for parameter estimates

[codeanticode]

Hi, I'm using POMP for MLE of the parameters in an SEIR epi model. In order to calculate the CI's for these parameters, I'm applying the Monte Carlo profile (MCAP) method (described in this paper). However, I'm getting a very flat likelihood surface, at least near the optimum parameters, for example:

This results in very side CIs... Do you have any general recommendations/comments about this? I can provide more details if needed, the code is here.

[kingaa]

At first glance, it seems the problem may lie with the profile likelihood calculation, rather than with the MCAP. Why do you think you see what, to a very high precision, is flatness?

[codeanticode]

Yes, it could be the case that's a problem with the likelihood calculation... or how the values are collected in my MCAP code. I will take a closer look and report back. At this point I don't have any good guesses, judging from the plot I posted (for other parameters, they look similar) the likelihood values do change along what seems to be a "cliff"...

[kingaa]

When searching, keep a lookout for discretization artifacts. In particular, it may happen that a parameter is effectively rounded off, so that the likelihood is actually perfectly flat across some range.

[codeanticode]

Ok so I looked closer into it, and I think I found some causes for for the problem of the "flat" loglik surface.

• I was using a distance of 1 to select the values of the parameter to include in the mcap calculation (here), irrespective of the overall range of the parameters. So if the parameter had a range of 0-1, then I was sampling the same points every time, which resulted in a constant surface. The solution to this would be to take a distance equal to (pmax - pmin) / (2 x N-1) where N is the number of subdivisions of the parameter range in the profile. Sounds good?
• Even after doing that, for several parameters I get something like this (format below is value -> max likelihood for that value):
  [1] "0 -> -2851.9764894194"
  [1] "0.214285714285714 -> -1082.06722108034"
  [1] "0.428571428571429 -> -268.769853681114"
  [1] "0.642857142857143 -> -239.769232190653"
  [1] "0.857142857142857 -> -237.077799880441"
  [1] "1.07142857142857 -> -240.499255719695"
  [1] "1.28571428571429 -> -233.518926901806"
  [1] "1.5 -> -233.961626798772"
  [1] "1.71428571428571 -> -234.642974260484"
  [1] "1.92857142857143 -> -231.092063045379"
  [1] "2.14285714285714 -> -240.042777122267"
  [1] "2.35714285714286 -> -237.715835794579"
  [1] "2.57142857142857 -> -244.934494544671"
  [1] "2.78571428571429 -> -251.464917699779"
  [1] "3 -> -257.279890410204"
  So the values are not being rounded off or anything but it's a situation where the maximum is attained quickly once the parameter reaches and remains fairly close to the maximum. So it is kind of a steep hill followed by a plateau going int he direction of 0 -> 3. So for this kind of surfaces the MCAP approximation would not be suitable?
• That last question makes me wonder if the quadratic approximation to the loglik surface is only valid in small neighborhoods in the parameter space, after all, even with all other parameters held constant, sampling over the entire range of the free parameter could yield a surface that's very much non-quadratic, which seems to be what's happening in my case...
• Another question... in the part of the code that generates the model for the profile likelihood calculation, I currently have a hard-coded constant for sub setting: loglik > max(loglik,na.rm=TRUE) - 20 (here). A student of mine at the time wrote that code and I believe it was adapted from a MCAP example that she found online (the JRSI paper only lists the code of the MCAP function), but I lost the link. I wonder if you might have a full example of MCAP calculation so I can compare with my implementation.

[ionides]

How about this R code to fix up the profile.
I think mcap3 below looks reasonable, but there is always room for more tweaking.

prof <- read.csv(header=T,text='
"theta", "loglik"
0, -2851.9764894194
0.214285714285714, -1082.06722108034
0.428571428571429, -268.769853681114
0.642857142857143, -239.769232190653
0.857142857142857, -237.077799880441
1.07142857142857, -240.499255719695
1.28571428571429, -233.518926901806
1.5, -233.961626798772
1.71428571428571, -234.642974260484
1.92857142857143, -231.092063045379
2.14285714285714, -240.042777122267
2.35714285714286, -237.715835794579
2.57142857142857, -244.934494544671
2.78571428571429, -251.464917699779
3, -257.279890410204
')

library(panelPomp)
mcap1 <- mcap(prof$loglik,prof$theta)

plot.profile <- function(mcap1,ylab,xline=2.5,yline=2.5,xlab="",quadratic=FALSE,...){
  if(missing(ylab)) ylab <- "log likelihood"
  ggplot() + geom_point(aes(x = mcap1$parameter, y = mcap1$lp)) + 
    geom_line(mapping = aes(x = mcap1$fit$parameter, y = mcap1$fit$smoothed),
              color = "red") +
    {if(quadratic) geom_line(mapping = aes(x = mcap1$fit$parameter, y = mcap1$fit$quadratic),
                             color = "blue",
                             lwd = 1.25)} +
    labs(x = xlab, y = ylab) +
    geom_vline(xintercept = mcap1$ci, color = "red") + 
    geom_hline(yintercept = max(mcap1$fit$smoothed, na.rm = T) - mcap1$delta, color = "red") +
    theme(panel.border = element_rect(colour = "black", fill=NA))
  }

library(ggplot2)
plot.profile(mcap1)

# the first two points fall off a big likelihood cliff
# this region is best avoided for an MCAP
# we can see what is going on more clearly once we rescale:

plot(loglik~theta,data=prof,ylim=c(-300,max(prof$loglik)))

# now repeat the analysis with them removed:
removed <- c(1,2)
mcap2 <- mcap(prof$loglik[-removed],prof$theta[-removed])
plot.profile(mcap2)

# with relatively few profile points, additional smoothing may be required
# let's increase the default smoothing parameter
mcap3 <- mcap(prof$loglik[-removed],prof$theta[-removed],lambda=1)
plot.profile(mcap3)

[ionides]

panelPomp is at https://github.com/cbreto/panelPomp

but for completeness, here is the mcap code, which is currently a fairly basic implementation.

##' Monte Carlo adjusted profile
##'
##' Given a collection of points maximizing the likelihood over a range
##' of fixed values of a focal parameter, this function constructs 
##' a profile likelihood confidence interval accommodating both
##' Monte Carlo error in the profile and statistical uncertainty present
##' in the likelihood function.
##' 
##' @param lp a vector of profile likelihood evaluations.
##' 
##' @param parameter the corresponding values of the focal parameter.
##'
##' @param confidence the required level of the confidence interval.
##'
##' @param lambda the loess parameter used to smooth the profile.
##'
##' @param Ngrid the number of points to evaluate the smoothed profile.
##'
##' @return mcap returns a list including the smoothed profile, 
##' a quadratic approximation, and the constructed confidence interval.
##' 
##' @author Edward L. Ionides
##'
##' @export
##'
mcap <- function(lp,parameter,confidence=0.95,lambda=0.75,Ngrid=1000){
  smooth_fit <- loess(lp ~ parameter,span=lambda)
  parameter_grid <- seq(min(parameter), max(parameter), length.out = Ngrid)
  smoothed_loglik <- predict(smooth_fit,newdata=parameter_grid)
  smooth_arg_max <- parameter_grid[which.max(smoothed_loglik)]
  dist <- abs(parameter-smooth_arg_max)
  included <- dist < sort(dist)[trunc(lambda*length(dist))]
  maxdist <- max(dist[included])
  weights <- rep(0,length(parameter))
  weights[included] <- (1-(dist[included]/maxdist)^3)^3
  quadratic_fit <- lm(lp ~ a + b, weights=weights,
    data = data.frame(lp=lp,b=parameter,a=-parameter^2)
  )
  b <- unname(coef(quadratic_fit)["b"] )
  a <- unname(coef(quadratic_fit)["a"] )
  m <- vcov(quadratic_fit)
  var_b <- m["b","b"]
  var_a <- m["a","a"]
  cov_ab <- m["a","b"]
  se_mc_squared <- (1 / (4 * a^2)) * (var_b - (2 * b/a) * cov_ab + (b^2 / a^2) * var_a)
  se_stat_squared <- 1/(2*a)
  se_total_squared <- se_mc_squared + se_stat_squared
  delta <- qchisq(confidence,df=1) * ( a * se_mc_squared + 0.5)
  loglik_diff <- max(smoothed_loglik) - smoothed_loglik
  ci <- range(parameter_grid[loglik_diff < delta])
  list(lp=lp,parameter=parameter,confidence=confidence,
    quadratic_fit=quadratic_fit, quadratic_max=b/(2*a),
    smooth_fit=smooth_fit,
    fit=data.frame(
      parameter=parameter_grid,
      smoothed=smoothed_loglik,
      quadratic=predict(quadratic_fit, 
        list(b = parameter_grid, a = -parameter_grid^2)
      )
    ),
    mle=smooth_arg_max, ci=ci, delta=delta,
    se_stat=sqrt(se_stat_squared), 
    se_mc=sqrt(se_mc_squared), 
    se=sqrt(se_total_squared)
  )
}

[codeanticode]

@ionides Thank you so much for your suggestion! Removing the values in the bottom of the likelihood cliff seem to solve the problem (together with the fix I mentioned earlier of using the length of each subdivision it the parameter range). The loglik values I posted were from a quick test run, but I will re-run it with the "production" parameters I will report back. Thanks again!

[codeanticode]

The MCAP calculation of the CIs is working pretty nicely! However, I have some questions regarding the quadratic fit when the maximum of the loglik is near the edge of the parameter range. Consider the following extreme case, where the parameter is bounded to the [0, 1] interval and the loglik maximum is very close to 0. Using a lambda smoothing of 1, the maximum of the quadratic function falls outside the range:

so the CI is undefined. Lowering lambda to 0.75 starts to look better:

but the CI is too wide. Only lowering a lot, to 0.3, I get a reasonable 95% CI of (0.02, 0.12):

I also thought of applying the logit transformation to extend the range to -inf, +inf. So here you have the corresponding plots for lambda = 1, 0.75, and 0.3:

lambda = 1

lambda = 0.75

lambda = 0.3

The final 95% CI applying the logit transform with 0.3 smoothing is a bit tighter at (0.02 0.09)

So I guess my question is what do you recommend to do in these cases? It looks like parameters where the loglik maximum is near the edges of the range cause trouble to the MCAP method, at least you need to tweak the lambda parameter manually until getting a good fit. Applying a variable transformation such as the logit for bounded variables seems to help as well.

But then there is the question (at least for me) of how to interpret the CI after these manipulations. Thank you!

[ionides]

Here, the Monte Carlo noise is small, so it is mostly an issue of using profile likelihood methods, rather than MCAP specifically. To build up a good profile likelihood, you need more evaluation points in the region with likelihood within 1.92 log units of the maximum. A discussion of profile likelihood in the context of ARMA models is given in (https://ionides.github.io/531w21/05/slides.pdf). When the likelihood is far from quadratic, or has boundaries, or is non-smooth, various statistical methods have difficulties. Profile likelihood also needs care in these situations, but I think the worked example in the above slides shows that it can nevertheless be a good choice compared to alternatives.

[codeanticode]

Thanks for pointing to this resource, I will look into it.

[kingaa]

mcap is now provided in the package (as of version 4).

[codeanticode]

That's great news. I eventually managed to use MCAP properly and got reasonable CIs for the parameters in my model. Thank you!