Exploring Differences between Two Decades of Mental Health Related Emergency Department Visits by Youth via Recurrent Events Analyses

Two collections of pediatric MHED (Mental Health Emergency Department) visit records, one from 2002–2010 and the other from 2010–2017, are available to us by two different data extractions from the population-based administrative data sources: the National Ambulatory Care Reporting System (NACRS) and the Alberta Health Care Insurance Plan (AHCIP). The MHED data include both basic demographic characteristics (e.g., Sex and Region) and the visit characteristics (e.g., Date and Diagnosis of the ED visit, and Age in integer years at the ED visit) on all the MHED visits made by residents under 18 years old of the province of Alberta, Canada. Viewing an ED visit as the event of interest, we formulate the MHED records into recurrent event data. Analyses of the recurrent event data are conducted to compare mental health related emergency department visits by children and youths of the two decades with respect to the frequency in age and how it is associated with risk factors/exposures such as Sex and geographic Region. In the rest of this paper, the two datasets are referred to as MHED data of “early” and “late” groups, short for the decades starting in 2000 and 2010, respectively.

Administrative health data in general contain rich information for investigating public health issues such as infectious diseases, disease surveillance, and health service use, cross-sectionally and longitudinally. These data can be particularly useful if the data are population-based (i.e., including all residents in a defined geographic region), collected over time, at individual-level, and accessible at low cost (Burgun et al., 2017). The data may include records of physician visits, laboratory investigations, health service use, diagnoses, and treatments (Gavrielov-Yusim and Friger, 2014). Disciplines like pharmacoepidemiology (Tamariz et al., 2012), epidemiologic surveillance (Ferrante et al., 2013) and health services research (Penfold et al., 2018) frequently use administrative health data. Administrative health data can be more cost-effective than resource-demanding studies to investigate, say, trends in infectious diseases (Gavrielov-Yusim and Friger, 2014). Recently, administrative health data have been used to study trends in the coronavirus disease 2019 (COVID-19) pandemic (Navaratnam et al., 2021), uptake of vaccinations (Shariff et al., 2022), and longer term outcomes from COVID-19 (Chevinsky et al., 2021).

Because of the special features of Canada’s health care system, administrative health databases are created to administer Canadian provincial and territorial health insurance plans; however, many restrictions and regulations apply to their access and use. In addition, these administrative health databases are developed and maintained to serve non-research purposes. Those features yield challenges to analyses of administrative health data in general. This paper presents strategies for dealing with the particular challenges confronted in our analyses of the MHED data. The ideas and methods are in fact applicable to many other administrative data analyses.

Each data extraction of the two MHED datasets was administered with calendar time. A scientifically meaningful time origin is usually subject-specific. We focus on the use of age as the individual time in the paper. Also, eligible subjects were limited to those aged less than 18 at the time of their ED visit and the data extraction provided only the age in integer years. Those aspects result in the first major challenge that the information availability varies from individual to individual. Specifically, the two MHED datasets include the subjects’ information over individual-specific intervals. That is, the MHED records of a subject are subject to double-censoring; the left-censoring and right-censoring times vary from subject to subject. The structure of doubly-censored recurrent event data (e.g. Zhang and Li, 1996; Cai and Sun, 2003) indicates a particular demand of loosely-structured models to be employed in data analyses such as the marginal models of recurrent events (e.g. Lawless and Nadeau, 1995; Lin et al., 2000). Since the range of age differences in the MHED records within individual subjects and among them is rather wide, we consider the extended Cox regression model with age-varying coefficients for the counting process of an individual’s cumulative MHED visits over age. Using the indicator of the late group (i.e. the 2010 group) as one of the covariates, the model permits exploring age-dependent differences between the two decades and other risk/exposure (covariate) effects on the MHED visit frequency. Analyses under the model include all the corresponding stratified analyses with the data from the two decades. We follow Hu and Rosychuk (2016) to adapt the local partial likelihood procedure of Cai and Sun (2003) and Tian et al. (2005) for estimating the model parameters with the recurrent event data.

