Intervention Time Series Analysis and Forecasting of Organ Donor Transplants in the US during the COVID-19 Era

Abstract

The COVID-19 pandemic has had a catastrophic effect on the healthcare system including organ transplants worldwide. The number of living donor transplants performed in the US was affected more significantly by the pandemic with a 22.6% decrease in counts from 2019 to 2020 due to concerns of unnecessarily exposing potential living donors and living donor recipients to possible COVID-19 infection. This paper examines donor transplant counts obtained from the United Network for Organ Sharing from January 2002 to August 2021 using an intervention time series model with March 2020 as the intervention event. Specifically, donor transplant counts are analyzed across the different organs, donor types, and some major individual sociocultural factors, which are potential conditions contributing to disparities in achieving donor transplant equity such as age, ethnicity, and gender. In addition, the kidney allocation policy implemented in March 2021 is introduced as a second intervention event for kidney donor transplants. Overall, forecasts generated by our methods are more accurate than those using seasonal autoregressive integrated moving average models without interventions and seasonal naive methods. The intervention time series model provides a forecast accuracy comparable to the exponential smoothing method.

1. Introduction

More than 107,000 patients in the United States are awaiting organ transplants. Every day, an estimated 33 people in the United States die while waiting for their transplant. For the past four years, on average, only 37,500 organs were transplanted annually, including around 31,000 organs from deceased donors and 6500 organs from living donors [1]. The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that caused the coronavirus disease (COVID-19) pandemic has, among many other things, devastatingly impacted solid organ transplantation (SOT) worldwide, in turn affecting potential donors, candidates, and recipients. While the number of deceased donor transplants remained unaffected, there was an overall decrease of 22.6% from 2019 to 2020 in the number of living donor transplants performed in the US [2]. Ref. [3] noted a strong temporal association between the increase in COVID-19 infections and a striking reduction in overall SOT procedures. Although deceased donation showed an increase of 10.1% and living donor transplants showed a 14.2% increase in 2021 over 2020, their totals were still lower compared to previous years [4]. Due to limited resources and preparedness for the pandemic, the healthcare system, including the organ transplant system, was disrupted. There were some predictive tools and platforms developed as a response to the pandemic. Ref. [5] introduced tools that used a machine learning method to predict the need for hospitalization using a large set of different COVID-19 patient types in southeastern Spain. The CIoTVID IoT platform was presented in [6] to help fight the propagation of the COVID-19 pandemic. To better understand the COVID-19 pandemic’s effect on SOT, it is important to statistically analyze, model, and forecast the profound impact of COVID-19 across various social factors that lead to inequities in donor transplants. The evaluation of the impact of any organ allocation policy that was introduced as an emergency response to the pandemic is also very essential. For example, on 15 March 2021, the organ procurement and transplantation network (OPTN) launched a new policy to match kidney transplant candidates with organs from deceased donors [7]. This kidney allocation policy replaced the 11 OPTN organ allocation regions with a more consistent measure of distance between the donor hospital and the transplant hospital for each candidate and aimed to increase equity in access to transplants. This policy was intended to reduce geographical disparity, pretransplant deaths, and maximize pediatric transplantation. Although the organ allocation system has been subjected to multiple changes with the introduction of new policies over time, our interest lay here in testing the effect of the most recent kidney allocation policy since it brought big changes to the organ allocation system.

Literature on the medical or clinical outcomes of solid organ transplantation is by now vast. However, there are limited studies available specific to analyzing donor transplant counts and forecasting them using a time series analysis in the COVID-era. Ref. [8] examined deceased donor and living donor kidney transplant trends using density mapping and linear regression with an interrupted time series analysis, which is also called intervention time series analysis (ITSA). They found that kidney transplants were affected significantly in recent months due to COVID-19. In [9], adult transplantation data from 1990 to 2019 were analyzed using an autoregressive integrated moving average model (ARIMA) in order to forecast the 2020 expected rates of organ transplants and waitlist registrations. A significant association between the pandemic and the deficit in kidney transplantation and waitlist registration existed. The ARIMA model with interventions could be an efficient method for policymakers to generate forecasts and to evaluate policy interventions or the effect of any major incident affecting donation rates, waitlist registrations, and transplant rates [10]. It is very flexible and easily accounts for seasonality, autocorrelation, underlying trends, and the different types of intervention effects [11,12]. Ref. [13] employed a quarterly interrupted time series regression analysis for estimating the impact of EU regulatory label changes for diclofenac in 2013 among people with cardiovascular disease in Denmark, the Netherlands, England, and Scotland. A review of the interrupted time series analysis in drug utilization research can be found in [14].

In this study, we model the effect of the COVID-19 pandemic as an intervention event on donor transplant time series in the US across different donor types, organs, and select major individual sociocultural factors such as age, ethnicity, and gender, contributing to disparities in achieving donor transplant equity [15]. Research in the past and data over the years have demonstrated that disparities occur at every stage of the transplant process and waitlist registrations with different factors contributing to each [16,17]. Accurate forecasting methods could help public health response with better preparedness and mitigation strategies in events such as the pandemic and any seasonal episodes. Additionally, the kidney allocation policy introduced in March 2021 is modeled as a second intervention event for kidney donor transplants. The forecasts of organ transplant counts from the intervention time series analysis model are then compared with those from the seasonal ARIMA (SARIMA) model without intervention, the exponential smoothing (ETS) method, and the seasonal naive (SNAIVE) method. Studies like ours could help improve the quality and response for maximizing donors in the organ transplant system.

The rest of this paper proceeds as follows. Section 2 introduces the intervention time series analysis model that underlies our work and briefly reviews the ETS and SNAIVE forecasting methods. The forecast accuracy measures used in the study as well as the Diebold–Mariano tests for equality of forecast accuracy are also presented in Section 2. Section 3 summarizes the model estimates and forecasts for various donor transplant time series. Concluding remarks and future works close the article in Section 4.

2. Materials and Methods

2.1. Data Sources

Monthly donor transplant data in the US across donor type, organs, age group were collected from the United Network for Organ Sharing (UNOS) [2,18] for the period from January 2002 to August 2021. For ethnicity and gender, monthly donor transplant data were collected from January 2002 to May 2021. Since the COVID-19 pandemic was declared a national emergency on 13 March 2020, the OPTN approved emergency policy action on 17 March 2020. Thus, all the data prior to March 2020 were considered preintervention with March 2020 being the intervention event for the study. For kidney transplants, the kidney allocation policy was considered as the second intervention event in March 2021.

2.2. Intervention Time Series Analysis

