A Semiparametric Model

We propose a semiparametric Poisson regression model (that of Wood 2017) ^{Reference Wood22} to estimate whether some deaths due to some causes increased after the passage of Hurricane Maria in Puerto Rico while controlling for the temporal variation in the mortality rate of the causes of death and the population size. Population displacement plays an important role, which can be directly influenced by economic, social, labor, academic, climatic variations, and other factors. ^{Reference Neyra23,Reference Acosta, Kishore and Irizarry24}

Additionally, index functions will be deployed in the model to evaluate if there are effects of changes in the behavior of certain leading causes of death, including some non-communicable diseases (NCDs) in different post-emergency periods.

Let D _t = number of deaths at time index t, N _t = population size at time t. We assume that D _t follows a Poisson distribution. For $i = 0, \ldots, I,$ we assign ${P_{i,t}}$ as an indicator of the period i for time t. Specifically, $i = 0$ represents the pre-emergency period, $i = 1$ is a given period after the emergency, and so on. For $j = 1, \ldots, J$ , we consider ${c_{j,t}}$ as an indicator of cause of death j at time t. On the other hand, we define $da{y_t} = $ day of the year, and $yea{r_t} = $ a categorical variable that captures the mortality trend through the years. We define the following generalized additive semiparametric Poisson regression model with interaction, ^{Reference Ruppert, Wand and Carroll25}

