A spatially explicit N-mixture model for the estimation of disease prevalence

Abstract

This article develops an approximate N-mixture model for infectious disease counts that accounts for under-reporting as well as spatial dependence induced by person-to-person spread of disease. We employ the model to estimate actual case counts in Oregon of chlamydia, an easily-treated but usually asymptomatic sexually transmitted disease. We describe a combined parametric bootstrap to account for uncertainty in parameter estimates as well as sampling variability in actual case counts. A simulation study illustrates that our method performs well in many scenarios when the model is correctly specified, and also gives reasonable results when the model is misspecified, and no spatial dependence exists.

Keywords

Case counts Chlamydia normal approximation Oregon parametric bootstrap

1 Introduction

Estimating the true prevalence of under-reported or under-ascertained diseases is difficult to accomplish without knowledge of local demographics or extensive studies (Gibbons et al.(2014) Gibbons, Mangen, Plass, Havelaar, Brooke, Kramarz, Peterson, Stuurman, Cassini, Fèvre, et al.). Additionally, understanding the dynamics of a disease is an important step in understanding changes in prevalence over time (Ostfeld et al.(2005) Ostfeld, Glass, and Keesing). With certain diseases, such as infectious diseases, we would expect there to be spatial dependence among the prevalence of diseases at locations near each other. We aim to provide a spatially explicit statistical method for surveilling imperfectly detected diseases using the reportable disease infrastructure already established through the Centers for Disease Control and the Oregon State Health Authority.

(Brintz et al.(2018) Brintz, Fuentes, and Madsen) established the use of the asymptotic N-mixture model as an appropriate model to make disease prevalence estimates for high-abundance imperfectly-detected or under-reported diseases but disregarded any spatial dependence of the disease prevalence. Other extensions of the N-mixture model have used factors such as weather, habitat, density dependence, to account for spatial patterns (Rossman et al.(2016) Rossman, Yackulic, Saunders, Reid, Davis, and Zipkin). In (Zhao et al.(2017) Zhao, Royle, and Boomer), we see the first spatially explicit N-mixture model which models survival, reproduction, emigration, and abundance while representing movement among adjacent habitat patches. They model the spatial relationship using an immigration random variable, suggesting that the number of immigrants to site $i$ is the proportion of those emigrating from adjacent site $j$ , weighted by the number of sites adjacent to site $j$ . While this structure makes sense for its application of wildlife, it may not appropriately model the spread of disease. In the context of infectious disease, we cannot assume that a discrete proportion of cases leaving one site will move to a neighbour, especially with vastly different populations of neighbouring sites. Therefore, we propose a spatially explicit model that includes covariance information between measured sites by using adjacent sites’ prevalence, rather than their abundances.

Our spatial model employs a neighbourhood structure to account for the spread of disease between adjacent sites, whereas the model of (Brintz et al.(2018) Brintz, Fuentes, and Madsen) had no mechanism to capture the spatial pattern of infection. As a further advancement, we augment the parametric bootstrapped confidence intervals to account for sampling variation in disease counts as well as uncertainty in the model's parameter estimates.

This article is organized as follows. Section 2 presents the Oregon chlamydia data which we use as a motivating example and subsequently analyse in Section 7. Section 3 briefly describes N-mixture models in general and an estimability diagnostic. The proposed spatial asymptotic N-mixture model is described in Section 4. In Section 5, we compare the spatial model to the non-spatial asymptotic model using the diagnostic given in Section 3. In Section 6, we use simulated data with spatial covariance to test the new model's performance and also compare its performance to a model without spatial covariance information. In Section 7, we include an analysis of the chlamydia data using the new spatial model. Section 8 gives some concluding remarks.

2 The chlamydia data