Intervention time series analysis is commonly used to measure the effect of an intervention event on the time series of interest [19]. The time of occurrence of the intervention event is assumed to be known and the intervention can be natural or man-made [20]. In general, the intervention event affects the time series by changing the mean function of the process. In this study, we investigated the effect of the COVID-19 pandemic on organ donor transplantation time series.

The time series model we consider in this section with a single intervention event that occurred at a known time T has the form:

$$Y_t=m_t+N_t$$

(1)

where $m_t$ denotes the intervention effect on time series $Y_t$, and $N_t$ represents the underlying time series without the effect of the intervention event. Here, the effect of the COVID-19 pandemic is described by the model

$$m_t=\left({\omega }_0+\frac{{\omega }_{1}B}{1-\delta B}\right)P_t^{\left(T\right)},0<\delta <1$$

(2)

where B is a backward shift operator, which is defined by $B{Y}_t=Y_{t-1}$. $P_t^{\left(T\right)}$ is a pulse function at time T that is used to model an intervention with temporary effect, given by

$$\begin{array}{c}\hfill P_t^{\left(T\right)}=\left\{\begin{array}{cc}1,\hfill & t=T\hfill \\ 0,\hfill & t\ne T.\hfill \end{array}\right.\end{array}$$

The model in Equation (2) can be rewritten as

$$m_t=\delta m_{t-1}+{\omega }_0{P}_t^{\left(T\right)}-{\omega }_0\delta P_{t-1}^{\left(T\right)}+{\omega }_1{P}_{t-1}^{\left(T\right)}.$$

(3)

Note that when $\delta =0$, it reduces to an MA(1) type filter:

$$m_t=({\omega }_0+{\omega }_{1}B)P_t^{\left(T\right)}.$$

In the case when ${\omega }_0=0$, it is called an AR(1) type filter:

$$m_t=\delta m_{t-1}+{\omega }_1{P}_{t-1}^{\left(T\right)}$$

with a delay of one time-unit before the intervention takes effect. Here, T denotes March 2020 and we assume that the initial condition $m_0=0$. Thus, $m_T={\omega }_0$, $m_{T+1}={\omega }_1$, and for $t>T+1$, $m_t={\omega }_1{\delta }^{t-T-1}$, which would represent the transient effect of the COVID-19 pandemic. Specifically, the COVID-19 pandemic caused an immediate change in the organ transplant counts with a magnitude ${\omega }_0$ in March 2020 and a strong decreasing effect with a magnitude ${\omega }_1$ in April 2020, followed by a gradual decay at the rate $\delta$ back to the original pre-COVID level with no permanent effect. Here, ${\omega }_0$, ${\omega }_1$, and $\delta$ are unknown parameters to be estimated.

With the intervention effect specified in Equation (2), only the mean level of the time series $\left\{Y_t\right\}$ in Equation (1) is affected by the intervention event. Phrased another way, the time series model for $\left\{N_t\right\}$ is the same before and after the intervention event. Therefore, the series $\{Y_t,t<T\}$, which is referred to as the preintervention data, can be used to identify the time series model for $\left\{N_t\right\}$. Here, $\left\{N_t\right\}$ is typically referred to as an unperturbed process. It is assumed to follow a multiplicative SARIMA for the monthly organ transplant count time series in this study. The seasonal period is $\nu =12$ for monthly data. The SARIMA model is denoted by SARIMA($p,d,q$) × ($P,D,Q$)${}_{\nu }$ and defined as

$${(1-B)}^d{(1-B^{\nu })}^D{N}_t=\frac{\theta \left(B\right)\Theta \left(B^{\nu }\right)}{\varphi \left(B\right)\Phi \left(B^{\nu }\right)}{e}_t,$$

(4)

where $\left\{e_t\right\}$ is white noise with mean zero and variance ${\sigma }_e^2$, and d and D are nonnegative integers. In this study, $d=1$ and/or $D=1$ are used to make $N_t$ stationary through differencing. The regular (nonseasonal) and seasonal operators in Equation (4) are specified as

$$\begin{array}{ccc}\hfill \varphi \left(B\right)& =& 1-{\varphi }_{1}B-{\varphi }_2{B}^2-\dots -{\varphi }_p{B}^p\hfill \\ \hfill \Phi \left(B^{\nu }\right)& =& 1-{\Phi }_1{B}^{\nu }-{\Phi }_2{B}^{2\nu }-\dots -{\Phi }_P{B}^{P\nu }\hfill \\ \hfill \theta \left(B\right)& =& 1+{\theta }_{1}B+{\theta }_2{B}^2+\dots +{\theta }_q{B}^q\hfill \\ \hfill \Theta \left(B^{\nu }\right)& =& 1+{\Theta }_1{B}^{\nu }+{\Theta }_2{B}^{2\nu }+\dots +{\Theta }_Q{B}^{Q\nu }\hfill \end{array}$$

where $\varphi \left(B\right)$ is the regular AR operator with p parameters of ${\varphi }_1,\dots ,{\varphi }_p$; $\theta \left(B\right)$ is the regular MA operator with q parameters of ${\theta }_1,\dots ,{\theta }_q$; $\Phi \left(B^{\nu }\right)$ is the seasonal AR operator with P parameters of ${\Phi }_1,\dots ,{\Phi }_P$; and $\Theta \left(B^{\nu }\right)$ is the seasonal MA operator with Q parameters of ${\Theta }_1,\dots ,{\Theta }_Q$. In total, there are $p+P+q+Q+1$ nonzero parameters in Equation (4). These parameters were estimated by the maximum likelihood method.

Besides the COVID-19 pandemic’s effect on the organ donor transplantation time series, the kidney donor series was modeled with an additional intervention event for the kidney allocation policy in March 2021. The two intervention events were modeled by

$$m_t=({\omega }_0+{\omega }_{1}B+{\omega }_2{B}^2)P_t^{\left(T\right)}+{\omega }_3{S}_t^{\left(T_k\right)},$$

(5)

where $T_k$ denotes March 2021. As above, T denotes March 2020. The COVID-19 pandemic effect here is specified by an MA(2)-type filter in the first term of Equation (5). ${\omega }_0$ represents the instantaneous pandemic effect in March 2020. ${\omega }_1$ and ${\omega }_2$ give the effect of the pandemic in April 2020 and May 2020, respectively. Next, the permanent intervention effect of the kidney allocation policy was modeled using a step function $S_t^{\left(T_k\right)}$, which was given by

$$\begin{array}{c}\hfill S_t^{\left(T_k\right)}=\left\{\begin{array}{c}1,t\ge T_k\hfill \\ 0,t<T_k.\hfill \end{array}\right.\end{array}$$

That is, $S_t^{\left(T_k\right)}=0$ during the preintervention period and $S_t^{\left(T_k\right)}=1$ throughout the postintervention period. ${\omega }_3$ in Equation (5) is the effect of the kidney allocation policy since March 2021. Notice that $S_t^{\left(T_k\right)}$ can be also expressed in terms of $P_t^{\left(T_k\right)}$:

$$S_t^{\left(T_k\right)}=\frac{1}{(1-B)}{P}_t^{\left(T_k\right)}.$$

2.3. Model Selection

The intervention time series analysis worked as follows. First, the SARIMA model in Equation (4) was fitted to the preintervention data. The best model was selected using a bias-corrected Akaike information criterion (AICc), which could be obtained as

$$\mathrm{AICc}=-2ln\widehat{L}+2\kappa +2\kappa \frac{(\kappa +1)}{n-\kappa -1},$$

(6)

where $\kappa$ is the number of estimated parameters in the model and $\widehat{L}$ is the maximum likelihood function value for the model. The smaller the AICc value, the better the model fit. Next, the selected best model for $\left\{N_t\right\}$ was applied to the whole series with the COVID-19 intervention effect in Equation (2).

2.4. Forecast Accuracy

In order to evaluate the forecast performance, we divided the organ transplant time series into two datasets: the training dataset ranging from January 2002 to August 2020, and the test dataset ranging from September 2020 to August 2021 for donor type, age group, and organ type, whereas for ethnicity and gender, the test dataset ranged from September 2020 to May 2021. Forecasts were then calculated at the selected forecast horizons (h) on a rolling forecasting origin.

The accuracy of the forecasts was measured by the mean absolute percentage error (MAPE) and root mean square error (RMSE). They were computed by averaging over the test sets:

$$\mathrm{MAPE}=\frac{{\displaystyle \sum _{i=1}^{n_f}}\left|\frac{{\widehat{Y}}_i-Y_i}{Y_i}\right|}{n_f}$$

(7)

and

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^{n_f}}{({\widehat{Y}}_i-Y_i)}^2}{n_f}}$$

