Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC

Produced using R version 3.5.3 and pomp2 version 2.0.13.1.

Introduction

This document shows how the results of He et al. (2010) can be reproduced. Specifically, it implements the model of that paper and computes the likelihood of the data under that model. The intentions are to:

Whilst He et al. (2010) examined measles in 20 towns, we will focus on London alone. The full 20 city dataset is available on the pomp website and the MLEs obtained by He et al. (2010) for the full 20 towns are included in the R code accompanying this document.

Preliminaries

To get started, we must install pomp if it is not already installed. The following commands install pomp from source (on github).

library(devtools)
install_github("kingaa/pomp")

If we were to run into trouble in the above, we would consult the package website.

Now we’ll load the packages we’ll need, and set the random seed, to allow reproducibility.

set.seed(594709947L)
library(ggplot2)
library(plyr)
library(tidyverse)
library(pomp2)
stopifnot(packageVersion("pomp2")>="2.0.9.1")

Data and covariates

Now we’ll load the data and covariates. The data are measles reports from 20 cities in England and Wales, digitized from the printed public records and graciously provided by Prof. Bryan Grenfell. We also have information on the population sizes and birth-rates in these cities; we’ll treat these variables as covariates.

load("twentycities.rda")
measles %>% 
  mutate(year=as.integer(format(date,"%Y"))) %>%
  filter(town==TOWN & year>=1950 & year<1964) %>%
  mutate(time=(julian(date,origin=as.Date("1950-01-01")))/365.25+1950) %>%
  filter(time>1950 & time<1964) %>%
  select(time,cases) -> dat
demog %>% filter(town==TOWN) %>%
  select(-town) -> demog

Let’s plot the data and covariates.

dat %>% ggplot(aes(x=time,y=cases))+geom_line()