Comment on "Under-reported data analysis with INAR-hidden Markov chains"

TL;DR

This paper critiques a previous model for disease reporting, showing it can be reformulated as a different INAR model with a geometric lag, affecting how reporting probabilities are estimated.

Contribution

It demonstrates that a Poisson INAR(1) model with underreporting is equivalent to a fully observed INAR(inf) model with geometric lag, highlighting the impact of lag structure assumptions.

Findings

Poisson INAR(1) with underreporting is equivalent to an INAR(inf) with geometric lag.

Estimated reporting probabilities depend on the lag structure.

Reformulation affects interpretation of disease reporting models.

Abstract

In Fernandez-Fontelo et al (Statis. Med. 2016, DOI 10.1002/sim.7026) hidden integer-valued autoregressive (INAR) processes are used to estimate reporting probabilities for various diseases. In this comment it is demonstrated that the Poisson INAR(1) model with time-homogeneous underreporting can be expressed equivalently as a completely observed INAR(inf) model with a geometric lag structure. This implies that estimated reporting probabilities depend on the assumed lag structure of the latent process.

Figures

Equations

$$X_t=\alpha \circ X_{t-1}+W_t$$

$${\tilde{X}}_t=\{\begin{array}{ll}{X}_t& \text{ with probability }1-\omega \\ q\circ X_t& \text{ with probability }\omega ,\end{array}$$

$${\tilde{X}}_t=q\circ X_t.$$

$$X_t$$

$$({\alpha }_1\circ X_t,\dots ,{\alpha }_p\circ X_t)\mid X_t$$

$$B_t\mid X_t$$

$$A_t^{(j)}$$

$${\alpha }_i\circ X_t$$

$${\alpha }_i=\beta {\gamma }^{i-1},i=1,2,\dots$$

$$X_t$$

$$B_t\mid X_t$$

$$A_t^{(j)}$$

$${\alpha }_i\circ X_t$$

$${\lambda }_Y=\frac{{\lambda }_X(1-{\gamma }_X)\frac{q_X}{q_Y}}{1-{\gamma }_X-\left(1-\frac{q_X}{q_Y}\right){\beta }_X};{\beta }_Y={\beta }_X\frac{q_X}{q_Y};{\gamma }_Y={\gamma }_X+\left(1-\frac{q_X}{q_Y}\right){\beta }_X.$$

$$X_t$$

$${\tilde{X}}_t$$

$$Y_t=0.17\circ Y_{t-1}+0.06\circ Y_{t-2}+0.02\circ Y_{t-3}+0.01\circ Y_{t-4}+\cdots +W_t(0.82)$$

$$X_t=\alpha \circ X_{t-1}+W_t;W_t\stackrel{\text{iid}}{\sim }\text{Pois}({\lambda }_X)$$

$${\tilde{X}}_t={\tilde{U}}_{t\cdot }+{\tilde{V}}_{t\cdot }.$$

$${\tilde{U}}_{t\cdot }=\sum _{i=0}^{\infty }{\tilde{U}}_{t,i}.$$

$${\tilde{V}}_{t\cdot }=\sum _{i=1}^{\infty }{\tilde{V}}_{t,i}.$$

$${\tilde{U}}_{t,\cdot }\stackrel{\text{iid}}{\sim }\text{Poisson}({\lambda }_Y);{\lambda }_Y=\sum _{i=0}^{\infty }\{{\alpha }_X(1-q_X){\}}^i{q}_X{\lambda }_X=\frac{q_X{\lambda }_X}{1-{\alpha }_X(1-q_X)}.$$

$${\beta }_Y={\alpha }_X{q}_X,{\gamma }_Y={\alpha }_X(1-q_X).$$

$${\tilde{B}}_t\mid {\tilde{X}}_t$$

$$A_t^{(j)}$$

$${\tilde{V}}_{t+i,i}$$

$${\tilde{X}}_t=\sum _{i=1}^{\infty }{\tilde{V}}_{t,i}+{\tilde{U}}_{t,\cdot }$$

$${\lambda }_Y$$

$${\beta }_Y$$

$${\gamma }_Y$$

Full text

Epidemiology, Biostatistics and Prevention Institute, University of Zurich,

Hirschengraben 84, 8001 Zurich, Switzerland

johannes.bracher@uzh.ch

Bibliographic note: This is the pre-reviewing version of a letter published in Statistics in Medicine 38(5), 893–898:

https://onlinelibrary.wiley.com/doi/full/10.1002/sim.8032

Fernández-Fontelo et al published a reply which raises some interesting further points:

https://onlinelibrary.wiley.com/doi/10.1002/sim.8033

After a 12 month embargo period the accepted version will be made available on arxiv.

1 Introduction

We read with great interest the article by Fernández-Fontelo et al [3] who discuss underreporting in count time series models, a problem which is encountered in many real-world settings. In case studies, reporting probabilities for several diseases are estimated under the assumption of a hidden INAR(1) process. In our comment we develop this idea further and point out an identifiability problem which arises when the INAR(1) assumption is relaxed. Specifically it is shown that it is not possible to distinguish between time-homogeneous underreporting and a geometric lag structure in Poisson INAR models. Estimation of reporting probabilities from time series data thus relies on the correct specification of the lag structure of the latent process.

The Poisson INAR(1) model [1] with parameters $\lambda>0$ and $0\leq\alpha<1$ is defined as

[TABLE]

where the $W_{t}$ independently follow a $\text{Poisson}(\lambda)$ distribution. The operator $\circ$ denotes binomial thinning, i.e. $\alpha\circ X_{t-1}=\sum_{i=1}^{X_{t-1}}Z_{i}$ with $Z_{i}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Bern}(\alpha)$. The thinning operations are assumed to be independent of each other and of $\{W_{t}\}$; further, the thinning operations at each $t$ and $W_{t}$ are assumed to be independent of $X_{t-1},X_{t-2},\dots$. Fernández-Fontelo et al [3] introduce an underreported version of this model. However, they assume that counts are only subject to underreporting with a certain probability $\omega$; otherwise reporting is assumed to be complete. The observed process is thus