(8)

where ${\widehat{Y}}_i$ is the forecast value for $Y_i$, and $n_f$ is the total number of forecasts. The lower values of MAPE and RMSE indicate better forecasts.

In this study, we compared our ITSA model with three forecasting methods for the monthly donor transplant time series: SARIMA without intervention, ETS, and SNAIVE.

2.5. Exponential Smoothing

An exponential smoothing method is an algorithm that generates point forecasts. The forecasts are weighted averages of past observations, with the weights decaying exponentially as the observations get older. Ref. [21] developed a fully automatic forecasting procedure based on exponential smoothing methods in a state-space framework.

For each exponential smoothing method, the underlying state-space model has either additive (A) errors or multiplicative (M) errors. These two models produce equivalent point forecasts but different prediction intervals. In addition, each method has a trend component and a seasonal component. The possibilities for the two components described in [22] are {none (N), additive (A), additive damped ($A_d$), multiplicative (M), and multiplicative damped ($M_d$)} for the trend component and {none (N), additive (A), and multiplicative (M) } for the seasonal component.

The triplet ($E,T,S$) refers to the error, trend, and seasonality. For example, the model ETS(A, N, N) has additive errors, no trend, and no seasonality, which describes the simple exponential smoothing method with additive errors. Similarly, the model ETS(A, A, N) has additive errors, additive trend, and no seasonality, which describes Holt’s linear method with additive errors. The Holt–Winters’ multiplicative method [23,24], denoted by ETS(M,A,M), has multiplicative errors, additive trend, and multiplicative seasonality. In total, there are 30 potential ETS models, each corresponding to a combination of different types of errors, trends, and seasonality. Note that six additive ETS models { (A,N,N), (A,A,N), (A,$A_d$,N), (A,N,A), (A,A,A), (A,$A_d$,A)} are equivalent to ARIMA models.

For seasonal time series forecasting in the study, we illustrated the ETS methods with seasonality using ETS(M,A,M). Details for other ETS models with seasonality are available in [21,22].

The point forecasts using ETS(M,A,M) are given by

$${\widehat{Y}}_{t+h|t}=({\ell }_t+h{b}_t)s_{t+h-\nu (k+1)},$$

(9)

where ${\widehat{Y}}_{t+h|t}$ denotes the h-step-ahead forecast based on all the observations up to time t. ${\ell }_t$ represents the level of the series at time t, $b_t$ is the trend at time t, and $s_t$ denotes the seasonal component of the series at time t. Here, $\nu$ is the seasonal period, and k is the integer part of $\frac{h-1}{\nu }$. ${\ell }_t$, $b_t$, and $s_t$ are modeled as follows:

$$\begin{array}{ccc}\hfill {\ell }_t& =& \alpha \frac{Y_t}{S_{t-\nu }}+(1-\alpha )({\ell }_{t-1}+b_{t-1})\hfill \\ \hfill b_t& =& \beta ({\ell }_t-{\ell }_{t-1})+(1-\beta )b_{t-1}\hfill \\ \hfill s_t& =& \gamma \frac{Y_t}{({\ell }_{t-1}+b_{t-1})}+(1-\gamma )s_{t-\nu }\hfill \end{array}$$

The smoothing parameters $\alpha$, $\beta$, and $\gamma$ need to be estimated from the data. The equivalent state-space model for ETS(M,A,M) has the form:

$$\begin{array}{ccc}\hfill \mathrm{Observation}\mathrm{equation}:Y_t& =& ({\ell }_{t-1}+b_{t-1})+s_{t-\nu }(1+{ϵ}_t)\hfill \\ \hfill \mathrm{State}\mathrm{equation}:{\ell }_t& =& ({\ell }_{t-1}+b_{t-1})(1+\alpha {ϵ}_t)\hfill \\ \hfill b_t& =& b_{t-1}(1+\alpha \beta {ϵ}_t)\hfill \\ \hfill s_t& =& s_{t-\nu }(1+\gamma {ϵ}_t)\hfill \end{array}$$

The errors $\left\{{ϵ}_t\right\}$ are independent and identically distributed Gaussian random variables with mean zero and constant variance ${\sigma }^2$. This state-space model generates point forecasts identical to those from Equation (9) but also generates prediction intervals and allows for a model selection.

In general, the automatic ETS forecasting procedure fits each of the 30 state-space models and then selects the one that fits the data best. In this study, the best ETS model is selected by the AICc.

