Characterising patterns in routinely reported longitudinal HIV data in South Africa using a Bayesian multiplicative interaction model
•,,,,,,,.
...
Abstract
Introduction We consider an analytical problem of characterising patterns and identifying discrepancies between database systems for longitudinal aggregated healthcare data involving multiple facilities.
Methods We used routinely collected data on the registered number of people living with HIV who initiated antiretroviral treatment (ART) in 69 South African facilities in 2019; reported in the Three Interlinked Electronic register (Tier.net) and the District Health Information System. A Bayesian multiplicative interaction model quantified the average time effect as realised through the heterogeneous facility-specific slopes and quantified discrepancies between the two database sources.
Results The estimated average trends showed a slight dip in June and a large dip in December. The estimated slopes identified clusters of facilities based on their ranges of fluctuations over time. The differences in average monthly ART initiations between the two database sources had a median of 1.6 (IQR 0.8–3.3), while 3 outlying facilities differed by at least 10 ART initiations between the 2 sources.
Conclusion Multiplicative interaction models are a powerful tool for quantifying average trends over time and for evaluating discrepancies between reporting systems for multiple facilities with heterogeneous time slopes. The Bayesian framework enables efficient estimation for a very large number of parameters.
What is already known on this topic
Large databases that capture longitudinal outcomes from routine clinical care in multiple facilities are a valuable resource for researchers and public health decision-makers. We consider an evaluation of HIV programme data on the number of people initiating antiretroviral treatment (ART) in multiple facilities in South Africa. These data are often analysed using traditional longitudinal models that assume a common trend for all reporting healthcare facilities.
What this study adds
We illustrate a Bayesian analytical method which offers an objective approach that accounts for the heterogenous facility by time interactions that contribute to an overall trend in these routinely reported longitudinal HIV data. This method is also extended to quantify discrepancies when data from the same facility are reported in different database sources. The model estimated average trends of the number of ART initiations, with dips in June and December, and it estimated heterogeneous slopes that help distinguish clusters of facilities based on their ranges of fluctuations over time. The differences in average monthly ART initiations between the two database sources had a median of 1.6 (IQR 0.8–3.3), while 3 outlying facilities differed by at least 10 ART initiations between the 2 sources.
How this study might affect research, practice or policy
This is a valuable analytical approach that should be considered as a standard for quantifying trends in longitudinal healthcare data from multiple facilities and for detecting reporting discrepancies.
Introduction
We consider an analytical problem of evaluating aggregated longitudinal healthcare data reported within the same data management system or healthcare facility. We also consider the task of quantifying differences in such data when they are reported in multiple database systems. Examples of these data include cascades of numbers of people screened for, diagnosed with and started on treatment for a particular disease, recorded over multiple timepoints. Accurate data are crucial for properly monitoring development, research, planning and decision-making in healthcare.1–3
Multiple databases that store the same or related healthcare indicators are often created by various stakeholders (eg, government vs non-governmental organisations) to serve specific purposes such as research, resource allocation and policy-making. There are multiple reasons for the discrepancies between data sources or for the observed fluctuations within longitudinal records. Fluctuations may reflect seasonal effects, changes in disease patterns in the catchment populations changes in treatment guidelines (eg, initiating people on treatment at different timepoints from time of diagnosis) or external shocks to the health system (eg, the COVID-19 pandemic3). Whereas, discrepancies may be due to human error when reporting or entering the data, rounding and poor reporting of missing data.4
Approaches for evaluating these types of data may involve simple arithmetic operations such as absolute and proportional differences between subsequent data points to study trends over time, or calculation of differences between data source averages. These can be informative, but they do not offer a summary of the overall time-related fluctuations while taking into account facility differences. Furthermore, when conducting such assessments for multiple facilities with repeated data points over time, the exercise can be overwhelming or even infeasible. The analytical method we propose here offers a summary of the overall trends while accounting for variations in facility records, and it is an objective approach that assumes that the true time effects are observed only through the contributions made by all facilities’ data.
The goal of this study is to use Bayesian multiplicative regression models to describe time effects in longitudinal healthcare data; determine differences in patterns between facilities; and identify discrepancies when data are reported in different database systems. The model assumes that the true time effects (eg, disease burden changes and seasonal effects) are unknown and are derived from the reported data. The model also assumes that neither of the data sources is superior to the other.
Patient and public involvement statement
Individual patients were not involved in the recruitment to and conduct of the study. We used facility-level aggregated data.
Materials and methods
Motivating example
The motivation for this analysis was a review of monthly clinic data on people living with HIV (PLHIV) abstracted from two data management systems in South Africa. These were collected as preliminary data to select sites for an ongoing cluster randomised trial, the FEDISA (“finish”) tuberculosis (TB) study, to increase TB preventative therapy (TPT) prescribing. The first data source is the South African Three Interlinked Electronic register (known as Tier.net)5 6 which is a national electronic medical register and electronic medical record that tracks laboratory results, treatments and intermittent illness and pregnancies of patients living with HIV, TB as well as maternal and child health indicators. The system is three tiered in that it is made up of the paper-based medical registers (tier 1), an electronic version of the paper-based registers (tier 2) and the full electronic medical record (tier 3) which contain information from tiers 1 and 2.4 The lower tiers are used in low-resourced areas and the goal is that, ultimately, all facilities will transition to tier 3 as resources and skilled personnel improve.
The second data source is the District Health Information System (DHIS) now called DHIS V.2.7 8 This is an open-source information management system that was developed at the University of Oslo and is used in 73 countries. In South Africa, it was first used in 1996 and it has now been used by all health districts nationally since 2001. In this system, aggregate level data are reported by the facilities to monitor progress in healthcare delivery, as well as target resource requirements. Some of the problems associated with data quality on DHIS V.2 have been poor reporting frequency and inadequate training for staff members responsible for reporting.9
Evaluations of the accuracy of data from these types of information systems can be done at the individual patient level by checking granular information such as treatment dates, treatment outcomes and lab test information.4 9–11 The current analysis focuses on evaluating the data at the aggregate facility level. We illustrate a method for using facility-level data to study time trends and discrepancies between database systems using monthly counts of PLHIV who initiated antiretroviral treatment (ART) reported in the year 2019 (ie, before the COVID pandemic interruptions). We consider 69 clinics across the North West and Free State provinces, from which 36 FEDISA TB study clinics were selected.
Multiplicative interaction regression models for longitudinal counts of ART initiations
The data have a hierarchical structure with monthly counts of ART initiations at the lowest level nested within each facility, and each facility’s data are recorded in the two database sources. For each facility, a maximum of 12 counts for the 2019 calendar year were abstracted from each database by the FEDISA study staff.
We used a hierarchical model to characterise the parameters at different levels of this structure and we implemented this model within the Bayesian framework where all parameters are considered random draws from some known distributions. We denote the outcome indicator, which is the count at month j in clinic k, as . The counts were highly skewed and they were log-transformed to reduce the variability in the data and to ensure that the predicted values remain above 0 when the normal parametric distribution is applied. Henceforth, the count outcome refers to the log-transformed values. We consider a model that characterises the changes in the clinic’s reported counts in each month as:
where
μ
is the intercept and reflects the kth clinic’s deviation from this intercept, and is the month effect, is the facility’s coefficient that reflects the amount of change in the count at the jth month; and is the random error term which is assumed to be distributed as an independent normal variable with mean 0 and variance . This model is reparametrised so that the clinic-specific mean is centred on .
The random error variance is allowed to vary between facilities so that ~ . This type of model has been used in various contexts including food or drink product rankings by human assessors, calibration of measuring instruments and sensitivity of plant varieties to different growing conditions.12–16 In these contexts, would be the food product or planting environment effects, while has been interpreted as the relative expansiveness of scale use by the assessors or the relative sensitivity of plant varieties to the environment. One of the goals in the context of food product tasting has been to estimate the assessor performance across experiments relative to other assessors in the absence of an (objective) measure of true stimulus intensity. In other words, there is a heterogenous structural relation between product and assessor, for example. Thus, the stimulus intensity is realised through the assessors’ scale expansiveness. In our context, we analogously use this model to quantify relative levels of month-to-month fluctuations in the number of ART initiations and study overall time trends as informed by these fluctuations, and to use facility-specific indicator averages estimated from this to quantify discrepancies between the data sources.
Thus, is modelled as a normally distributed random variable with mean 0 and variance . This parametrisation is referred to as centreing and it enables the interpretation of as the facility’s mean count. The ’s are constrained to have mean 1 for model identifiability, and they quantify the facility’s relative agreement with the overall trend observed in the same period.
To model data from multiple data sources, the multiplicative interaction model was extended to accommodate this layer of the hierarchy (denoted by subscript l) to become
In this form, each facility’s intercept, slope and random error vary by source, and we assume that the time effect is the same as the data are contemporaneous.
The dataset also included self-reported information by the facilities on the number of staff nurses who prescribe ART and TPT in the clinics, and there were differences in the number of ART participants by province. Thus, the model was extended to include the effects of these covariates which were categorised as province 1 and 2; and number of nurses up to four versus more than four, respectively. Thus, the models become
Bayesian model implementation and prior values
The advantages of using the Bayesian framework here include the flexibility to place constraints on the slope parameters; avoiding unreasonable error variance estimates, which is a problem that has been previously observed for these types of models12–14; and to enable inclusion of expert knowledge into the analysis. The Bayesian model estimation involves combining prior information about parameters to be estimated with the likelihood of the observed data to obtain a plausible distribution (called the posterior distribution) of the parameters conditional on the data. Likely parameter values are then repeatedly sampled from this distribution to draw inference. We implemented these models using the R Jags package.17 This uses Monte Carlo Markov Chain sampling (MCMC)18 to draw a series of estimates from the posterior distributions of the unknown parameters. Any missing count outcomes are assumed to occur at random and are imputed with the average estimated by the MCMC algorithm. The first set of samples drawn is discarded (ie, considered burn-in) and the rest are used to provide final posterior means and standard errors of the parameters, and to generate 100(1−α)% hightest posterior density intervals (HPDI) for the parameters. In simple terms, the HPDIs are the range of values within which 100(1−α)% of the posterior probability density lies. In our implementation, we generated 9000 MCMC samples after a burn-in of 1000 samples. The hierarchical model assumes the missing data are missing at random and imputes the average for those facilities which may have missing indicators for some of the months.
Table 1A lists the prior distributions of the model parameters and table 1B lists the corresponding prior values used for both Models 1b and 2b. The prior values of the parameters are chosen based on our experience conducting cluster-randomised trials involving facilities that provide TB and HIV services in South Africa. For mathematical convenience, the prior distributions are chosen to be conjugate with the likelihood distributions. Thus, the log of ART initiation counts, the mean parameters and slopes were modelled with normal distributions with means and variances as listed in table 1. The variances were assumed to follow an inverse gamma (IG) distribution with df ‘d*’ that reflect the uncertainty around the parameter values, and the sample variances ‘s*’, sigma*−2~gamma(d*/2, d*s*2/2) which is equivalent to d*s2* sigma*−2~X2(d*).19 Based on prior experience with the data from facilities in South Africa, we assume that across clusters there is a median of about 20 ART initiations per facility per month with a skewed distribution that can range from 0 for small facilities to more than 50 in large facilities. We conducted a robustness check by running these analyses using non-informative prior distributions.
Table 1
|
(A) Prior distributions. (B) Prior values
Results
Characterising trends within a single data source
Our method does not assume that one data source is superior to the other and thus, for the sake of neutrality, the results are reported without specifically naming the source. Figure 1 shows plots of reported ART initiations for the two data sources. These indicate some clear divergence in a few facilities, for example, 17, 24 and 30. Figure 2 shows monthly reported counts and superimposed predicted posterior means (dashed lines) from the model for facilities reported in source 1 for 16 of the facilities (the rest of the facilities are shown in online supplemental figure 1). These indicate good predictive properties of the model. The traceplots of the slope () estimates are shown in online supplemental figure 2 indicating a consistent pattern throughout the MCMC iterations. There were 7 missing data points (4 for facility 15 and 3 for facility 56). These plots show differences in patterns and slopes of counts over time between facilities, making the multiplicative model a reasonable approach. There are roughly common time trends where ART initiations in December tend to be low. Figure 3 shows the posterior means of the estimated latent time effects and the corresponding 95% highest posterior density region. This shows relatively even counts for all months except a slight dip in June and a large dip in December. Figure 4 shows estimated facility-specific slopes versus within-facility variance . This shows a cluster of facilities which are in agreement with the estimated time trends as informed by the contribution of all facilities’ data, and these have relatively low estimated within-facility variance. In contrast, for example, facility 15 has 4 months’ worth of missing data, high variability and a flat trend.
Monthly reported counts (solid lines) and predicted posterior means (dashed lines) for 16 facilities reported in source 1. The labels above the plots indicate the facility number, the estimated mean counts and their SD, and the facility-specific slopes . ART, antiretroviral treatment.
Facility-specific slopes versus within facility variance .
Characterising trends and evaluating discrepancies between data sources
Figure 1 of the reported ART initiations charting two data sources shows that 7 out of 69 (10%) of the facilities show some noticeable divergence in the trend of this indicator. In contrast with source 1 which had only 7 (0.8%) missing data points, source 2 had a total of 44 (5%) missing data points across 8 facilities, with facilities 51, 57 and 58 missing 9 data points each, facility 35 missing 7 and the rest missing 4 or less data points. The missing values were reported by source 2 when source 1 was reporting very small monthly numbers (0, 1 and 2) of ART initiations for those clinics. This suggests zeros may have been entered as missing values, and this would warrant further investigation to reconcile the two sources.
Figure 5 shows differences in posterior mean intercepts from the two sources and this shows similar trends except for 7 (10)% facilities that diverged some. The model results showed a median difference in the average number of reported ART initiations of 1.6 (IQR=0.8–3.3) between the two data sources, and particularly divergent facilities 17, 24 and 69 which differed by an average of 9.5, 10 and 19, respectively, between the two sources. The model results accurately reflect the observed data as shown in figure 6 where for facilities 69 and 24, the reported data trends for the two sources clearly diverge while for clinics 6 and 42 the counts from the different sources closely match. Figure 7 shows the estimated facility-specific slopes from source 1 versus source 2 with source 2 showing lower estimates overall. The estimated time effects from the combined data sources analysed together are similar to that of the source 1 data alone (online supplemental figure 3).
Posterior mean estimates of the multiplicative slope from the two sources, text labelling refers to clinic numbers.
Sensitivity analysis
We conducted a robustness check in which we defined non-informative priors for the beta2, beta 3~N(0,0.001) and all the variance parameters set to IG(1,1) as well as IG (0.1,0.1) and this led to the same posterior inference.
Discussion
We have demonstrated an analytical approach to describe time effects in longitudinal healthcare data; determine differences in patterns between facilities; and quantify discrepancies when data are captured in more than one database. Multiplicative interaction models enabled us to characterise the heterogeneous changes of facility-level number of ART initiations over time. The estimated averages that account for these heterogeneous changes were then used to assess and quantify discrepancies between the two data sources. This model assumed no known true time effects and that these are realised through the magnitude of the changes in facility-reported data points.
With respect to ART initiation counts used for our illustration, the time effect of ART initiations remained fairly constant over the year with a slight increase in May and a dip in June, as well as a large dip in December. It is common to observe reduced care-seeking or diagnoses of diseases during the Christmas holiday season in South Africa. The slight May–June changes may be an indication of changes in care-seeking behaviours during this time when there is a considerable drop in temperatures as the winter season sets in, as well as the start of the schools’ winter break. The two data sources reported very similar trends and the model results indicate a median difference of less than two average number of ART initiations for the majority of the facilities. Three facilities that diverged by 10 or more average number of ART initiations warrant closer investigations.
Our analytical approach enables one to study trends for different facilities (or other cluster designations) in more detail than the graphical evaluation of longitudinal healthcare data. For example, Benade et al3 showed the effect of the COVID-19 pandemic on trends in ART initiations in South Africa using DHIS V.2 data. If there is a desire to further quantify the heterogeneous impact of the pandemic in different geographical areas or facilities, the multiplicative interaction model would be an ideal tool for such an analysis. The model can also be used for prediction of future realisations of the indicators.
We implemented this analysis in the Bayesian framework due to the complexity of the model structure and the flexibility to extend the model to predict future outcomes. If one was seeking a predictive model, model comparisons could be conducted using deviance information criterion or its alternatives as discussed by Spieglehalter et al.20 We evaluated model robustness of posterior estimates to changes in prior values of the variance parameters and found that these changes did not alter the posterior inference. Nonyane and Theobald, Brockhoff and Skovgaard and Smith et al12–14 noted challenges in implementing this type of model in the frequentist framework including zero estimates of residual variance estimates when they are modelled to differ between individuals, as well as unreliable maximum likelihood estimates for the individual slopes. Nabugoomu et al21 described a residual maximum likelihood estimation in models with location by variety interaction terms to estimate the means and sensitivity coefficients, with locations treated as random effects. Pødenphant and Kristensen developed an R package (Mumm) for a mixed model with an interaction between the random error term (for cluster) and a fixed term.22 Other approaches to evaluate trends in longitudinal data and possible discrepancies may be considered. An example is the latent variable modelling, which identifies latent clusters of longitudinal data.23 This, however, does not enable one to evaluate the heterogeneous slopes for each individual unit of analysis’ such as facilities in our data, but instead, group trends and slopes are quantified.
The accuracy and reliability of longitudinal healthcare data sources is crucial for monitoring programmes and making policy decisions as well as conducting research. Our purpose was not to evaluate the accuracy of the reported indicator values but rather to demonstrate a method for assessing trends and discrepancies when these data are reported by different sources. Evaluations of the accuracy have been done for these databases using various approaches. For example, in one of the earlier studies, Garieb et al9 conducted a study on a small sample of DHIS-reporting clinics and found problems with missing and out-of-range values which they attributed to human error. Over the years, improvements to the system and trainings have been conducted. Discrepancies between Tier.net or DHIS V.2 and other data sources have also been reported. For example, Etoori et al10 found that, in comparison to local HIV surveillance data, Tier.net misclassified 36% of patient outcomes of losses to follow-up, mortality and transfers out.
Strengths and limitations
The Bayesian analytical framework for this application is a powerful tool that enables one to incorporate prior expert knowledge of the unknown parameters to the analysis, as well as to estimate multiple parameters efficiently. Over the last 20 years, easy-to-use open-source software packages for Bayesian analysis have become available, making the Bayesian implementation feasible. We note, however, that obtaining prior information can also be a limitation. We applied this model to a normally distributed outcome, but it is straightforward to generalise to other distributions of the indicators as well.
Conclusion
For the number of ART initiations over time, estimated overall trends over the course of 2019 remained steady with a reduction during the December holiday period. The two data sources show similar reporting for most of the facilities with an estimated median difference of number of ART initiations of 1.6 (IQR: 0.8–3.3). Multiplicative interaction models are a powerful tool for quantifying average trends over time and for evaluating discrepancies between reporting systems. The Bayesian framework enables efficient parameter estimation in this context with a very large number of parameters. We recommend that researchers consider this analytical approach to evaluate trends in routine clinical data and for identifying deviations that warrant additional attention.
