Objectives

To modify and run a maximum likelihood analysis for a nonlinear POMP model of sufficient size and complexity to provide a foundation for the final project.

All the following questions relate to the case study in Chapter 12 of the notes, using iterated filtering to maximize the likelihood for the boarding school influenza model represented by the pomp object bsflu2.

Carry out the following exercises, and write an Rmd file presenting your code and explanations. Use stew or bake to carry out the computations. Scale your computations to a reasonable runtime given the computational resources you have available. Submit to Canvas a zip file with the Rmd file and additional files containing the R objects cached by stew or bake.

Optionally, you can carry out the final version of the computations using the Great Lakes cluster. You can develop and debug your code on a different machine and then run a longer version, with a larger Monte Carlo effort, on Great Lakes. This will be a useful skill for the final project, but not essential.

Question 8.1. Assessing and improving algorithmic parameters.

Use the diagnostic plots on slide 57 of Chapter 12 to form a hypothesis on how you might be able to improve the choice of the algorithmic parameters (i.e., the arguments to the call to mif2 that relate to the operation of the algorithm and are not part of the model). Compare the diagnostic plots with and without your proposed modification, to assess the success of your hypothesis.

First we recall the results from chapter 12 notes

summary(results_global$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -75.56  -75.43  -75.01  -75.01  -74.73  -74.36
plot(mifs_global)

We can see that the effective sample size at time 12 is low compared to other time periods. In this case, we will increase the number of particles for this time period to 70,000, as opposed to 20,000 for all time periods, which is the default for run_level 2.

Np_new = rep(bsflu_Np, length(time(bsflu2, t0=TRUE)))
Np_new[c(12, 13)] = 70000

mifs_global = box_eval("box-eval-8_1-%d.rda", Np_new)
results_global = lik_global_eval("lik_global_eval-8_1-%d.rda")

summary(results_global$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -76.35  -75.05  -74.86  -74.57  -74.26  -71.67
plot(mifs_global)

We see that the effective sample size overall is higher and we achieved a maximum likelihood that is also much higher at -71.67.

Question 8.2. Finding sharp peaks in the likelihood surface.

Even in the small, 3 parameter, boarding school influenza example, it takes a considerable amount of computation to find the global maximum (with values of \(\beta\) around 0.004) starting from uniform draws in the specified parameter box. The problem is that, on the scale on which “uniform” is defined, the peak around \(\beta\approx 0.004\) is very narrow. Propose and implement a more favorable way to draw starting parameters for the global search, which is less dependent on the scale. Your solution may involve taking logarithms, since this converts scale factors to additive factors: ranges that are uniform on a logarithmic scale therefore have good scale invariance properties.

Rather than drawing samples of \(\beta \sim Unif(0.001,0.1)\), we will draw samples \(\beta \sim e^V\) where \(V \sim Unif(\log(0.001), \log(0.01))\).

title = "box-eval-8_2-%d.rda"
stew(file=sprintf(title,run_level),{
  t_global <- system.time({
    mifs_global <- foreach(i=1:bsflu_Nglobal,.combine=c) %dopar% {
      mif2(
        mifs_local[[1]],
        params=c(Beta=exp(runif(1,log(bsflu_box[1,1]),log(bsflu_box[1,2])))
                 ,apply(bsflu_box[2:3,],1,function(x) runif(1,x[1],x[2])),
                 bsflu_fixed_params),
        Np=bsflu_Np,
        Nmif=bsflu_Nmif,
        cooling.fraction.50=bsflu_cooling.fraction.50
      )}
  })
},seed=1270401374,kind="L'Ecuyer")


results_global = lik_global_eval("lik_global_eval-8_2-%d.rda")

summary(results_global$logLik,digits=5)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -76.87  -75.24  -74.69  -74.78  -73.85  -73.24
plot(mifs_global)

Comparing the maximum likelihood we obtained with this new sampling to notes 12, we find that we have improved the likelihood from -74.36 to -73.24. We also notice that the range of starting values for \(\beta\) is narrower than the range obtained from the previous sampling scheme.

Question 8.3. Construct a profile likelihood.

How strong is the evidence about the specific value of the contact rate, \(\beta\), given the bsflu2 model and data? Use mif2 to construct a profile likelihood and corresponding approximate confidence interval for this parameter.

Below we construct the profile likelihood for \(\beta\).

# Construct parameter box 
ranges = apply(results_global[,3:5], 2,range)
starts = profileDesign(
  Beta = exp(seq(log(bsflu_box[1,1]), log(bsflu_box[1,2]), length = 50)),
  lower = ranges[1,-1],
  upper = ranges[2,-1],
  nprof = 50
)

title = "beta_prof.rda"
stew(file=title,{
  t_global <- system.time({
    beta_prof = foreach(i=1:nrow(ranges), 
                     .combine=rbind, .packages=c("pomp")) %dopar% {
               
               mif2(
                 mifs_local[[1]],
                 params= c(unlist(starts[i,]), bsflu_fixed_params),
                 rw.sd=rw.sd(mu_IR=bsflu_rw.sd, 
                             rho = bsflu_rw.sd),
                 Np=1000,
                 Nmif=20
                 ) %>%
                 mif2() -> mf
               
               replicate(5, 
                         mf %>% pfilter() %>% logLik()
               ) %>%
                 logmeanexp(se=TRUE) -> ll
               
               data.frame(as.list(coef(mf)),loglik=ll[1],loglik.se=ll[2])
               }
  })
},seed=1270401374,kind="L'Ecuyer")

# Construct CI for Beta
l_prof = max(beta_prof$loglik)
l_hat = l_prof - 0.5 *qnorm(0.975)^2
ind = which(beta_prof$loglik >= l_hat)
lower_CI = beta_prof$Beta[min(ind)]
upper_CI = beta_prof$Beta[max(ind)]

beta_prof %>%
  group_by(Beta) %>%
  filter(rank(-loglik)<=2) %>%
  ungroup() %>%
  ggplot(aes(x=Beta,y=loglik,
             ymin=loglik-2*loglik.se,ymax=loglik+2*loglik.se))+
  geom_point()+
  geom_errorbar()+
  geom_smooth(method="loess",span=0.2)+
  scale_x_log10() + 
  geom_hline(aes(yintercept=l_prof), linetype="dashed") +
  geom_hline(aes(yintercept=l_hat), linetype="dashed") +
  geom_vline(aes(xintercept=lower_CI), linetype="dashed") +
  geom_vline(aes(xintercept=upper_CI), linetype="dashed") + 
  ggtitle(TeX("Profile Likelihood for $\\beta$"))

We find that the approximate 95% confidence interval for \(\beta\) is [0.00339, 0.00625].

Question 8.4. This feedback question is worth credit.

  1. Explain which parts of your responses above made use of a source, meaning anything or anyone you consulted (including classmates or office hours) to help you write or check your answers. All sources are permitted. To encourage responsible use of these sources while maintaining class integrity, we require a response to this question, even if this may occasionally just say that you worked out everything entirely by yourself. See the syllabus for additional information on grading.
  2. How long did this homework take? Report on any technical difficulties that arose.

(a) For most of the homework I relied on Lecture 12 notes. For constructing the profile likelihood in 8.3, I found sample code from kingaa’s website.

(b) The homework took me about 12 hours to complete. Running the code and debugging the results took the longest time.