2.6. Seasonal Naive

A seasonal naive model is a seasonal benchmark model for time series forecasting with seasonality. It works well for highly seasonal data. The h-step-ahead point forecast using data up to time t can be obtained by

$${\widehat{Y}}_{t+h|t}=Y_{t+h-\nu (k+1)}.$$

(10)

$\nu$ and k in Equation (10) are defined the same as in Equation (9). The forecast using SNAIVE is set to be the last observed value from the same season. For example, the forecast for next April is equal to the last observed April value in the series. Thus, it depends only on one observation.

2.7. Diebold–Mariano tests

The forecast accuracy measures mentioned in Section 2.4 are commonly used to compare forecast performance among two or more forecast methods. However, the difference in the MAPE or RMSE between the two forecast methods might just be due to random chance. Ref. [25] proposed the Diebold–Mariano (DM) test for assessing the equality of forecast accuracy.

Let $u_t$ be the h-step-ahead forecast error at time t. The loss associated with forecast error is denoted by $L\left(u_t\right)$. Now, define the loss differential at time t between forecasts one and two as $d_t=L\left(u_{t,1}\right)-L\left(u_{t,2}\right)$. Thus, the null hypothesis of equal forecast accuracy for two forecasts is $E\left(d_t\right)=0$. That is, the expected value of the loss differential series is 0. Note that the forecast errors are typically autocorrelated. The DM test assumes that the loss differentials are weakly stationary and have short memory. The DM test statistic is then given by

$$S=\frac{\overline{d}}{\sqrt{Var\left(\overline{d}\right)}}$$

(11)

where $\overline{d}$ is the sample average over the $n_f$ loss differentials,

$$\overline{d}=\frac{{\sum }_{t=1}^{n_f}{d}_t}{n_f}$$

(12)

and the variance of $\overline{d}$ can be calculated using

$$Var\left(\overline{d}\right)=\frac{{\gamma }_0+2{\sum }_{k=1}^{h-1}{\gamma }_k}{n_f}.$$

(13)

Here, ${\gamma }_k$ is the autocovariance of $\left\{d_t\right\}$ at lag k. It can be estimated by

$${\widehat{\gamma }}_k=\frac{{\sum }_{t=k+1}^{n_f}(d_t-\overline{d})(d_{t-k}-\overline{d})}{n_f}$$

(14)

Under the null hypothesis, S is asymptotically normally distributed with mean zero and unit variance. Thus, the DM test can be easily performed in practice.

For the finite sample, if $n_f$ is small, [25] reported that the DM test had a large type I error rate. That is, the DM test was more likely to incorrectly reject the null hypothesis of equal forecast accuracy when it is true for small $n_f$. With some modifications, [26] developed a modified DM test statistic

$$S^{*}=\sqrt{\frac{n_f+1-2h+\frac{1}{n_f}h(h-1)}{n_f}}S,$$

(15)

which improved the problem of over-size for the original DM test. The $S^{*}$ statistic here is compared with critical values from the Student’s t distribution with $(n_f-1)$ degrees of freedom. In this study, the modified DM test was used in Section 3 with a mean squared error loss function.

3. Results

Disclaimer: The data reported here were supplied by the United Network for Organ Sharing as the contractor for the Organ Procurement and Transplantation Network. The interpretation and reporting of these data are the responsibility of the author(s) and in no way should be seen as an official policy of or interpretation by the OPTN or the U.S. Government.

3.1. Donor Type

Donor transplant data consist of living and deceased donor type transplants. It is essential to maximize the number of living donors since it significantly decreases the wait time for patients in need of a transplant. It also helps people waiting for deceased donor organs by lowering the number of people on the waiting list. Living organ transplantation has many advantages over transplantation from a deceased donor [27]. Particularly for kidneys, the national wait time for a deceased donor kidney transplant can be as long as six years, and living donations significantly decrease the wait time for patients in need of a kidney transplant since the living kidney donor is often a close relative of the person getting the transplant, and there is a better chance of a good genetic match and less chance of rejection. As a result, living donor kidneys tend to last twice as long as deceased donor kidneys. Moreover, lower doses of immunosuppressive drugs may be used with fewer side effects. Figure 1 plots monthly counts of all donor transplants from 2002 to 2021 by donor type. The sudden fall in donor counts in March and subsequently in April 2020 was caused by the COVID-19 pandemic. It seems that the counts still did not go back completely to their prepandemic levels until August 2021.

A SARIMA$(1,1,2)\times {(0,0,2)}_{12}$ model with a period of 12 and the intervention effect in Equation (2) was fit to the living donor series up until August 2020. The monthly living donor transplant counts along with a model fit can be viewed in the top panel of Figure 2. The fitted model captured key features of the donor transplant counts fairly well, especially the effects of COVID-19, which are shown in the bottom panel of Figure 2. Notice that the COVID-19 pandemic caused a significant sudden change in the living donor transplant counts with a magnitude of −252.74 in March 2020. This further dropped to 509.71 in April 2020 followed by a gradual rebound at the rate of 33% to its pre-COVID level.

Figure 3 displays the living donor transplant counts from January 2019 to August 2021 and the forecasts from September 2020 to August 2021 using the ITSA model, SARIMA without intervention, ETS(M,N,A), and SNAIVE. Forecasts from ITSA and ETS performed similarly and were more accurate compared to those of SARIMA without intervention and SNAIVE. Living donor counts for the first six months were better forecasted using ETS and the latter ones using ITSA.

The accuracy measures are summarized in

Table 1. Forecast accuracy measures by donor type at various forecast horizons.

	h	MAPE (%)		RMSE
		Living	Deceased	Living	Deceased
ITSA	1	12.16	5.46	72.78	191.79
	2	16.76	6.08	93.09	210.46
	4	18.95	7.64	104.88	232.60
	8	9.58	3.97	65.78	183.95
SARIMA	1	13.49	6.27	113.19	235.95
	2	21.53	7.45	171.14	262.09
	4	26.73	8.22	195.57	274.17
	8	37.20	8.08	244.44	303.53
ETS	1	10.28	4.94	68.39	178.38
	2	10.57	6.09	69.86	204.10
	4	10.48	6.60	71.19	217.95
	8	11.45	7.39	76.62	240.55
SNAIVE	1	26.94	10.06	192.28	398.95
	2	29.35	10.81	200.83	416.41
	4	32.01	11.03	217.98	440.09
	8	29.25	12.63	240.27	513.26

