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

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

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}\]

  1. 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.
  2. \(\hat\lambda^{[2]}\) is the log mean exp of the total log likelihood for all units.
  3. \(\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
CC-BY_NC

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