The second major challenge confronted when analyzing the MHED data is the inherent truncation issue of administrative health data. Administrative health data usually provide information only on people who seek health services due to a particular disease or condition, which is a subgroup of the general population. To provide administrators and policy-makers with an overall view of the disease burden and trends of using health services, interests are usually in the population as a whole on its measures such as rates of using health services. Examinations of only health services users for a particular disease or condition during a limited time period form a truncated sample of the whole population and examinations of their data only lead to biased conclusions in general. Most well-established approaches in analysis of truncated data focus on the situations where observations on the time to an event are subject to length-biased sampling; see, for example, Wang (1996), Shen et al. (2009), and Qin et al. (2011). There is relatively less work concerning zero-truncated recurrent event data.

In our application, individuals in Alberta who are younger than 18 years old in the two decades, no matter if they are included in the MHED cohort or not, may be framed as a random sample of the general population. By the features of Canadian medical insurance system with population-based administrative databases, eligible subjects who are not included in the pediatric MHED cohort can be safely assumed not to have MHED visits during the two time periods. Augmented by that fact, each of the two MHED datasets can be viewed as a collection of doubly-censored recurrent event records where the covariates of all the subjects not in the MHED cohort are missing. There is considerably rich literature available on statistical analysis of right-censored recurrent event data (Cook and Lawless, 2007). Most well-established statistical tools are not directly applicable, however, to analysis of the MHED data when targeting the general population since covariates of the individuals not in the cohort are unavailable. Hu and Lawless (1996) considered a similar situation where the awareness of a study unit and the associated observation period is unavailable until the study unit experiences at least one event of interest. They proposed to handle the truncated data using information on the size of the total sample including study units not having events over the observation period and considered robust nonparametric estimators of the rate and mean function for the process of recurrent events. We adapt this idea to tackle the truncation issue. Particularly, we aim to overcome the challenge of unknown covariates, which is posed by the zero-truncated data.

Some approaches have recently been proposed to utilize aggregated data or summary statistics such as disease incidence rate and survival probability to correct bias when analyzing a smaller-scale sample from the target population; see, for example, Zheng et al. (2021), Gao and Chan (2021), and Li et al. (2022). Those publications consider the situations where either the covariate effects are constant over time or effects are assumed to be the same within the source cohort and the target population. We are interested in exploring time/age-varying effects of the covariates on the aforementioned MHED visits. We analyze the MHED data, the recurrent event data, in a combination of publicly available population-based aggregated data to provide insights into the MHED visit patterns of the general population in addition to targeting on the pediatric MHED cohort, the subpopulation whose MHED visit records are extracted and available to us.

More practical challenges were encountered when the MHED data were analyzed. Privacy laws prevented patient birthdates to be extracted from the AHCIP during the early time, and thus the 2002–2010 data do not include the birthdates of the subjects. This data feature partly prompted Hu and Rosychuk (2016) to develop a marginal regression model to handle recurrent events with coarsened censoring times, which resulted because the start date for an individual (i.e., birthdate) was not known exactly but known to be contained in an interval (i.e., based on the age at ED visit). We adapt their approach in the analyses of this paper. Further, the two MHED datasets were extracted at different times, and the data extractions did not allow unique personal identifiers. Thus, the subjects with MHED visits in both the early and late datasets cannot be identified and linked. Despite being aware of the fact, we start with assuming that all of the individuals in the MHED data are independent with each other in the analyses. We then conducted a simulation study to examine the consequences if the assumption is violated.

The contributions of this paper are twofold. We present innovative strategies for processing doubly-censored recurrent event data, and overcoming the truncation induced by the data collection. In addition, through exploring the pediatric MHED visit records, we provide new insights into how the mental health of children and youth changes over time. We organize the rest of this paper as follows. The notation and modeling are introduced together with discussions on useful model specifications in Section 2. Section 3 presents analyses of the aforementioned MHED data from 2002–2010 and 2010–2017 concerning the pediatric MHED cohort. Section 4 shows strategies and analyses of the MHED data supplemented by aggregated population information concerning the general population. In Section 5, we report simulation studies conducted to numerically investigate two potential issues in the data analyses of Sections 3 and 4. Some final remarks are given in Section 6.

