Detecting contaminated birthdates using generalized additive models

Distribution of birthdates in a database

For the purpose of this paper, we consider only databases that identify one or more groups of individuals—such as patients or clinicians. Many clinical or administrative databases fall into this category including electronic health records (EHR), emergency department information systems, databases of controlled drug prescriptions, and medical claims databases. With a clear identification of individuals in such a database, each date d defines a set of N(d) individuals born on that day.

In an ideal situation where the actual numbers of births and deaths for every day are available, for example in a region with a well-maintained birth and death registry and limited population migration, the actual number of births in the population for every date can be estimated. But with any given database, it is useful to distinguish 1) the general population P₁ that includes everyone living in a region at a certain time period, 2) the “at risk” population P₂ that includes everyone who in theory could be included in the database (e.g., males to a prostate cancer database), and 3) the group of people P₃ who are actually in the database. Individuals in a database P₃ can be regarded as a sample of the “at risk” population P₂. Consequently the number N(d) can be regarded as realization of a birthdate distribution defined on P₂.

Here k is the number of individuals born on day d and λ(d) = ∫_t ∈ dμ(t)dt is the aggregated intensity for date d and equals the expected value of N(d).

In certain situations with larger data variance, a negative binomial distribution can be used instead. That requires one more over-dispersion parameter to fit. In this application, it is difficult to assess over-dispersion; hence we prefer the simpler Poisson model assumption.

Let μ _d be the mean of μ(t) on day d and l be the length of a day. Then λ(d) = μ _d ⋅ l. If μ(t) changes slowly, sequence < λ(d) > can be regarded approximately as the result of sampling μ(t) daily and then multiplying the sample with the constant l. The sequence < λ(d) > is determined by the birthdate distribution of the general population P₁ and the representation of P₁ in the database.

The birthdate distribution of the general population P₁ is determined on a large scale by the age profile of the population and on a small scale by seasonal and weekly fluctuations in child births. The age profile of the population (as in [10]) is the aggregated result of changes in birth and death rates. For example, the post-WWII baby boom has contributed to the aging population structure in US, Canada, and Australia [11, 12]. As such changes are often driven by long lasting demographic forces such as economic development, war, and progress of medical science [13, 14], they manifest as slow and smooth variation in the sequence < λ(d)>. In contrast, seasonal and weekly variations in child births (as in [15]) act on a shorter temporal scale. More recently, most child births occur in a hospital environment and scheduled Caesarean section or labour induction are more frequent. These factors lead to more births on weekdays than weekends, which produces cyclic dips in the sequence < λ(d) > [16, 17].

How the general population is represented in a patient database is mostly determined by the nature of the disease. For example, a database of prostate cancer patients should contain only males; a database of skin cancer patients probably contains more Caucasians than patients of other races. In terms of the birthdate distribution, one common consideration is that a disease may pose higher risk to certain age groups. For example, individuals born before year 1930 may be disproportionately represented in a database of chronic non-cancer pain. Like the population age profile, the nature of the disease acts on a larger temporal scale and it should not affect the smoothness of the sequence < λ(d)>. For example, a person born on March 2nd, 1980 and a person born on March 3rd, 1980 should have very similar chances of being included in a database, assuming all other conditions are equal.

Finally for most hospitals, a patient database grows out of paper records. As a hospital serves more patients, its database may cover a larger portion of the at-risk population P₂. However, it is reasonable to assume that any such change in the size of the database will affect different age groups proportionally, and is independent of the shape of < λ(d) >.

To summarize, the birthdate distribution of a database can be modelled using a function λ(d) that is generally smooth but contains weekly dips due to weekend reduction in births.

Tell-tale indicators of birthdate contamination: discrepancy between the expected and observed counts

When birthdates in a database are systematically contaminated, for example when a missing value is consistently replaced with the zero date-time, incorrect birthdates may be repeatedly generated. In such cases, the frequency of an incorrect birthdate d will be higher than other dates in the database. That is, the observed number of patients with birthdate d, denoted N(d), will be larger than the expected number λ(d). Therefore identifying incorrect birthdates can be achieved through identifying the date d whose person count N(d) is well above the expected count λ(d). Of course this observation is not new. For example, the difference between N(d) and λ(d) is called within deviation by Dasu and Johnson [18]. To the authors’ knowledge, however, no previous attempt has been made to statistically model this discrepancy between N(d) and λ(d) for error detection.

In many applications, only a small number of the most likely erroneous birthdates will be manually checked. If required, the expected number of birthdate errors in database can be estimated from the results of manually checking small samples of the database [19].

A generalized additive model for birthdate distribution

The expected birthdate distribution λ(d) can be estimated from the data < N(d) > by assuming smoothness of the function λ(d), with proper handling of the reduction in births on weekends. In the previous section, we have seen that explicit modelling of the mechanisms shaping the age profile would be very complex. Nevertheless, both long-term and seasonal variations in births/deaths can be recovered by smoothing of counts N(d) in a relatively short temporal window. The reduction in weekend births can be modelled by a multiplicative factor that only affects Saturdays, Sundays, and public holidays. As the reduction is a gradual development driven by increasing hospital births, elective Caesarean section, and labour induction, the multiplicative factor can also be modelled by a smooth function. Finally, it is worth noting that government policies to encourage fertility can also lead to more births on a particular day [20], but such events are very rare and should be explicitly modelled on a case-by-case basis.

