Bayesian Dynamic Generalized Additive Model for Mortality during COVID-19 Pandemic
Wei Zhang Antonietta Mira Ernst C. Wit

Nodeless superconductivity and topological nodal states in molybdenum carbide
    Toni Shiroka
    
=============================================================================





	


§ ABSTRACT 
While COVID-19 has resulted in a significant increase in global mortality rates, the impact of the pandemic on mortality from other causes remains uncertain. To gain insight into the broader effects of COVID-19 on various causes of death, we analyze an Italian dataset that includes monthly mortality counts for different causes from January 2015 to December 2020. Our approach involves a generalized additive model enhanced with correlated random effects. The generalized additive model component effectively captures non-linear relationships between various covariates and mortality rates, while the random effects are multivariate time series observations recorded in various locations, and they embody information on the dependence structure present among geographical locations and different causes of mortality. Adopting a Bayesian framework, we impose suitable priors on the model parameters. For efficient posterior computation, we employ variational inference, specifically for fixed effect coefficients and random effects, Gaussian variational approximation is assumed, which streamlines the analysis process. The optimisation is performed using a coordinate ascent variational inference algorithm and several computational strategies are implemented along the way to address the issues arising from the high dimensional nature of the data, providing accelerated and stabilised parameter estimation and statistical inference.


Keywords: generalized additive model, state space model, variational inference 



§ INTRODUCTION



The COVID-19 pandemic has had profound impacts across the globe. Researches focus on various aspects such as health disparities linked to racial and socio-economic factors and the healthcare system's adaptation in terms of testing, contact tracing, and vaccine rollouts <cit.>. A key area of investigation is excess mortality, which provides an overarching view of the pandemic's impact on human health. This includes factors such as government lockdown measures and disruptions to non-COVID healthcare services <cit.>. While excess mortality offers a general perspective, examining the consequences of specific causes of death as a result of these factors during the pandemic is crucial for developing more targeted future mitigation strategies. For example, the pandemic has contributed to increases in deaths from chronic conditions as observed by <cit.>. There also has been a rise in accidental deaths, homicides or suicides <cit.>. Contributing to this research field, we analyze Italian monthly death counts recorded in 21 Italian regions from 2015 to 2020, categorized according to the International Classification of Diseases, 10th Revision (ICD-10) <cit.>, to understand the effects on specific human mortality during the pandemic better.

We apply the generalized additive model (GAM) <cit.> to study non-linear relationships between response, cause-specific mortality rate, and continuous covariates, including lockdown intensity level and age, in terms of smooth functions. GAMs have demonstrated considerable potential to study COVID-19 mortality <cit.>. In the more general context of spatial-temporal data modeling, GAMs are well-suited for analyzing changes over time and across different geographic regions, crucial in understanding the differentiating consequences. One way of employing GAMs to model the change of spatial pattern over time is via a tensor product smoother <cit.>. Alternatively, <cit.> introduce dynamic GAMs where extra dynamic spatial random effects are incorporated into the mean of response variable, offering a solution to forecasting discrete time series while estimating relevant nonlinear predictor associations that conventional generalized linear models (GLM) plus spatial temporal random effects are not able to account for <cit.>. In our study, we follow <cit.>'s approach and assume that the mortality rate of specific cause in a region for a certain month combines a GAM regression component and random effects. The random effects are further assumed to be correlated as opposed to being independent, facilitating statistical inference on dependence structure between geographical regions and various causes of death <cit.>. We adopt a Bayesian approach and assign suitable priors on model parameters to better infer correlation structure from the data.


The model framework is a special case of non-Gaussian state space models where the observation equation is Poisson distribution. Inference with state space models includes filtering to estimate the current state given past observations and smoothing to estimate past states given all observations up to the current time <cit.>. Filtering and smoothing is exact in linear Gaussian state space models where Kalman Filter can be applied. When the model is non-Gaussian, particle filter approximates the states using a set of weighted particles <cit.>. The posterior samples in <cit.> are drawn either in the Gibbs sampling software JAGS <cit.> or with the Hamiltonian Monte Carlo in Stan <cit.>. However, the approximation for the target deteriorates as the dimension increases <cit.>. Due to the high dimensionality inherent in the data and model, these sampling algorithms are impractical. To address the issue, we use a variational inference algorithm for fast approximation <cit.>; more specifically, the joint variational density of fixed effect coefficients and random effects takes the Gaussian form. Various methods exist for optimizing the mean and variance of this approximating Gaussian distribution. For instance, <cit.> propose to parameterize the density in terms of its mean and a lower triangular scale matrix whereas <cit.> incorporate sparsity in the Cholesky factors of the precision matrix. Both developed stochastic gradient methods for optimization. Our approach instead involves Newton's method and fixed point iteration as in <cit.> to achieve faster convergence in fewer iterations. This aligns with recent advancements with recent development in Gaussian variational approximations for high-dimensional state space models <cit.>, but distincts in that we do not define structure of the proposed variational approximation thanks to our chosen optimization technique.
  
The rest of the paper is organized as follows. In Section <ref>, we formulate the model and specify the priors imposed on parameters. In Section <ref>, we demonstrate how to derive the ELBO in this model setup and present the variational algorithm for posterior inferences. We then apply the model and the algorithm to the Italian monthly mortality data in Section <ref>. Finally, Section <ref> gives some concluding remarks and points to future work.



§ BAYESIAN DYNAMIC GAM



Let Y_n,t be the mortality count of instance n at time t, n=1,…,N, t=1,…,T. For each n, covariates such as ordinal a_n for age categorical g_n for gender, l_n for L geographical locations and k_n for K different causes of death are available. Additionally, we are interested in the lockdown effect on outcome variable Y_n,t, therefore we include a stringency index that quantifies the intensity of government restriction policies for each geographical location over time and we denote it by r_n,t, which depends on n through l_n. We assume that Y_n,t is Poisson random variable with rate equal to

    ϵ_n,texp[𝐱^β_nβ + f^r(r_n,t) + f^a(a_n) +  f^k,r(k_n, r_n,t) + f^k,a(k_n, a_n) + f^g,a(g_n, a_n) + z^*_n,t]

where ϵ_n,t is the offset, 𝐱^β_nβ are parametric terms which include gender effect. f^r(r_n,t) and f^a(a_n) are natural cubic splines for government intervention effect and age effect respectively and they take the form 

    f^r(r) = ∑_j=1^J u^r_j(r)β^r_j,    f^a(a) = ∑_j=1^J u^a_j(a)β^a_j,

where J is the number of knots, u^r_j, u^a_j are natural cubic spline bases, β^r_j and β^a_j are coefficients to be estimated. The smoothness penalization terms associated with the bases are defined as

    λ^r_1(β^r)' S^r_1β^r + λ^r_2(β^r)' S^r_2β^r,   λ^a_1(β^a)' S^a_1β^a + λ^a_2(β^a)' S^a_2β^a.