This section describes the formulation of the MHED data and the modeling considered in the data analyses.

Denote by $\mathcal{O}$ the whole sample of children and youths who were under 18 years old and residents of Alberta, Canada during 2002–2017; $\mathcal{O}_{1}$, the sub-sample that only includes the subjects who made MHED visits during the time, i.e., the MHED cohort. That yields $\mathcal{O}=\mathcal{O}_{1}\bigcup\mathcal{O}_{0}$, where $\mathcal{O}_{0}$ is the sub-sample of the eligible subjects who had no MHED visit during the time. Moreover, denote the early and late groups (i.e., the 2000 decade and the 2010 decade groups) by $\mathcal{O}_{E}$ and $\mathcal{O}_{L}$ and, consequently, the MHED cohorts are $\mathcal{O}_{E1}$ and $\mathcal{O}_{L1}$. Thus $\mathcal{O}=\mathcal{O}_{E}\bigcup\mathcal{O}_{L}$ and $\mathcal{O}_{1}=\mathcal{O}_{E1}\bigcup\mathcal{O}_{L1}$.

Let $N(a)$ represent the number of MHED visits made by a subject up to age $a$, and $\mathcal{H}(a)=\{N(t):0\leq t<a\}$ be the history information prior to age $a$. We use $Z$ to denote the covariates of interest, which in this paper are the two categorical variables Sex and Region. Further, let $X$ be the indicator of whether or not a subject belongs to the late subgroup $\mathcal{O}_{L}$, and $[W^{E}_{L},W^{E}_{R}]$ and $[W^{L}_{L},W^{L}_{R}]$ be the data extraction windows (in calendar time) of the early and late MHED data, respectively. The left and right censoring times (in age) are then, respectively, $C^{E}_{L}=\max(0,W^{E}_{L}-B)$ and $C^{E}_{R}=\min(18,W^{E}_{R}-B)$ for a subject from the early subgroup $\mathcal{O}_{E}$ with birthdate $B$; $C^{L}_{L}=\max(0,W^{L}_{L}-B)$ and $C^{L}_{R}=\min(18,W^{L}_{R}-B)$, for $\mathcal{O}_{L}$. The above variable names with subscript $i$ are the corresponding values associated with subject $i$.

Given that the extraction windows are known, the MHED data can be expressed as follows:

where $\mbox{\sc Data}_{E1,i}=\Big{\{}N_{i}(a)-N_{i}(C_{Li}):~{}a\in(C^{E}_{Li},C^{E}_{Ri}]\Big{\}}\bigcup\big{\{}Z_{i}\big{\}}$ and $\mbox{\sc Data}_{L1,i}=\Big{\{}N_{i}(a)-N_{i}(C_{Li}):~{}a\in(C^{L}_{Li},C^{L}_{Ri}]\Big{\}}\bigcup\big{\{}B_{i},Z_{i}\big{\}}$ are the available data of subject $i$ from the early cohort $\mathcal{O}_{E1}$ and the late cohort $\mathcal{O}_{L1}$, respectively. Note that $C^{E}_{Li}$ and $C^{E}_{Ri}$ are unavailable for $i\in\mathcal{O}_{E1}$ due to the missing birthdate $B_{i}$ in our application.

In this paper, we consider statistical inference on the target population, denoted by $\mathcal{P}^{*}=\mathcal{P}^{*}_{E}\bigcup\mathcal{P}^{*}_{L}$ in the following two settings: the MHED cohort of the two decades $\mathcal{O}_{1}=\mathcal{O}_{E1}\bigcup\mathcal{O}_{L1}$ is a random sample of $\mathcal{P}^{*}=\mathcal{P}_{1}$, or the whole group of interest $\mathcal{O}=\mathcal{O}_{E}\bigcup\mathcal{O}_{L}$ is a random sample of $\mathcal{P}^{*}=\mathcal{P}=\mathcal{P}_{1}\cup\mathcal{P}_{0}$. In the rest of this paper, $\mathcal{P}^{*}$ together with its subpopulations $\mathcal{P}_{E}^{*}$ and $\mathcal{P}_{L}^{*}$ is specified in the context.

