A Test for Differential Ascertainment in Case-Control Studies with Application to Child Maltreatment

TL;DR

This paper introduces a statistical test to detect differential ascertainment in case-control studies, addressing bias in odds ratio estimation caused by multiple data sources, with an application to child maltreatment death data.

Contribution

It presents a novel method to test for differential ascertainment and demonstrates how to correct bias in odds ratio calculations in case-control studies with multiple data sources.

Findings

Detected differential ascertainment by race in child maltreatment death data.

Showed that ignoring differential ascertainment leads to biased odds ratios.

Provided a correction method to improve estimate accuracy.

Abstract

We propose a method to test for the presence of differential ascertainment in case-control studies, when data are collected by multiple sources. We show that, when differential ascertainment is present, the use of only the observed cases leads to severe bias in the computation of the odds ratio. We can alleviate the effect of such bias using the estimates that our method of testing for differential ascertainment naturally provides. We apply it to a dataset obtained from the National Violent Death Reporting System, with the goal of checking for the presence of differential ascertainment by race in the count of deaths caused by child maltreatment.

Tables

Table 1. Table 1: Complete contingency table M e superscript 𝑀 𝑒 M^{e} for the count of individuals with exposure e 𝑒 e in the three-list setting.

		List 1
		Y		N
		List 2
		Y	N	Y	N
List 3	Y	$M_{111}^e$	$M_{101}^e$	$M_{011}^e$	$M_{001}^e$
	N	$M_{110}^e$	$M_{100}^e$	$M_{010}^e$	$M_{000}^e$

Table 2. Table 2: The point estimate and 95 % percent 95 95\% confidence interval for each parameter of the model, when the differential ascertainment parameter is set to be θ ∈ { − 1 , − 0.5 , 0 , 0.5 , 1 } 𝜃 1 0.5 0 0.5 1 \theta\in\{-1,-0.5,0,0.5,1\} . We used α i = − 0.5 subscript 𝛼 𝑖 0.5 \alpha_{i}=-0.5 and α i ​ j = 0.2 subscript 𝛼 𝑖 𝑗 0.2 \alpha_{ij}=0.2 , while the mean counts are γ E = 500 , γ U = 1000 formulae-sequence subscript 𝛾 𝐸 500 subscript 𝛾 𝑈 1000 \gamma_{E}=500,\gamma_{U}=1000 . Each population is generated 200 times and ascertained with 3 lists.

	$\theta$
	-1	-0.5	0	0.5	1
${\hat{\alpha }}_1$	-0.518 (-0.76, -0.28)	-0.507 (-0.74, -0.29)	-0.504 (-0.72, -0.32)	-0.498 (-0.67, -0.32)	-0.504 (-0.70, -0.32)
${\hat{\alpha }}_2$	-0.522 (-0.84, -0.19)	-0.505 (-0.82, -0.16)	-0.509 (-0.84, -0.22)	-0.508 (-0.78, -0.24)	-0.507 (-0.78, -0.22)
${\hat{\alpha }}_3$	-0.510 (-0.94, -0.08)	-0.504 (-0.91, -0.08)	-0.509 (-0.88, -0.11)	-0.490 (-0.81, -0.13)	-0.490 (-0.88, -0.10)
${\hat{\alpha }}_{12}$	0.221 (-0.16, 0.59)	0.203 (-0.14, 0.59)	0.209 (-0.09, 0.55)	0.199 (-0.11, 0.51)	0.213 (-0.08, 0.49)
${\hat{\alpha }}_{13}$	0.203 (-0.13, 0.56)	0.193 (-0.13, 0.55)	0.205 (-0.09, 0.52)	0.187 (-0.08, 0.50)	0.183 (-0.09, 0.50)
${\hat{\alpha }}_{23}$	0.210 (-0.10, 0.58)	0.214 (-0.19, 0.63)	0.202 (-0.12, 0.50)	0.196 (-0.06, 0.52)	0.201 (-0.04, 0.49)
${\hat{\gamma }}_E$	512 (382, 698)	511 (413, 620)	503 (430, 576)	500 (453, 565)	499 (457, 552)
${\hat{\gamma }}_U$	1013 (886, 1169)	1007 (889, 1139)	1005 (897, 1159)	998 (898, 1112)	1002 (902, 1131)
$\hat{\theta }$	-1.010 (-1.33, -0.71)	-0.521 (-0.76, -0.29)	0.006 (-0.18, 0.18)	0.504 (0.34, 0.66)	0.999 (0.81, 1.19)

Table 3. Table 3: Contingency table of exposed (white) and unexposed (black) cases. The cell at the bottom right is missing because it corresponds to the count of individuals that are not ascertained in any of the three lists.

		DC
		Y		N
		LE
		Y	N	Y	N
CME	Y	$207,\text{166}$	$53,\text{35}$	$139,\text{118}$	$58,\text{35}$
	N	$23,\text{24}$	$15,\text{11}$	$44,\text{40}$	?

Table 4. Table 4: MLE estimates of the parameters using the incomplete data in Table 3 , and allowing a different capture propensity parameter θ j subscript 𝜃 𝑗 \theta_{j} for each list.