for both living and deceased donor transplant time series at various forecast horizons. For example, the RMSE values of one-step-ahead forecasts were 72.78 and 68.39 for ITSA and ETS, respectively, which were much lower than the corresponding RMSE value of 192.28 from SNAIVE and a value of 113.19 from the SARIMA model forecasts without intervention. The estimated ITSA model parameters for both donor types are provided in

Table A1. ITSA model parameter estimates for donor types.

	Donor Type
Coefficients	Living	Deceased
${\varphi }_1$	−0.66
${\theta }_1$	−0.07	−0.80
${\theta }_2$	−0.75
${\Theta }_1$	0.60	−1.50
${\Theta }_2$	0.25	0.56
${\omega }_0$	−252.74	−559.54
$\delta$	0.33	0.15
${\omega }_1$	−509.71	−787.22

in the Appendix A.

Table 2. Results of the modified Diebold–Mariano test comparing ITSA to ETS, SARIMA without intervention, and SNAIVE by donor type.

Donor Type	h	ETS	SARIMA	SNAIVE
Living	1	0.32	−1.11	−1.66
	4	3.52	−1.47	−2.25
Deceased	1	0.46	−1.12	1.40
	4	0.14	−0.86	−1.19

Bold indicates a statistical difference between the two forecasting methods at the 10% level of significance. Negative (positive) values indicate that ITSA performed better (worse) than ETS, SARIMA without intervention, or SNAIVE.

provides the results of the modified Diebold–Mariano test for equality of forecast accuracy comparing ITSA to ETS, SARIMA without intervention, and SNAIVE across forecast horizons one and four. A bold value indicates that the two forecasting methods were significantly different at the 10% level of significance. If the test statistic value is negative, then it shows that ITSA outperformed its competing method. For example, the DM test statistic was $-2.25$ for the four-step-ahead forecasts of living donor counts comparing ITSA and SNAIVE, suggesting that ITSA performed significantly better compared to SNAIVE. Similarly, ETS outperformed ITSA for this example with $S^{*}=3.52$.

To further evaluate the forecasting performance of our ITSA model, we collected weekly metrics from the OPTN Metrics website [2] ranging from September 2021 to March 2022. Monthly living donor transplant counts were then obtained by summing the weekly donor counts for each month during this period. Here, the new calculated monthly living donor counts from September 2021 to March 2022 were treated as the test dataset and the donor counts from January 2002 to August 2021 as the training dataset.

The fitted living donor counts with the training data are plotted in the top panel of Figure 4. Evidently, the ITSA model depicted the effects of the COVID-19 pandemic efficiently. The pandemic effects shown in the bottom panel of Figure 4 indicate that living donor transplants regained their losses after about two years after March 2020.

Forecasts for the test data are displayed in Figure 5 and matched up with the data reasonably well. In addition, the forecasts were not quite as volatile as the actual observations.

3.2. Organ

Heart, kidney, liver, and lung donor transplant counts were analyzed in our study. Figure 6 plots the monthly transplant counts of the four organs from both living and deceased donors. Apparently, the number of kidney donor transplants was affected most significantly by the COVID-19 pandemic, which is our focus here.

A SARIMA$(0,1,2)\times {(0,1,1)}_{12}$ model with two intervention effects, namely, the occurrence of COVID-19 in March 2020 and the introduction of the kidney allocation policy in March 2021 in Equation (5) was fit to the kidney donor series from January 2002 to August 2021. The model fit plotted in the top panel of Figure 7 follows the observed kidney transplant trend very closely. The estimated effects of both COVID-19 and kidney allocation policy change are shown in the bottom panel of Figure 7.

As expected, the COVID-19 pandemic caused a significant sudden change in kidney donor transplant counts. The number of kidney donor transplants dropped by a magnitude of 438.85 in March 2020, 869.88 in April, and 317.91 in May. The effect of the kidney allocation policy appeared positive with a magnitude of 28.74, suggesting an increase in the number of kidney donor transplants due to the new kidney allocation policy. Given that only five months were included in the model fitting after the policy took effect, it was not surprising that its effect was insignificant.

In addition, for heart, liver, and lung donor transplant counts, we fit an ITSA model with only the COVID-19 pandemic as the intervention event. The estimated model parameters for all the organs focused on in this study are provided in

Table A2. ITSA model parameter estimates for organs.

	Organ
Coefficients	Heart	Liver	Lung	Kidney
${\theta }_1$	−0.87	−0.78	−0.81	−0.59
${\theta }_2$				0.04
${\Theta }_1$	−0.74	−0.81	0.27	−0.87
${\omega }_0$	−42.21	−107.91	−50.18	−438.85
${\delta }_1$	0.04	0.1	0.21
${\omega }_1$	−88.75	−163.80	−78.21	−869.88
${\omega }_2$				−317.91
${\omega }_3$				28.74

in the Appendix A.

Next, we compared the forecasts from September 2021 to April 2022 using the ITSA model with intervention events to those from SARIMA without intervention, ETS(M,A,M), and SNAIVE. All forecasts are plotted in Figure 8. Our method performed well, especially for the latter months of 2021.

Table 3. Forecast accuracy measures by organ type at various forecast horizons.

	h	MAPE (%)				RMSE
		Heart	Liver	Lung	Kidney	Heart	Liver	Lung	Kidney
ITSA	1	3.69	5.50	11.68	9.97	13.94	55.40	27.19	181.45
	2	3.45	6.37	12.13	9.73	13.64	64.03	28.57	178.25
	4	3.63	7.67	14.16	11.23	14.23	79.37	31.07	182.16
	8	3.13	9.36	7.23	-	11.27	88.98	16.97	-
SARIMA	1	4.64	6.11	11.73	9.46	18.06	60.76	30.39	166.30
	2	4.92	7.28	13.03	10.28	19.03	71.05	33.79	170.12
	4	4.94	8.79	15.33	12.37	19.50	89.14	35.92	194.96
	8	5.36	11.77	11.92	-	23.02	110.43	31.76	-
ETS	1	4.59	5.83	8.78	9.26	16.72	57.05	24.22	165.18
	2	5.09	6.95	10.30	9.06	18.60	67.71	26.64	164.43
	4	5.58	8.84	10.78	11.35	21.00	91.24	26.25	180.73
	8	6.99	13.26	9.52	-	25.47	122.95	24.38	-
SNAIVE	1	9.33	11.62	15.47	20.36	45.48	128.04	43.15	352.98
	2	10.06	12.66	15.32	17.84	47.48	133.73	43.88	310.89
	4	0.08	13.99	18.02	15.35	45.23	145.76	48.37	269.19
	8	9.26	16.63	14.66	-	52.95	162.84	49.30	-

