Worked Solutions to the Exercises

Exercise 8.1

A PanelPOMP model with all parameters unit-specific is identical to a collection of independent POMP models. There may not be anything to be gained by using the PanelPOMP framework in this case, except perhaps that the panelPomp package can help with the book-keeping for working with the collection of panels.

One can expect an overhead from using panelPomp software when the simpler pomp software is sufficient.

PanelPOMP models with all unit-specific parameters can be used to test panelPomp methodology against pomp methodology since the latter also applies in this special case.

Exercise 8.2

We can first check the list of help topics:

library(help=panelPomp)

as                      Coercing 'panelPomp' objects as 'list',
                        'pompList' or 'data.frame'
contacts                Contacts model
get_dim                 Get single column or row without dropping names
mif2                    PIF: Panel iterated filtering
panelGompertz           Panel Gompertz model
panelGompertzLikelihood
                        Likelihood for a panel Gompertz model via a
                        Kalman filter
panelPomp               Constructing 'panelPomp' objects
panelPomp-package       Inference for PanelPOMPs (Panel Partially
                        Observed Markov Processes)
panelPomp_methods       Manipulating 'panelPomp' objects
panelRandomWalk         Panel random walk model
panel_loglik            Handling of loglikelihood replicates
panel_logmeanexp        Log-mean-exp for panels
params                  Convert to and from a 'panelPomp' object
                        'pParams' slot format and a one-row
                        'data.frame'
pfilter                 Particle filtering for panel data
plot                    panelPomp plotting facilities
simulate                Simulations of a panel of partially observed
                        Markov process

And then check out the promising ones:

?panelPomp_methods

Exercise 8.3

Since a PanelPOMP is a collection of independent POMP models, they can each be filtered separately.
The pfilter method for class panelPomp does exactly this, after extracting the parameter vector for each unit from the shared and unit-specific parameters belonging to the panelPomp object.

contacts <- contacts()
pf <- pfilter(contacts,Np=100)
class(pf)

[1] "pfilterd.ppomp"
attr(,"package")
[1] "panelPomp"

class(unitobjects(pf)[[1]])

[1] "pfilterd_pomp"
attr(,"package")
[1] "pomp"

Exercise 8.4

It turns out that the SMC algorithm implemented by pfilter() gives an unbiased Monte Carlo estimate of the likelihood.
For inferential purposes, we usually work with the log likelihood.
Due to Jensen’s inequality, SMC has a negative bias as an estimator of the log likelihood, i.e., it systematically underestimates the log likelihood.
Usually, the higher the Monte Carlo variance on the likelihood, the larger this bias.
Thus, lower Monte Carlo variance on the log likelihood is commonly associated with higher estimated log likelihood.
Heuristically, products propagate error rapidly. Averaging Monte Carlo replicates over units before taking a product reduces this error propagation. It does not lead to bias in the likelihood estimate since independence over units and Monte Carlo replicates insures that the expected product of the averages is the expected average of the products.

Some notation for a more formal investigation

Let $\hat\lambda_u^{(k)}$ be the $k$ th replication of the Monte Carlo log likelihood evaluation for unit $u$ .
Let $\hat L_u^{(k)}=\exp\big\{\hat\lambda_u^{(k)}\big\}$ be the corresponding likelihood.
Let $\hat\lambda^{(k)}=\sum_{u=1}^U \lambda_{u}^{(k)}$ be an estimate the log likelihood of the entire data based on replication $k$ .
Let $\hat L^{(k)}=\exp\big\{\hat\lambda^{(k)}\big\}$ be the corresponding estimate of the likelihood.
Different possible estimates of the actual log likelihood $\lambda=\sum_{u=1}^U \lambda_u$ are

$\begin{eqnarray} \hat\lambda^{[1]} &=& \frac{1}{K}\sum_{k=1}^K \hat\lambda^{(k)} \\ \hat\lambda^{[2]} &=& \log \left( \frac{1}{K}\sum_{k=1}^K \exp \big\{\hat \lambda^{(k)} \big\} \right) \\ \hat\lambda^{[3]} &=& \sum_{u=1}^U\frac{1}{K}\sum_{k=1}^K \hat\lambda^{(k)}_u \\ \hat\lambda^{[4]} &=& \sum_{u=1}^U \log \left( \frac{1}{K}\sum_{k=1}^K \exp\big\{\hat \lambda^{(k)}_u \big\} \right) \end{eqnarray}$

Check that $\hat\lambda^{[1]}$ and $\hat\lambda^{[3]}$ are equal. However, they are inconsistent, since $\hat\lambda^{(k)}_u$ is a biased estimate of $\lambda_u$ and the bias does not disappear when we take an average over replicates.
$\hat\lambda^{[2]}$ is the log mean exp of the total log likelihood for all units.
$\hat\lambda^{[4]}$ is the sum of the log mean exp for each unit separately.

To compare variances, it is convenient to move back to the likelihood scale:

$\begin{eqnarray} \hat L^{[2]} &=& \frac{1}{K}\sum_{k=1}^K \prod_{u=1}^U L^{(k)} \\ \hat L^{[4]} &=& \prod_{u=1}^U \frac{1}{K}\sum_{k=1}^K \hat L^{(k)}_u \end{eqnarray}$

Bretó, Ionides, and King (2020) showed that $\hat L^{[4]}$ is smaller than $\hat L^{[2]}$ .
In a limit where $U$ and $K$ both grow together, $\hat L^{[4]}$ can be stable while $\hat L^{[2]}$ increases to infinity.

Produced with R version 4.3.3, pomp version 5.6, and panelPomp version 1.1.0.1.

References

Bretó C, Ionides EL, King AA (2020). “Panel Data Analysis via Mechanistic Models.” J Am Stat Assoc, 115(531), pre–published online. https://doi.org/10.1080/01621459.2019.1604367.

Top of this document
Previous page
Back to the lesson
Course homepage
Acknowledgments

Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.