The **R** codes in this script are available for download. This work is licensed under the Creative Commons attribution-noncommercial license. Please share and remix noncommercially, mentioning its origin.

This document shows how the results of He, Ionides, and King (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:

- facilitate new data analyses based on variants of this model by showing how it is implemented in
**pomp** - establish baseline likelihoods for these data against which future analyses can be compared
- illustrate the construction of a complex model in
**pomp**

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.

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(tidyverse)
library(pomp)
```

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
|> filter(town==TOWN) |>
demog select(-town) -> demog
```

Let’s plot the data and covariates.

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

```
|>
demog pivot_longer(-year) |>
ggplot(aes(x=year,y=value))+geom_point()+
facet_wrap(~name,ncol=1,scales="free_y")
```

Now we’ll smooth the covariates, storing the results in a new data frame. Note that, for reasons discussed in He et al. (2010), we lag the birthrate by 4 yr.

```
|>
demog reframe(
time=seq(from=min(year),to=max(year),by=1/12),
pop=predict(smooth.spline(x=year,y=pop),x=time)$y,
birthrate=predict(smooth.spline(x=year+0.5,y=births),x=time-4)$y
covar ) ->
```

View the covariates:

```
plot(pop~time,data=covar,type='l')
points(pop~year,data=demog)
```