reports the MAPE and RMSE values for kidney, heart, liver, and lung at different forecast horizons. Overall, the model forecasts for heart and liver were better using ITSA, and ETS performed better for lung. For kidney, ITSA, SARIMA without intervention, and ETS performed similarly and were all better than SNAIVE.

The DM test results are summarized in

Table 4. Results of the modified Diebold–Mariano test comparing ITSA to ETS, SARIMA without intervention, and SNAIVE by organ type.

Organ	h	ETS	SARIMA	SNAIVE
Heart	1	−0.82	−1.18	−1.63
	4	−0.83	−0.79	−1.47
Liver	1	−0.36	−1.61	−2.21
	4	−1.08	-1.27	−2.21
Lung	1	1.04	−0.94	−1.38
	4	0.79	−0.84	−1.41
Kidney	1	0.88	0.90	−2.17
	4	0.08	-0.26	−1.44

Bold indicates a statistical difference between the two forecasting methods at the 10% level of significance. Negative (positive) values indicate that ITSA performed better (worse) than ETS, SARIMA without intervention, or SNAIVE.

for organs. Overall, ITSA, ETS, and SARIMA without intervention generated comparable forecasts across all organs. For liver and kidney, ITSA performed significantly better than SNAIVE.

3.3. Age Group

Donor age (including both living and deceased donors) was grouped into five intervals: 65+, 50–64, 35–49, 18–34, and ≤17. Figure 9 displays the monthly counts of all donor transplants performed in the US across the different age groups from January 2002 to August 2021. A significant drop in donor transplant counts around March 2020 was evident in all age groups, especially for ages over 35 years. It became noticeable early in the COVID-19 pandemic that older age was a risk factor for becoming severely ill with COVID-19. However, the virus’s impact on older adults went beyond a higher risk for serious infection. It may have been due to limited access to care for all health conditions, as well as considerable social and economic hardships [28]. Hence, in this subsection, our focus is on the 65+-year-old age group.

With the intervention effect in Equation (2), a SARIMA $(0,1,2)\times {(0,1,1)}_{12}$ was fit to the 65+-year-old donor transplant counts. The model fit and the estimated COVID-19 effects are plotted in Figure 10. The COVID-19 pandemic caused a sudden change in the 65+-year-old donor transplant counts with a magnitude of −174.72 in March 2020. This further reduced to −321.13 in April 2020 followed by a faster rebound at a rate of 13% to its pre-COVID level. For all age groups, the estimated ITSA model parameters are provided in

Table A3. ITSA model parameter estimates for age groups.

	Age Group
Coefficients	65+	50–64	35–49	18–34	$\le 17$
${\varphi }_1$			−0.16
${\theta }_1$	−0.64	−0.82	−0.81	−0.83	−0.90
${\theta }_2$	−0.11				−0.06
${\Phi }_1$			−0.37	−0.61	0.37
${\Phi }_2$				−0.26	0.15
${\Theta }_1$	−0.86	−0.98
${\omega }_0$	−174.72	−239.65	−112.85	−60.55	−29.93
$\delta$	0.13	0.06	0.23	0.45	0.30
${\omega }_1$	−321.13	−448.44	−265.88	−114.28	−77.29

in the Appendix A.

The forecasts generated for the 65+-year-old donor transplants from September 2020 to August 2021 are shown in Figure 11. ETS(M,A,M) was selected and performed better compared to other methods for the 65+-year-old age group. In addition,

Table 5. Forecast accuracy measures by age group at various forecast horizons.

	h	MAPE (%)					RMSE
		65+	50–64	35–49	18–34	≤17	65+	50–64	35–49	18–34	≤17
ITSA	1	4.12	6.81	5.69	6.71	6.91	41.27	104.08	56.11	29.30	16.88
	2	4.74	8.00	6.49	7.45	7.06	46.62	116.91	60.94	31.64	16.98
	4	5.69	8.77	7.42	8.41	7.91	52.99	126.29	64.49	34.03	18.11
	8	6.17	6.50	3.83	3.55	5.18	58.04	98.80	36.31	15.52	15.76
SARIMA	1	4.88	7.19	7.49	7.91	8.41	45.59	112.28	59.09	36.88	20.14
	2	5.51	8.07	9.22	8.91	9.33	50.98	125.23	96.80	41.30	22.07
	4	5.11	8.77	10.71	10.42	10.99	50.88	132.05	104.98	45.86	24.30
	8	5.96	9.01	10.85	4.89	10.48	52.15	141.89	122.05	17.28	24.21
ETS	1	5.22	6.88	6.52	6.33	8.14	45.92	109.60	56.17	28.95	18.66
	2	5.55	8.17	7.02	7.05	7.50	49.15	128.15	63.40	32.14	16.59
	4	4.97	8.73	7.79	8.85	8.42	46.86	138.20	73.26	37.63	18.64
	8	5.73	9.11	9.62	8.79	6.04	54.53	146.68	89.77	39.18	14.87
SNAIVE	1	11.76	12.22	11.68	12.85	0.07	134.41	232.17	127.27	62.29	35.69
	2	12.66	13.34	12.22	13.53	15.35	140.32	242.49	132.27	64.80	37.07
	4	13.94	15.34	13.18	15.61	0.07	152.87	266.35	143.51	71.32	11.61
	8	14.78	15.58	14.85	16.29	20.42	179.89	298.91	170.29	80.71	46.56

summarizes the forecast accuracy measures for all age groups at various forecast horizons. The MAPE and RMSE values for ITSA were consistently lower than those for SARIMA without intervention, SNAIVE, and ETS for the 50–64-year-old age group.

The DM test results for age groups provided in

Table 6. Results of the modified Diebold–Mariano test comparing ITSA to ETS, SARIMA without intervention, and SNAIVE by age group.

Age Group	h	ETS	SARIMA	SNAIVE
65+	1	−0.67	−0.78	−1.40
	4	−0.37	−0.31	−1.87
50−64	1	−0.63	−1.05	−1.59
	4	−0.37	−0.31	−1.87
35−49	1	0.00	−1.12	−1.46
	4	−0.26	−1.73	−1.79
18−34	1	0.07	−1.04	−1.70
	4	−0.20	−0.90	−1.58
≤17	1	−0.64	−1.70	1.82
	4	−0.35	−2.06	−2.10

Bold indicates a statistical difference between the two forecasting methods at the 10% level of significance. Negative (positive) values indicate that ITSA performed better (worse) than ETS, SARIMA without intervention, or SNAIVE.