Here S^r_1, S^r_1, S^a_1, S^a_2 contain known coefficients. β^r=(β^r_1, …, β^r_J)' and β^a=(β^a_1, …, β^a_J)' are vectors of spline coefficients. From a Bayesian perspective, this is equivalent to imposing the following multivariate normal priors on β_r and β_a

    β^r ∼𝒩(0, (λ^r_1S^r_1+λ^r_2 S^r_2)^-1) 
    β^a ∼𝒩(0, (λ^a_1 S^a_1+λ^a_2 S^a_2)^-1).

λ^r_1, λ^r_2, λ^a_1 and λ^a_2 control the level of roughness.

As for the two-way interaction terms, f^k,r(k, r) models interactions between cause of death and stringency index in a non-parametric manner. f^k,a(k, a) accounts for interactions between causes of death and age while f^g,a(g, a) captures interactions between gender and age. The three interactions terms are constructed in the following way. For each level k of causes of death or each g of gender, a unique spline is specified <cit.>. We formulate f^k,r(k, r) as an example,  f^k,a(k, a) and f^g,a(g, a) are modeled in a similar way. For each k other than the baseline, f^k,r(k, r) is assumed to be natural cubic spline such that

    f^k,r(k, r) = ∑_j=1^J u^k,r_k,j(r)β^k,r_k,j.

The J(K-1) dimensional vector β^k,r=(β^k,r_2,1,…,β^k,r_2,J,β^k,r_3,1…,β^k,r_K,J)' is jointly penalised by S^k,r_1 and S^k,r_2 with smoothing parameters λ^k,r_1 and λ^k,r_2 and the penalisation term is λ^k,r_1(β^k,r)' S^r_1β^k,r + λ^k,r_2(β^k,r)' S^r_2β^k,r,
which, in terms of Bayesian prior, translates to

    β^k,r∼𝒩(0, (λ^k,r_1S^k,r_1+λ^k,r_2 S^k,r_2)^-1).


The last term z^*_n,t embeds the spatial-temporal structure in the data; in fact, z^*_n,t depends on n through l_n and k_n, therefore in total, we have LK time series of length T. z^*_n,t can be modeled with great flexibility. We make the following latent state assumption. Let 𝐳^*_t=(z^*_1,t, …, z^*_N,t)' and 𝐳_t=(z_1,t, …, z_LK,t)' such that

    𝐳^*_t = (I_LK⊗1)𝐳_t,

where 1 stands for a vector whose elements are all equal to 1. The dimension of 1 is determined by the number of age groups and gender. The assumption implies that residual mortality rates z^*_n,t of the same causes of death in the same region are identical regardless of age and gender, which is reasonable since we have already taken into account age and gender effect in GAM component. We further assume that marginally the latent state 𝐳_t∼𝒩(μ,Σ) and P(𝐳_t - μ) follows a simple autoregressive (AR) process such that 

    P(𝐳_t - μ) = Φ P(𝐳_t-1 - μ) + ϵ_t,   ϵ_t∼𝒩(0, I_LK-ΦΦ'),

where P is an upper triangular matrix such that P'P=Ω=Σ^-1, i.e. P are the Cholesky factors of the precision matrix Ω.
Φ contains multivariate autoregressive coefficients. For simplicity, we require Φ to be diagonal with diagonal entries ϕ=c(ϕ_1, …, ϕ_LK)' such that -1<ϕ_1, …, ϕ_LK<1. The pre-multiplying of 𝐳_t by P is to ensure that marginally 𝐳_t ∼𝒩(0, Σ), t=0, …, T where Σ incorporates dependence structure of geographical locations as well as mortality causes. Reorganizing 𝐳_t as matrices Z_t of size L× K, we assume that 

    Z_t ∼MN(μ, Σ^k, Σ^l),

which is a matrix normal distribution equivalent to 

    𝐳_t ∼𝒩(μ, Σ^k ⊗Σ^l),

where ⊗ stands for the Kronecker product between Σ^k and Σ^l, the covariance matrices of the causes of death dependence and regional dependence respectively. We impose the following priors on the model parameters. Firstly, μ is assigned a normal prior 𝒩(0, σ^2_μ I_LK). Ω^k = (Σ^k)^-1 is the precision matrix that underlines the conditional independence of mortality causes. For simplicity, we stick to the usual Wishart distribution

    Ω^k ∼Wishart(δ^k, θ^k I_K).





The same reasoning applies to the precision matrix Ω^l=(Σ^l)^-1 and it follows that

    Ω^l ∼Wishart(δ^l, θ^l I_L).









We complete the model with a hierarchy of prior specification on the linear regression coefficients, on penalization parameters that control smoothness of splines and on the autoregressive coefficients, on the cause specific mean. Denoting by 𝐱_n the vector containing 𝐱^β_n and all bases in GAMs for n, the hierarchical model can be summarized as

    Y_n,t∼Poisson[ ϵ_n,texp(𝐱_nβ^* + z^*_n,t)],   β^* = (β', (β^r)', (β^a)', (β^k,r)', (β^k,a)', (β^g,a)')', 
    β∼𝒩(0, σ^2_β I),   β^r ∼𝒩(0, (λ^r_1S^r_1+λ^r_2 S^r_2)^-1),   β^a ∼𝒩(0, (λ^a_1 S^a_1+λ^a_2 S^a_2)^-1), 
    β^k,r∼𝒩(0, (λ^k,r_1S^k,r_1+λ^k,r_2 S^k,r_2)^-1),   β^k,a∼𝒩(0, (λ^k,a_1S^k,a_1+λ^k,a_2 S^k,a_2)^-1), 
    β^g,a∼𝒩(0, (λ^g,a_1S^g,a_1+λ^g,a_2 S^g,a_2)^-1), 
    λ^r_1, λ^r_2, λ^a_1, λ^a_2, λ^k,r_1, λ^k,r_2, λ^k,a_1, λ^k,a_2, λ^g,a_1, λ^g,a_2 ∼Gamma(α_λ, β_λ), 
    𝐳^*_t = (I_LK⊗1)𝐳_t, 
    
           P(𝐳_t-μ) = Φ P(𝐳_t-1-μ) + ϵ_t,   ϵ_t∼𝒩(0, I_LK-ΦΦ'),   μ∼𝒩(0,σ^2_μ I_LK), 
    Φ = diag(ϕ_1,…,ϕ_LK),   ϕ_1+1/2, …, ϕ_LK+1/2i.i,d∼Beta(α_ϕ, β_ϕ), 
       P'P = Ω = Ω^k ⊗Ω^l,  
               
            Ω^k ∼Wishart(δ^k, θ^k I_K),   Ω^l ∼Wishart(δ^l, θ^lI_L).

























§ POSTERIOR INFERENCE VIA VARIATIONAL APPROXIMATION



To make posterior inference, one may devise suitable Markov Chain Monte Carlo (MCMC) algorithms to obtain posterior samples. However, the algorithm may fail to deliver desirable convergent output within reasonable time when it is difficult to explore the geometry of the target distribution due to the sheer dimension of the data. Therefore we resort to variational inference approach for fast approximation. The target posterior distribution is 

    p(.    β^*,λ,𝐳_0:T, .ϕ, μ, Ω^k, Ω^l |𝐲) 
       ∝ p(𝐲|β^*, 𝐳_0,T)p(β^* |λ)p(λ)p(𝐳_0:T|ϕ, μ, Ω^k, Ω^l)p(μ)p(ϕ)p(Ω^k)p(Ω^l),

and the goal is to find the optimal q^* (β^*,λ,𝐳_0:T, μ,ϕ, Ω^k, Ω^l ) from a pre-specified family of distributions such that the Kullback–Leibler (KL) divergence, defined as

    KL[q (·) || p(·|𝐲)] = E_q[logq(·)/p(·|𝐲)],

with expectation taken with respect to q(·), is minimized. This is equivalent to maximizing the evidence lower bound (ELBO) given by

    !ELBO[ p(·|𝐲)] 
    
        =   E_q [logp(𝐲|β^*, 𝐳_0,T)p(β^* |λ)p(𝐳_0:T|ϕ, μ, Ω^k, Ω^l)p(λ)p(μ)p(ϕ)p(Ω^k)p(Ω^l)/q (β^*,𝐳_0:T, λ,μ,ϕ, Ω^k, Ω^l )].

We further assume that the family of candidate  approximation q(·) can be factorized as

    q (β^*,𝐳_0:T, λ,μ,ϕ, Ω^k, Ω^l ) = q (β^*,𝐳_0:T)q(λ)q(μ)q(ϕ)q(Ω^k)q(Ω^l),

where q (β^*,𝐳_0:T) is multivariate normal distribution with mean 𝐦 and covariance matrix M. This is essentially the variational Gaussian approximation (VGA) that has been widely implemented in literature and its theoretical properties when applied to Poisson data are studied 
by <cit.>. M can be full or it can be sparse block diagonal matrix so that there is no correlation between β^* and 𝐳_0:T, which means that q (β^*,𝐳_0:T) can be further factorized as the product of two multivariate normal densities, q (β^*) and q (𝐳_0:T). The sparse version of M greatly reduces algorithm complexity. However, in our application, since the dimension of 𝐳_0:T is overwhelmingly larger than β^*, the block diagonal assumption does not produce much gain; therefore we stick to the full matrix characterisation of M. Furthermore, we choose the following variational densities implied by mean field approximation: for each element in λ, it is a point mass at λ^q,r_1, λ^q,r_2, λ^q,a_1, λ^q,a_2, λ^q,k,r_1, λ^q,k,r_2, λ^q,k,a_1, λ^q,k,a_2, λ^q,g,a_1, λ^q,g,a_2. Together they are λ^q. The Dirac measure choice as a variational density avoids evaluating expectation with respect to otherwise non-trivial variational densities. For q(μ), we assume independent normal approximation densities q(μ)=𝒩(μ^q, diag[(σ^q)^2]). For q(ϕ), even though it is analytically possible to compute the expectation when assuming that (ϕ_1+1)/2i.i,d∼Beta(α^q,ϕ_1, β^q,ϕ_1), …, (ϕ_LK+1)/2 i.i,d∼Beta(α^q,ϕ_LK, β^q,ϕ_LK), the optimization scheme we adapt diverges as only the means of beta variational distributions is identifiable, therefore we use Dirac measures on q(ϕ) so they are ϕ^q. Finally, we set  q(Ω^k) = Wishart(δ^q,k, D^q,k) and q(Ω^l) = Wishart(δ^q,l, D^q,l).



 §.§ ELBO calculation

With (<ref>), the ELBO can be written as the sum of expectations

    ELBO[ p(·|𝐲)] 
    
        =   E_q [logp(𝐲|β^*, 𝐳_0,T)p(β^* |λ)p(𝐳_0:T|ϕ, μ, Ω^k, Ω^l)/q(β^*,𝐳_0:T)] + E_q(λ)[logp(λ))/q(λ))] 
        + E_q(μ) [logp(μ)/q(μ)] + E_q(ϕ) [logp(ϕ)/q(ϕ)] + E_q(Ω^k) [logp(Ω^k)/q(Ω^k)] + E_q(Ω^l ) [logp(Ω^l )/q(Ω^l )].