${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_{12}$	${\hat{\alpha }}_{13}$	${\hat{\alpha }}_{23}$	${\hat{\gamma }}_E$	${\hat{\gamma }}_U$	${\hat{\theta }}_1$	${\hat{\theta }}_2$	${\hat{\theta }}_3$
0.056	0.867	0.175	0.575	0.954	0.863	580	459	-0.002	-0.241	0.152

Table 5. Table 5: MLE estimates of the parameters using the incomplete data in Table 3 , fitting the model in ( 9 ) under the null hypothesis H 0 : θ 1 = θ 2 = θ 3 = θ : subscript 𝐻 0 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 𝜃 H_{0}:\theta_{1}=\theta_{2}=\theta_{3}=\theta .

${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_{12}$	${\hat{\alpha }}_{13}$	${\hat{\alpha }}_{23}$	${\hat{\gamma }}_E$	${\hat{\gamma }}_U$	$\hat{\theta }$
0.074	0.750	0.289	0.572	0.951	0.848	580	459	-0.033

Table 6. Table 6: MLE estimates of the parameters using the incomplete data in Table 3 . In the first row we used the null hypothesis H 0 : θ = 0 : subscript 𝐻 0 𝜃 0 H_{0}:\theta=0 , while in the second row we used H 0 : γ E / γ U = 1.256 : subscript 𝐻 0 subscript 𝛾 𝐸 subscript 𝛾 𝑈 1.256 H_{0}:\gamma_{E}/\gamma_{U}=1.256 so γ ^ E subscript ^ 𝛾 𝐸 \hat{\gamma}_{E} is not directly estimated bu just set to be 1.256 ⋅ γ ^ U ⋅ 1.256 subscript ^ 𝛾 𝑈 1.256\cdot\hat{\gamma}_{U} .

${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_{12}$	${\hat{\alpha }}_{13}$	${\hat{\alpha }}_{23}$	${\hat{\gamma }}_E$	${\hat{\gamma }}_U$	$\hat{\theta }$
0.055	0.730	0.266	0.574	0.953	0.852	579	461
0.055	0.730	0.266	0.574	0.953	0.852		460	0.004

Table 7. Table 7: MLE estimates of the parameters in the model. We used a saturated log-linear model to estimate the count of the missing cells. In the second line the estimates are under the null θ = 0 𝜃 0 \theta=0 .

${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_{12}$	${\hat{\alpha }}_{13}$	${\hat{\alpha }}_{23}$	$\hat{\theta }$
0.104	0.773	0.305	0.583	0.961	0.861	-0.092
0.052	0.720	0.247	0.584	0.961	0.868

Table 8. Table 8: MLE estimates of the parameters in the extended model. In the first row we have the estimates of the parameters when using capture-recapture to first estimate the population size (and it is very similar to the first row in Table 7 ). In the second row we use the Poisson assumption on the total count of cases, and the estimates are similar to the ones in Table 5 .

${\hat{\alpha }}_1$	${\hat{\alpha }}_2$	${\hat{\alpha }}_3$	${\hat{\alpha }}_{12}$	${\hat{\alpha }}_{13}$	${\hat{\alpha }}_{23}$	${\hat{\gamma }}_W$	${\hat{\gamma }}_B$	$\hat{\mu }$	$\hat{\sigma }$
0.118	0.721	0.345	0.733	0.971	0.819			-0.110	0.0009
0.074	0.750	0.288	0.572	0.952	0.848	580	459	-0.033	0.0043

Equations

$$\frac{\mathbb{P}(e=E|Y=1)}{\mathbb{P}(e=U|Y=1)}\cdot \frac{\mathbb{P}(e=U|Y=0)}{\mathbb{P}(e=E|Y=0)}=\frac{\mathbb{P}(Y=1|e=E)}{\mathbb{P}(Y=0|e=E)}\cdot \frac{\mathbb{P}(Y=0|e=U)}{\mathbb{P}(Y=1|e=U)}.$$

$$OR=\frac{\frac{N^E}{N^E+N^U}\cdot \frac{T^U-N^U}{T^U+T^E-N^U-N^E}}{\frac{N^U}{N^E+N^U}\cdot \frac{T^E-N^E}{T^U+T^E-N^U-N^E}}=\frac{N^E\cdot (T^U-N^U)}{N^U\cdot (T^E-N^E)}.$$

$$RR=\frac{\mathbb{P}(Y=1|e=E)}{\mathbb{P}(Y=1|e=U)}=\frac{\mathbb{P}(e=E|Y=1)\cdot \mathbb{P}(e=U)}{\mathbb{P}(e=U|Y=1)\cdot \mathbb{P}(e=E)}=\frac{N^E\cdot T^U}{N^U\cdot T^E};$$

$$X^{e}_{ij}=\left\{\begin{array}[]{ll}1&\mbox{if individual $i$ is ascertained by list $j$}\\ 0&\mbox{otherwise}\end{array}\right.$$

$$\mathbb{P}(X_{ij}^e=1)=\frac{e^{{\theta }_{ij}^e+{\alpha }_j}}{1+e^{{\theta }_{ij}^e+{\alpha }_j}}$$

$$p_{xyz}^E:=\mathbb{P}\left(X_{i1}^E=x,X_{i2}^E=y,X_{i3}^E=z\right)=\frac{e^{x\cdot ({\theta }_1+{\alpha }_1)}}{1+e^{{\theta }_1+{\alpha }_1}}\cdot \frac{e^{y\cdot ({\theta }_2+{\alpha }_2)}}{1+e^{{\theta }_2+{\alpha }_2}}\cdot \frac{e^{z\cdot ({\theta }_3+{\alpha }_3)}}{1+e^{{\theta }_3+{\alpha }_3}}$$

$$p_{xyz}^U:=\mathbb{P}\left(X_{i1}^U=x,X_{i2}^U=y,X_{i3}^U=z\right)=\frac{e^{x\cdot {\alpha }_1}}{1+e^{{\alpha }_1}}\cdot \frac{e^{y\cdot {\alpha }_2}}{1+e^{{\alpha }_2}}\cdot \frac{e^{z\cdot {\alpha }_3}}{1+e^{{\alpha }_3}}$$

$$\mathcal{L}(M^e|N^e)\sim \text{Multinomial}(M^e|N^e,p_{xyz}^e)$$

$$\mathcal{L}(M^E,M^U|N^E,N^U)=\mathcal{L}(M^E|N^E)\cdot \mathcal{L}(M^U|N^U).$$

$$\mathcal{L}(M_{obs}^e)$$

$$=\sum _{m=0}^{\infty }\mathcal{L}(M_{obs}^e,M_{000}^e=m,N^e=N_{obs}^e+m)$$

$$=\sum _{m=0}^{\infty }\mathcal{L}(M_{obs}^e,M_{000}^e=m|N^e=N_{obs}^e+m)\cdot \mathbb{P}(N^e=N_{obs}^e+m)$$

$$\mathcal{L}(M_{obs}^e)$$

$$=\frac{{\prod }_{x,y,z}{(p_{xyz}^e)}^{M_{xyz}^e}\cdot {\gamma }_e^{N_{obs}^e}\cdot e^{-{\gamma }_e}}{{\prod }_{x,y,z}{M}_{xyz}^e!}\cdot \sum _{m=0}^{\infty }\frac{({\gamma }_e\cdot p_{000}^e{)}^m}{m!}$$

$$=\frac{{\prod }_{x,y,z}{(p_{xyz}^e)}^{M_{xyz}^e}\cdot {\gamma }_e^{N_{obs}^e}\cdot e^{{\gamma }_e(p_{000}^e-1)}}{{\prod }_{x,y,z}{M}_{xyz}^e!}$$

$$\mathcal{L}(M_{obs}^E,M_{obs}^U)=\mathcal{L}(M_{obs}^E)\cdot \mathcal{L}(M_{obs}^U).$$

$$p_1^E:=$$

$$\Rightarrow$$

$$p_{2,x}^E:=$$

$$\Rightarrow$$

$$p_{3,x,y}^E:=$$

$$\Rightarrow$$

$$p_{xyz}^E$$

$$p_{xyz}^U$$

$$\log\left(\mathcal{L}(M^{e}_{obs})\right)\propto\sum_{x,y,z}M_{xyz}^{e}\cdot\log\big{(}\gamma_{e}\cdot p_{xyz}^{e}\big{)}+\gamma_{e}(p_{000}^{e}-1).$$

$$\Lambda =-2({\ell }_0-\ell )$$

$$\mathbb{E}\left[M_{000}^E\right]$$

$$\mathbb{E}\left[M_{000}^U\right]$$

$$\theta <0\iff p_{000}^E>p_{000}^U\iff N^E/N^U>\mathbb{E}[N_{obs}^E]/\mathbb{E}[N_{obs}^U].$$

$$\mathcal{L}(M^e|N^e)={\int }_{-\infty }^{\infty }\mathcal{L}(M^e|N^e,\theta )\cdot {\varphi }_e(\theta )d\theta$$

$$\mathcal{L}(M_{obs}^e)$$

$$={\int }_{-\infty }^{\infty }\frac{{\prod }_{x,y,z}{(p_{xyz}^e)}^{M_{xyz}^e}\cdot {\gamma }_e^{N_{obs}^e}\cdot e^{{\gamma }_e(p_{000}^e-1)}}{{\prod }_{x,y,z}{M}_{xyz}^e!}\cdot {\varphi }_e(\theta )d\theta$$

$$\ell (M_{obs}^E,M_{obs}^U)=\sum _{x,y,z}\left(M_{xyz}^E\log ({\gamma }_E\cdot p_{xyz}^E)+M_{xyz}^U\log ({\gamma }_U\cdot p_{xyz}^U)\right)+{\gamma }_E(p_{000}^E-1)+{\gamma }_U(p_{000}^U-1)$$

$$\left(\theta =0,{\alpha }_1=1,{\gamma }_E=\frac{1+e}{2},{\gamma }_U=\frac{1+e}{2}\right)\text{and}\left({\theta }^{\prime }=0,{\alpha }_1^{\prime }=0,{\gamma }_E^{\prime }=e,{\gamma }_U^{\prime }=e\right)$$

$$\left(\theta =1,{\alpha }_1=0,{\gamma }_E=1+e,{\gamma }_U=2e\right)\text{and}\left({\theta }^{\prime }=0,{\alpha }_1^{\prime }=1,{\gamma }_E^{\prime }=1+e,{\gamma }_U^{\prime }=1+e\right)$$

$${\gamma }_e\cdot p_{ij}^e={\gamma }_e^{\prime }\cdot {p^{\prime }}_{ij}^e\text{for any exposure and any }i,j\text{s.t. }i+j>0.$$

Full text

A Test for Differential Ascertainment in Case-Control Studies with Application to Child Maltreatment

Matteo Sordello and Dylan S. Small

Department of Statistics, Wharton School, University of Pennsylvania

We propose a method to test for the presence of differential ascertainment in case-control studies, when data are collected by multiple sources. We show that, when differential ascertainment is present, the use of only the observed cases leads to severe bias in the computation of the odds ratio. We can alleviate the effect of such bias using the estimates that our method of testing for differential ascertainment naturally provides. We apply it to a dataset obtained from the National Violent Death Reporting System, with the goal of checking for the presence of differential ascertainment by race in the count of deaths caused by child maltreatment.

00footnotetext: MSC 2010 subject classifications: 62F03, 62P25

1 The Problem of Differential Ascertainment in Case-Control Studies

In case-control studies, one wants to understand the association between some risk factor and a rare outcome. The odds of a risk factor are compared among cases (subjects with the outcome) vs. controls (subjects without the outcome). When it is possible to get a complete count of the cases or to obtain a random sample of the cases, then estimating the prevalence of the risk factor among the cases is straightforward. But for some settings, it is not possible to get a complete count or take a random sample of the cases. This problem does not arise strictly in case-control studies, but in any situation in which one wants to count the number of individuals that belong to a certain group. In such circumstances we say that we are in the presence of imperfect ascertainment. It can happen in a multitude of different contexts, for example when counting the number of people addicted to a certain drug (Domingo-Salvany et al.,, 1995), civilians that died during a war (Ball et al.,, 2002; Price et al.,, 2014) or people that suffer from diabetes (International Working Group for Disease Monitoring and Forecasting,, 1995). A further complication might arise if the group of people under consideration is divided into two subpopulations, according to a variable. In such cases, the ascertainment can be not only imperfect, but possibly easier in one of these two subpopulations, and harder in the other. For example, one can think that there is a threshold on the variable age under which it is easier to ascertain if a person is a drug addict, or that gender is a relevant factor in the ascertainment of civilian deaths or of individuals suffering from diabetes. In practical terms, this means that the ratio of the observed counts in the two subpopulations is different from the real underlying ratio. Here we say that we are in the presence of differential ascertainment, a situation that is harder to detect than simple imperfect ascertainment.

In this work, we focus on a specific situation where we suspect that differential ascertainment could be involved, and develop the machinery to perform a rigorous statistical analysis. The question we want to answer concerns deaths caused by maltreatment of children aged 10 or less. Notice that here we are in the presence of imperfect ascertainment, since it is not always possible to assess whether or not a death has been caused by child maltreatment (Herman-Giddens et al.,, 1999). We are interested in knowing whether the deaths caused by child maltreatment are differentially ascertained according to the variable race. There are two main reasons we are interested in this question. The first one is that most of the research done in this area has the goal of providing guidelines that can be used to develop and evaluate intervention strategies, ultimately to reduce the number of deaths caused by child maltreatment. If the real proportions in the population are different from the ones observed, then the resources put forward to solve the problem may be poorly allocated and less effective than they should be. The second reason is that differential ascertainment can lead case-control studies to provide false conclusions. To explain this, we introduce some notation for case-control studies (for a thorough discussion of case-control studies see Keogh and Cox, (2014)). For a case-control study, let $e\in\{E,U\}$ be the bivariate exposure, $E$ stands for exposed and $U$ for unexposed, and $Y\in\{0,1\}$ be the outcome of interest, where $Y=1$ is the indicator of the cases and $Y=0$ of the controls. We would like to infer the distribution of $Y|e$, and in particular to compute the risk ratio (RR), which is defined as $\mathbb{P}(Y=1|e=E)/\mathbb{P}(Y=1|e=U)$. Unfortunately, given that it is structured retrospectively, a case-control study can only provide information on the distribution of $e|Y$. The following equality comes to our aid, creating a link between these distributions.

[TABLE]

The left hand side is called the odds ratio (OR), and it is a fundamental quantity to assess the association between exposure and outcome. Let $T^{E}$ be the total number of exposed individuals and $T^{U}$ the total number of unexposed individuals in the whole population of interest (so counting both cases and controls). We assume that $T^{E}$ and $T^{U}$ are known. Now let $N^{E}$ be the number of exposed cases and $N^{U}$ the number of unexposed cases. The odds ratio is then

[TABLE]

We assume throughout the paper that it is possible to obtain a random sample from the two groups of exposed and unexposed controls, or simply that the ratio $(T^{U}-N^{U})/(T^{E}-N^{E})$ is known or can be estimated. Plugging the estimate of the ratio $(T^{U}-N^{U})/(T^{E}-N^{E})$ into the odds ratio (1) reduces the problem of estimating OR to estimating the ratio between the exposed and unexposed cases, $N^{E}/N^{U}$. The distribution of OR will depend on the randomness in the estimates for $N^{E}/N^{U}$ and $(T^{U}-N^{U})/(T^{E}-N^{E})$ combined.

If there is imperfect ascertainment of the cases, the real values of $N^{E}$ and $N^{U}$ are unknown. If one makes the assumption that the real proportion of exposed to unexposed cases is consistently estimated by the observed proportion, then plugging the observed proportion into (1) provides a consistent estimate of the OR. However, if there is differential ascertainment, we know by definition that the use of the observed ratio instead of the real $N^{E}/N^{U}$ will lead to some bias in the computation of OR, which can be severe, as we will show in Section 2.2.

Another situation where we would care about the count of exposed and unexposed cases is in a case-cohort study (Prentice,, 1986), in which the cases have imperfect ascertainment but the controls have perfect ascertainment. In case-cohort studies, a census of cases in a cohort is taken and then a sub-cohort of the cohort is sampled (the sample includes both cases and control). For a case-cohort study with perfect ascertainment of cases, one can estimate the risk ratio by observing that

[TABLE]

the quantity $N^{E}/N^{U}$ is obtained from the census of cases and $T^{U}/T^{E}$ is estimated from the sub-cohort which is a random sample of the population in the cohort. If there is imperfect ascertainment of cases, then assuming that the cases are only a small fraction of the population in the cohort, the exposure ratio in the population, $T^{U}/T^{E}$, can still be estimated from the sub-cohort and again our attention focuses on finding an estimate for $N^{E}/N^{U}$.

Problems of this nature, where the total size of a population is unknown, have been deeply studied in the capture-recapture literature. Bishop et al., (1975) and White, (1982) wrote extensively on incomplete tables and how to estimate the size of a closed population. Cormack, (1989) considered various extensions, including the situation where the target population is open. Braeye et al., (2016) compared different capture-recapture estimators in situations where their assumptions are violated, while Seber et al., (2000) and Lee et al., (2001) discussed the problem of matches between lists that are not recorded correctly. One commonly used model for the catchability of an individual in the capture-recapture literature is the Rasch model (Rasch,, 1960; Agresti,, 1994). Here we consider a simple modification of the original Rasch model. The model results in $N^{E}$ and $N^{U}$ being realizations of a Poisson random variables with means $\gamma_{E}$ and $\gamma_{U}$ respectively. We estimate those means, and introduce the parameters $\theta_{j}$ representing the amount of differential ascertainment in each list. We use these $\theta_{j}$’s to test if the differential ascertainment is significantly different from zero. If it is significant, then we use $\hat{\gamma}_{E}$ and $\hat{\gamma}_{U}$ in the computation of OR or RR. If it is not significant, then the observed counts of exposed and unexposed cases can be used instead of the estimated means.

2 The Model

As anticipated before, we consider a bivariate exposure $e\in\{E,U\}$ and a total number of cases $N^{e}$ in each group. We first introduce the classical version of the Rasch model (Rasch,, 1960), which relies on two important assumptions.

Assumption 2.1** (Individual Independence).**

The ascertainment of each individual in any list is independent from the ascertainment of any other individual.

Assumption 2.2** (List Independence).**

For any individual, the ascertainment happens in a list independently from any other list.

Assumption 2.2 is not completely reasonable, since we can imagine that being ascertained in one or more lists makes an individual more or less prone to be ascertained again. This mechanism is explained in Bishop et al., (1975) under the name of social visibility/invisibility or trap fascination/avoidance, depending on the context. In Section 2.1 we will get rid of this assumption and explain how to include dependence among the lists. We will sometimes refer to the concept of ascertainment as capture, given that our framework is essentially the same as the capture-recapture one.

The number of lists used is denoted by $J$. Conditionally on the value $N^{e}$, the ascertainment matrix $X^{e}$ has dimension $N^{e}\times J$, and, for $i=1,...,N^{e}$ and $j=1,...,J$, each entry takes values

[TABLE]

These Bernoulli random variables take value $1$ with probability

[TABLE]

The parameters $\theta^{e}_{ij}$ represent the individual propensity for an individual with exposure $e$ to be captured from list $j$, and $\alpha=(\alpha_{1},...,\alpha_{J})$ are the capture strengths of each list, which are assumed to be independent from the exposure. This model naturally creates, for each exposure $e$, a contingency table $M^{e}$ with $2^{J}$ entries, each entry being the number of individuals that are ascertained in a specific combination of lists. We denote $M^{e}_{k_{1},...,k_{J}}$, where $k_{1},...,k_{J}\in\{0,1\}$, to be the count of individuals that are ascertained in all the lists with corresponding index that takes value $1$, and not ascertained in the lists with index [math]. For example, when $J=3$ (the setting we will use from now on, since it is what we will be dealing with in the NVDRS dataset), $M^{e}_{101}$ is the count of individuals that appear in list $1$ and list $3$, but not in list $2$. In Table 1 we show a contingency table for the three lists setting.

We consider two simplified versions of this model. In one we assume that the value of the individual capture propensity is zero for all the unexposed cases, and only depends on the list, but not the specific individual, for the exposed. In this framework we set $\theta^{U}_{ij}=0$ and $\theta^{E}_{ij}=\theta_{j}$. Then we further simplify the model by only allowing a single parameter $\theta$ for all the exposed cases and all the lists. While the former setting allows more flexibility in the model, the second one provides a clear answer to the question of the presence of differential ascertainment in the data (see Remark 2.1). We will explain how to test if it is necessary to consider different values for each $\theta_{j}$. The following description of the model will use the notation where the capture propensity is list dependent for the exposed cases.

In the Appendix we introduce another variation, where we allow each individual to have a different propensity to be captured. We make the assumptions that exposed cases have individual propensity $\theta_{ij}^{E}\sim N(\mu,\sigma^{2})$, and unexposed individuals have $\theta_{ij}^{U}\sim N(0,\sigma^{2})$. The magnitude of the newly introduced parameter $\sigma$ can tell us if there is significant variability among the individuals. We will see, however, that in the study of the NVDRS data we don’t really need to introduce this variability, since our estimate of $\sigma$ is going to be extremely small and all the other parameters very close to the ones estimated with the initial setting.

Using Assumption 2.2, we compute the probability for each exposed individual to appear in a specific cell identified by subscripts $xyz$, which is just the product of the single lists’ probabilities

[TABLE]

while for unexposed cases we have

[TABLE]

Thanks to Assumption 2.1, the likelihood of a complete table $M^{e}$ with total count $N^{e}$ is given by a multinomial distribution

[TABLE]

where the probabilities are given by (3) or (4). A consequence of Assumption 2.1 is that the contingency tables associated with different exposures are independent. Then the total likelihood of the two complete contingency tables we are considering is

[TABLE]

Notice that up to now we were able to condition on the total number of cases $N^{e}$, since we were dealing with contingency tables where no cell was missing. In real applications, however, the contingency tables are not complete, and the cells $M^{e}_{000}$ are missing, since they represent the elements that are not ascertained in any list. This means that $N^{e}$ is now unknown, and the definition of the likelihood is not as immediate as before. We then add the assumption that, for any $e$, the count $N^{e}\sim\text{Poisson}(\gamma_{e})$. If we let $M_{obs}^{e}$ be the contingency table without the unobserved cell and $N_{obs}^{e}$ the corresponding count of all the observed cases, then the likelihood of the observed table is

[TABLE]

where $\left\{M^{e}_{obs},M^{e}_{000}=m|N^{e}=N^{e}_{obs}+m\right\}$ is a complete contingency table with $m$ observations in the bottom-right cell, hence distributed according to a multinomial distribution with probabilities given in (3) or (4), and $N^{e}$ is Poisson with mean $\gamma_{e}$. In its extended form, the likelihood of one incomplete contingency table is then

[TABLE]

where the products are only on the $x,y,z$ such that $x+y+z>0$. The likelihood of the two incomplete contingency tables is, again using Assumption 2.1,

[TABLE]

2.1 Dependence among lists

While the independence among the individuals is an important assumption, we can extend this model to include dependence among the lists. We consider the dynamical Rasch model (Verhelst and Glas,, 1993), which eliminates the independence assumption among the lists with a simple modification of the probabilities in (2). Instead of having them fixed, we generate them sequentially based on the previous realizations of the ascertainment variables. These are the same for each individual, and are defined, for those that are exposed, as

[TABLE]

and if the number of lists is greater we can keep generating $\mathbb{P}(X^{E}_{ik}=1)$ and $X^{E}_{ik}$ based on the previous values of $X^{E}_{ij},\ j=1,...,k-1$. Here we decided to consider only the two-list interactions, given by the parameters $\alpha_{12},\alpha_{13}$ and $\alpha_{23}$. It is also possible to use, as the number of lists grows, higher-order interactions. For unexposed individuals the procedure is analogous, just setting as before $\theta_{j}=0$ for $j=1,2,3$. Then we have, with the notation introduced before

[TABLE]

The model for two incomplete and independent contingency tables is then the same as before, only with the probabilities in (10) substituted for (3) and (4). The identifiability of the parameters in this model is discussed in the Appendix. To fit the model, one simply maximizes the log-likelihood obtained in (9), where

[TABLE]

A generalized likelihood ratio test can then be performed to verify if it is important to use different values for the $\theta_{j}$’s in different lists. The null hypothesis is $H_{0}:\theta_{1}=\theta_{2}=\theta_{3}$, and it is tested against the alternative that at least one $\theta_{j}$ is different from the others. The test statistic that we can use is

[TABLE]

where $\ell_{0}$ is the log-likelihood evaluated in the maximum likelihood estimator computed under $H_{0}$ (so when all $\theta_{j}$’s are equal), and $\ell$ is the same estimator found on the model with no constraints. The statistic $\Lambda$ is compared against the quantiles of a $\chi^{2}$ distribution with $2$ degrees of freedom. If no significant difference is found between the $\theta_{j}$’s, then the simpler model that only include a single value for $\theta$ can be used.

The next step after obtaining the estimates for the parameters of the model is to check if the estimated ratio $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ is substantially different from the observed one, $r:=N^{E}_{obs}/N^{U}_{obs}$. Since it is hard to understand if the difference is significant, there are two things that we can do, very similar in practice. Find a confidence interval for the ratio $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ or, if we are using the model where a single value for $\theta$ is allowed, perform hypothesis testing on it. In both cases we can use a bootstrap procedure.

• •

To test the null hypothesis $H_{0}:\theta=0$, we need to obtain the empirical distribution for the estimate $\hat{\theta}$ under $H_{0}$. We first estimate the parameters of the model (9) under the null, and then generate with those parameters new simulated populations; first the new total counts of cases using a Poisson distribution, then the new complete contingency tables using a Multinomial distribution. After excluding the unobserved cells $M^{e}_{000}$ from the tables, we fit the same model again and obtain the estimate for the parameter $\theta$. This empirical distribution is used to compute the relevant quantiles and check if the point estimate for $\hat{\theta}$ obtained with the real data is significantly far from zero.

• •

To get the confidence interval for $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ under the null hypothesis $H_{0}:\gamma_{E}/\gamma_{U}=r$, we fit again the model (9) on the real data imposing the constraint $\gamma_{E}=r\cdot\gamma_{U}$. Then, as before, create the new incomplete contingency tables, fit the model and compute and store the ratio of interest.

Remark 2.1**.**

Interpretation of the parameter $\theta$ as responsible for differential ascertainment. Let’s consider $N^{E}$ and $N^{U}$ fixed. We want to compute $N^{E}/N^{U}$ in order to estimate the odds ratio, but since the cells $M^{E}_{000}$ and $M^{U}_{000}$ are unobserved we cannot have a precise count of cases. However from (10), when $\theta_{1}=\theta_{2}=\theta_{3}=\theta$ we know that

[TABLE]

We also know that $N^{E}=\mathbb{E}\left[N_{obs}^{E}+M^{E}_{000}\right]$ and $N^{U}=\mathbb{E}\left[N_{obs}^{U}+M^{U}_{000}\right]$. This tells us that $\mathbb{E}\left[N_{obs}^{E}\right]/\mathbb{E}\left[N_{obs}^{U}\right]=N^{E}/N^{U}$ if and only if $\mathbb{E}\left[M^{E}_{000}\right]/\mathbb{E}\left[M^{U}_{000}\right]=N^{E}/N^{U}$, which is true only if $\theta=0$, thanks to the monotonicity of the denominator of $\mathbb{E}\left[M^{E}_{000}\right]$ with respect to $\theta$. The largest is $|\theta|$, the larger the difference will be between the ratio of the observed counts and the real ratio $N^{E}/N^{U}$. Notice also that

[TABLE]

If different values for the $\theta_{j}$’s are allowed, this interpretation is no longer possible. In fact, if one sets for example $\theta_{1}=\alpha_{2}-\alpha_{1}$, $\theta_{2}=\alpha_{1}-\alpha_{2}$ and $\theta_{3}=0$, then $p_{000}^{E}=p_{000}^{U}$ even though the lists ascertain the cases with different probabilities for exposed and unexposed individuals. $\square$

2.2 The Effect of Differential Ascertainment in Estimating OR

We present here a simulation that shows that the use of only the observed counts can lead to severe bias in the estimate of the odds ratio. We set $\gamma_{E}=500$ and $\gamma_{U}=1000$ to be the average count of cases among the exposed and unexposed population, and consider a single differential ascertainment parameter $\theta\in[-1,1]$ equal for all lists. We generate the two counts $N^{E}$ and $N^{U}$ from a Poisson distribution with mean $\gamma_{E}$ and $\gamma_{U}$ respectively. At every iteration, the ascertainment of the cases happens with the probabilities defined in (10). Here we only use the main effects, set to be $\alpha_{i}=-0.5$, since the interactions do not influence the counts in the unobserved cells. We consider the situation where the ascertainment is performed with $J\in\{1,3,5\}$ lists. The average of the observed ratio of exposed over unexposed individuals is recorded in Figure 1, together with the $95\%$ confidence bands and the true underlying value $0.5$. The first immediate realization that we have is that, when we created differential ascertainment by setting $\theta\neq 0$, the observed ratio can be extremely biased. It is also clear that the bias increases with $|\theta|$ in each direction (even though non symmetrically), getting far from the true ratio $\gamma_{E}/\gamma_{U}=0.5$. However, such bias gets smaller the more lists we use. This is because, when $J$ grows, even if the individual ascertainment probabilities for each list are small, the chance that at least one list will capture an individual grows, so the count of individuals that are missed by all the lists gets lower.

This analysis tells us that, in a situation where we don’t have a large number of different sources of information available, and we are unsure about the presence of differential ascertainment, it can be risky to just trust the observed ratio. When only one source of information is available, there is unfortunately not much we can do. But we show in the next section how the model just introduced can recognise if this problem if present (estimating the underlying value of $\theta$) and mitigate it using the estimated true counts.

2.2.1 How the model performs when there is differential ascertainment

We now look at the data generated in Section 2.2, and see if our model is able to detect the presence of differential ascertainment. We focus on the three-list setting, the one relevant for our data analysis in Section 4. The parameters used are the same as before. Figure 2 shows the point estimate and $95\%$ confidence bands for the estimated ratio $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ and for the estimated parameter $\hat{\theta}$ when the differential ascertainment is $\theta\in[-1,1]$. Compared to Figure 1, we see that the estimated ratio is always centered around the true ratio $0.5$, and that the width of the confidence bands is not extremely sensitive to the value of $\theta$. In the right part of the plot, instead, we see that the estimate $\hat{\theta}$ is correctly centered around the true value of $\theta$, and that again the confidence bands are not very wide. Figure 2 also allows us to see that, if we use the estimated ratio to test for the presence of differential ascertainment, we would reject the null when the $95\%$ confidence bands do not cover the average observed value, which happens when $|\theta|\gtrsim 0.4$. When using the estimated parameter $\hat{\theta}$, instead, the rejection happens already for $|\theta|\gtrsim 0.25$. This means that, when we are interested in knowing if differential ascertainment is present in our data, it is better to use the estimated $\hat{\theta}$ directly instead of the ratio of exposed and unexposed cases. In Table 2 we see the average estimates for all the parameters of the model out of $200$ simulations, together with a $95\%$ confidence interval.

One can think that a simpler alternative to this method is to directly estimate $N^{E}$ and $N^{U}$, for example using a log-linear model. Using a bootstrap procedure, similar to what we explained at the end of Section 2.1, we can then obtain a confidence interval for the ratio $N^{E}/N^{U}$ to attach to the point estimate. The literature on this topic is very rich. For a detailed explanation, see Bishop et al., (1975) Moreover, Pledger, (2000) defines a unified log-linear framework allowing for unequal capture probabilities among the population (see also Chao,, 1987; Pollock et al.,, 1990; Coull and Agresti,, 1999) and Tilling and Sterne, (1999) include the use of covariates in log-linear models. This procedure, however, even though it is valid in providing a reasonable confidence interval for the ratio of interest, does not yield directly reliable statistical evidence on the presence or not of differential ascertainment. If we use a saturated log-linear model on the data generated in Section 2.2, we can see in Figure 3 that the width of the confidence bands heavily depends on the differential ascertainment parameter. Moreover we notice that often, even when differential ascertainment is present, the confidence interval obtained in this way will cover the observed value $N^{E}_{obs}/N^{U}_{obs}$, so that we cannot reject the null hypothesis that no differential ascertainment is present. Even though this procedure is not able to tell us if we are facing differential ascertainment, we can still use it to first estimate the count of individuals in the missing cells, and then apply the complete model in (6). We implement this idea in the Appendix, and see that the results are consistent with the ones obtained with our method.

3 Description of the NVDRS Dataset

To address our question about differential ascertainment of death due to child maltreatment by race, we used a dataset obtained from the National Violent Death Reporting System (NVDRS). Here we provide a brief description of the NVDRS, but a more detailed explanation can be found in Blair et al., (2016) The NVDRS is a surveillance system created in 2002 from a pilot program called the National Violent Injury Statistics System. In 2010, 19 States were participating in the data collection; 14 additional States joined the system in 2014. The goal of the system is to monitor the occurrence of homicides, suicides, unintentional firearm deaths, deaths of undetermined intent, and deaths from legal intervention (excluding legal executions) in the US. The system uses different data sources and merges them together to obtain a more complete description of the event. The primary data sources used for the system are death certificates, coroner/medical examiner reports (including toxicology reports) and law enforcement reports. Coroner/medical examiner reports and law enforcement reports contain narratives that can provide detailed information regarding the circumstances and characteristics of incidents. The manner of death is captured by International Classification of Diseases, Tenth Revision (ICD-10) codes. Other relevant information collected by the system are demographics (for victims and suspects), method of injury (eg, sharp instrument), place of injury, information on the victim’s life stressors, the victim-suspect relationship (eg, spouse or partner), the presence of intimate partner violence (IPV), toxicology (for victims), weapon information and whether other crimes occurred that are related to the incident (eg, robbery followed by homicide).

As already mentioned, we focus our attention on the deaths of children aged 10 or less, and consider as cases the deaths that were caused by child maltreatment. We are interested in determining if race played a role in the ascertainment of the cases. The timeframe of interest is 2010 to 2015. We follow Leeb et al., (2008) (pg. 11) in considering child maltreatment to be “any act or series of acts of commission or omission by a parent or other caregiver that results in harm, potential for harm, or threat of harm to a child.”

The dataset consists of 1983 deaths of children from ages 0 to 10 years old and 62 variables that describe the deaths. Among these variables, there is one that describes the police report, one for the coroner report and others linked to the death certificate, for example the code for cause of death, the manner of death and some descriptive variables. From these variables, we constructed three separate lists, all of them trying to count the total number of deaths caused by child maltreatment. The death certificate (DC) combines information from the ICD-10 codes, the death causes and the caregiver being identified as a suspect. For the police report (LE), and coroner report (CME), there are in the dataset two variables that define which deaths were considered to be caused by child maltreatment according to the police report and coroner’s report respectively. Unfortunately, these variables cannot be trusted to capture all of the cases. A further analysis, described in the Appendix, was needed to incorporate all the cases that were described in the narratives but not automatically reported in the lists.

We want to point out that we are not the first ones to propose a study related to case-controls on the NVDRS data. Koziol-Mclain et. al., (2006) analysed the risk factors for femicide-suicide in a case-control study that made use of data from $11$ cities in the NVDRS dataset.

4 Analysis of the NVDRS Data

After the manual scrutiny of the three lists, the total number of cases observed is 968. As anticipated, we define the exposure based on the variable race, and consider exposed the cases where the victim is white and unexposed the cases where the victim is black. The observed ratio of white and black cases is $N^{E}_{obs}/N^{U}_{obs}=539/429=1.256$. Table 3 is the contingency table for the observed cases, divided by whites and blacks (written in boldface). We first want to test if there is evidence that the individual ascertainment parameters $\theta_{j}$ are list specific, or if instead we can simply assume that there exist a single parameter $\theta$, possibly responsible for differential ascertainment.

To do so, we fit the incomplete model (9) to these data using the probabilities in (10). The order of the lists is the one presented in Table 1, so DC represents list 1, LE list 2 and CME list 3. With this notation, we get the parameter estimates in Table 4 and a value for the log-likelihood of $\ell=3457.38$.

There seems to be a fairly large difference in the $\theta_{j}$’s, especially between $\theta_{2}$ and $\theta_{3}$. To test if this difference is significant, we fit on the same data the constrained model, which assumes $H_{0}:\theta_{1}=\theta_{2}=\theta_{3}$. Here, we get the parameter estimates in Table 5, and a value for the log-likelihood of $\ell_{0}=3455.60$. The generalised likelihood ratio statistic is then $\Lambda=2\cdot(3457.38-3455.60)=3.56$, which has p-value of $0.169$ on a $\chi^{2}$ distribution with $2$ degrees of freedom. Given this result, we cannot reject the null hypothesis that all the $\theta_{j}$’s have the same value. We will then continue our analysis assuming from now on that there is only one parameter $\theta$. We want to test if there is any evidence of differential ascertainment, i.e., if such $\theta$ is significantly different from [math], and if the estimated ratio is different from the observed one.

The estimated values $\hat{\gamma}_{E}$ and $\hat{\gamma}_{U}$ in Table 5 tell us that the unobserved values in the missing cell are $580-539=41$ white cases and $459-429=30$ black cases. Our estimate of the true underlying ratio of exposed and unexposed cases is $\hat{\gamma}_{E}/\hat{\gamma}_{U}=1.264$, which is close to the observed one (yet slightly bigger). This is consistent with the fact that $\hat{\theta}$ is very close to [math] and negative, so we get the result anticipated in Remark 2.1.

A confidence interval for $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ under the null hypothesis $H_{0}:\gamma_{E}/\gamma_{U}=1.256$ can be obtained in the way explained is Section 2.1, by first fitting the incomplete model (9) under the constraint $\gamma_{E}=1.256\cdot\gamma_{U}$ (the maximum likelihood estimators are reported in Table 6) and then creating new populations. When fitting the model on these new populations, we get a $95\%$ confidence interval of $[1.097,1.430]$, which means that $H_{0}$ is not rejected. The empirical distribution for the ratio is reported in the right panel of Figure 4, together with the estimated ratio.

In the left panel of Figure 4, instead, we have the empirical distribution of $\hat{\theta}$ under the null hypothesis $H_{0}:\theta=0$. It is obtained fitting the model on the NVDRS data using the probabilities $p^{U}_{xyz}$ in (10) for both exposed and unexposed cases (the estimates are reported in Table 6). Using those values, we generate new incomplete populations, and fit the incomplete model to obtain a $95\%$ confidence interval for $\theta$ under the null of $[-0.199,0.193]$. Notice that, keeping all other parameters fixed, this range of values for $\theta$ can make the probability $p_{000}^{E}$ range in $[0.046,0.093]$.

From this analysis we can conclude that there is no evidence of differential ascertainment by race in the NVDRS data that we analysed.

5 Discussion

In this work, we developed a method to test for the presence of differential ascertainment when data are collected in multiple lists. We showed through a simulation the dangerous effect that differential ascertainment can have on estimating the odds-ratio in observational studies, and used the proposed method to detect when it is present and correct its effect. We applied it to a dataset obtained from the National Violent Death Reporting System, where the cases of child maltreatment that led to death of children from 0 to 10 years old are recorded. When considering race as an exposure, the observed ratio of the cases is $1.256$. We showed that the estimated ratio is very close to the observed one, and used the confidence intervals for $\hat{\theta}$ and for $\hat{\gamma}_{E}/\hat{\gamma}_{U}$ to test whether there is any significant presence of differential ascertainment.

In our analysis we found no evidence of differential ascertainment by race in the NVDRS data. The point estimate for the differential ascertainment parameter $\theta$ is $-0.033$ which is contained in the $95\%$ confidence interval for the parameter under the null $\theta=0$. The confidence interval tells us that plausible values for the amount of differential ascertainment are approximately $-0.2$ to $0.2$. In the extreme cases, the probability of not ascertainment of an element with exposure in each of the three list can go from $0.046$ to $0.093$. This tells us that, although we cannot guarantee that differential ascertainment is not present in our data, we have no evidence of a significant impact. Notice that in this work we focused on the situation in which only the cases can be affected by differential ascertainment. However, this model can also be used in the situation where the controls might suffer from the same complication. In that case, one needs to keep into consideration the randomness in the estimate of the ratios for both cases and controls, to get a confidence interval for the odds ratio (1).

Apart from the direct application of the method on multiple lists from the whole population of interest, one can think about using it as an exploratory analysis on a smaller sample from the whole population, when it may be expensive to collect data from multiple lists rather than just one. In this case, an initial test for differential ascertainment in the small sample can give an idea of the possible bias that we would get if we were to only ascertain the data using a single list.

Acknowledgements

The authors are grateful to Dr. Joanne Klevens for useful discussions and suggestions on the NVDRS data.

Appendix A Modification of the NVDRS dataset

Reading the police and coroner narratives, it gets clear that the completeness of the variables LE_DeathAbuse and CME_DeathAbuse cannot be trusted. In multiple circumstances the narrative explicitly classifies an individual as a victim of child abuse or maltreatment, while the corresponding variable does not report it as such. For this reason we had to manually check all the narratives of the individuals that were not initially classified as cases. Below are the details of how the cases were defined in each of the three lists.

• DC:

For creating the death certificate list, we didn’t have access to a full narrative, but only to a group of variables. Following a procedure similar to the one described in Klevens and Leeb, (2010) we coded a death as being caused by child maltreatment if (1) the ICD-10 code was Y06 or Y07 (where Y06 refers to neglect and abandonment and Y07 to other maltreatment syndromes), (2) the ICD-10 code was one of assault, i.e., X85, X86, X87, X88, X89, X90, X91, X92, X93, X94, X95, X96, X97, X98, X99, Y00, Y01, Y02, Y03, Y04, Y08 or Y09 and the variable DeathCause contained one of the following words: abandon, head trauma, shaken, abusive or abuse, and (3) the ICD-10 code was one of assault and the suspect was the caregiver.

This list captured 298 white cases and 236 black cases, for a total of 534.

• LE:

We started considering as cases the individuals for which LE_DeathAbuse = 1. Moreover we also read the narratives of the victims that were initially coded as controls, to account for the inefficiency of the variable LE_DeathAbuse. Here we followed the definition of abuse contained in Leeb et al., (2008) and in particular we focused on the following aspects:

• –

The perpetrator has to be a caregiver of the victim, i.e., a person that “at the time of the maltreatment is in a permanent or temporary custodial role.”

• –

When considering acts of commission (Child Abuse) the action should be intentional, but the harm to the child need not be the intended consequence. The acts of omission (Child Neglect) are more rare because NVDRS tends to not capture them, but also in this case the harm to the child need not be the intended consequence.

Initially there were 495 cases captured by the variable LE_DeathAbuse regarding individuals of race black or white. We manually added another 266, bringing the total number to 761. Of these, 413 are white and 348 are black.

• CME:

We used the same rule as for the police narrative, initially using the variable CME_DeathAbuse and then reading the remaining narratives.

Here we initially had 534 cases captured by the variable CME_DeathAbuse, and we manually added another 277. Of the 811 total cases, 457 are white and 354 are black.

Appendix B Data Analysis Estimating the Missing Counts First

As anticipated in Section 2.2.1, we can use a log-linear model to first estimate the total counts $N^{E}$ and $N^{U}$, and then apply the complete model in (6). We decided to adopt the saturated model, to impose no constraints on the parameters. It is also possible here to use the model selection procedure proposed in Ball et al., (2002), in case one wants to employ a more parsimonious model. The number of unobserved exposed cases is estimated to be $47$, while for the unobserved unexposed cases is $26$. The estimated ratio is then $586/455=1.288$. To get the $95\%$ confidence interval for the ratio $N^{E}/N^{U}$ using this log-linear model, we resample the observed data using a multinomial distribution with the empirical probability for each cell. Then estimate the count of the missing cells using the saturated (or selected) log-linear model. The confidence interval that we would obtain is $[1.16,1.43]$.

If we want to also study the presence of differential ascertainment, instead, we can assume that the total counts of cases is $N^{E}=586$ exposed and $N^{U}=455$ unexposed, which will be considered fixed. We can now fit the complete model (6), in which the only source of randomness comes from the multinomial distributions $M^{e}|N^{e}$ of the cases in the two contingency tables. Again we expect the estimate for $\theta$ to be negative and small in magnitude, as per Remark 2.1. The estimates that we get from this model are reported in the first line of Table 7.

We now want to find the empirical distribution of $\hat{\theta}$ under $H:\theta=0$. We start by fitting the model with $\theta=0$ on the complete tables, with the numbers $47$ and $26$ in the bottom right cells of the exposed and unexposed contingency tables respectively. The estimates we got are reported in the second line of Table 7. Using them, we generate $1000$ new populations from two multinomial distributions with probabilities given by (10). We then fit the complete model on these populations, and get the empirical distribution for $\hat{\theta}$ which is plotted in Figure 5, together with the quantiles that define a $90\%$ and $95\%$ confidence interval. Let $\alpha=0.05$, the relevant quantiles are $q_{\alpha/2}=-0.157,\ q_{\alpha}=-0.127,\ q_{1-\alpha}=0.129$ and $q_{1-\alpha/2}=0.149$.

Appendix C An Extension of the Model

As anticipated in Section 2, we can consider an extension of the model that keeps into account the individual differences in the propensity to be captured. Instead of simply setting $\theta_{i}^{U}=0$ and $\theta_{j}^{E}=\theta$ for any $i=1,...,N^{E}$, $j=1,...,N^{U}$, we allow these values to vary, and assume they are normally distributed, $\theta_{i}^{U}\sim N(0,\sigma^{2})$ and $\theta_{i}^{E}\sim N(\mu,\sigma^{2})$. Each value is generated independently, and the parameter $\mu$ now represents the mean amount of differential ascertainment.

We start by examining the situation in which the two contingency tables are complete. In our main example on NVDRS data, this means that we apply the log-linear model described before, and add the values 85 and 63 in the empty cells of the exposed and unexposed contingency tables respectively. We can consider the total counts $N^{e}$ as fixed, and the likelihood of each contingency table is then:

[TABLE]

Here $\mathcal{L}(M^{e}|N^{e},\theta)$ is the multinomial likelihood where the total count is $N^{e}$ and the probabilities are given in (10). The density $\phi_{e}(\theta)$ is the normal density with mean $\mu$ for the exposed, and [math] for the unexposed cases. When fitting this likelihood on the NVDRS data, again using Assumption 2.1, we get the estimates reported in the first row of Table 8.

We happily notice that these values are extremely similar to the ones in Table 7, and that $\hat{\sigma}$ is very small. This suggests that the original model with $\theta_{i}^{U}=0$ and $\theta_{i}^{E}$ fixed is appropriate for the data we are studying. To confirm this theory, we consider the natural extension of our original model, where the contingency tables are incomplete and the individual capture strengths are normally distributed. We assume as before that $N^{e}\sim\text{Poisson}(\gamma_{e})$. Starting from (2), and treating the complete contingency table as in (12), we have

[TABLE]

Fitting this likelihood to our data, we get the estimates in the second line of Table 8. Again it looks like allowing variability in the individual capture strength is not necessary, since the estimates are nearly the same as the ones in Table 5, and $\hat{\sigma}$ is very small.

Appendix D Identifiability of the Parameters

Here we discuss the problem of identifiability in our model, and show that when the number of lists is $J\geq 2$ all the parameters are identifiable. Identifiability means that if we have two distinct set of parameters, then the corresponding probability distributions of our model are different. To prove identifiability one usually assumes that two set of parameters induce the same distribution and proves that this implies that the parameters are actually the same in the two set. We denote the second set of parameters with a prime symbol ’, and the probabilities defined with those parameters by ${p^{\prime}}_{ij}^{e}$.

The likelihood for each incomplete contingency table when the number of lists is $J=3$ is in (2). When the number of lists is different, only the sum needs to be changed to take into account all the observed entries of the contingency tables. Considering together exposed and unexposed cases as in (9), and ignoring the denominator of the likelihood (since it does not depend on the parameters), we get that the log-likelihood of the model is proportional to

[TABLE]

When $J=1$ (there is only one list) we do not have identifiability. In fact for the two sets of parameters

[TABLE]

the log-likelihood is $(M_{1}^{E}+M_{1}^{U})\cdot\log\left(\frac{e}{2}\right)-e$, while for

[TABLE]

the log-likelihood is $M_{1}^{E}+M_{1}^{U}-2e$. This proves that we can get the same likelihood with distinct sets of parameters, hence the model is not identifiable.

When the number of lists is $J=2$, we show that we already have identifiability. In this case, in fact, we only add the parameter $\alpha_{2}$ (since the higher order interaction is always excluded). We see from (14) that the conditions for having the same distribution with the two set of parameters is that

[TABLE]

Here we have that, for the exposed contingency table, the probabilities are

[TABLE]

[TABLE]

and same for the unexposed, without the parameter $\theta$. If we combine together the equation involving $p_{11}^{U}$ and the one with $p_{10}^{U}$, we get that $\alpha_{2}=\alpha_{2}^{\prime}$. Similarly using the equation with $p_{11}^{U}$ and the one with $p_{01}^{U}$ we have $\alpha_{1}=\alpha_{1}^{\prime}$. It is simple now to use the equations for the exposed cases and show that $\theta=\theta^{\prime}$ and finally that $\gamma_{E}=\gamma_{E}^{\prime}$ and $\gamma_{U}=\gamma_{U}^{\prime}$. This proves that the model is identifiable.

When the number of lists is larger, we have $2^{J}-1$ expressions $\gamma_{U}\cdot p_{i_{1},...,i_{J}}^{U}=\gamma_{U}^{\prime}\cdot{p^{\prime}}_{i_{1},...,i_{J}}^{U}$, from which we get identifiability of the $2^{J}-2$ lists parameters ($J$ for the individual lists, $\binom{J}{2}$ for the two-list interaction and so on). From three of the remaining $2^{J}-1$ expressions $\gamma_{E}\cdot p_{i_{1},...,i_{J}}^{E}=\gamma_{E}^{\prime}\cdot{p^{\prime}}_{i_{1},...,i_{J}}^{E}$ we get identifiability of $\theta,\gamma_{E}$ and $\gamma_{U}$.

Remark D.1**.**

The same result on identifiability holds if we allow different $\theta_{j}$ coefficients for each list. In the situation with only $2$ lists we now have $6$ observations and $6$ parameters, and similarly as before it is possible to get identifiability of $(\alpha_{1},\alpha_{2},\theta_{1},\theta_{2},\gamma_{E},\gamma_{U})$. When we extend this result to $J>2$ lists we can use the $2^{J}-1$ expressions $\gamma_{U}\cdot p_{i_{1},...,i_{J}}^{U}=\gamma_{U}^{\prime}\cdot{p^{\prime}}_{i_{1},...,i_{J}}^{U}$ to we get identifiability of the $\alpha$ parameters and the remaining $2^{J}-1$ expressions $\gamma_{E}\cdot p_{i_{1},...,i_{J}}^{E}=\gamma_{E}^{\prime}\cdot{p^{\prime}}_{i_{1},...,i_{J}}^{E}$ to solve for $\theta_{1},...,\theta_{J},\gamma_{E}$ and $\gamma_{U}$. $\square$