show that forecasts using ITSA were in general more accurate than those using the other three methods for almost all the age groups. However, this difference was not very significant. For the age group of ≤17, the DM statistic ($-2.10$) was significant, implying that ITSA outperformed SNAIVE for the four-step-ahead forecasts.

3.4. Ethnicity

This section analyzes donor transplant counts from different ethnic groups. Figure 12 shows the monthly donor counts from White, Black, Asian, Hispanic, Multiracial, American Ind/Alaska Native, Native Hawaiian/other Pacific Islander. According to the analysis by the Institute for Health Metrics and Evaluation [29], Hispanic and Latino Americans were more likely to die from COVID-19 than non-Hispanic White Americans. The elevated risk in Hispanic, Latino, and non-Hispanic Black Americans meant that their risk of death at age 50 was equal to the risk of death for a 65-year-old non-Hispanic White American. This disparity could be attributed to several potential sources, including Hispanic and Latino populations’ increased work exposure during the pandemic [30]. Ethnic minorities on the transplant waitlist can get a better chance to find a good donor match if there is a more diverse donor registry since rare immune system markers match better for people from a similar ethnicity [31]. In 2019, Hispanics made up 18.4% of the national population. The number of organ transplants performed on Hispanics in 2020 was about 30% of the number of Hispanics currently waiting for a transplant. The number of transplants performed on Whites was 48.8% of the number currently waiting. While 20.5% of the total candidates currently waiting for transplants were Hispanic, they comprised 14.6% of organ donors in 2020. Almost 69% of organs recovered from Hispanic patients in 2020 were from deceased donors [32].

A SARIMA$(0,1,1)\times {(0,1,1)}_{12}$ model with the COVID-19 intervention effect in Equation (2) was fit to the Hispanic donor transplant time series from January 2002 to August 2020. Monthly Hispanic donor transplant counts along with the model fit are plotted in the top panel of Figure 13. The bottom panel of Figure 13 shows the COVID-19 effects. The model fit the data quite well and the counts almost rebounded back gradually to their prepandemic levels by the summer of 2020. Specifically, the COVID-19 pandemic caused a sudden drop in Hispanic donor transplant counts by 104.34 in March 2020. This further dropped to −250.93 in April 2020 followed by a fast rebound at the rate of 19% to its pre-COVID level. The estimated ITSA model parameters for all the major ethnic categories focused on in our study are provided in

Table A4. ITSA model parameter estimates for ethnicity.

	Ethnicity
Coefficients	White	Black	Hispanic	Asian
${\theta }_1$	−0.83	−0.99	−0.82	−0.89
${\varphi }_1$		−0.06
${\Theta }_1$	−0.78	−0.90	−0.85	0.14
${\omega }_0$	−386.83	−149.97	−104.34	−27.18
$\delta$	0.22	0.46	0.19	0.41
${\omega }_1$	−768.97	−196.54	−250.93	−31.22

in the Appendix A.

The forecasts of Hispanic donor transplants from September 2022 to May 2021 generated using ITSA, SARIMA without intervention, ETS(M,Ad,A), and SNAIVE are shown in Figure 14. The overall forecasting ability of the ITSA model was good. ETS provided more accurate forecasts than the other three methods for the Hispanic donor transplant series. The MAPE and RMSE values of the major ethnic groups are reported in

Table 7. Forecast accuracy measures by ethnicity at various forecast horizons.

	h	MAPE (%)				RMSE
		Hispanic	White	African American	Asian	Hispanic	White	African American	Asian
ITSA	1	12.05	4.86	6.82	7.02	74.91	130.51	43.53	40.39
	2	14.11	6.02	9.33	7.34	84.40	151.80	53.46	39.54
	4	16.41	6.47	15.24	8.51	89.45	151.35	76.45	41.60
	8	6.96	5.19	14.05	3.07	122.36	135.65	69.42	17.23
SARIMA	1	13.69	5.08	7.73	7.64	87.82	171.76	48.80	43.04
	2	16.01	6.50	10.20	7.35	100.48	196.25	59.48	38.34
	4	19.56	6.72	16.18	5.54	109.03	203.83	80.63	31.65
	8	16.69	12.21	14.87	5.37	122.36	314.61	74.62	31.27
ETS	1	9.67	5.29	7.22	7.22	70.14	158.23	46.69	46.69
	2	12.10	5.78	6.72	6.72	80.78	179.67	40.45	40.45
	4	13.66	5.66	5.46	5.46	84.57	188.96	30.03	30.03
	8	14.65	11.49	2.85	2.85	106.19	285.78	20.09	20.09
SNAIVE	1	21.27	12.74	12.92	12.92	154.02	413.60	79.90	79.90
	2	22.54	14.19	12.02	12.02	162.16	438.60	75.12	75.12
	4	28.75	17.05	11.22	11.22	186.80	500.79	75.23	75.23
	8	33.69	27.99	20.51	20.51	259.75	738.082	106.04	106.04

for various forecast horizons. These measures were lower for the White population when using ITSA.

Table 8. Results of the modified Diebold–Mariano test comparing ITSA to ETS, SARIMA without intervention, and SNAIVE by ethnicity.

Ethnicity	h	ETS	SARIMA	SNAIVE
Hispanic	1	0.60	−1.49	−1.49
	4	0.18	−0.88	−1.07
White	1	−1.12	−1.41	−1.56
	4	-0.46	-0.56	-1.26
African American	1	−0.29	−1.26	−1.58
	4	1.13	−1.64	0.02
Asian	1	−0.83	−0.63	−1.97
	4	1.30	0.94	−0.84

Bold indicates a statistical difference between the two forecasting methods at the 10% level of significance. Negative (positive) values indicate that ITSA performed better (worse) than ETS, SARIMA without intervention, or SNAIVE.

provides the results of the modified Diebold–Mariano test for equality of forecast accuracy comparing ITSA to ETS, SARIMA without intervention, and SNAIVE across forecast horizons one and four by ethnicity. In general, ITSA generated better forecasts across the White, African American, and Asian population, when compared to the other three methods for the one-step-ahead forecasts as the test statistics were negative. In particular, ITSA forecasts were significantly better compared to SNAIVE for the Asian population. Other than that, there seemed to be no significant difference among all four methods in forecasting donor transplant counts by ethnicity.

3.5. Gender