The first term on the right hand side of the equation bears most importance as it connects likelihood with prior and it is also the most computationally heavy part to optimize due to the dimension of 𝐳_0,T. Denote the mean and covariance matrix of joint multivariate normal prior on β^* and 𝐳_0:T by 𝐦_0 and M_0.

M_0 is block diagonal, whose upper diagonal block corresponds to the covariance of β^* and lower diagonal block is the vector autoregressive matrix that derives from the assumed 𝐳_t dynamics. The lower diagonal block is


    (I_T⊗ P^-1)[  I_LK     Φ   Φ^2     ⋯ Φ^T-1;     Φ  I_LK     Φ     ⋯ Φ^T-2;   Φ^2     Φ  I_LK     ⋯ Φ^T-3;     ⋮     ⋮     ⋮     ⋱     ⋮; Φ^T-1 Φ^T-2 Φ^T-3     ⋯  I_LK;       ](I_T⊗(P^-1)').

It is more convenient to use the precision matrix as the term repeatedly appears in the object function that we aim to optimize. The precision matrix can be expressed as (I_T⊗ P')R(I_T⊗ P) with

    !R= [             (I_LK - Φ^2)^-1           -(I_LK - Φ^2)^-1Φ                           0                           ⋯                           0;           -(I_LK - Φ^2)^-1Φ (I_LK - Φ^2)^-1(I_LK + Φ^2)           -(I_LK - Φ^2)^-1Φ                           ⋯                           0;                           0           -(I_LK - Φ^2)^-1Φ (I_LK - Φ^2)^-1(I_LK + Φ^2)                           ⋯                           0;                           ⋮                           ⋮                           ⋮                           ⋱                           ⋮;                           0                           0                           0                           ⋯             (I_LK - Φ^2)^-1;                             ].











The first terms is hence

    E_q [logp(𝐲|β^*, 𝐳_0,T)p(β^* |λ)p(𝐳_0:T|μ, ϕ, Ω^k, Ω^l)/q(β^*,𝐳_0:T)] 
    
    =    𝐲'X𝐦-∑_i=1^NTϵ_iexp(𝐱'_i𝐦+1/2𝐱'_i M 𝐱_i) -1/2[𝐦-E_q(𝐦_0)]'E_q(M_0^-1)[𝐦-E_q(𝐦_0)] 
        -1/2tr[E_q(M_0^-1)M] + 1/2log|M| - 1/2tr[E_q(M_0^-1))COV_q(𝐦_0)] + 1/2E_q(log|M_0^-1|) 
        + constant

Here, X is a two-block design matrix. The left block matrix consists of all regressors while the right block is the kronecker product I_LKT⊗1. The length of vector 1 relies on specific data configuration; in our case, it equals the product between number of age groups and gender. 𝐱_i is the column vector of i-th row of the design matrix X. The expected value E_q(𝐦_0)=(0', (1_T⊗μ^q)')' and E_q(M_0^-1) is

    !E_q(M_0^-1)=[                            1/σ_β^2I                                   0                                   ⋯                                   0                                   0;                                   0         λ_1^q,r S_1^r+λ_2^q,r S_2^r                                   ⋯                                   0                                   0;                                   ⋮                                   ⋮                                   ⋱                                   0                                   0;                                   0                                   0                                   0 λ_1^q,g,a S_1^g,a+λ_2^q,g,a S_2^g,a                                   0;                                   0                                   0                                   0                                   0        E_q[(I_T⊗ P')E_q(R)(I_T⊗ P)] ],

where 

    E_q( R ) =[ E^1_q( R ) E^2_q( R )          0          ⋯          0; E^2_q( R ) E^3_q( R ) E^2_q( R )          ⋯          0;          0 E^2_q( R ) E^3_q( R )          ⋯          0;          ⋮          ⋮          ⋮          ⋱          ⋮;          0          0          0          ⋯ E^1_q( R ) ],

with diagonal matrices

    E^1_q( R ) = diag(1/1-(ϕ^q)^2),   E^2_q( R ) = diag(-ϕ^q/1-(ϕ^q)^2),   E^3_q( R ) = diag(1+(ϕ^q)^2/1-(ϕ^q)^2).






Here we vectorize the function by writing its arguments in terms of vectors ϕ^q. The expectation is therefore

    !E_q[(I_T⊗ P')E_q(R)(I_T⊗ P)] = [ E_q[P'E^1_q( R )P] E_q[P'E^2_q( R )P]                  0                  ⋯                  0; E_q[P'E^2_q( R )P] E_q[P'E^3_q( R )P] E_q[P'E^2_q( R )P]                  ⋯                  0;                  0 E_q[P'E^2_q( R )P] E_q[P'E^3_q( R )P]                  ⋯                  0;                  ⋮                  ⋮                  ⋮                  ⋱                  ⋮;                  0                  0                  0                  ⋯ E_q[P'E^1_q( R )P] ].

Recall that P'P=Ω = Ω^k ⊗Ω^l, when the variational densities of precision matrices Ω^k and Ω^l are Wishart(δ^q,k, D^q,k) and Wishart(δ^q,l, D^q,l), the derived Cholesky upper triangular matrix P can be written as P = (A^kV^q,k) ⊗(A^lV^q,l) with (V^q,kA^k)'A^kV^q,k = Ω^k, (V^q,lA^l)'A^lV^q,l = Ω^l. V^q,k, V^q,l are Cholesky factors of D^q,k and D^q,l and 

    A^k = [   c^k_1 n^k_1,2 n^k_1,3       ⋯ n^k_1,K;       0   c^k_2 n^k_2,3       ⋯ n^k_2,K;       0       0   c^k_3       ⋯ n^k_3,K;       ⋮       ⋮       ⋮       ⋱       ⋮;       0       0       0       ⋯   c^k_K ],    
    A^l = [   c^l_1 n^l_1,2 n^l_1,3       ⋯ n^l_1,K;       0   c^l_2 n^l_2,3       ⋯ n^l_2,K;       0       0   c^l_3       ⋯ n^l_3,K;       ⋮       ⋮       ⋮       ⋱       ⋮;       0       0       0       ⋯   c^l_K ],

with c^k_i∼χ^2_δ^q,k-i+1, c^l_i∼χ^2_δ^q,l-i+1, n^k_i,j, n^l_i,j∼𝒩(0,1) independently. This is known as the Bartlett decomposition <cit.>. Thanks to the decomposition and diagonal assumption on Φ, we are able to work out the expectation of the quadratic forms. In the following, we demonstrate the details of deriving E_q[P'E^1_q( R )P]. E_q[P'E^2_q( R )P] and E_q[P'E^3_q( R )P] have similar structures.

    E_q[P'E^1_q( R )P]     = E_q[(V^q,k⊗ V^q,l)'(A^k⊗ A^l)'E^1_q( R )(A^k⊗ A^l)(V^q,k⊗ V^q,l)] 
        = (V^q,k⊗ V^q,l)'E_q[(A^k⊗ A^l)'E^1_q( R )(A^k⊗ A^l)](V^q,k⊗ V^q,l).

Now the problem simplifies to take the expectation of (A^k⊗ A^l)'E^1_q( R )(A^k⊗ A^l) with




















    A^k⊗ A^l = [  c^k_1 A^l n^k_1,2A^l n^k_1,3A^l          ⋯ n^k_1,KA^l;          0  c^k_2 A^l n^k_2,3A^l          ⋯ n^k_2,KA^l;          0          0  c^k_3 A^l          ⋯ n^k_3,KA^l;          ⋮          ⋮          ⋮          ⋱          ⋮;          0          0          0          ⋯  c^k_K A^l ].

Since c^k_i,n_i,j^k are independent, the expected value is a block diagonal matrix with diagonal entries E[(c_1^k)^2]E{(A^l)'[E^1_q( R )]_1:LA^l}, E[(n_1,2^k)^2]E{(A^l)'[E^1_q( R )]_1:LA^l} +E[(c_2^k)^2]E{(A^l)'[E^1_q( R )]_(L+1):2LA^l}, …, E[(n_1,K^k)^2]E{(A^l)'[E^1_q( R )]_1:LA^l} +⋯+ 

E[(c_K^k)^2]E{(A^l)'[E^1_q( R )]_((K-1)L+1):KLA^l}










where [E^1_q( R )]_1:L denotes the diagonal blocks indexed from 1 to L and so on. By the same reasoning, E{(A^l)'[E^1_q( R )]_1:LA^l} is a diagonal matrix whose diagonal entries are E[(c_1^l)^2][E^1_q( R )]_1, E[(n_1,2^l)^2][E^1_q( R )]_1+E[(c_2^l)^2][E^1_q( R )]_2, …, E[(n_1,L^l)^2][E^1_q( R )]_1+⋯+E[(c_L^l)^2][E^1_q( R )]_L. The same structure also holds for the remaining matrix expectations. The last two quantities in the expectation E_q [logp(𝐲|β^*, 𝐳_0,T)p(β^* |λ)p(𝐳_0:T|μ, ϕ, Ω^k, Ω^l)/q(β^*,𝐳_0:T)] are

    COV_q(𝐦_0)= [                  0                  0;                  0 I_T ⊗diag[(σ^q)^2] ],

and

    E_q(log|M_0^-1|) =    
        log|λ_1^q,rS_1^r+λ_2^q,rS_2^r| + 
    log|λ_1^q,aS_1^a+λ_2^aS_2^q,a| + log|λ_1^q,k,rS_1^k,r+λ_2^q,k,rS_2^k,r| 
        + log|λ_1^q,k,aS_1^k,a+λ_2^q,k,aS_2^k,a| + log|λ_1^q,g,aS_1^g,a+λ_2^q,g,aS_2^g,a| 
    
       +(T-1)∑_i=1^LKlog(1/1-(ϕ_i^q)^2)+LT[log|D^q,k|+∑_k=1^K ψ(δ^q,k-k+1/2)] 
        + KT[log|D^q,l|+∑_l=1^L ψ(δ^q,l-l+1/2)] + constant,

with ψ(·) representing the digamma function.
The remaining terms in the ELBO are

    E_q(λ) [logp(λ)/q(λ)] =    (α - 1)logλ^q - βλ^q + constant 
    E_q(μ) [logp(μ)/q(μ)] =     -1/2σ_μ^2μ^q'μ^q - 1/2σ_μ^2∑(σ^q)^2
        + 1/2∑log[(σ^q)^2] + constant 
    
        E_q(ϕ) [logp(ϕ)/q(ϕ)] =    ∑_i=1^LK (α_ϕ-1)log(1+ϕ_i^2) + (β_ϕ-1)log(1-ϕ_i^2) + constant 
    E_q(Ω^k) [logp(Ω^k)/q(Ω^k)] =    δ^k/2log|D^q,k| + δ^k-δ^q,k/2∑_i=1^K ψ(δ^q,k-i+1/2) - δ^q,k/2θ^ktr(D^q,k) 
        +δ^q,kK/2+logΓ_K(δ^q,k/2) + constant 
    E_q(Ω^l) [logp(Ω^l)/q(Ω^l)] =    δ^l/2log|D^q,l| + δ^l-δ^q,l/2∑_i=1^L ψ(δ^q,l-i+1/2) - δ^q,l/2θ^ltr(D^q,l) 
        +δ^q,lL/2+logΓ_L(δ^q,l/2) + constant

Now that we have formulated the ELBO, a coordinate ascent variational inference (CAVI) algorithm is devised and employed to maximize the target; that is we optimize the ELBO with respect to 𝐦, M, λ^q, ϕ^q, δ^q,k, D^q,k, δ^q,l, D^q,l sequentially while keeping the other parameters fixed. This is outlined in Algorithm <ref>.
               



 §.§ optimizing with respect to m and M



The first order conditions of optimal 𝐦, M are

    ∂ELBO/∂𝐦 = X'    𝐲 - ∑_i=1^NTϵ_i exp(𝐱'_i𝐦+ 1/2𝐱'_i M𝐱_i)𝐱_i - E_q(M_0^-1)[𝐦-E_q(𝐦_0)] 
    ∂ELBO/∂ M = 1/2{   -X'diag[ϵ_1exp(𝐱'_1𝐦+1/2𝐱'_1 M 𝐱_1), …, 
       ϵ_NTexp(𝐱'_NT𝐦+1/2𝐱'_NT M 𝐱_NT)]X 
        -E_q(M_0^-1)+M^-1}

We follow the numeric algorithm proposed by <cit.> and update 𝐦 using Newton's method, which requires the computation of Hessian matrix

    H^ELBO_𝐦 = -X'diag   [ϵ_1exp(𝐱'_1𝐦+1/2𝐱'_1 M 𝐱_1), …, .
       . ϵ_NTexp(𝐱'_NT𝐦+1/2𝐱'_NT M 𝐱_NT)]X  - E_q(M_0^-1),

and the new 𝐦^(h+1) at iteration h+1 is updated based on the previous value 𝐦^(h)

    𝐦^(h+1) = 𝐦^(h) -(H^ELBO_𝐦)^-1∂ELBO/∂𝐦.

To obtain the optimal M, a fixed-point method is employed and at iteration h+1, it updates M according to

    M^(h+1) = g(M^(h)) = (X'diag   [ϵ_1exp(𝐱'_1𝐦+1/2𝐱'_1 M^(h)𝐱_1), …, .
       . ϵ_NTexp(𝐱'_NT𝐦+1/2𝐱'_NT M^(h)𝐱_NT)]X +E_q(M_0^-1))^-1,

where M^(h) stands for the value at iteration h.

In practice, it is possible that the fixed-point iteration algorithm takes many iterations to converge or it oscillates between two values. Since to compute g(M^(h)) is the most expensive step in the whole algorithm due to the matrix inversion operation of the potentially high dimensional matrix, we wish to avoid such situations from happening, therefore we employ the Anderson acceleration technique as outlined in <cit.>. Essentially, the Anderson acceleration algorithm uses a combination of S previous iterates to form the next iterate, improving convergence rate. The pseudo code of the algorithm is given in Algorithm <ref>.












 §.§ optimizing with respect to lambda



We optimize λ^q jointly through Newton's method as well. The steps to obtain first order conditions and Hessian matrices with respect to λ_1^q,r and λ_2^q,r are illustrated as examples and the rest of the smoothing parameters have similar formulas. The first order conditions of λ^q,r_1, λ^q,r_2 are

    ∂ELBO/∂λ_1^q,r =     -1/2(𝐦^r)'S_1^r𝐦^r-1/2tr(S_1^rM^r)+1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_1^r]+ α-1/λ_1^q,r - β, 
    ∂ELBO/∂λ_2^q,r =     -1/2(𝐦^r)'S_2^r𝐦^r-1/2tr(S_2^rM^r) +1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_2^r]+ α-1/λ_2^q,r - β,

where 𝐦^r and M^r take corresponding entries in 𝐦 and M that are associated with the spline of stringency index. The Hessian matrix block associated with λ^q,r_1, λ^q,r_2 is

    !H_λ_1^q,r,λ_2^q,r^ELBO = [  -1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_1^r(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_1^r] - α-1/(λ_1^q,r)^2                    -1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_1^r(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_2^r];                    -1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_1^r(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_2^r] -1/2tr[(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_2^r(λ_1^q,rS_1^r+λ_2^q,rS_2^r)^-1S_2^r]  - α-1/(λ_2^q,r)^2 ].

Similarly we can derive first order conditions and Hessian matrices with respect to λ_1^q,a, λ_2^q,a, λ_1^q,a, λ_2^q,a, λ_1^q,k,r, λ_2^q,k,r, λ_1^q,k,a, λ_2^q,k,a, λ_1^q,g,a, λ_2^q,g,a.
The joint updating of smoothing parameters λ^q at iteration h+1 is then

    (λ^q)^(h+1) = (λ^q)^(h) - (H^ELBO_λ^q)^-1∂ELBO/∂λ^q.

Here H^ELBO_λ^q is a block diagonal matrix with block entries H_λ_1^q,r,λ_2^q,r^ELBO, H_λ_1^q,a,λ_2^q,a^ELBO,
H_λ_1^q,k,r,λ_2^q,k,r^ELBO,
H_λ_1^q,k,a,λ_2^q,k,a^ELBO and
H_λ_1^q,g,a,λ_2^q,g,a^ELBO. Note that λ^q is subject to positive constraints, therefore in each step, we apply a projected Newton's method with backtracking line search to find the optimal solution <cit.>.



 §.§ optimizing with respect to mu and sigma


First order conditions with respect to μ_i^q, the i-th element of μ^q is

    ∂ELBO/∂μ_i^q = E_q(M_0^-1)[𝐦-E_q(𝐦_0)]∂E_q(𝐦_0)/∂μ_i^q-1/σ^2_μμ_i^q.

The second derivatives are

    ∂^2 ELBO/∂μ_i^q∂μ_j^q = (∂E_q(𝐦_0)/∂μ_i^q)'E_q(M_0^-1)∂E_q(𝐦_0)/∂μ_j^q-1/σ^2_μ

when i=j and 

    ∂^2 ELBO/∂μ_i^q∂μ_j^q = (∂E_q(𝐦_0)/∂μ_i^q)'E_q(M_0^-1)∂E_q(𝐦_0)/∂μ_j^q

when i≠ j. Newton's method is then directly applicable to obtain optimal μ^q.

First order conditions to update (σ^q_i)^2, the i-th element of (σ^q)^2 are

    ∂ELBO/∂(σ^q_i)^2 = - 1/2tr[E_q(M_0^-1) COV_q(𝐦_0)/∂(σ^q_i)^2] - 1/2σ_μ^2
        + 1/2(σ^q_i)^2.

(σ^q)^2 can be updated by setting first order conditions equal to 0 and then solve the linear system. The positivity constraint is automatically satisfied here. 



 §.§ optimizing with respect to phi



To update a specific ϕ^q_i in the l-th position of sequence 1,…,L and k-th position of sequence 1,…,K, first note that the term containing ϕ^q_i in the expectation of quantity [(A^k)'⊗(A^l)']E^1_q( R )(A^k⊗ A^l) is 

    1/(1-(ϕ^q_i)^2)diag[(0,…,0,δ^q,k-k+1,1,…,1)'⊗(0,…,0,δ^q,l-l+1,1,…,1)'],

with l-th entry of vector (0,…,0,δ^q,k-k+1,1,…,1)' equal to δ^q,k-k+1 and k-th entry of vector (0,…,0,δ^q,l-l+1,1,…,1)' equal to δ^q,l-l+1. Similarly, in the expectation of [(A^k)'⊗(A^l)']E^2_q( R )(A^k⊗ A^l) and [(A^k)'⊗(A^l)']E^3_q( R )(A^k⊗ A^l), we change 1/(1-(ϕ^q_i)^2) to -ϕ^q_i/(1-(ϕ^q_i)^2) and (1+(ϕ^q_i)^2)/(1-(ϕ^q_i)^2) accordingly. The first order condition with respect to ϕ^q_i is therefore

    ∂ELBO/∂ϕ^q_i = -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l) [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]'  
       (I_T ⊗ V^q,k⊗ V^q,l)'∂E^*(ϕ,δ^q,l, δ^q,k)/∂ϕ^q_i}
    
        -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l)M(I_T ⊗ V^q,k⊗ V^q,l)'∂E^*(ϕ,δ^q,l, δ^q,k)/∂ϕ^q_i}
    
        +    T-1/2(1/1-ϕ^q_i - 1/1+ϕ^q_i) + α_ϕ-1/1+ϕ^q_i - β_ϕ-1/1-ϕ^q_i

with E^*(ϕ,δ^q,l, δ^q,k) = E_q[(I_T ⊗ A^k⊗ A^l)'E_q( R )(I_T ⊗ A^k⊗ A^l) ].

We can further derive the second order condition

    ∂^2 ELBO/∂ (ϕ^q_i)^2 = -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l) [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' 
       (I_T ⊗ V^q,k⊗ V^q,l)'∂^2 E^*(ϕ,δ^q,l, δ^q,k)/∂ (ϕ^q_i)^2}
    
        -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l)M(I_T ⊗ V^q,k⊗ V^q,l)'∂^2 E^*(ϕ,δ^q,l, δ^q,k)/∂ (ϕ^q_i)^2}
    
        +    T-1/2(1/(1-ϕ^q_i)^2 + 1/(1+ϕ^q_i)^2) - α_ϕ-1/(1+ϕ^q_i)^2 - β_ϕ-1/(1-ϕ^q_i)^2

Newton's method is then applied to update ϕ^q_i for i=1,…,LK sequentially. If the number LK is large in the application, ϕ^q_i can also be updated in parallel as the second order condition involves only ϕ^q_i.



 §.§ optimizing with respect to deltal and Vl



δ^q,l and V^q,l are updated jointly. First order conditions are

    ∂ELBO/∂δ^q,l = -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l) [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' 
       (I_T ⊗ V^q,k⊗ V^q,l)' ∂E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,l}
    
        -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l)M(I_T ⊗ V^q,k⊗ V^q,l)'∂E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,l}
    
        +KT/4   ∑_l=1^L ψ_1(δ^q,l-l+1/2) + δ^l - δ^q,l/4∑_l=1^Lψ_1(δ^q,l-l+1/2)-1/2θ^ltr(D^q,l) + L/2

for degree of freedom δ^q,l,

    ∂ELBO/∂ V^q,l_i,i = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,i}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,i} 
        + KT+δ^l/V^q,l_i,i - δ^q,lV^q,l_i,i/θ^l

for diagonal elements V^q,l_i,i, and

    ∂ELBO/∂ V^q,l_i,j = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}-δ^q,lV^q,l_i,j/θ^l,

for off-diagonal elements V^q,l_i,j with i<j.
In the Hessian matrix, the second derivatives are

    ∂^2 ELBO/(∂δ^q,l)^2 =    KT/8∑_l=1^L ψ_2(δ^q,l-l+1/2) - 1/4∑_l=1^L ψ_1(δ^q,l-l+1/2)+δ^l - δ^q,l/8∑_l=1^L ψ_2(δ^q,l-l+1/2),


    ∂^2 ELBO/(∂ V^q,l_i,i)^2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i,i
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,i}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i,iE^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,i}-KT+δ^l/(V^q,l_i,i)^2 - δ^q,l/θ^l,
 and

    ∂^2 ELBO/(∂ V^q,l_i,j)^2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i,j
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i,jE^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}- δ^q,l/θ^l.