[TABLE]

where $0<q\leq 1$ is a reporting probability and, given $X_{t}$, $\tilde{X}_{t}$ is independent of the past. The assumption that reporting is 100% complete during some periods seems strong in many contexts. And indeed, in two case studies (weekly number of human papillomavirus cases in Girona, Spain; annual deaths from mesothelioma in Great Britain) the estimates $\hat{\omega}$ are close to 1 with confidence intervals including this value. This indicates that $\omega=1$, i.e. time-homogeneous underreporting, is an important special case. In the following the focus will thus be on models where, instead of (2), the reporting process is

[TABLE]

2 Interplay of underreporting and geometric lags in INAR models

We now generalize the INAR(1) from (1) to an INAR($\infty$) model with geometric lags and show how closely this aspect is related to underreporting. The starting point is the INAR($p$) model introduced by Alzaid and Al-Osh [2] which is defined as

[TABLE]

with $\sum_{i=1}^{p}\alpha_{i}<1$.

This INAR($p$) model does not allow for $p=\infty$ as the multinomial distribution requires a finite number of categories. We therefore reformulate the multinomial distribution in (5) as (compare e.g. [5], p. 33)

[TABLE]

where $I$ is the indicator function. In the interpretation given in Weiß [7] (p. 46), $B_{t}$ is the number of individuals (among the $X_{t}$ from time $t$) which will be renewed during $t+1,\dots,t+p$. The variables $A_{t}^{(j)},j=1,\dots,B_{t}$ are the respective waiting times. This formulation extends more easily to the case $p=\infty$, specifically consider a geometric lag structure

[TABLE]

with $0<\beta<1-\gamma<1$. As $\sum_{i=1}^{\infty}\alpha_{i}=\beta/(1-\gamma)$ we get $\alpha_{i}/\left(\sum_{k=1}^{\infty}\alpha_{k}\right)=(1-\gamma)\gamma^{i-1}$, i.e. the $A_{t}^{(j)}$ follow a geometric distribution with parameter $1-\gamma$ and support $\{1,2,\dots\}$. The INAR($\infty$) process $\{X_{t}\}$ with parameters $\beta$ and $\gamma$ can thus be defined as

[TABLE]