We assume that, for age $a\in(0,18)$ and $i\in\mathcal{P}^{*}$,

where $\Lambda_{0}(a)=\int_{0}^{a}\lambda_{0}(u)du$ with $\lambda_{0}(\cdot)$ the baseline rate function. This model is an extension of the proportional rates/means model discussed by, for example, Pepe and Cai (1993), Lawless and Nadeau (1995), and Lin et al. (2000) in their marginal analysis of recurrent events.

The model includes many useful specifications. For example, it is the Poisson process model with the conditional intensity function in the same form as (2) if assuming the counting process $\big{\{}N(a):a\geq 0\big{\}}$ is Poisson. The Poisson process model reduces to the Andersen-Gill model when all the regression coefficients are time-independent (Andersen and Gill, 1982). Table 1 displays relevant semiparametric marginal models derived from (2).

We note that the model is equivalent to the following stratified model. With $\beta_{E}(a)=\beta(a)$, $d\Lambda_{0E}(a)=d\Lambda_{0}(a)$, $\beta_{L}(a)=\beta(a)+\gamma(a)$, and $d\Lambda_{0L}(a)=\exp\big{\{}\alpha(a)\big{\}}d\Lambda_{0}(a)$,

This representation provides a meaningful interpretation of the regression coefficients of (2) in our application. More discussion is given in the following sections. When interests are in knowing about each of the subpopulations, one can conduct a stratified analysis without specifying how the two subpopulations are correlated.

This section focuses on the MHED cohort. The MHED data from the two decades are analyzed with the target population $\mathcal{P}^{*}=\mathcal{P}_{1}$ including $\mathcal{O}_{1}=\mathcal{O}_{E1}\bigcup\mathcal{O}_{L1}$ as a random sample, where $\mathcal{O}_{1}$ is defined in Section 2 as all the subjects in the MHED cohort from the early and late decades. We begin with a descriptive analysis of the MHED data. A comparison of the MHED visits between the two decades is conducted through analyses of the MHED data under the model in (2) and its specifications given in Table 1.

The aforementioned data provide information only on the individuals in the pediatric MHED cohort, who were Alberta residents, under 18 years of age and experienced at least one MHED visit during 2002–2017. Therefore, the analysis reported in Section 3 cannot reveal the frequency and associations with the factors/exposures of the ED visits concerning the whole province. We thus conducted alternative analyses of the MHED data following the approach proposed by Hu and Lawless (1996) to adjust for the truncation. The key idea is to make use of the fact that eligible subjects who are not included in the MHED cohort must not have had any MHED visit during the two time periods. This assumption is plausible due to the single-payer government medical insurance system in Canada and that population-based administrative databases are used in the administration of health services. The MHED data augmented by that information can then be framed into a collection of doubly-censored recurrent event data where all of the covariate data of the subjects not included in the MHED cohort are missing. The following presents an analysis of the MHED data integrated with publicly available population census information.

We conducted simulation studies to examine the finite-sample performance of our approach and to verify the findings of the MHED data analyses. Two settings were considered in the simulation: Setting 1 mimicked the real data example, generating recurrent events associated with two eras; Setting 2 simulated a situation where subjects were divided into two generations based on their birthdates and our approach is applied to analyze the generated data. Setting 2 allows us to explore the performance of our approach with a misspecified model. In each setting, we conducted analyses of the generated data concerning with the target population $\mathcal{P^{*}}$ being $\mathcal{P}_{1}$, the population of subjects with at least one event during a given period, and $\mathcal{P}$, the whole population of all individuals no matter if included in $\mathcal{P}_{1}$. When $\mathcal{P^{*}}=\mathcal{P}_{1}$, we checked how analysis outcomes are affected when the independence assumption for the two groups does not hold. When the target population is $\mathcal{P^{*}}=\mathcal{P}=\mathcal{P}_{1}\bigcup\mathcal{P}_{0}$, we examined how well the approximation in the approach performs using the population-level census information.