Chlamydia is a sexually transmitted infection caused by the bacterium Chlamydia trachomatis. Chlamydia is imperfectly detected because most people who have it do not have symptoms. Since the infection can be cured if detected, a statistical model should allow for the reduction of cases over time. The infectious nature suggests that neighbouring population centres may have similar prevalence, so a model that accounts for spatial dependence and allows an increase in cases through sexual transmission is appropriate. (Brintz et al.(2018) Brintz, Fuentes, and Madsen) analysed chlamydia data for 2007–2016 obtained from the Oregon Health Authority (http://www.oregon.gov/oha/) using an approximate N-mixture model that accounted for non-detection of the disease but assumed no spatial dependence in case counts. In this article, we extend the model to account for spatial dependence and apply it to Oregon chlamydia data from the years 2010–2018. The spatial model uses annual population sizes by county which we were able obtain beginning in 2010.

Table 1 illustrates the structure of the data, which includes annual population estimates and observed counts of reported chlamydia cases. The observed chlamydia counts of the most populous counties are two to three orders of magnitude larger than those of the least populous counties.

Table 1

The chlamydia data includes yearly population estimates and observed disease cases for all 36 of Oregon's counties. The table shows data for the five most populous and the five least populous counties

County	Population				Observed Cases
	2010	2011	…	2018	2010	2011	…	2018
Multnomah	737 291	749 699	…	811 880	3296	3992	…	5459
Washington	531 645	540 963	…	597 695	1390	1642	…	2404
Clackamas	376 790	379 485	…	416 075	945	974	…	1394
Lane	351 923	354 005	…	379 611	1276	1297	…	1844
Marion	315 951	318 049	…	346 868	1395	1497	…	1887
⋮	⋮	⋮		⋮	⋮	⋮		⋮
Grant	7464	7409	…	7176	14	11	…	24
Wallowa	7012	6978	…	7081	9	3	…	9
Gilliam	1882	1951	…	1894	3	4	…	3
Sherman	1779	1750	…	1708	4	5	…	4
Wheeler	1447	1416	…	1366	3	0	…	3

Figure 1

A heat map of the 2016 reported chlamydia prevalence in all Oregon counties

Figure 1 shows a map of Oregon's 36 counties, shaded by reported chlamydia prevalence in 2016, and suggests spatial dependence, since counties close to each other have similar prevalence. We will utilize county prevalence in the spatial component of our proposed model.

3 N-mixture models

(Royle (2004)) introduced N-mixture models for estimating the size $N$ of a population given imperfect detection of individuals. The model assumes counts $n_{it}$ of distinct individuals are observed at $i = 1, \dots, R$ sites during $t = 1, \dots, T$ sampling occasions, where the populations in the $R$ sites are assumed closed with respect to births, deaths, and migration during the $T$ occasions. Under the closure assumption,

n_{i t} \sim Binomial (N_{i}, p)

with $N_{i}$ the size of the population at site $i$ and $p$ the unknown detection probability. The parameters of interest are the site population sizes $N_{i}$ and the total population size $N = \sum_{i = 1}^{R} N_{i}$ . The objective of (Royle (2004)) was a model with good performance given sparse data (small $R$ and $T$ , small $p$ , and perhaps many $n_{it} = 0$ ), and his innovation was to model the $N_{i}$ as iid random variables with density $f (N_{i}; θ)$ , typically Poisson or negative binomial, then integrate them out of the joint likelihood to obtain

L (p, θ) = \prod_{i = 1}^{R} \{\sum_{N_{i} = max_{i} n_{it}}^{\infty} (\prod_{t = 1}^{T} Bin (n_{it}; N_{i}, p)) f (N_{i}; θ)\},

(3.1)

where $θ$ represents the parameter vector of $f (N_{i}; θ)$ , and Bin (.; N_i,p) represents the binomial likelihood with index $N_{i}$ and success probability $p$ . The integrated likelihood (3.1) is approximated by summing over $N_{i}$ to a large finite bound $K$ , then numerically maximized to obtain maximum likelihood estimates $\hat{θ}$ and $\hat{p}$ . Since the sites’ population sizes $N_{1}, \dots, N_{R}$ are iid, and since their common expectation can be expressed in terms of $θ$ , we can estimate $E (N)$ as $\hat{E} (N_{i}) \cdot R$ , where $\hat{E} (N_{i})$ is $E (N_{i})$ with $θ = \hat{θ}$ .

(Dail and Madsen (2011)) relaxed the population closure assumption. The goal of their model was to estimate $N_{it}$ , the population size at site $i$ and time $t$ , by modelling the joint distribution of $N_{i 1}, \dots, N_{iT}$ to account for changes in population size at site $i$ over time, while retaining the assumption of independence between sites. In particular, the exact likelihood for the generalized N-mixture model is

L (p, θ | {n_{it}}) = \prod_{i = 1}^{R} \{\sum_{N_{i 1} = n_{i 1}}^{\infty} \dots \sum_{N_{iT} = n_{iT}}^{\infty} [\prod_{t = 1}^{T} Bin (n_{it}; N_{it}, p)] f (N_{i 1}, \dots, N_{iT}; θ)\},

(3.2)

with $p$ as above and prior density for the abundance parameters modelled under a Markov assumption as $f ({N_{i 1}, \dots, N_{iT}}; θ) = f (N_{i 1}; θ) \prod_{t = 2}^{T} f (N_{it} | N_{it - 1}; θ)$ . Parameter vector $θ$ includes parameters of the initial abundance distribution $f (N_{i 1}; θ) = Poisson (λ)$ and population dynamics parameters $ω$ , the survival probability, and $γ$ , the arrival rate. Specifically, denoting the number of population units ‘surviving’ from time $t - 1$ to time $t$ as $S_{it}$ , and the gain in population units between time $t - 1$ and time $t$ as $G_{it}$ , let

S_{i t} |N_{i t - 1} \sim Bin (N_{i t - 1}, ω)

(3.3)

\begin{array}{l} G_{i t} |N_{i t - 1} \sim Poisson (γ) \\ N_{i t} |N_{i t - 1} = S_{i t} |N_{i t - 1} + G_{i t} |N_{i t - 1} \end{array}

(3.4)

and as in (Royle (2004)), approximation of the summations to infinity are made using a large finite $K$ . This approximation of the likelihood is computationally intractable when the observed abundances $n_{it}$ are large because the chosen value of $K$ must also be sufficiently large relative to the observations, so the generalized N-mixture model's utility is limited to situations where the $n_{it}$ are small.

(Brintz et al.(2018) Brintz, Fuentes, and Madsen) adapted the generalized N-mixture model to model counts of disease cases, using the Oregon chlamydia data for 2007–2016 mentioned Section 2 as a case study. The chlamydia data includes several observed counts exceeding 1 000, much too large for the generalized N-mixture model. Therefore, they proposed a multivariate normal (MVN) approximation of the generalized N-mixture likelihood (3.2) based on a normal approximation to the binomial and to the Poisson, thus avoiding the infinite sums.

The sites are assumed to be independent, so the MVN approximated likelihood is

L (θ | \{n_{i t}\}) \approx \prod_{i = 1}^{R} M V N (μ_{i}, Σ_{i}),

(3.5)

where MVN (⋅,⋅) denotes the multivariate normal density function, μ_i is the T × 1 mean vector for site i, and Σ_t is the T × T covariance matrix for site i. Parameter vector θ includes parameters for the distribution of initial prevalence N_i1/pop_i, where ${pop}_{i}$ is the population of site $i$ , and population dynamics parameters to model the conditional distribution of N_it given N_it−1. For notational convenience, we now include detection probability p as an element of θ. The elements of μ_t and Σ_t are expressed in terms of the elements of θ.

In addition to avoiding the infinite sums of (3.2), the MVN approximation (3.4) has a closed-form expression for the Fisher information matrix $I$ of the log likelihood. The inverse of $I$ provides the asymptotic variance-covariance matrix of the maximum likelihood parameter estimators. (Brintz et al.(2018) Brintz, Fuentes, and Madsen) proposed using the determinant of $I$ as a diagnostic for parameter estimability, which is difficult for some combinations of parameters within the N-mixture model framework due to the lack of information about both $p$ and $N_{i}$ . For example, if the true detection probability $p$ is close to 0, then the data may not contain sufficient information to decide if $p$ is small or if the $N_{i}$ are small. Similarly, if there are little or no actual population dynamics, then we should not expect to successfully estimate the population dynamics parameters. After describing the spatial model in Section 4, we compare the spatial and non-spatial models using this estimability diagnostic.

4 The spatial N-mixture model

The asymptotic N-mixture model of (3.4) assumes independence between sites, so the observed data vector ${n_{it}}$ , $t = 1, \dots, T$ , $i = 1, \dots, R$ , has a $TR \times TR$ block diagonal covariance matrix where the $T \times T$ blocks $Σ_{i}$ represent dependence in time within a site but the diagonal structure reflects independence between sites. In the spatial model, we wish to model spatial dependence between sites in addition to time dependence within a site, so the spatial model's covariance matrix will not have a block-diagonal structure.

Our focus is on estimating disease prevalence, so we now assume that, as is typical in disease surveillance, observations exist for all sites at all time periods within the scope of the study. In the chlamydia data, the sites are Oregon's counties, and the time periods are the years from 2010 to 2018. In this context, $n_{it}$ and $N_{it}$ represent the observed and actual number of chlamydia cases in county $i$ during year $t$ respectively. The ‘survivor’ random variable $S_{it}$ of (3.3) should be reinterpreted as the number of people in county $i$ who were infected at time $t - 1$ and remain infected at time $t$ , or who were cured and reinfected between time $t$ and time $t - 1$ . The ‘gain’ random variable $G_{it}$ of (3.3) represents the number of new cases in county $i$ between time $t - 1$ and time $t$ , and consists of both new infections and people who are already infected and move into county $i$ between time $t - 1$ and time $t$ .

Given that parameter identifiability is a general concern with N-mixture models, we model spatial dependence without including extra parameters that might be confounded with parameters from the non-spatial components of the model. In order to accomplish this, we incorporate spatial dependence based on adjacency, much like the conditional auto-regressive (CAR) model of (Abellan et al.(2008) Abellan, Richardson, and Best).

In particular, we construct an $R \times R$ adjacency matrix $A$ consisting of 0s and 1s, such that $A_{ij} = A_{ji} = 1$ if sites $i$ and $j$ are neighbours, and $A_{ij} = 0$ otherwise. For the chlamydia data, we define two counties to be neighbours if they share a border. We also specify that site $i$ is its own neighbour, so $A_{ii} = 1$ . Below, we will refer to the size of site $i$ ’s neighbourhood as $A_{i \cdot} \equiv \sum_{s = 1}^{R} A_{is}$ .

For sexually transmitted diseases such as chlamydia, an adjacency matrix defined in this way is sensible because of the connection between the disease and the sex trade. Our model does not specify how adjacency is defined, so other criteria may be used to define an appropriate neighbourhood structure for other diseases.

Preliminary modelling attempts of the chlamydia data suggested that using an equal weighted neighbourhood was not appropriate due to the possibility of adjacent counties having very different population sizes. With equal weighting, the prevalence in very large or small populations disproportionately affects their neighbours. Therefore, we adjust our model to use a weighted prevalence model as described below.

As above $i = 1, \dots, R$ indexes sites and $t = 1, \dots, T$ indexes sampling occasion. It will be convenient to define notation ${N_{α, t_{0}}}$ to represent the set of $N_{it}$ for all sites $i = 1, \dots R$ and all times $t = 1, \dots, t_{0}$ where $α$ refers to ‘all’ sites.

Components of the spatial model are

n_{i t} | \{N_{i t}\} \sim independent Binomial (N_{i t}, p)

(4.1)

N_{i 1} \sim independent Poisson (β p o p_{i 1})

(4.2)

S_{i t} | \{N_{i t}\} \sim independent Binomial (N_{i, t - 1}, ω), t > 1

(4.3)

where ${pop}_{it}$ denotes the population in county $i$ at time $t$ , $p$ represents the constant detection rate, and $S_{it}$ represents those infected in site $i$ at time $t - 1$ and at time $t$ (cf. equation (3.3)). The parameter $β$ represents the initial prevalence of the disease, and the parameter $ω$ represents the probability of a case persisting until the next time period. Therefore, in the disease surveillance context, we refer to $ω$ as the persistence probability rather than the survival probability.

The recruitment process (cf. equation (3.3)) will now incorporate the previous time period's prevalence of all of site $i$ ’s neighbours so that

\begin{matrix} G_{it} | {N_{α, t - 1}} & \sim independent Poisson (γ {\tilde{N}}_{it - 1}), t > 1 \end{matrix}

(4.4)

where

{\tilde{N}}_{it - 1} = {pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} N_{rt - 1} / {pop}_{rt - 1}}{A_{i \cdot}}