This ensures, like the multinomial distribution in (5), that $(\alpha_{i}\circ X_{t})\mid X_{t}\sim\text{Bin}(X_{t},\alpha_{i})$ while $\sum_{i=1}^{\infty}\alpha_{i}\circ X_{t}\leq X_{t}$. Again for each $t$ it is assumed that $W_{t}$ and $\alpha_{i}\circ X_{t}\mid X_{t},i=1,\dots$ are independent of the past history, i.e. $X_{t-k}$ and $\alpha_{j}\circ X_{t-k}$ for $j,k=1,2,\dots$. Such geometric lag structures are closely linked to underreporting as the following can be shown (proofs in the Appendix):

1. (A)

Consider a Poisson INAR(1) process $\{X_{t}\}$ with parameters $\lambda_{X},\alpha_{X}$ and the underreported process $\{\tilde{X}_{t}=q_{X}\circ X_{t}\}$. There is a Poisson INAR($\infty$) process $\{Y_{t}\}$ of type (10)–(13) which is is equivalent to $\{\tilde{X}_{t}\}$; its parameters are $\lambda_{Y}=\lambda_{X}q_{X}/\{1-\alpha_{X}(1-q_{X})\};\beta_{Y}=\alpha_{X}q_{X};\gamma_{Y}=\alpha_{X}(1-q_{X})$. 2. (B)

Consider a Poisson INAR($\infty$) process $\{X_{t}\}$ of type (10)–(13) with parameters $\lambda_{X},\beta_{X}$ and $\gamma_{X}$. The underreported process $\{\tilde{Y}_{t}=q_{Y}\circ Y_{t}\},q_{Y}=\beta_{X}/(\beta_{X}+\gamma_{X})$ where $Y_{t}$ is an INAR(1) process with parameters $\lambda_{Y}=\lambda_{X}(\beta_{X}+\gamma_{X})(1-\gamma_{X})/\beta_{X};\alpha_{Y}=\beta_{X}+\gamma_{X}$ is equivalent to $\{X_{t}\}$. 3. (C)

More generally, consider a Poisson INAR($\infty$) process $\{X_{t}\}$ with parameters $\lambda_{X},\beta_{X},\gamma_{X}$ and $\{\tilde{X}_{t}=q_{X}\circ X_{t}\}$. For any $q_{Y}\in[q_{X}\beta_{X}/(\beta_{X}+\gamma_{X}),1]$ there is a Poisson INAR($\infty$) process $\{Y_{t}\}$ so that $\{\tilde{Y}_{t}=q_{Y}\circ Y_{t}\}$ is equivalent to $\{\tilde{X}_{t}\}$. The parameters of $\{Y_{t}\}$ are

[TABLE]

For each underreported INAR($\infty$) process and, as a special case, each underreported INAR(1) process, there are thus many different INAR($\infty$) representations. Each of them features a different combination of reporting probability $q_{X}$ and decay factor $\gamma_{X}$ for the autoregressive parameters. In Fernández-Fontelo et al [3] identifiability is ensured by the assumption that the latent process is indeed INAR(1), i.e. $\gamma_{X}=0$. Without substantial prior knowledge to justify a precise value of $\gamma_{X}$, however, the reporting probability $q_{X}$ cannot be estimated.

3 Application to human papillomavirus cases in Girona, Spain

For illustration of the argument we revisit the analysis of reported human papillomavirus cases in Girona, Spain (2010–2014) from Fernández-Fontelo et al [3]. The authors assume the latent model (1) with reporting process (2), but as $\hat{\omega}=0.92$ with a confidence interval from 0.78 to 1.07 they suggest that the simpler version (3) could be used as well. We will thus simplifyingly pretend that their parameter estimates come from this simpler model (even though the model then cannot accomodate overdispersion; this could be addressed via a different immigration distribution). The estimated data generating process (eq. (22)–(23) in [3]) is then

[TABLE]

The most interesting result from a public health perspective is the estimated reporting probability of 0.33 as it directly translates to an estimate of the unobserved disease burden. Using statement (C) from Section 2, however, we can pick any $q_{Y}\in[0.33,1]$ and obtain the parameters $\hat{\beta}_{Y},\hat{\gamma}_{Y},\hat{\lambda}_{Y}$ of an INAR($\infty$) process $\{Y_{t}\}$ so that $\{\tilde{Y}_{t}=q_{Y}\circ Y_{t}\}$ is equivalent to $\{\tilde{X}_{t}\}$. In Figure 1these parameters are diplayed as functions of $q_{Y}$.