Off-diagonal entries in the Hessian matrix are 

    ∂^2 ELBO/∂δ^q,l∂ V^q,l_i,j = -tr{   E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,l[𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]'  
       (I_T ⊗ V^q,k⊗ V^q,l)'∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}
    
        -tr{   E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,lM (I_T ⊗ V^q,k⊗ V^q,l)'∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,j}- V^q,l_i,j/θ^l,
 with i≤ j, and 

    ∂^2 ELBO/∂ V^q,l_i_1,j_1∂ V^q,l_i_2,j_2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' ∂(I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i_2,j_2
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i_1,j_1}
    
        -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂^2 (I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i_1,j_1∂ V^q,l_i_2,j_2}
    
        -tr{   M ∂(I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,l_i_2,j_2E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i_1,j_1}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂^2 (I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i_1,j_1∂ V^q,l_i_2,j_2}

for i_1≤ j_1, i_2≤ j_2, i_1≠ i_2,j_1≠ j_2.

By definition δ^q,l>L-1 and the diagonal entries are real positive numbers. The constrained optimization problem is again solved using projected method with backtracking line search algorithm. The Newton's method has the merit to converge faster than linear rate to the optimal solution when the current value is in the neighborhood of the solution, however, when the current value is far away from the optimal point, the performance is less satisfying. Therefore, we perform both Newton's method and Cauchy's method in each iteration. The Cauchy's method simply uses the gradient to construct the updating rule and it has superior convergence radius compared to Newton's method.






















































 §.§ optimizing with respect to deltak and Vk


The procedure to update δ^q,k and V^q,k is essentially identical to the one outlined in Section <ref>. The optimization is performed using a hybrid of Cauchy and Newton's methods while the constraints are satisfied with projected method with backtracking line search. Below, we formulate first order conditions and the Hessian matrix.


    ∂ELBO/∂δ^q,k = -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l) [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' 
       (I_T ⊗ V^q,k⊗ V^q,l)' ∂E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,k}
    
        -1/2tr{   (I_T ⊗ V^q,k⊗ V^q,l)M(I_T ⊗ V^q,k⊗ V^q,l)'∂E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,k}
    
        +LT/4   ∑_k=1^K ψ_1(δ^q,k-k+1/2) + δ^k - δ^q,k/4∑_k=1^Kψ_1(δ^q,k-k+1/2)-1/2θ^ktr(D^q,k) + K/2

for degree of freedom δ^q,k,

    ∂ELBO/∂ V^q,k_i,i = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,k)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,l_i,i}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,i} 
        + LT+δ^k/V^q,k_i,i - δ^q,kV^q,k_i,i/θ^k