Gender disparities have always been prevalent in healthcare and continue to exist even in organ donation. Although organs in and of themselves are gender neutral and can be exchanged between the sexes, women account for up to two-thirds of all organ donations. Ref. [33] found no clear reasons why women were more willing to risk undergoing surgery than men and that their gender disparity was not mirrored in the demand for donated organs. More men than women are recipients, and women are less likely to complete the necessary steps to receive donated organs. Many research studies into diseases and treatments are skewed with a higher number of male participants, and conducting more research involving women can ensure women are properly represented in medical research. This will help alleviate gender bias [34]. Ref. [35] studied the impact of sex mismatch on transplant outcome and stated that it still remained a matter of debate. Their research found that sex and gender interaction could affect graft success and recommended that they should be taken into account in the management of organ-transplanted patients. In a national survey, both sexes had similar reasons for becoming and not becoming organ donor. However, compared with men, women were more willing to donate their organs to family members and strangers. Systematic changes and strategies to build trust in the medical system were also identified as means to change the stance of non-organ-donors from both sexes [36]. Figure 15 depicts donor transplant counts across male and female genders from January 2002 to May 2021 and the focus here was the female gender.

We fit a SARIMA$(0,1,1)\times {(1,1,0)}_{12}$ with COVID-19 intervention in Equation (2) to the female donor transplant series from January 2002 to August 2020. The fitted donor counts and the estimated COVID-19 effects are plotted in Figure 16. As the above, the COVID-19 pandemic caused a sudden change in the female donor transplant counts with a magnitude of −208.36 in March 2020. This further dropped to a magnitude of 784.12 in April 2020 followed by a rebound at the rate of 30% to its pre-COVID level with no permanent effect. In addition, the ITSA model parameter estimates for both males and females are provided in

Table A5. ITSA model parameter estimates for gender.

	Gender
Coefficients	Female	Male
${\varphi }_1$		−0.16
${\theta }_1$	−0.87	−0.79
${\Phi }_1$	−0.46
${\Theta }_1$		−0.74
${\Theta }_2$		−0.02
${\omega }_0$	−208.36	−326.80
$\delta$	0.30	0.02
${\omega }_1$	−784.12	−613.78

in the Appendix A.

The forecasts of female donor counts from September 2020 to May 2021 are plotted in Figure 17 using ITSA, SARIMA without intervention, ETS(M,N,M), and SNAIVE. The forecasts using SARIMA without intervention were more accurate for the first six months and ITSA gave better forecasts for the latter months in 2021 for female donor counts. The forecast accuracy measures are reported in

Table 9. Forecast accuracy measures by gender at various forecast horizons.

	h	MAPE (%)		RMSE
		Male	Female	Male	Female
ITSA	1	6.97	6.00	151.43	96.07
	2	7.57	7.20	163.41	110.79
	4	8.46	7.16	220.10	177.23
	8	2.41	2.77	53.03	45.15
SARIMA	1	6.67	9.31	158.27	217.90
	2	8.59	9.31	184.78	217.90
	4	9.25	13.91	189.85	289.89
	8	9.84	26.74	227.79	472.06
ETS	1	5.81	6.53	137.09	119.81
	2	7.29	6.49	160.64	127.23
	4	7.18	8.10	153.22	142.24
	8	8.57	14.41	185.94	237.84
SNAIVE	1	13.91	13.58	347.11	310.31
	2	15.64	14.84	368.17	328.66
	4	17.46	18.98	407.49	378.96
	8	23.72	34.09	554.77	593.45

for various forecast horizons. Both MAPE and RMSE values were generally lower for the ITSA model than SARIMA without intervention and SNAIVE. ETS seemed to perform equally to ITSA.

The DM test results for gender provided in

Table 10. Results of the modified Diebold–Mariano test comparing ITSA to ETS, SARIMA without intervention, and SNAIVE by gender.

Gender	h	ETS	SARIMA	SNAIVE
Female	1	−0.72	−1.20	−1.27
	4	−0.26	−0.85	−0.98
Male	1	0.53	−0.30	−1.73
	4	0.22	−0.15	−1.15

show that forecasts using ITSA are normally more accurate than those using SARIMA without intervention and SNAIVE for both males and females based on the negative values of the test statistic. However, this difference was not statistically significant.

4. Conclusions and Future Work

The occurrence of the COVID-19 pandemic provided an important research opportunity to investigate its nature and dynamics and our findings could help public health professionals and organizations mitigate risk when unanticipated situations arise in the future. The effects of COVID-19 on donor transplants could be well modeled with an intervention time series analysis. For all donor transplant series in the US, the significant drop caused by the COVID-19 pandemic in March 2020 is less than that in April 2020. This is because the OPTN Executive Committee issued a COVID-19 emergency policy that came into effect starting 17 March 2020 and donor transplant counts were affected thereafter.

The ITSA model forecast better when compared to the general SARIMA model without interventions and SNAIVE and performed similarly to ETS. The Diebold–Mariano tests did not show as many significant differences among the forecast methods compared in the study as we had expected. Since the sample sizes in each of the DM tests were small, the DM tests need to be used with caution.

This study can easily be extended to forecast waitlist registrations and organ transplant counts with different intervention events. It can also be utilized by policymakers to examine the impact of policies and practices to improve equity metrics in the organ transplant system.

Note that there are many other factors across the individual biological/healthcare system/community that contribute to inequities among donor transplant candidates. They were not addressed since the scope of this study spanned only across individual sociocultural factors. In a future study, the combined effects of all factors mentioned above on donor transplant count time series can be evaluated in a mixed multi-way time series analysis of covariance (ANCOVA) with intervention effects. For example, if we want to study the combined effects of gender and ethnicity on the living kidney donor count time series, a two-way time series ANCOVA model with an intervention effect can be developed as follows:

$$Y_{ijt}={\mu }_t+m_t+{\alpha }_{it}+{\beta }_{jt}+{ϵ}_{ijt},$$

(16)

where $Y_{ijt}$ is the observed donor count at time t for the ith level of gender and jth level of ethnicity. ${\mu }_t$ includes trend and/or seasonality components, $m_t$ is the intervention effect, ${\alpha }_{it}$ is the gender effect, and ${\beta }_{jt}$ is the ethnicity effect. The error term $\left\{{ϵ}_{ijt}\right\}$ is assumed to be mean zero stationary time series. This model format is similar to the one in [37].

The parameter estimation and forecasts for the proposed time series ANCOVA model are quite challenging. Recently, a hybrid model, which combines both statistical and machine learning components has been gaining popularity in time series forecasting as it outperforms either pure traditional statistical methods or pure machine learning methods [38,39]. The hybrid method of the proposed time series ANCOVA model and deep neural networks seems an appealing approach to time series forecasting.