(4.5)

for $t > 1$ . Equations (4.2) and (4.3) model the mean ‘recruitment’ of cases to site $i$ at time $t$ as a parameter $γ$ times the average prevalence $N_{rt - 1} / {pop}_{rt - 1}$ from site $i$ ’s neighbours at time $t - 1$ , scaled by the population of site $i$ at time $t - 1$ . In this context, the parameter $γ$ may be thought of as describing the infection rate rather than the rate of recruitment. With this formulation, we assume an increased prevalence at a site's neighbours will have a positive effect on the number of new cases. Additionally, since the model assumes neighbours only affect each other's prevalence, a more populous county will not have an unrealistic effect on a lesser populous neighbour's new cases. We assume $S_{it}$ and $G_{it}$ are conditionally independent given ${N_{α, t - 1}}$ .

For simplicity, our implementation of the model assumes no spatial dependence within the first time period. The model could be augmented to include an initial spatial dependence structure, but here we focus on the spatial dependence induced by disease transmission and movement of infected individuals after a baseline time point.

Modelling $p$ , $β$ , $ω$ , and $γ$ as constant for all counties is another simplifying assumption. Further refinements of the model could describe one or more of these parameters in terms of covariates.

Our model also assumes negligible ‘edge effects’, that is, we suppose that Oregon's counties bordering other states are relatively unaffected by their out-of-state neighbours. In practice, one could include population and prevalence data from these neighbours in the analysis.

To complete the specification of the spatial model, we specify the conditional distribution of $N_{it}$ given the past as

N_{it} | {N_{α, t - 1}} = S_{it} | N_{it - 1} + G_{it} | {N_{α, t - 1}}, t > 1 .

(4.6)

In order to implement this model, we approximate the joint likelihood of the observed data ${n_{it}}$ as MVN with mean and covariance determined by the model defined in equations (4.1) through (4.4). Observe that the model specifies the ${n_{it}}$ are conditionally independent given the ${N_{it}}$ so that, by the law of total covariance,