for diagonal elements V^q,k_i,i, and

    ∂ELBO/∂ V^q,k_i,j = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}-δ^q,kV^q,k_i,j/θ^k,

for off-diagonal elements V^q,k_i,j with i<j.
In the Hessian matrix, the second derivatives are

    ∂^2 ELBO/(∂δ^q,k)^2 =    LT/8∑_k=1^K ψ_2(δ^q,k-k+1/2) - 1/4∑_k=1^K ψ_1(δ^q,k-k+1/2)+δ^k - δ^q,k/8∑_k=1^K ψ_2(δ^q,k-k+1/2),


    ∂^2 ELBO/(∂ V^q,k_i,i)^2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i,i
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,i}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i,iE^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,i}-LT+δ^k/(V^q,k_i,i)^2 - δ^q,k/θ^k,
 and

    ∂^2 ELBO/(∂ V^q,k_i,j)^2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i,j
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i,jE^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}- δ^q,k/θ^k.

Off-diagonal entries in the Hessian matrix are 

    ∂^2 ELBO/∂δ^q,k∂ V^q,k_i,j = -tr{   E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,k[𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]'  
       (I_T ⊗ V^q,k⊗ V^q,l)'∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}
    
        -tr{   E^*(ϕ,δ^q,l, δ^q,k)/∂δ^q,kM (I_T ⊗ V^q,k⊗ V^q,l)'∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i,j}- V^q,k_i,j/θ^k,
 with i≤ j, and 

    ∂^2 ELBO/∂ V^q,k_i_1,j_1∂ V^q,k_i_2,j_2 = -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' ∂(I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i_2,j_2
       E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i_1,j_1}
    
        -tr{   [𝐦-E_q(𝐦_0)] [𝐦-E_q(𝐦_0)]' (I_T ⊗ V^q,k⊗ V^q,l)' 
       E^*(ϕ,δ^q,l, δ^q,k)∂^2 (I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i_1,j_1∂ V^q,k_i_2,j_2}
    
        -tr{   M ∂(I_T ⊗ V^q,k⊗ V^q,l)'/∂ V^q,k_i_2,j_2E^*(ϕ,δ^q,l, δ^q,k)∂(I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i_1,j_1}
    
        -tr{   M (I_T ⊗ V^q,k⊗ V^q,l)'E^*(ϕ,δ^q,l, δ^q,k)∂^2 (I_T ⊗ V^q,k⊗ V^q,l)/∂ V^q,k_i_1,j_1∂ V^q,k_i_2,j_2}

for i_1≤ j_1, i_2≤ j_2, i_1≠ i_2,j_1≠ j_2.






























































































§ ITALIAN MORTALITY ANALYSIS



To analyze non-linear relationship between mortality rates and covariates as well as evolving mortality trends of COVID-19 and other causes before and during the pandemic, we applied our method to the provisional monthly mortality data from Italy. The data, spanning from January 2015 to December 2020 (a total of T=72 months), includes declarations of 18 different causes of death reported by physicians for all deaths in Italy. These causes of death are listed in Table <ref>. Additionally, the death counts are compiled across 420 categories, derived from combining 10 age groups, 2 genders, and L=21 Italian regions, which together with mortality causes, gives 7,560 samples. In total, there are 544,320 observations. A detailed description of this mortality data is available at https://www.istat.it/it/archivio/240401. In our analysis, we set K=17 and exclude cause No.14 complications of pregnancy, childbirth and the puerperium. The reasons for discarding this category are the following. First, death counts are almost all 0's (30,165 out of 30,240); second, estimation with this category is unstable and the convergence rate of the CAVI algorithm drastically slows down; third, excluding one category does not affect the inference on other causes of death. Consequently, we examine N=7,140 time series Y_n,t for n=1,…,N, t=1,…,T. In terms of non-Gaussian state space models, the dimension of observations is 7,140 while the dimension of latent states is 357. Due to high-dimensional structure of the data and the model, conventional sampling algorithms are not applicable, which motivates to approximate the posterior distribution in a variational approach.



We have briefly described covariates included in the model in Section <ref> when formulating the GAM component. Now we elaborate more on one key variable, the Italian Stringency Index (ISI) developed by <cit.>, analogous to the Oxford Stringency Index (OSI) by <cit.>, which measures non-pharmaceutical interventions implemented in Italy to combat COVID-19. The index tracks both national and regional intervention intensity and it is ideal for studying mortality rate in regional level. We focus on non-linear interactions f^k,r(k, r) between ISI and different causes of death, acknowledging the potential varying impacts of the pandemic on other mortality causes as highlighted in existing studies. As for the offsets ϵ_n,t, we consider days per month, monthly aggregated COVID-19 cases by region, and regional population figures for other causes of death. Notably, for external causes like trauma and poisoning, we incorporate a mobility index from the Google COVID-19 Community Mobility Reports (Google LLC "Google COVID-19 Community Mobility Reports". https://www.google.com/covid19/mobility/), to model mortality rate changes per mobility unit. Parameters in the prior distributions are 
α_λ = 1, β_λ = 1000, α_ϕ = 10, β_ϕ = 10,
δ^k = K = 17, δ^l = L, θ^k = K-2 = 15, θ^l = L-2 = 19, σ^2_β = 10,
σ^2_μ = 1. To check whether the proposed CAVI method achieves convergence, we start the algorithm from different initialization. The trajectories of ELBO are shown in Figure <ref>; an optimal ELBO value of 6,564,738 is reached within reasonable steps despite that the initial values are randomly generated.



We first discuss the non-linear relationship between mortality rates and ISI. Figure <ref> shows that GAMs predict similar smoothing patterns for the following mortality causes, endocrine, nutritional and metabolic diseases in Figure <ref>, psychic and behavioral disorders in Figure <ref>, diseases of the nervous system and sense organs in Figure <ref>, diseases of the circulatory system in Figure <ref>, diseases of the respiratory system in Figure <ref>, and external causes of trauma and poisoning in Figure <ref>. When the Italian government implemented loose intervention measures, the mortality rates went down as the intervention became more strict until intervention reached certain median-to-low level (ISI between 20 and 40), the mortality rates took a reverse trend and went up as more strict policies were carried out. In general, early and mild interventions like social distancing and improved hygiene practices not only reduced the spread of COVID-19 but also other infectious diseases, which could benefit individuals with a wide range of health conditions. The initial stages of the pandemic also led to increased health awareness and changes to more regular lifestyle. However, these benefits subsided as healthcare resources were increasingly diverted to treat COVID-19 patients, individuals with other health conditions might have faced delayed or reduced access to necessary health care. More specifically, for some individuals with psychic and behavioral disorders, when the intervention measures were mild, staying at home might have initially reduced stressors such as workplace pressures or social anxiety, potentially benefiting mental health, however, the potential benefits were replaced by increased isolation and loneliness, heightened anxiety and stress as well as substance abuse when situation deteriorated and higher lockdown levels were enforced. As for diseases of the circulatory system and diseases of the respiratory system, initial reduction of mortality rates at low ISI levels could also be attributed to factors like reduction in air pollution, reduced exposure to allergens, less physical strain from commuting as a result of early intervention procedures. Lastly, we comment on external causes of trauma and poisoning. In the early stages of the pandemic, there was rising public vigilance and a focus on safety, which might have extended beyond COVID-19 precautions to general safety practices, potentially reducing accidents. Meanwhile, people were adapting to new routines which might have included safer practices at home, reducing the risk of domestic accidents or poisonings. When the stringency measures prolonged, people started to be impacted by mental health issues as well as social and economic hardship, which could potentially contribute to the upward trend.



Another common non-linear pattern of predicted mortality rates by ISI shared among some mortality causes is an overall downward trend. The mortality causes exhibiting such behaviors are some infectious and parasitic diseases in Figure <ref>, tumors in Figure <ref>, diseases of the digestive system in Figure <ref>, diseases of the skin and sub-cutaneous tissue in Figure <ref> and congenital malformations and chromosomal anomalies in Figure <ref>. For these causes, benefits of social distancing due to high stringency levels dominates factors such as disrupted health care access, delayed diagnoses and treatments. We also note that credible intervals are wider when ISI is close to 100, suggesting distinct mortality patterns appearing in different regions and age groups.

We selected to display three representative mortality patterns between age and causes of death in Figure <ref>. Almost all mortality causes show the same feature that the mortality increases with age, similar to the pattern of some infectious and parasitic diseases shown in Figure <ref>. The mortality rate of COVID-19, however, peaked around slightly over 80 years old. Although there have been studies pointing out that age is a positive predictor of the morality rate <cit.>, their stratification usually sets 80 years old and above as a group. The age groups in the data set that we analyse are more refined with 4 separate age groups for age 80+. Our result offers new insights on the relationship between COVID-19 mortality and age.



As for temporal evolving z_n,t^* in the Poisson rates, we focus on three aspects. The approximated posterior distribution of μ indicates the average mortality levels of each cause of death in each region unexplained by covariates in the model. From Figure <ref>, we can see that both tumors in Figure <ref> and diseases of the circulatory system in Figure <ref> contribute most to death counts whereas some morbid conditions that originate in the perinatal period shown in Figure <ref> have the lowest mortality rate on average. Figure <ref> also reveals certain geographical disparity. For instance, Piedmont, Liguria, Apulia and Abruzzo are the regions mostly affected by COVID-19 according to Figure <ref>, however this conclusion is in contradiction with the observation that regions such as Lombardy, Veneto, and Campania experienced many cases and deaths. Recall that μ measures unexplained mortality rate per 1,000 cases since we include COVID-19 cases in offset ϵ_n,t, the discrepancy could be attributed to factors including healthcare system quality, timing of the outbreak, testing capacity and so on. For example, regions with less comprehensive testing may report fewer mild or asymptomatic cases, resulting in a higher mortality rate, as the total cases are underestimated. Nevertheless, to answer this question, we need  access to more covariates. Certain regional gap is more in line with conventional knowledge. Figure <ref> shows that endocrine, nutritional and metabolic diseases are more deadly in the south as southern regions have higher rates of obesity and diabetes. Temperature is another risk factor behind the phenomenon.





We now discuss the other two aspects associated with z^*_n,t, the dependence structure captured by (V^q,l)'V^q,l for the scale matrix D^q,l of the variational Wishart density and the correlation captured by (V^q,k)'V^q,k for the scale matrix D^q,k. Since D^q,l contains spatial information that controls the dynamic of z^*_n,t, we compute the partial correlation matrix, that is the normalized inverse, of D^q,l and then set the entries in the inverse matrix whose absolute value are below 0.1 to 0 while replace the remaining entries with 1. Figure <ref> maps the resulting adjacency matrix into Italy where red edges between two regions stand for 1's in the adjacency matrix, indicating a direct statistical relationship between them that is not explained by their relationships with other regions. There edges form four clusters; northern Italian regions are connected, and  Lombardy appears to be the central hub for disease spread with many edges connected to it, possibly due to higher mobility or being a major transportation center. The remaining three clusters lie in the center and the southern Italy. The behavior indicates that regions inside each cluster might experience similar trends in mortality rates due to migration patterns, economic ties, climate conditions, or population behaviors. On the other hand, the absence of an edge between two regions implies conditional independence, meaning that once the influence of all other regions is controlled for, there is no direct statistical relationship between these two regions. This is the case with most regions.

Lastly, we turn to the correlation structure among various mortality causes. Figure <ref> shows strong positive correlation inside the following 5 causes of death, endocrine, nutritional and metabolic diseases, psychic and behavioral disorders, diseases of the nervous system and sense organs, diseases of the circulatory system, and diseases of the respiratory system. In terms of the relationship between COVID-19 and other causes of death, we only observe weak negative relationship between COVID-19 and symptoms, signs, abnormal results and ill-defined causes, which suggests that improved detection of COVID-19 death likely reduced the number of deaths classified under ill-defined causes. The lack of correlation between between COVID-19 and other causes of death indicates that after controlling government intervention intensity during the pandemic, COVID-19 itself does not have strong influence on moralities in other death categories. 





§ SUMMARY AND FUTURE WORK



We develop an effiecient variational inference procedure for Bayesian dynamic GAMs where the dependent variable follows a Poisson distribution. The GAM component in the Poisson rate captures non-linear relationship between the outcome and covariates while the dynamics in the rate has a state space representation. The latent states in the model are assumed to be stationary with covariance matrix defined using a Kronecker product that disentangles one large covariance matrix into smaller covariance matrices. High-dimensional setting is considered, which prohibits the use of exact sampling algorithms, therefore we adopt variational algorithm to achieve fast convergence. The approach is applied to Italian mortality data to facilitate our understanding of mortality pattern changes during the COVID-19 outbreak.

Even though we concentrate on explanatory power of the model, it is potentially suitable for forecasting as the state space component models dynamics in mortality patterns. Another extension makes more sophisticated assumptions on the covariance matrix of the latent states. For instance, it is interesting to impose hierarchical G-Wishart distribution prior instead of Wishart prior on the precision matrices so that conditional dependence structure can be directly inferred from the model. Alternatively, the covariance matrix can be time-varying, allowing for more flexibility to account for changes in the dependence structure, especially before the COVID-19 pandemic and after the pandemic. We leave these discussion to future work.



apalike