So if we relax the INAR(1) assumption from the original article and assume a latent INAR($\infty$) process, the data are equally compatible with a whole range of reporting probabilities. The most extreme re-formulation of $\{\tilde{X}_{t}\}$ from equations (14), (15) is a completely observed INAR($\infty$) process $\{Y_{t}\}$ with

[TABLE]

where the lagged terms after $p=4$ are cut off as the respective parameters become very small. In certain settings generation time distributions of infectious diseases may give us an idea about an appropriate lag structure and thus value of $\gamma_{Y}$, making the model identifiable. The period of communicability of HPV is unknown, but likely to be at least as long as the persistence of lesions [4]. Development of lesions is assumed to take 2–3 months in most cases. For weekly data a more spread-out lag structure may thus be more appropriate than an INAR(1) specification.

Appendix: Derivation of statements (A), (B), (C) from Section 2

Consider a Poisson INAR(1) process $\{X_{t}\}$ with

[TABLE]

and the underreported process $\{\tilde{X}_{t}=q_{X}\circ X_{t}\}$. It is now shown that, as stated in (A), $\{\tilde{X}_{t}\}$ is equivalent to an INAR($\infty$) process $\{Y_{t}\}$ of type (10)–(13) with parameters $\lambda_{Y},\beta_{Y},\gamma_{Y}$. The argument is easiest understood when expressed in terms of the survival interpretation of $\{\tilde{X}_{t}\}$:

• (i)

New individuals can be born at each time step $t$; their number follows a Poisson distribution with rate $\lambda_{X}$.

• (ii)

Individuals already present at $t$ have a probability $\alpha_{X}$ of still being alive at time $t+1$.

• (iii)

Alive individuals are observed with probability $q_{X}$ at each time step.

All births, deaths and observation events are assumed to be independent. $X_{t}$ is the number of individuals alive at time $t$, $\tilde{X}_{t}$ is the number of those who are observed. Now denote by $\tilde{U}_{t\cdot}$ the number of individuals observed in $t$ which have not been observed previously, and by $\tilde{V}_{t\cdot}$ the number of individuals observed in $t$ which have already been observed at a previous time point so that

[TABLE]

The term $\tilde{U}_{t\cdot}$ can be further decomposed by when the individuals were born; denoting by $\tilde{U}_{t,i},i=0,1,\dots$ the number of individuals first observed in $t$ and born in $t-i$ one gets

[TABLE]

Similarly, $\tilde{V}_{t\cdot}$ is decomposed by when the individuals were last observed; denoting by $\tilde{V}_{t,i},i=1,2,\dots$ the individuals last observed in $t-i$ and observed again in $t$ this leads to

[TABLE]

Figure 2 illustrates the definition of $\tilde{U}_{t,i}$ and $\tilde{V}_{t,i}$ with a simple example.

Obviously an individual born in $t$ can only be observed for the first time in at most one out of $\{t,t+1,\dots\}$, i.e. $W_{t}$ is split up into $U_{t,0},U_{t+1,1},\dots$ and a part of individuals which is never observed. The probability that an individual born in $t$ is first observed in $t+i,i=0,1,\dots$ is $\alpha_{X}^{i}(1-q_{X})^{i}q_{X}$ (the individual has to survive $i$ times, stay unobserved $i$ times and finally be observed once). The splitting property of the Poisson distribution ([6], 53) and the independence between the $W_{t},t=0,\pm 1,\dots$ then imply that all $\tilde{U}_{t,i},i=0,1,\dots;t=0,\pm 1,\dots$ independently follow Poisson distributions. Their rates depend on $i$, specifically $\tilde{U}_{t,i}\sim\text{Poisson}(\alpha_{X}^{i}(1-q_{X})^{i}q_{X}\lambda_{X})$. Consequently, the $\tilde{U}_{t,\cdot}$, too, are independent of each other. As sums of independent Poisson random variables they likewise follow a Poisson distribution ([6], 14):

[TABLE]