\begin{matrix} cov (n_{it}, n_{ju}) & = & E [cov (n_{it}, n_{ju} | {N_{it}}] + cov [E (n_{it} | {N_{it}}), E (n_{ju} | {N_{it}})] \\ = & p^{2} cov (N_{it}, N_{ju}) . \end{matrix}

(4.7)

If $Σ$ is the variance-covariance matrix of the ${n_{it}}$ , the normal approximation to the model in equations (4.1) through (4.4) can be expressed as

L (β, p, ω, γ | {n_{it}}) \approx MVN (μ, Σ),

(4.8)

where $μ$ is the vector of means with elements $μ_{it} = E (n_{it})$ .

To write the mean parameter $μ_{it}$ for site $i$ at time $t$ in terms of the model parameters $p$ , $β$ , $ω$ , and $γ$ , we start with $t = 1$ and calculate

μ_{i 1} = E (n_{i 1}) = E [E (n_{i 1} | {N_{it}})] = p β {pop}_{i 1}

(4.9)

by the law of total expectation and equations (4.1) and (4.1). For $t > 1$ , with similar reasoning and given constant persistence (formerly survival) and infection (formerly recruitment) rates $ω$ and $γ$ from (4.1) and (4.2), the expected observed count $E (n_{it})$ is defined recursively as

μ_{it} = p (ω μ_{it - 1} + γ {\tilde{μ}}_{it - 1}),

(4.10)

where ${\tilde{μ}}_{it - 1}$ is defined analogously to ${\tilde{N}}_{it - 1}$ in equation (4.3), replacing each $N_{rt - 1}$ with $μ_{rt - 1}$ .

By ordering the $n_{it}$ first by time $t$ and then by site $i$ within $t$ , $Σ$ can be broken up into $T^{2}$ $R \times R$ sub-matrices where the $R \times R$ matrices along the main diagonal represent the spatial covariance of the observed abundances between sites at time $t$ , for $t = 1, \dots, T$ . The off-diagonal blocks represent the covariance matrices among different sites at different time lags.

Because the ${n_{it}}$ are conditionally independent given the ${N_{it}}$ , all spatial and temporal dependence in the ${n_{it}}$ is due to dependence among the ${N_{it}}$ . Let $σ_{ij}^{tu}$ denote cov (N_it,N_ju) where $t$ and $u$ index the $T$ time periods and $i$ and $j$ index the $R$ sites. From (4.5), a single element of $Σ$ is $σ_{ij}^{tu} / p^{2}$ . For example, $σ_{1 1}^{1 2} / p^{2} = cov (n_{11}, n_{12}) = Σ_{1, R + 1}$ , the element occupying the first row and $(R + 1)$ st column of $Σ$ .

We define the $σ_{ij}^{tu}$ recursively as follows. Appendix A.1.1 gives details of the derivations.

For $t = u = 1$ , the diagonal elements $σ_{ii}^{11}$ of the $R \times R$ spatial covariance matrix for time $t = 1$ are all equal to $β {pop}_{i 1}$ and the off-diagonal elements $σ_{ij}^{11}, i \neq j$ are all 0, since we are assuming all spatial similarity at the first time period is due only to population size. This independence assumption is not required, and a dependence model could be incorporated at time $t = 1$ . With the independence assumption,

σ_{ij}^{11} = \{\begin{matrix} β {pop}_{i 1} & i = j \\ 0 & i \neq j \end{matrix} .

(4.11)

For $t > 1, u = 1$ , $σ_{ij}^{t 1}$ can be written using the law of total covariance as

cov (N_{it}, N_{j 1}) = ω σ_{i, j}^{t - 1, 1} + γ {pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} σ_{r, j}^{t - 1, 1} / {pop}_{rt - 1}}{A_{i \cdot}} .

For $t > 1$ and $u > 1$ , again as a result of the law of total covariance,

\begin{matrix} σ_{ij}^{tu} & = ω^{2} σ_{ij}^{t - 1 u - 1} + ω γ ({pop}_{j u - 1} \sum_{r = 1}^{R} \frac{A_{jr} σ_{ir}^{t - 1 u - 1}}{A_{j \cdot} {pop}_{r u - 1}} + {pop}_{i t - 1} \sum_{r = 1}^{R} \frac{A_{ir} σ_{rj}^{t - 1 u - 1}}{A_{i \cdot} {pop}_{r t - 1}}) \\ + γ^{2} {pop}_{i t - 1} {pop}_{j u - 1} \sum_{r} \sum_{s} \frac{A_{ir} A_{js} σ_{rs}^{t - 1 u - 1}}{A_{i \cdot} A_{j \cdot} {pop}_{rt - 1} {pop}_{su - 1}} . \end{matrix}

(4.12)

Finally, for a diagonal element of $Σ$ , $t = u$ and $i = j$ , and

\begin{matrix} σ_{i i}^{tt} & = ω^{2} σ_{i i}^{t - 1 t - 1} + ω γ (2 {pop}_{i t - 1} \sum_{r = 1}^{R} \frac{A_{ir} σ_{ir}^{t - 1 t - 1}}{A_{i \cdot} {pop}_{r t - 1}}) \\ + γ^{2} {pop}_{i t - 1}^{2} \sum_{r = 1}^{R} \sum_{s = 1}^{R} \frac{A_{ir} A_{is} σ_{rs}^{t - 1 t - 1}}{{pop}_{r t - 1} {pop}_{s t - 1} A_{i}^{2}} \\ + ω (1 - ω) μ_{it - 1} + γ {\tilde{μ}}_{it - 1} . \end{matrix}

(4.13)

The resulting mean vector and covariance matrix are of the form

μ = [\begin{matrix} μ_{1 1} \\ μ_{2 1} \\ ⋮ \\ μ_{R 1} \\ μ_{1 2} \\ ⋮ \\ μ_{R 2} \\ ⋮ \\ μ_{R T} \end{matrix}] Σ = p^{2} [\begin{matrix} σ_{1 1}^{1 1} & \dots & σ_{1 R}^{1 1} & σ_{1 1}^{1 2} & \dots & σ_{1 R}^{1 2} & \dots & σ_{1 1}^{1 T} & \dots & σ_{1 R}^{1 T} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ σ_{R 1}^{1 1} & \dots & σ_{R R}^{1 1} & σ_{R 1}^{1 2} & \dots & σ_{R R}^{1 2} & \dots & σ_{R 1}^{1 T} & \dots & σ_{R R}^{1 T} \\ σ_{1 1}^{2 1} & \dots & σ_{1 R}^{2 1} & σ_{1 1}^{2 2} & \dots & σ_{1 R}^{2 2} & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ σ_{R 1}^{2 1} & \dots & σ_{R R}^{2 1} & σ_{R 1}^{2 2} & \dots & σ_{R R}^{2 2} & ⋮ & ⋮ \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ σ_{1 1}^{T 1} & \dots & σ_{1 R}^{T 1} & \dots & \dots & \dots & \dots & σ_{1 1}^{T T} & \dots & σ_{1 R}^{T T} \\ ⋮ & ⋮ \\ σ_{R 1}^{T 1} & \dots & σ_{R R}^{T 1} & \dots & \dots & \dots & \dots & σ_{R 1}^{T T} & \dots & σ_{R R}^{T T} \end{matrix}]

Equations (4.7) through (4.11), give the parameters of the approximate likelihood (4.6) in terms of parameter $[β, p, ω, γ]$ . We maximize $log [MVN (μ, Σ)]$ with respect to $θ = [β, p, ω, γ]$ to obtain maximum likelihood estimates (MLEs) $\hat{θ} = [\hat{β}, \hat{p}, \hat{ω}, \hat{γ}]$ .

We obtain yearly estimates of total cases using the approach in (Dail and Madsen (2011)) as follows. We estimate $N_{i 1}$ , the case count in the initial time period for each site $i$ , using the estimated prevalence coefficient $\hat{β}$ and the known populations ${pop}_{i 1}$ , and then step through the other time periods using estimates $\hat{ω}$ and $\hat{γ}$ such that

{\hat{N}}_{i 1} = \hat{β} p o p_{i 1}

(4.14)

{\hat{N}}_{i t} = \hat{ω} {\hat{N}}_{i t - 1} + \hat{γ} {\hat{\tilde{N}}}_{i t - 1}, t > 1

(4.15)

\hat{N_{\cdot t}} = \sum_{i = 1}^{R} {\hat{N}}_{i t}

(4.16)

where ${\hat{\tilde{N}}}_{it - 1}$ is the expression in equation (4.3) with the $N_{rt - 1}$ replaced with estimates ${\hat{N}}_{rt - 1}$ . Note that ${\hat{N}}_{it}$ estimates $E (N_{it})$ and ${\hat{N}}_{\cdot t}$ estimates $E (N_{\cdot t})$ .

We determine approximate 95% confidence intervals for the $N_{it}$ and $N_{\cdot t}$ , not for their expectations, via a combined parametric bootstrap as follows. From the theory of maximum likelihood, the asymptotic distribution of $\hat{θ}$ is $N [θ, I^{- 1} (θ)]$ , where $I (θ)$ denotes the Fisher information matrix as a function of the model parameters $θ$ . For the MVN model, $I (θ)$ is available in closed form, so we estimate the asymptotic variance-covariance matrix of $\hat{θ}$ as $[I (\hat{θ})]^{- 1}$ . Thus, we may approximate the sampling distribution of $\hat{θ}$ as $N (\hat{θ}, [I (\hat{θ})]^{- 1})$ , and generate 10,000 samples of $[\hat{p}, \hat{β}, \hat{ω}, \hat{γ}]$ by numerically simulating random draws from this distribution.

Similarly, we approximate the sampling distribution of the $N_{it}$ by a parametric bootstrap. Given the $b$ th bootstrapped realization of $[\hat{p}, \hat{β}, \hat{ω}, \hat{γ}]$ , which we denote $θ^{(b)} = [{\hat{p}}^{(b)}, {\hat{β}}^{(b)}, {\hat{ω}}^{(b)}, {\hat{γ}}^{(b)}]$ , simulate $N_{it}^{(b)}$ for $i = 1, \dots, R$ and $t = 1, \dots, T$ according to the model defined in (4.1) through (4.4) with parameters equal to the elements of $θ^{(b)}$ . Calculating (4.12) from the $N_{it}^{(b)}$ yields $N_{\cdot}^{(b)}$ . Thus we obtain $b = 1, \dots 10, 000$ bootstrapped estimates of $N_{it}$ and $N_{\cdot t}$ for $i = 1, \dots, R$ and $t = 1, \dots, T$ , which we take as parametric bootstrap approximations of the sampling distributions of the ${\hat{N}}_{it}$ and ${\hat{N}}_{\cdot t}$ . Lower and upper $1 - α$ confidence bounds for $N_{it}$ or $N_{\cdot t}$ are taken as the $α / 2$ and $1 - α / 2$ quantiles from the appropriate bootstrapped distribution.

5 Estimability

We now compare the spatial model with the non-spatial using the Fisher information-based diagnostic described in Section 3. The asymptotic variance-covariance matrix of the the MLE of $θ$ is given by the the inverse information matrix $[I (θ)]^{- 1}$ . To determine under what settings the MLEs of the parameters of interest become highly correlated enough to impede estimation, we can calculate the determinant of $I (θ)$ for a variety of $θ$ . We consider that a perfectly negative or positive correlation between two parameters’ MLEs would suggest that they are completely confounded with one another, and this would be reflected as $det [I (θ)] = 0$ .

The $mn$ th entry of $I$ in the MVN (μ, ∑) model is

I_{m n} = \frac{\partial μ'}{\partial θ_{m}} Σ^{- 1} \frac{\partial μ}{\partial θ_{n}} + \frac{1}{2} t r (Σ^{- 1} \frac{\partial Σ}{\partial θ_{m}} Σ^{- 1} \frac{\partial Σ}{\partial θ_{n}}),

(5.1)

where $μ$ and $Σ$ are defined above and $[θ_{1} θ_{2} θ_{3} θ_{4}] = [p β ω γ]$ (Mardia and Marshall (1984)). In comparing the spatial and non-spatial models, we note that detection probability $p$ and persistence rate $ω$ have the same interpretation in both models. The infection parameter $γ$ must be interpreted differently, as it serves as the parameter which determines the strength of the spatial relationship among neighbours in the spatial model whereas in the non-spatial model, gains into site $i$ may come from anywhere, not primarily neighbours.

One improvement of the spatial model over the non-spatial model is that the $ω$ and $γ$ parameters are no longer completely confounded for estimates of $ω$ close to zero. As an example, we plug in point estimates for the chlamydia data from (Brintz et al.(2018) Brintz, Fuentes, and Madsen), that is, $β = 0.00392, p = 0.668, ω = 0, γ = 1.072$ , and we find that the Fisher information-based asymptotic correlation between MLEs of $ω$ and $γ$ is $-$ 0.9992. While these parameters are highly negatively correlated, they are no longer completely confounded with each other and their asymptotic variance can be estimated using the spatial model.

It should be noted that (5.1) is based on the true parameters. If parameter confounding occurs for a particular dataset, it will be difficult to obtain parameter estimates to substitute into (5.1). It may be possible to instead use historical information to set bounds for the parameters, then investigate identifiability from (5.1) by examining a limited number of plausible combinations of parameter values.

6 Simulations

We conducted two simulation studies. The first explores the performance of the spatially explicit model, and the second investigates its robustness compared to the non-spatial model. In the first study, we used the Oregon adjacency matrix and county populations for the nine years 2010 to 2018 (Table 1) to simulate abundance data from 24 scenarios constructed using $β \in (0.005, 0.05)$ , $p \in (0.4, 0.7, 0.9)$ , $ω \in (0.5, 0.8)$ , and $γ \in (0.3, 0.6)$ .

We generated 1 000 simulated datasets for each scenario. For each simulated dataset, we optimized the log of approximate likelihood (4.6) with respect to $θ = [β, p, ω, γ]$ using the optim function of (R Core Team (2016)) with method L-BFGS-B and box constraints $0.00001 \leq β \leq 0.99$ , $0.01 \leq p \leq 0.99$ , $0.01 \leq ω \leq 0.99$ , $0.01 \leq γ \leq 1$ . Preliminary explorations suggested the optimization was sensitive to the initial value of $p$ but not of $β$ , $ω$ , or $γ$ . Therefore, we started the optimization search at three different initial points with $p_{init} \in {0.1, 0.5, 0.9}$ , $ω_{init} = γ_{init} = 0.5$ , and $β_{init} = \frac{1}{R} \sum_{i = 1}^{R} n_{i 1} / {pop}_{i 1}$ , the observed prevalence in year $i = 1$ averaged over the $R = 36$ counties. When the three optimizations resulted in different final values, we set $\hat{θ}$ to the result that maximized the log-likelihood. We deemed a simulation successful when optim returned a convergence code of 0 and all final values were properly within their box constraints. When a simulation was unsuccessful, we discarded it. In practice, one would look further for starting values that lead to a successful optimization, though this was not feasible within our simulation study.

We calculated the absolute relative error of the prevalence estimates for county $i$ in year $t$ as

{ARE}_{it} = \frac{1}{n_{s}} \sum_{k = 1}^{n_{s}} \frac{| N_{it}^{(k)} - {\hat{N}}_{it}^{(k)} |}{E (N_{it}^{(k)})},

(6.1)

where $N_{it}^{(k)}$ is the realization of $N_{it}$ in the $k$ th simulation, and ${\hat{N}}_{it}^{(k)}$ is the estimate of $N_{it}^{(k)}$ . The denominator of each term in (6.1) is $E (N_{it}^{(k)})$ rather than $N_{it}^{(k)}$ because the latter is occasionally simulated as 0. We then summarized the performance of all prevalence estimates for a given simulation scenario as

ARE = \frac{1}{RT} \sum_{i = 1}^{R} \sum_{t = 1}^{T} {ARE}_{it} .

(6.2)

We calculated a 90% confidence interval for $N_{\cdot 9}$ , the total abundance for the last time period, according to the combined parametric bootstrap procedure described in Section 4, and recorded its success and width. The combinations of the settings chosen reflect the observed slow increase in chlamydia prevalence over time on the Oregon Health Authority website.

Table 2

Each row is based on the successful simulations out of an attempted 1 000 simulated datasets using the given $β$ , $p$ , $ω$ and $γ$ values. ARE is the summarized absolute relative error from equation (6.2) averaged over the simulations, Cov% is the proportion of 90% confidence intervals for $N_{\cdot 9}$ that included the true value, and Width is the average width of these confidence intervals. $n_{s}$ is the number of simulations out of 1 000 that successfully converged to a value properly inside the box constraints

$β$	$p$	$ω$	$γ$	ARE	Cov%	Width	$n_{s}$
0.005	0.4	0.5	0.3	0.0980	87.7	3 243	991
0.005	0.4	0.5	0.6	0.1138	89.8	40 738	997
0.005	0.4	0.8	0.3	0.0444	87.8	28 517	978
0.005	0.4	0.8	0.6	0.0333	87.8	20 4513	956
0.005	0.7	0.5	0.3	0.0204	90.5	1 573	994
0.005	0.7	0.5	0.6	0.0265	88.7	19 873	997
0.005	0.7	0.8	0.3	0.0106	91.1	12 417	995
0.005	0.7	0.8	0.6	0.0107	89.9	96 128	989
0.005	0.9	0.5	0.3	0.0361	96.9	1 122	834
0.005	0.9	0.5	0.6	0.0322	96.0	12 975	860
0.005	0.9	0.8	0.3	0.0083	93.6	7 292	977
0.005	0.9	0.8	0.6	0.0081	92.5	55 899	961
0.050	0.4	0.5	0.3	0.0645	90.1	29 600	995
0.050	0.4	0.5	0.6	0.0724	90.4	377 475	991
0.050	0.4	0.8	0.3	0.0285	88.7	266 551	986
0.050	0.4	0.8	0.6	0.0260	91.1	193 4590	984
0.050	0.7	0.5	0.3	0.0338	90.3	15 079	988
0.050	0.7	0.5	0.6	0.0256	89.8	189 804	988
0.050	0.7	0.8	0.3	0.0013	91.3	112 203	982
0.050	0.7	0.8	0.6	0.0127	89.5	889 888	984
0.050	0.9	0.5	0.3	0.0301	96.1	9 633	878
0.050	0.9	0.5	0.6	0.0202	97.3	116 096	875
0.050	0.9	0.8	0.3	0.0120	87.1	61 990	919
0.050	0.9	0.8	0.6	0.0041	88.5	471 553	956

Table A.1 in the Appendix gives estimates of the four parameters $β$ , $p$ , $ω$ , and $γ$ , averaged over the $n_{s}$ simulations for each of the 24 scenarios. The results there suggest that all four estimators are approximately unbiased.

Table 2 shows simulation results for the abundance estimates. ARE is the summarized absolute relative error (6.2), averaged over the $n_{s}$ successful simulations for each scenario. Higher true detection level $p$ was associated with smaller ARE as was larger $ω$ . These settings lead to larger observed counts $n_{it}$ . Larger $γ$ was associated with a slightly increased ARE in several scenarios. This setting is where the spatial dependence is stronger. The level of $β$ had little effect on ARE except when $p = 0.4$ and $ω = 0.05$ , where ARE was slightly smaller for $β = 0.05$ . In this setting, the $n_{it}$ will be much larger than when $β = 0.005$ .

Turning to confidence interval success percentage, we note that the nominal 90% coverage was no less than 87.1%. When $p = 0.9$ and $ω = 0.5$ , 90% nominal coverage exceeds 95% for all combinations of $γ$ and $β$ . For these settings, detection is nearly perfect and persistence is so low that the disease vanishes from some of the smaller counties by time $t = 9$ . In all other scenarios, observed coverage rates are no more than 3.6% from the nominal 90%.

Average confidence interval widths were larger for higher settings of initial prevalence $β$ , infection rate $γ$ , and persistence $ω$ . These settings all lead to larger observed counts $n_{it}$ , and therefore to larger estimated counts $N_{it}$ .

The number of successful simulations out of the 1 000 attempted is at least 834, and usually over 950. The lowest values of $n_{s}$ are when $p = 0.9$ and $ω = 0.5$ . These are the four simulation scenarios that yield the highest confidence interval success rate.

The second simulation study explored the same parameter settings but simulated data from the spatial model described in Section 4 as well as from the non-spatial model N-mixture model given in Section 3. Both models were fit to both types of simulations to explore robustness to model misspecification. We simulated 1 000 datasets for each simulation type and combination of parameters.

Table 3

Each row is based on simulated data with $R = 36$ sites and $T = 9$ time periods using the given $β$ , $p$ , $ω$ , and $γ$ values. The two columns labelled Spatial Data represent 1 000 datasets simulated under the spatial model described in 4. The two columns labelled Non-spatial Data represent 1 000 datasets simulated under the generalized non-spatial M-mixture model of 3. Both the spatial model and the non-spatial model are fit to each dataset, and the ARE (6.2) is recorded. Omitted from the table are 10 scenarios where there were fewer than 750 successful optimizations out of 1 000 for any one of the four model fits

				Non-spatial data		Spatial data
				Spatial	Non-spatial	Spatial	Non-spatial
$β$	$p$	$ω$	$γ$	ARE	ARE	ARE	ARE
0.005	0.4	0.5	0.3	0.124	0.0794	0.098	1.47
0.005	0.4	0.5	0.6	0.234	0.0459	0.114	2.36
0.005	0.4	0.8	0.3	0.181	0.115	0.0444	1.7
0.005	0.7	0.5	0.3	0.0277	0.0113	0.0204	0.597
0.005	0.7	0.5	0.6	0.0088	0.0183	0.0265	1.44
0.005	0.7	0.8	0.3	0.0182	0.0066	0.0106	0.212
0.005	0.9	0.5	0.3	0.0836	0.0115	0.0361	0.337
0.005	0.9	0.5	0.6	0.0974	0.011	0.0322	0.722
0.005	0.9	0.8	0.3	0.0366	0.0053	0.0083	0.0905
0.05	0.4	0.5	0.3	0.111	0.0503	0.0645	0.826
0.05	0.7	0.5	0.3	0.0512	0.0153	0.0338	0.414
0.05	0.7	0.8	0.3	0.04	0.0027	0.00132	0.0278
0.05	0.9	0.5	0.3	0.103	0.0196	0.0301	0.27
0.05	0.9	0.8	0.3	0.0874	0.01	0.012	0.009

Table 3 gives results of the robustness simulations for the 14 scenarios where all four of the model fits involved more than 750 successful optimizations out of the 1 000 attempted. The ARE (6.2) using the correct model is smaller than the ARE from the misspecified model in all scenarios except one. The exception occurs when $β = 0.005$ , $p = 0.7$ , $ω = 0.5$ , and $γ = 0.6$ , and the spatial model is fit to non-spatial data. The difference in ARE between the misspecified and correct models is greater for the spatial data than the non-spatial data in 12 of the 14 scenarios shown, indicating that the spatial model is more robust than the non-spatial model.

7 Chlamydia results

We fit the chlamydia data to the spatially explicit model. Point estimates, asymptotic variances, and approximate 95% confidence intervals are shown in Table 4.

Table 4

Output of R's optim function. MLE are the estimates, $\hat{Var} (MLE)$ are the diagonal elements of the inverse of the asymptotic information matrix evaluated at the MLE. Approximate 95% confidence limits (CL) are calculated as $MLE \pm 1.96 \sqrt{\hat{Var} (MLE)}$

parameter	MLE	$\hat{Var} (MLE)$	CL
$β$	0.0053	0.14e-7	(0.0046, 0.0060)
$p$	0.6694	0.00221	(0.5774, 0.7614)
$ω$	0.9164	0.00103	(0.8535, 0.9793)
$γ$	0.1693	0.00101	(0.1070, 0.2316)

The estimated initial prevalence was $\hat{β} = 0.0053$ with an approximate 95% confidence interval of $(0.0051, 0.0055)$ . The estimated detection probability was $\hat{p} = 0.6692$ with approximate 95% confidence interval $(0.5773, 0.7611)$ suggesting a substantially less-than-perfect rate for detecting chlamydia cases in Oregon. The results also include non-zero estimates for both persistence parameter $ω$ and infection parameter $γ$ , allowing a potentially more interpretable model than the non-spatial model which estimated $ω = 0$ . Given that the estimated detection rate is well below 100%, including those with latent symptoms or those in hidden populations, we would not expect all cases of chlamydia to be cured or to leave the county. With the spatial model, we estimated $\hat{ω} = 0.9165$ with approximate 95% confidence interval $(0.8536, 0.9794)$ suggesting a high level of disease persistence and a lack of relocation from year to year for those with chlamydia. Note that the model doesn't distinguish between individuals that are not cured and those that are cured and reinfected within a single year. We estimated $\hat{γ} = 0.1691$ with approximate 95% confidence interval $(0.1068, 0.2314)$ . The non-zero estimate of $γ$ suggests that there is a spatial effect of infection from adjacent counties.

Figure 2

Both observed total counts $n_{\cdot t}$ and estimated total counts $\hat{N_{\cdot t}}$ of Chlamydia in Oregon for 2010–2018. Vertical bars represent 95% confidence intervals. The graph includes predicted total counts and 95% prediction interval for 2019

Figure 2 illustrates observed total counts $n_{\cdot t}$ and estimated total counts ${\hat{N}}_{\cdot t}$ for $t = 1 \dots 9$ , corresponding to calendar years 2010–2018. 95% confidence intervals for the $N_{\cdot t}$ are calculated using the combined parametric bootstrap described in Section 4. Since the model described in (4.1) through (4.4) specifies the $N_{it}$ depend on ${pop}_{it - 1}$ , the population sizes at time $t - 1$ , we can carry the combined parametric bootstrap forward one time step and obtain a point and interval prediction for $N_{\cdot 10}$ , the total case count for 2019.

Figure 2 shows a steeper increasing trend of chlamydia cases than was found with the non-spatial model and similar precision in the confidence intervals. However, the confidence intervals shown here estimate $N_{\cdot t}$ whereas those in (Brintz et al.(2018) Brintz, Fuentes, and Madsen) estimated $E (N_{\cdot t})$ . For 2019, we predict a total chlamydia case count of $43, 017$ with 95% prediction interval of $(36, 931, 49, 290)$ .

Instead of estimating state-wide totals, we can report county estimates ${\hat{N}}_{it}$ for all counties $i$ and any given year $t$ , as these are available within the parametric bootstrap described in Section 4. Table 5 gives point estimates and 95% confidence bounds of case counts for 2018. Note that if estimates of county-wise prevalence are preferred, these are immediate from Table 5 by dividing by the given county populations.

Table 5

County-wise observed and estimated case counts for 2018. Estimates and confidence limits are derived from the parametric bootstrap described in Section 4. Convert to point estimates and confidence limits for prevalence by dividing by the given county population

County	Population	Observed	Estimated	95% CL
Multnomah	811 880	5 459	7 651	(6 696, 8 574)
Washington	597 695	2 404	5 558	(4 863, 6 227)
Clackamas	416 075	1 394	3 896	(3 405, 4 373)
Lane	379 611	1 844	3 607	(3 155, 4 050)
Marion	346 868	1 887	3 244	(2 829, 3 646)
Jackson	219 564	994	2 116	(1 842, 2 383)
Deschutes	191 996	591	1 699	(1 472, 1 919)
Linn	127 335	472	1 188	(1 029, 1 345)
Douglas	110 283	328	1 092	(942, 1 238)
Yamhill	107 002	362	1 007	(869, 1 143)
Benton	92 101	453	875	(753, 994)
Josephine	87 393	294	852	(729, 971)
Polk	85 234	359	782	(673, 888)
Umatilla	77 516	304	788	(663, 910)
Klamath	67 653	347	665	(567, 763)
Coos	64 389	165	646	(545, 740)
Columbia	52 377	160	498	(421, 573)
Lincoln	49 388	128	468	(397, 539)
Clatsop	39 764	141	377	(316, 438)
Malheur	30 725	188	319	(261, 377)
Tillamook	26 787	53	256	(212, 300)
Wasco	26 505	119	259	(209, 310)
Union	26 461	85	266	(217, 315)
Hood River	23 428	99	227	(186, 268)
Jefferson	24 192	165	224	(183, 266)
Curry	22 813	50	228	(187, 270)
Crook	23 867	103	218	(175, 262)
Baker	16 006	33	166	(132, 200)
Morrow	11 372	52	116	(87, 146)
Lake	7 879	18	79	(59, 100)
Harney	7 329	18	74	(55, 94)
Grant	7 176	20	75	(56, 95)
Wallowa	7 081	9	71	(52, 91)
Gilliam	1 894	3	20	(11, 29)
Sherman	1 708	4	18	(9, 28)
Wheeler	1 366	3	14	(7, 22)

8 Discussion

We have developed a spatial extension to the asymptotic N-mixture model as a natural next step for modelling infectious disease data. The human to human connection suggests that a high prevalence in one location could encourage a high prevalence in neighbouring locations. As such, we induced the spatial dependence using a neighbourhood model. By basing the spatial dependence on prevalence by population rather than total case count over a year, we avoid letting highly populous counties disproportionately affect their more sparsely populated neighbours. The model assumes that the recruitment of cases into one county is based on the average prevalence of its neighbours (including itself) at previous time period. In contrast, the non-spatial model only conditions case recruitment on each county's own previous time period's prevalence.

We have shown through simulation that the spatial asymptotic approximation of the N-mixture model performs well in common scenarios and outperforms the non-spatial asymptotic approximation model with spatially dependent data. The $ω$ and $γ$ parameters are no longer confounded for values of $ω$ close to 0, and the interpretation of these parameters differ as compared to the non-spatial model. The spatial model also shows robustness in cases where the data is not spatially dependent.

In the Oregon Health Authority chlamydia data, we see somewhat similar parameter point estimates between the spatial and non-spatial models for both the initial prevalence parameter $β$ and the detection parameter $p$ . However, with an estimated detection rate $\hat{p} \approx 0.67$ , we would expect the persistence of the disease not to be $ω = 0$ as suggested by the non-spatial model, as some of the cases would be hidden or undiagnosed and not be treated. We also might not expect a persistence rate over 85% as suggested by the spatial model, although with the model's lack of distinction between uncured cases and reinfections, an unexpectedly large $\hat{ω}$ is perhaps not so surprising.

Future refinements of the spatial model include modelling the four parameters as functions of covariates, though we note that these must be at the same spatial scale as the observed counts. We may also consider redefining spatial adjacency so to better model the spatial pattern of infection. For example, in Oregon, there is more movement from north to south along the Interstate 5 corridor than there is in the east-west direction. This anisotropy could be modelled by extending neighbourhoods in the north-south direction.

9 Computing

The simulations and data analysis described in Sections 6 and 7 are coded in R (R Core Team (2016)). Source code and the Oregon chlamydia data are available at http://www.statmod.org/smij/archive.html.

A.1 Appendix

A.1.1 Details of covariance calculations of the MVN model

We define the $σ_{ij}^{tu} = cov (N_{it}, N_{ju})$ recursively by considering four cases:

$t = u = 1$

$t > 1, u = 1$ (and by symmetry this covers the case $t = 1, u > 1$ )

$t > 1, u > 1$ and $i \neq j$ or $t \neq u$ (non-diagonal element of $Σ$ )

$t = u$ and $i = j$ (diagonal element of $Σ$ )

In all cases, we calculate cov (N_it,N_ju) using the law of total covariance by conditioning on ${N_{α, v - 1}}$ with $v = max (t, u)$ :

cov (N_{it}, N_{ju}) = E [cov (N_{it}, N_{ju} | {N_{α, v - 1}})] + cov [E (N_{it} | {N_{α, v - 1}}), E (N_{ju} | {N_{α, v - 1}})] .

(A.1)

If $i \neq j$ or $t \neq u$ the first term of (A.1) is 0 because one of the two $N$ ’s is among the information we're conditioning on. If $i = j$ and $t = u$ , then the first term of (A.1) becomes

E [var (N_{it} | {N_{α, t - 1}})] \neq 0

which is why we have to consider the fourth case.

Case 1: Suppose $t = u = 1$ .

The diagonal elements $σ_{ii}^{11}$ of the $R \times R$ spatial covariance matrix are all equal to $β {pop}_{i 1}$ and the off-diagonal elements $σ_{ij}^{11}, i \neq j$ are all 0, since we are assuming all spatial similarity at the first time period is due only to population size. Thus,

σ_{ij}^{11} = \{\begin{matrix} β {pop}_{i 1} & i = j \\ 0 & i \neq j \end{matrix} .

(A.2)

Case 2: Suppose $t > 1, u = 1$ .

$σ_{ij}^{t 1}$ can be written using the law of total covariance as

\begin{matrix} cov (N_{it}, N_{j 1}) & = & E [cov (N_{it}, N_{j 1} | {N_{α, t - 1}})] + cov [E (N_{it} | {N_{α, t - 1}}), E (N_{j 1} | {N_{α, t - 1}})] \\ = & cov (ω N_{it - 1} + γ {\tilde{N}}_{it - 1}, N_{j 1}) \\ = & ω σ_{i, j}^{t - 1, 1} + γ {pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} σ_{r, j}^{t - 1, 1} / {pop}_{rt - 1}}{A_{i \cdot}} . \end{matrix}

(A.3)

The first term on the right-hand side of the first line of (A.3) is 0 because $N_{j 1}$ is among the information we're conditioning on, so cov (N_it, N_j1|{N_{α, 1}})=0.

Case 3: Suppose $t > 1$ and $u > 1$ where $i \neq j$ or $t \neq u$ .

Let $v = max {t, u}$ .

\begin{matrix} cov (N_{it}, N_{ju}) & = & E [cov (N_{it}, N_{ju} | {N_{α, v - 1}})] + cov [E (N_{it} | {N_{α, v - 1}}), E (N_{ju} | {N_{α, v - 1}})] \\ = & 0 + cov [E (S_{it} + G_{it} | {N_{α, v - 1}}), E (S_{ju} + G_{ju} | {N_{α, v - 1}}) \\ = & cov (ω N_{it - 1} + γ {\tilde{N}}_{it - 1}, ω N_{ju - 1} + γ {\tilde{N}}_{ju - 1}) \\ = & ω^{2} cov (N_{it - 1}, N_{ju - 1}) \\ + γ ω [cov ({\tilde{N}}_{it - 1}, N_{ju - 1}) + cov (N_{it - 1}, {\tilde{N}}_{ju - 1})] \\ + γ^{2} cov ({\tilde{N}}_{it - 1}, {\tilde{N}}_{ju - 1}) \\ = & ω^{2} σ_{ij}^{t - 1 u - 1} \\ + γ ω cov ({pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} N_{rt - 1} / {pop}_{rt - 1}}{A_{i \cdot}}, N_{ju - 1}) \\ + γ ω cov (N_{it - 1}, {pop}_{ju - 1} \sum_{r = 1}^{R} \frac{A_{jr} N_{ru - 1} / {pop}_{ru - 1}}{A_{j \cdot}}) \\ + γ^{2} cov ({pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} N_{rt - 1} / {pop}_{rt - 1}}{A_{i \cdot}}, {pop}_{ju - 1} \sum_{r = 1}^{R} \frac{A_{jr} N_{ru - 1} / {pop}_{ru - 1}}{A_{j \cdot}}) \\ = & ω^{2} σ_{ij}^{t - 1 u - 1} \\ + γ ω {pop}_{it - 1} \sum_{r = 1}^{R} \frac{A_{ir} σ_{rj}^{t - 1 u - 1} / {pop}_{rt - 1}}{A_{i \cdot}} \\ + γ ω {pop}_{ju - 1} \sum_{r = 1}^{R} \frac{A_{jr} σ_{ir}^{t - 1 u - 1} / {pop}_{ru - 1}}{A_{j \cdot}} \\ + γ^{2} {pop}_{it - 1} {pop}_{ju - 1} \sum_{r = 1}^{R} \sum_{s = 1}^{R} \frac{A_{ir} A_{js} σ_{rs}^{t - 1 u - 1}}{{pop}_{rt - 1} {pop}_{su - 1} A_{i \cdot} A_{j \cdot}} . \end{matrix}

(A.4)

Case 4: Suppose $t = u$ and $i = j$ . Here, $v = t$ .

For a diagonal element of

Σ

, we can use the calculation in (A.4) since cov (N_it, N_it) = var (N_it). However, since

E [cov (N_{it}, N_{it} | {N_{α, v - 1}})] \neq 0

, unlike the second line of equation (A.4), the following must be added to the expression on the last line of (A.4).

Table A.1

Average parameter estimates over the $n_{s}$ successful simulations for each of the 24 scenarios

$β$	$p$	$ω$	$γ$	$\hat{β}$	$\hat{p}$	$\hat{ω}$	$\hat{γ}$
0.005	0.4	0.5	0.3	0.005	0.398	0.497	0.303
0.005	0.4	0.5	0.6	0.006	0.391	0.512	0.588
0.005	0.4	0.8	0.3	0.005	0.402	0.787	0.313
0.005	0.4	0.8	0.6	0.005	0.407	0.779	0.621
0.005	0.7	0.5	0.3	0.005	0.700	0.503	0.297
0.005	0.7	0.5	0.6	0.005	0.697	0.510	0.590
0.005	0.7	0.8	0.3	0.005	0.699	0.804	0.297
0.005	0.7	0.8	0.6	0.005	0.700	0.798	0.602
0.005	0.9	0.5	0.3	0.005	0.873	0.514	0.286
0.005	0.9	0.5	0.6	0.005	0.879	0.516	0.584
0.005	0.9	0.8	0.3	0.005	0.895	0.806	0.295
0.005	0.9	0.8	0.6	0.005	0.894	0.806	0.594
0.050	0.4	0.5	0.3	0.053	0.400	0.496	0.304
0.050	0.4	0.5	0.6	0.054	0.397	0.499	0.601
0.050	0.4	0.8	0.3	0.051	0.403	0.795	0.305
0.050	0.4	0.8	0.6	0.051	0.403	0.793	0.607
0.050	0.7	0.5	0.3	0.052	0.690	0.502	0.298
0.050	0.7	0.5	0.6	0.051	0.695	0.503	0.597
0.050	0.7	0.8	0.3	0.050	0.704	0.794	0.306
0.050	0.7	0.8	0.6	0.051	0.697	0.795	0.605
0.050	0.9	0.5	0.3	0.052	0.878	0.510	0.290
0.050	0.9	0.5	0.6	0.051	0.887	0.504	0.596
0.050	0.9	0.8	0.3	0.051	0.892	0.803	0.297
0.050	0.9	0.8	0.6	0.050	0.899	0.798	0.602

\begin{matrix} E [cov (N_{it}, N_{it} | {N_{α, v - 1}})] & = & E [var (N_{it} | {N_{α, v - 1}})] \\ = & E [var (S_{it} + G_{it} | {N_{α, v - 1}})] \\ = & E [var (S_{it} | {N_{α, v - 1}}) + var (G_{it} | {N_{α, v - 1}})] \\ = & E [ω (1 - ω) N_{it - 1} + γ {\tilde{N}}_{it - 1}] \\ = & ω (1 - ω) μ_{it - 1} + γ {\tilde{μ}}_{it - 1} . \end{matrix}

(A.5)

There is no covariance term in the third line because cov (S_it,G_it|{N_α,v-1})=0.

A.1.2 Simulation results

Table A.1 gives averages of point estimates for parameters $β$ , $p$ , $ω$ , and $γ$ over the $n_{s}$ successful simulations out of 1 000 for each of the 24 scenarios described in Section 6. In all 24 cases, the estimates appear unbiased.

Footnotes

Acknowledgments

The authors disclosed receipt of the following financial support for the research, authorship and/or publication of this article:

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship and/or publication of this article:

This investigation was supported by the University of Utah Population Health Research (PHR) Foundation, with funding in part from the National Center for Research Resources and the National Center for Advancing Translational Sciences, National Institutes of Health, through Grant UL1TR002538.

References

Abellan

, Richardson

and Best

(2008) Use of space-time models to investigate the stability of patterns of disease. Environmental Health Perspectives , 116, 1111–19.

Brintz

, Fuentes

and Madsen

(2018) An asymptotic approximation to the N-mixture model for the estimation of disease prevalence. Biometrics , 74, 1512–18.

Dail

and Madsen

(2011) Models for estimating abundance from repeated counts of an open metapopulation. Biometrics , 67, 577–87.

Gibbons

, Mangen

M-JJ

, PlaßD Havelaar

, Brooke

, Kramarz

, Peterson

, Stuurman

, Cassini

, Fevre

and Kretzschmar

MEE

(2014). Measuring underreporting and under-ascertainment in infectious disease datasets: A comparison of methods. BMC Public Health , 14, 147.

Mardia

and Marshall

(1984) Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika , 71, 135–46.

Oregon Health Authority (2017). Oregon STD statistics . URL http://www.oregon.gov/oha/PH/DiseasesConditions/CommunicableDisease/DiseaseSurveil (last accessed 1 August 2017).

Ostfeld

, Glass

, and Keesing

(2005) Spatial epidemiology: An emerging (or re-emerging) discipline. Trends in Ecology & Evolution , 20, 328–36.

R Core Team (2016) R: A Language and Environment for Statistical Computing . Vienna: R Foundation for Statistical Computing. URL https://www.R-project.org/9last accessed 18 May 2021).

Rossman

, Yackulic

, Saunders

, Reid

, Davis

and Zipkin

(2016) Dynamic N-occupancy models: Estimating demographic rates and local abundance from detection-nondetection data. Ecology , 97, 3300–3307.

10.

Royle

(2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics , 60, 108–15.

11.

Zhao

, Royle

and Boomer

(2017). Spatially explicit dynamic N-mixture models. Population Ecology , 59, 293–300.