Here the log link function is used because λ(d) can only be positive; I(⋅) is the indicator function and s₁(d) and s₂(d) are smooth functions of d. The function s₁(d) models both long-term and seasonal variations in the birthdate distribution; s₂(d) models the gradual change in the reduced weekend births. Public holidays are similar to weekends; it only requires extra bookkeeping in the modelling process.

Here d is the dimension of the function domain (in our case d = 1), m is the order the smoothness penalty, ϕ _j are $\mathit{M}=\left(\begin{array}{c}\hfill \mathit{m}+\mathit{d}+1\hfill \\ \hfill \mathit{d}\hfill \end{array}\right)$ basis functions that spans the null space of a penalty energy function, and η _md is a radial-basis function. δ = (δ₁, δ₂, …, δ _n ) ^T and α = (α₁, α₂, …, α _M ) ^T are parameters with the constraint $\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{\delta }}_{\mathit{i}}{\mathit{\varphi }}_{\mathit{j}}\left({\mathit{x}}_{\mathit{i}}\right)=0$ for each j. Further details of thin-plate splines can be found in [23]. Just like other kernel methods, the computation of thin-plate splines has a complexity of O(n³). For ease of computation, a lower rank approximation of δ is used. Such reduced-rank thin-plate splines [23] are used in our application. The rank of thin-plate splines controls the smoothness of the function, and is selected by generalized cross-validation criterion (GCV) [24]. To capture seasonal variations in birthdate distribution, the maximum rank for the candidate splines should be set to a sufficiently large number. We recommend a maximum rank of 2 m, where m is the number of years covered by the birthdates.

Dasu and Johnson [18] generated a histogram similar to Figure 1 to help detect inadvertent censoring caused by a default date. Graphical features of the histogram such as spikes and V-shaped valleys were identified as indicators for missed or censored data. However, the presence and location of such graphical features are determined by visual inspection. In contrast, here with a GAM, the deviation of particular counts in the histogram can be quantified for probability based decisions.

A public health application

We use a health administration database to demonstrate the fitting of GAM and the identification of outliers. A drug regulatory authority maintains a database of drug dispensing records, which covers around 40,000 individuals with birthdates ranging from the 1890s to 2010s. The data is integrated from separate data collections from more than 3000 off-site pharmacies. As each pharmacy may choose different software to manage their data, quality assurance at the data collection stage is limited. For example, one system replaces the null date-time value with a default date, which changes whenever the system is updated. This has resulted in a typical downstream database with multiple sources of contamination. Our task was to identify those birthdates that may have been incorrectly introduced by the faulty data management software. As the data were transcribed from doctor’s prescriptions filled at off-site (community) pharmacies, they were geographically and functionally independent of the data custodians. These handwritten paper scripts were archived after data entry. Geographically locating, accessing and then physically combing through these handwritten documents is theoretically possible (with generous resources), but was impractical for the purposes of this study.

Figure 1 shows the number N(d) of individuals for each birthdate d. For evaluation purposes, most incorrect birthdates have already been manually identified by a domain expert. There were 58 errors in total. The authors have no knowledge of the criteria used by the domain expert to identify the errors and how complete and accurate they are. Nevertheless, Figure 1 shows at least three potential sources of database errors. First on the far left, errors 1899-12-30 and 1900-01-01 were likely to have been introduced by incorrect handling of the zero date-time in software systems. The error 1911-01-01 may also have resulted from the confusion of year 1911 and year 2011. Second in the middle, errors 1950-01-01, 1960-01-01, 1963-02-14, 1970-01-01, 1979-01-01, 1980-01-01 may result from incorrect self-reported birthdates. Finally on the right of the figure, a group of incorrect birth dates were introduced by the pharmacy management system that mixed birthdates with drug dispensing dates.

Outliers detected through generalized additive model

We fit an additive model in the form of Equation (2) to the birthdate counts in Figure 1. The gam function from the mgcv R-package [23] was used to estimate the intensity function λ(d). The smoothness was determined through GCV with an upper bound of 200 for the degree of freedom. Figure 2 shows the estimated intensity function and the outliers above the point-wise 99.99th percentile.

The detection resulted in 8 false positives and 7 false negatives. The false positives can be classified into two categories: 5 in the middle between years 1957 to 1980; 3 on the far right between years 2006 and 2007. From visual inspection, the 5 false positives in the middle are indeed outliers. We conjecture that these birthdates were caused by multiple identities of some patients. When data from multiple pharmacies were linked, patient identities were consolidated through names and addresses. Mismatch can happen during the consolidation process and the patient’s birthdate would be overrepresented. The 3 false positives on the right are also curious. They seem to suggest that there are dozens of very young (~6 years old) patients on controlled drugs. This warrants further investigation.