Similarly, an individual observed in $t$ can only be observed next in at most one out of $\{t+1,t+2,\dots\}$. The $\tilde{X}_{t}$ observed individuals are thus split up into $\tilde{V}_{t+1,1},\tilde{V}_{t+2,2},\dots$ and a part of individuals which is never observed again. The probability that an individual is observed next in $t+i,i=1,\dots$ is $\alpha_{X}^{i}(1-q_{X})^{i-1}q_{X}$ (the individual has to survive $i$ times, stay unobserved $i-1$ times and finally be observed once) or $\beta_{Y}\gamma_{Y}^{i-1}$ with

[TABLE]

The probability that the individual will be observed again at all is $\sum_{i=1}^{\infty}\beta_{Y}\gamma_{Y}^{i-1}=\beta_{Y}/(1-\gamma_{Y})$. Consequently, under the condition that the $j$-th of the $X_{t}$ individuals will be observed again, the waiting time $A_{t}^{(j)}$ until this occurs has probability mass function $\text{Prob}(A_{t}^{(j)}=i)=\beta_{Y}\gamma_{Y}^{i-1}/\{\beta_{Y}/(1-\gamma_{Y})\}=\gamma_{Y}^{i-1}(1-\gamma_{Y}),i=1,\dots$. It thus follows a geometric distribution with parameter $1-\gamma_{Y}$. Denoting the number of individuals observed in $t$ which will be observed again by $\tilde{B}_{t}$ we can thus write

[TABLE]

Combining this with equations (16), (18) and (19) to

[TABLE]

one can see that $\{\tilde{X}_{t}\}$ indeed follows the form (10)–(13) with parameters $\lambda_{Y},\beta_{Y},\gamma_{Y}$ (the $\tilde{V}_{t,i}$ in equation (24) correspond to the $\alpha_{Y,i}\circ X_{t-i}$ from equation (10) and $\tilde{U}_{t\cdot}$ corresponds to $W_{t}$). Also it is clear that both $\tilde{U}_{t,\cdot}$ and $\tilde{V}_{t+i,i}\mid\tilde{X}_{t}$ are independent of the past history of the observed process $\{\tilde{X}_{t}\}$, i.e. $\tilde{X}_{t-k}$ and $\alpha_{j}\circ\tilde{X}_{t-k},i,j,k=1,2,\dots$

Statement (B) follows directly as it is easily verified using (A) that the given $\{\tilde{Y}_{t}\}$ and $\{X_{t}\}$ are equivalent. The restriction $0<\beta_{X}<1-\gamma_{X}<1$ ensures that $\lambda_{Y}>0$ and $0<\alpha_{Y},q_{Y}\leq 1$ so that $\{Y_{t}\}$ and $\{\tilde{Y}_{t}\}$ are well-defined.

Statement (C) follows from statements (A) and (B) in the following way. Consider an INAR($\infty$) process $\{X_{t}\}$ and $\{\tilde{X}_{t}=q_{X}\circ X_{t}\}$. Using (B) an INAR(1) process $\{Z_{t}\}$ with parameters $\lambda_{Z}=\lambda_{X}(\beta_{X}+\gamma_{X})(1-\gamma_{X})/\beta_{X},\alpha_{Z}=\beta_{X}+\gamma_{X}$ can be constructed so that $\{\tilde{Z}_{t}=q_{Z}\circ Z_{t}\};q_{Z}=\beta_{X}/(\beta_{X}+\gamma_{X})$ is equivalent to $\{X_{t}\}$. Thus $\{\tilde{X}_{t}\}$ is in turn equivalent to $\{Z^{*}=(q_{Z}q_{X})\circ Z_{t}\}$. One can now choose any $q_{Y}\in[q_{X}q_{Z},1]$ and define $\{Y_{t}=(q_{X}q_{Z}/q_{Y})\circ Z_{t}\},\{\tilde{Y}_{t}=q_{Y}\circ Y_{t}\}$. The process $\{\tilde{Y}_{t}\}$ is then obviously equivalent to $\{Z^{*}_{t}\}$ and thus $\{\tilde{X}_{t}\}$. The proof is complete as statement (A) implies that $\{Y_{t}\}$ has a representation as an INAR($\infty$) process with parameters

[TABLE]

Acknlowledgements

The author would like to thank Leonhard Held, Małgorzata Roos and Christian H. Weiß for helpful comments.