$$\begin{gathered}\log \left( {{\mu _{i,j,t}}} \right) =\log \left( {{N_t}} \right) + {\beta _0} + {{\rm year}_t} + {\boldsymbol P}{'_{}}{\boldsymbol\alpha} + {\boldsymbol c}{'_t}\gamma + {\left( {{{\boldsymbol P}_t} \otimes {{\boldsymbol c}_t}} \right)^{'}}{\boldsymbol\omega} \\ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!+f\left( {da{y_t}} \right) + {\varepsilon _{i,j,t}},\end{gathered}$$

where ${\mu _{i,j,t}} = E\left( {{D_t}{\rm{|}}t,{\rm{\;}}{{\boldsymbol P'}_t},{\rm{\;}}{{\boldsymbol c'}_t},{\rm{\;}}{N_t},{\rm{\;}}da{y_t},{\rm{\;}}yea{r_t}} \right).$ We define ${{\rm{\boldsymbol P'}}_t} = \left( {{P_{0,t}}, \ldots, {P_{I,t}}} \right)$ and ${\boldsymbol \alpha} = \left( {{\alpha _0}, \ldots, {\alpha _I}} \right){\rm{'}}$ denotes the coefficients of the model for the I pre and post-emergency periods. Similarly, ${{{\boldsymbol c'}}_t} = \left( {{c_{1,t}}, \ldots, {c_{J,t}}} \right)$ and $\gamma = \left( {{\gamma _1},{\gamma _2} \ldots, {\gamma _J}} \right){\rm{'}}$ denotes the coefficients of the model for the J causes of death. While $\left( {\boldsymbol P}_t \otimes {\boldsymbol c}_t \right)^{\rm{'}} = \left( P_{0,t}{c}_{1,t}, \ldots, {P}_{I,t}{ c}_{J,t} \right)$ is the Kronecker product between vectors ${{\rm{\boldsymbol P'}}_t}$ and ${{\rm{\boldsymbol c'}}_t}$ . ${P_{i,t}}{c_{j,t}}$ represents the interaction between the I pre- or post-emergency periods and the J causes of death at time t and $\omega = \left( {{\omega _{11}}, \ldots, {\omega _{MQ}}} \right){\rm{'}}$ represents the coefficients of the interaction. On the other hand, the smooth function $f\left( {da{y_t}} \right)$ captures within-year variation in deaths while ${\beta _0}$ is the contrast of the smooth function and the explanatory variables, and ${\varepsilon _{i,j,t}}$ contains the variability not explained by the other terms. For post-emergency periods, ${P_{i,t}} = 1$ and the interaction coefficient make it possible to adapt the effect of the post-emergency period for cause j. The coefficients of the model (1) are estimated by penalized likelihood maximization, considering the population displacement caused by a state of emergency while the smoothing penalty parameters are determined by the restricted maximum likelihood (REML). Model (1) can be used to estimate excess death of cause of death j during period i post-emergency of index time t through the difference between the estimation of the model with ${P_{i,t}} = 1$ , ${c_{j,t}} = 1$ , and ${P_{i,t}}{c_{j,t}} = 1$ versus the estimated model with ${P_{i,t}} = 0$ , ${c_{j,t}} = 1$ , and ${P_{i,t}}{c_{j,t}} = 0$ . For $i \ge 1$ and $j \ge 1$ , we obtain:

$${\hat \varphi _{i,j,t}} = {\hat E}\left( {{D_t}{\rm{|}}t,{P_{i,t}} = 1,\;{c_{j,t}} = 1,\;{P_{i,t}}{c_{j,t}} = 1,\;N_t^*,\;da{y_t},\;yea{r_t}} \right),$$

$${\hat \chi _{i,j,t}} = \hat E\left( {{D_t}{\rm{|}}t,{P_{i,t}} = 0,\;{c_{j,t}} = 1,\;{P_{i,t}}{c_{j,t}} = 0,\;N_t^*,\;da{y_t},\;yea{r_t}} \right).$$

${\hat b_0}$ , ${\hat a_i}$ , ${\hat{d}_j}$ , and ${\hat e_{i,j}}$ , estimate ${\beta _0}$ , ${\alpha _i}$ , ${\gamma _j}$ , and ω _ij , respectively. Also, notice that ${\hat \varphi _{i,j,t}}$ uses $N_t^*$ , a population size unaltered by post-emergency displacement because ${\hat a_i}$ and ${\hat e_{i,j}}$ are already considering the impact of population displacement. Using it again for ${\hat \varphi _{i,j,t}}$ would suggest accounting for the effect twice. This way,

(1) $$\eqalign{& & {{\hat \varphi} _{i,j,t}- {{\hat \chi }}_{i,j,t} = \exp (\log (N_t^*) + \hat b_0 + {{\hat a}_i} + \hat{d}_j} + {{\hat e}_{ij}} + \hat f({\rm{da}}{{\rm{y}}_t}) + \widehat {year}_t) \cr & & \hskip 8pt- \exp (\log (N_t^*) + {{\hat b}_0} + {\hat{d}_j}+ \hat f({\rm{da}}{{\rm{y}}_t})+ \widehat {year}_t) \cr & & \hskip 8pt = N_t^*{e^{\hat b_0+ \hat{d}_j + \hat f(da{y_t}) + \widehat {year}_t}}(e^{\hat a_i + \hat e_{ij}} - 1). \cr}$$

The expression ${e^{{{\hat b}_0} + \hat{d}_{j} + \hat f\left( {da{y_t}} \right) + {{\widehat {year}}_t}}}N_t^*$ in Equation (1) represents the number of deaths typically seen, and $\left( {{e^{{{\hat a}_i} + {{\hat e}_{ij}}}} - 1} \right)$ represents percentage changes during the post-emergency periods adjusting for the interaction between the J causes of death and the post-emergency periods. Equation (1) is the maximum likelihood estimator for the excess deaths expected from the J causes of death in post-emergency period i. Using Equation (1), the accumulated excess of the J causes of death is,

(2) $$\mathop \sum \limits_{t = q}^r \left( {{{\hat \varphi }_{i,j,t}} - {{\hat \chi }_{i,j,t}}} \right)$$

for any period that begins at the index q and ends at r.

In general terms, model (1) allows statistical inference of the duration of the effect in the I post-emergency periods, the effect on the J causes of death during a state of emergency, and the effect caused by the interaction between the I post-emergency periods and the J causes of death in the population displacement.

Data

Our data correspond to daily deaths reported by the Puerto Rico Demographic and Vital Statistics Registry between January 1, 2014, to December 31, 2018, and includes a total of 148,490 records. This work is exempt from requiring human subjects board review. Death certificates are processed using Automatic Classification of Medical Entry (ACME) software. The input to ACME is the multiple cause-of-death International Classification of Diseases (ICD-10) codes assigned to each entity (eg, disease condition, accident, or injury) listed on death certificates. ACME then applies the World Health Organization (WHO) rules and guidelines to the ICD codes and selects an underlying cause of death. ^{Reference Lu26}