The 7 false negatives include 2 in the middle (1979-01-01 and 1980-01-01), 4 on the right (all in year 2001), and 1 on the right (2011-01-01). The false negative 2011-01-01 has too few patients to justify a statistical decision. The 4 false negatives in year 2001 coincide with a large number of incorrect birthdates that have been correctly identified. We believe that after the correctly identified errors have been cleaned, these 4 false negatives will likely be found by rerunning the identification procedure.

The estimated Poisson intensity function λ(d) (the blue line in Figure 2) reflects both the general age profile of the patient population and short-term variations in demographics. Figure 3 shows the segment of the curve for weekdays between year 1925 and 1970. The curve suggests a seasonal pattern of child birth. The dip around 1933–1934 suggests the effect of the ‘Great Depression’ on child birth in Australia. The elevated counts at the centre echo the baby boom in Australia following the second world war, which Salt defined to be the period between 1946 and 1961 [12]. The latter two features of the estimated function λ(d) are consistent with the official statistics [26].

The effect of reduced weekend birth (term s₂(d) in Equation (2)) is shown in Figure 4. It is consistent with the trend that weekend births have been significantly reduced due to elective Caesarean section. It also suggests that elective Caesarean section has been gaining popularity in the past 50 years.

Comparison with outlier detection based on a standard time-series model

As the sequence < N(d) > is also a time series, an autoregressive integrated moving average (ARIMA) model is a natural alternative to a generalized additive model. An ARIMA(p, n, q) model contains p autoregressive terms, n nonseasonal differences, and q moving average terms. It has the following mathematical expression: $\left(1-{\sum }_{\mathit{i}=1}^{\mathit{p}}{\mathit{\varphi }}_{\mathit{i}}{\mathit{L}}^{\mathit{i}}\right){\left(1-\mathit{L}\right)}^{\mathit{n}}\mathit{N}\left(\mathit{d}\right)=\left(1+{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathit{\theta }}_{\mathit{i}}{\mathit{L}}^{\mathit{i}}\right){\mathit{\epsilon }}_{\mathit{d}}$, where L is the lag operator, ϕ _i and θ ⁱ are autoregressive and moving average parameters, respectively. When a time series is season with m points per season, a seasonal ARIMA (SARIMA) model can be used. In a SARIMA(p, n, q)(P, N, Q) _m model, additional P autoregressive terms, N differences, and Q moving average terms are used to model the seasonality of the time series. More details of ARIMA models can be found in most time series text books (e.g., [27]).

Although ARIMA models may be more familiar to many people, they are not ideal for modelling the birthdate distribution. First, an ARIMA model assumes normality of data, which is not always appropriate as N(d) is a count. In particular, outlier detection based on normal percentiles may not be as accurate compared with Poisson percentiles. Next, it is not easy to incorporate other covariates in an ARIMA model and assess their effects on the counts. Most importantly, from a linear time series model like ARIMA, it is impossible to infer point-wise tail probabilities (see the green line in Figure 2). If an over-all tail probability estimate is used, incorrect birthdates at the two ends of the age range will be missed (see the third advantage of the GAM described previously).

We fit a seasonal ARIMA (SARIMA) model on the sequence < N(d) > adjusting for the weekly cycle. It results in a high-order model ARIMA(0, 1, 5)x(0, 0, 2) ₇ . By looking above the 99.99^th percentile in the residuals, we identified 69 outlier birthdates (See Figure 5). We compared these birthdates with the birthdates identified by the domain experts. A total of 36 birthdate errors were correctly matched, with 11 false negatives and 22 false positives, mostly in the period between year 1940 and 1980, in which many patients were born. Because no point-wise percentile can be inferred, this method tends to miss the birthdate errors near the left end of the birthdate range, where the data is sparse. For example, the date 1900-01-01 was identified by the GAM model (see Figure 2), but was missed by the ARIMA model. A comparison of sensitivity and specificity of the two methods are shown in Table 1.

The GAM-based method showed better performance in terms of sensitivity and specificity. The area under the curve (AUC) for SARIMA-based outlier detection was 0.991; AUC for GAM-based outlier detection was 0.998 (p-value 0.28).

Age profile of young oxycodone user

Finally, we show how identification of incorrect birthdates can prevent serious misinterpretation of data.

Oxycodone is an opioid analgesic used for pain management. A patient using oxycodone for an extended period is likely to develop dependence to the medication. Therefore in many places, including the state of Queensland in Australia, an oxycodone treatment episode longer than two months requires a report to the regulatory authority. Here we use the drug dispensing database to understand the age distribution of long-term users of oxycodone.

The age distribution for younger patients is shown in Figure 6^a. The figure shows unusual clusters of patients of age 5 or 10. These clusters are worrying as the main reason for Oxycodone prescription is chronic non-cancer pain, which is rare among children.

We applied the GAM-based identification method to find birthdate errors. The 99.99th percentile was used as the cut-off. The records with identified birthdate errors were then removed from the analysis. The corrected age distribution based on the cleaned subset is highlighted in red. Although the new distribution still shows several young patients whose presence in the database is worth further investigation, it is more consistent with the common understanding that young people rarely have chronic pain condition that warrants large quantities of opioids.