Bayesian inference and impact of parameter prior specification in flexible multilevel nonlinear models in the context of infectious disease modeling

1.   Introduction

Multilevel nonlinear models, also known as nonlinear mixed-effects models (NLMMs), offer a sophisticated approach to data analysis by accounting for both within- and between-subject variability. These models are particularly well-suited to leverage the richness of repeated measurements. Over the past few decades, there has been growing interest in applying NLMMs to infectious disease modeling to gain a more comprehensive understanding of disease dynamics and to improve control strategies. This interest is particularly driven by the need to account for substantial heterogeneity across countries when conducting in-depth analyses of micro-infection dynamics.

Estimation methods may lack robustness if not carefully designed. To enhance the application of these models in infectious disease analyses, Bayesian methods have been proposed as a valuable framework to incorporate a priori information more rigorously. This approach flexibly addresses sparse and unbalanced longitudinal data from individual subjects ^[1]. However, a common concern with Bayesian methods is their sensitivity to various aspects of the modeling process, including prior distributions and initial values ^[2]. A sensitivity analysis serves as a fundamental tool to investigate the model uncertainty. Prior distributions play a critical role in determining the extent to which subjective opinions influence the model estimation. As ^[3] noted, the effect of a prior distribution depends on factors such as its degree of informativeness and the sample size. This study builds on recent developments in Bayesian flexible multilevel nonlinear models (FMNLMs) for infectious disease analyses ^[3,4]. Specifically, it examines how prior informativeness and accuracy, along with the sample or cluster size, affect the parameter estimates in FMNLMs. The focus is on understanding whether strong but inaccurate subjective beliefs about epidemic parameters can significantly impact estimation outcomes.

The literature on the specific impact of different prior distributions in infectious disease modeling remains limited. This investigation aims to provide insights into the consequences of using informative prior distributions in this specific context. We focus on two flexible models, namely the scale mixtures of skew-normal nonlinear mixed model (SMSN-NLMM) and the semi-nonparametric nonlinear mixed model (SNP-NLMM), to examine how varying levels of prior informativeness can influence the model convergence, parameter estimation, and computation time through a simulation study. Starting from completely diffuse prior distributions, we introduce reasonable modifications to the prior assumptions, recompute the posterior quantities of interest, and assess whether these changes significantly affect the interpretations or conclusions.

In this framework, we assume that the random terms (random effects and residuals) follow either a scale mixtures of skew-normal (SMSN) distribution ^[5,6,7] or a semi-nonparametric (SNP) distribution ^[8,9]. Using specific epidemic data characteristics, such as sample size, cluster size, or reporting frequency, we evaluate the empirical performance of these models in recovering the epidemic parameters under varying levels of prior informativeness.

This paper is structured as follows: in Section 2, we introduce a case study of infectious disease dynamic models and describe the real dataset used in this analysis; Section 3 provides a concise summary of the multivariate flexible multilevel nonlinear models (SMSN-NLMM and SNP-NLMM) and their associated Bayesian estimation procedures; Section 4 presents a simulation study to evaluate the methodology; in Section 5, we apply the proposed methodology to the COVID-19 dataset described in Section 2, thus showcasing its utility and presenting the results; and finally, Section 6 presents a discussion of the findings and their implications before a short summary of the key findings in the conclusion.

2.   The data and infectious disease dynamic models: A case study of the COVID-19 pandemic in Francophone West Africa

Infectious disease outbreak data, such as reported cases and deaths, are collected within countries or regions. Countries that experience outbreaks earlier tend to accumulate more data compared to those in the same region that experience outbreaks later. To prepare for potential outbreaks, countries that have not yet reached the peak often leverage existing data from neighboring countries to model the average patterns and predict the trajectory of the epidemic. Additionally, the level of surveillance and preparedness against outbreaks varies across countries, thus leading to a diverse range of responses that can influence the epidemiological parameters. Sub-notification is a significant challenge in accurately understanding the true number of infections in the population due to the varying testing capacities between countries.

Figure 1 presents example cases, showing the reported number of daily COVID-19 cases across different countries. It displays daily reported cases from the first identified cases until October 1, 2020, in seven West African countries (Benin, Burkina Faso, Ivory Coast, Guinea, Mali, Niger, and Senegal). These countries were selected for comparison and discussion, as previously chosen by ^[10] to examine the epidemic dynamics in the region. The authors emphasized that these countries demonstrate the varying dynamics of the pandemic, despite implementing different government measures within the same geographical area.

The datasets can be accessed via the link https://github.com/CSSEGISandData/COVID-19/ provided by Johns Hopkins University in collaboration with the Center for Systems Science and Engineering (see ^[11]).

It is worth noting that, in contrast to countries in regions such as America or Europe, where cases have been reported daily since the first identified case, the Francophone West African countries presented here experienced several consecutive days with zero detected cases, despite being at different stages of the pandemic curve. The different stages of disease spread are crucial information that should be incorporated into the model. Moreover, infectious disease data often exhibit skewness, outliers, or heavy-tailed behavior ^[2,12,13]. For instance, Figure 2a shows the density curve of residuals from reported cases, thus revealing significant skewness in the reported responses. It is important to note that these reported responses have been subjected to linear transformations, as proposed in ^[13], $y_{ij} = z_{ij}/k_z$, where $z_{ij}$ is the number of reported cases at the $i^{th}$ country and the $j^{th}$ day since the first reported case, and $k_z$ is chosen to be the sample standard deviation from the country data, which is the smaller one in the observed data. Non-normality is also confirmed by the Q-Q plot of the fitted residuals in Figure 2b. To characterize skewness with heavy tails, which often appear in such epidemic data, we will develop a Bayesian flexible parametric (SMSN-NLMM) or semi-nonparametric (SNP-NLMM) multilevel nonlinear models in conjunction with Model 3.1 under the assumption of SMSN and SNP distributions. That is, we assume that both $\mathit{\boldsymbol{\varepsilon_{i}}}$ and $b_i$ either follow a multivariate SMSN distribution as proposed by ^[5,6,13] and ^[7] or a SNP as introduced by ^[8] and ^[9]. These methodologies jointly accommodate the different stages of the diseases while borrowing information from the different time series and considering the epidemic data structure to provide a more robust and reliable fit and prediction.

3.   Model formulation and Bayesian approach

3.1.   Nonlinear mixed-effect models (NLMMs)

NLMMs can be applied across various fields beyond infectious disease modeling. In this section, we focus on their application to infectious disease data. Specifically, we treat subjects as either countries or regions within a country.

Let us consider $n$ independent countries (or regions of a country), each with $m_{i}$ measurements from the $i^{th}$ country (or region of a country). Let $y_{ij}$ denote the observed response and $x_{ij}$ represent the vector of within-subject covariates for the $j^{th}$ measurements ($j = 1, ..., m_{i}$) from the $i^{th}$ country ($i = 1, ..., n$). The response vector for the $i^{th}$ country is denoted as $\mathit{\boldsymbol{y}}_{i} = (y_{i1}, ..., y_{im_i})^{\top}$. Next, let us define a scalar-valued differentiable nonlinear function $\mu(x_{ij}|\phi_{i})$ that describes the average trajectory of the response. The general nonlinear mixed-effects model is expressed as follows ^[14] :

where $\mathit{\boldsymbol{\varepsilon}}_{i}$ is a $n_{i}-$dimensional vector of independent within-country errors. The mixed effects parameter vector $\mathit{\boldsymbol{\phi}}_{i}$ is modeled as follows ^[14] :

In Eq (3.2), $\mathit{\boldsymbol{A}}_{i}$ is a design matrix of size (r, p) for the fixed effects that possibly depends on elements of ${\mathit{\boldsymbol{x}}}_{i}$, and $\mathit{\boldsymbol{B}}_{i}$ is a design matrix of size (r, q) for the random effects that allow the incorporation, for instance, of "time-varying" covariates in the random effects with r, which represents the number of subject-specific parameters included in the nonlinear function. $\mathit{\boldsymbol{\beta}}$ is a $p$-dimensional locator vector of fixed-effects. $\mathit{\boldsymbol{b}}_{i}$ is a $q$-dimensional vector of random-effects (assumed to be mutually independent across subjects and independent of the within-country errors $\mathit{\boldsymbol{\varepsilon}}_{i}$) associated with the $i^{th}$ country.

The random-effects vector $\mathit{\boldsymbol{b}}_{i}$ has the density function $g_b(.|\mathbf{\theta})$ and the conditional distribution of $\mathbf{y}_{i}$ given $\mathit{\boldsymbol{b}}_{i}$, and the covariate vector $x_{ij}$ has the density function $f_{Y}(.|\mathit{\boldsymbol{b}}_{i}, \mathit{\boldsymbol{\theta}})$. The estimation of the parameter vector $\mathit{\boldsymbol{\theta}}$ (which contains $\mathit{\boldsymbol{\beta}}$ and distribution parameters) of NLMMs is commonly achieved through maximization of the marginal likelihood $L(.|\mathbf{y})$ of the observed data $\mathbf{Y} = (\mathit{\boldsymbol{y}}_{1}^{\top}, ..., \mathit{\boldsymbol{y}}_{n}^{\top})^{\top}$ as follows:

Both $g_b(.|\mathit{\boldsymbol{\theta}})$ and $f_{Y}(.|\mathit{\boldsymbol{b}}_{i}, \mathit{\boldsymbol{\theta}})$ are often assumed to be multivariate normal densities; however, the observed data may not support these assumptions.

3.2.   Scale mixture of skew-normal nonlinear mixed-effects model (SMSN-NLMM) specification

The parametric flexible multilevel nonlinear model is the nonlinear mixed effects model in which the random terms (random effects and errors) follow a scale mixture of skew-normal (SMSN) distributions. It's expressed as follows: $\mathit{\boldsymbol{y}}_{i} = \mathit{\boldsymbol{\mu}}_{i}(\mathit{\boldsymbol{\phi}}_{i}) +\mathit{\boldsymbol{\epsilon}}_{i}\text{; }\quad\mathit{\boldsymbol{\phi}}_{i} = \mathit{\boldsymbol{A}}_{i}\mathit{\boldsymbol{\beta}} + \mathit{\boldsymbol{B}}_{i}\mathit{\boldsymbol{b}}_{i}$ with

where $\eta$ and $\mathit{\boldsymbol{\eta}}_b$ are location shifts used to zero means of residual errors and random effects, respectively, $\omega^{2}$ and $\mathit{\boldsymbol{\Omega}}_b$ are the scale parameters that know the degrees of freedom ($\nu$ and $\nu_b$) and the shape parameters ($\lambda$ and $\mathit{\boldsymbol{\lambda}}_b$) to effectively capture skewness and heavy-tailed behavior. Meanwhile, under a mixed effects design with $n$ subjects and $m_i$ measurements for the $i^{th}$ subject, the SMSN-NLMM is given for the $j^{th}$ outcome $Y_{ij}$ ($i = 1, \cdots, n$; $j = 1, \cdots, m_i$) as follows:

with $E\left(Y_{ij} \mid \mathbf{x}_{ij}, {\boldsymbol{\phi}_i} \right) = \mu(x_{ij}\mid \mathit{\boldsymbol{\phi_i}}) = \mu_{ij}$; $Var\left(Y_{ij} \mid {\boldsymbol{\phi}_i} \right) = \sigma^2$; $E(\mathit{\boldsymbol{b_i}}) = 0$; $Var(\mathit{\boldsymbol{b_i}}) = \mathit{\boldsymbol{\Sigma_b}}$. Let $\kappa$ be the mixture variable associated to response $y_{ij}$ with the cumulative distribution function $H _{\kappa}(\mathit{\boldsymbol{\nu}})$ and the density function $h_{\kappa}(\mathit{\boldsymbol{\nu}})$ and $U$, the mixture variable associated to random effect vector $b_i$ with cumulative distribution function $H _{u}(\mathit{\boldsymbol{\nu}}_b)$, and the density function $h_{u}(\mathit{\boldsymbol{\nu}}_b)$, whereas $\kappa_{et} = E\{\kappa^{-t/2}\} < \infty$ ($t = 1, 2$) and $U_{t} = E\{U^{-t/2}\} < \infty$ ($t = 1, 2$). Then, the expectation of the response, and the residual variance, and the variance-covariance matrix of random effects are expressed as follows ^[7]: $\xi_{ij} = \mu_{ij}-c\kappa_{e1}\delta_e$; $\omega^2 = \kappa_{e2}^{-1}(\sigma^2+c^2\kappa_{e1}^{2}\delta_e^2)$; $\mathit{\boldsymbol{\xi}}_b = -cU_{1}\delta_b$ and $\mathit{\boldsymbol{\Omega}}_b = U_{2}^{-1}(\Sigma_b+c^2U_{1}^{2}\delta_b\delta_b\top)$ with $c = \sqrt{2/\pi}$.

A key feature of SMSN-NLMM is that it can be formulated within a flexible hierarchical representation that is useful to easily write implementable STAN or BUGS codes and for analytical derivations. While referring to ^[15] and ^[16], it follows from (3.4) that

where $\delta_e = \omega \lambda (1 + \lambda^2)^{-1/2}$, ${\mathit{\boldsymbol{\delta}}}_{b} = (1 + {\boldsymbol{\lambda}}_ {b}^{\top} {\boldsymbol{\lambda}}_ {b})^{-1/2} \mathit{\boldsymbol{\Omega}}_b^{1/2} {\boldsymbol{\lambda}}_b$, $\bar{\omega}^2 = \omega^2 - \delta_e^2$, $\bar{\boldsymbol{\Omega}}_b = {\boldsymbol{\Omega}_b} - {\mathit{\boldsymbol{\delta}}}_{b} {\mathit{\boldsymbol{\delta}}}_{b}^{\top}$. $t$ and $s$ are positive standard half normal random variables (absolute value of a standard normal random variables) denoted $H\mathcal{N}(0, 1)$.

For instance, candidate shape parameters are $\delta_e$ and $\mathit{\boldsymbol{\delta}}_{b}$. Meanwhile, the parameters to be estimated are the fixed effects $\mathit{\boldsymbol{\beta}}^{\top}$, the shape parameters $\delta_e$ and $\mathit{\boldsymbol{\delta}}_{b}$, the degrees of freedom $\nu$ and $\nu_b$, and the variance components $\sigma^2$ and $\Sigma_b$. The parameters $\nu$ and $\nu_b$ of the mixture variables control the shape of the SMSN distributions with $\delta_e$ and $\mathit{\boldsymbol{\delta}}_b$, especially the skewness and the excess of kurtosis. Taking $\mathit{\boldsymbol{Y}}$ as the complete data (observed data with a missing data), the complete likelihood function associated with $\mathit{\boldsymbol{Y}}$ under the hierarchical representation can be expressed using Bayes Theorem as follows:

where the notation $f_{x}(x)$ refers to the density function of $x$. $\mathit{\boldsymbol{\theta}} = (\mathit{\boldsymbol{\beta}}^{\top}, \sigma^{2}, \mathit{\boldsymbol{\varrho}}^\top, \mathit{\boldsymbol{\rho}}^\top, \delta_e, \mathit{\boldsymbol{\delta}}_{b}, \nu, \nu_b)^{\top}$ takes $\mathit{\boldsymbol{\varrho}}^\top$, $\mathit{\boldsymbol{\rho}}^\top$ as a random effects dispersion parameters matrix ($\mathit{\boldsymbol{\Sigma}}_b$ is supposed to be unstructured) diagonal elements and off diagonal elements, respectively. $\mathit{\boldsymbol{\kappa}}_{i} = (\kappa_{i1}, ..., \kappa_{im_i})^{\top}$ and $\mathit{\boldsymbol{T}}_{i} = (t_{i1}, ..., t_{im_i})^{\top}$.

3.3.   Semi non-parametric nonlinear mixed-effects model (SNP-NLMM) specification

SNP-NLMM is the nonlinear mixed effects model in which the random terms (random effects and errors) follow SNP distributions. The SNP-NLMM is expressed as follows: $\mathit{\boldsymbol{y}}_{i} = \mathit{\boldsymbol{\mu}}_{i}(\mathit{\boldsymbol{\phi}}_{i}) +\mathit{\boldsymbol{\epsilon}}_{i} \text{; }\quad \mathit{\boldsymbol{\phi}}_{i} = \mathit{\boldsymbol{A}}_{i}\mathit{\boldsymbol{\beta}} + \mathit{\boldsymbol{B}}_{i}\mathit{\boldsymbol{b}}_{i}$ with

where the subscript $i$ represents the subject index, $\eta$ and $\mathit{\boldsymbol{\eta}}_b$ are the location shifts used to zero means of error and random effects, respectively, a shape vector $\mathit{\boldsymbol{\psi}} \in [-\pi/2, \pi/2]$, and a scale matrix $\mathit{\boldsymbol{\Sigma}} > 0$ (positive definite) with Cholesky factor $\mathit{\boldsymbol{R}}$ $(\mathit{\boldsymbol{RR}}^{T} = \mathit{\boldsymbol{\Sigma}})$ and integer valued expansion order $K$. Here, we assume that the residual and the random effects are mutually independent across $i$, and write the following:

where $R$ is the residual dispersion parameter, $\mathit{\boldsymbol{R}}_b$ is a ($q \times q$) upper triangular scale matrix for random effects, and $Z_{ij} \sim SSNP(\psi, K)$ and $\mathit{\boldsymbol{Z}}_{b_i} \sim \mathcal{SSNP}_q(\mathit{\boldsymbol{\psi}}_b, K)$ are standard SNP variables.

$SSNP_{q}(\mathit{\boldsymbol{\psi}}, K)$ defined in the sequel denotes a standard SNP variable with a null location vector and the identity scale $(SNP_q(\mathit{\boldsymbol{0}}; \mathit{\boldsymbol{I}}_q; \mathit{\boldsymbol{\lambda}}; K))$. A standard $SSNP_{q} (\psi, K)$ variable $\mathit{\boldsymbol{z}}$ has a density of the following form:

where $P_{K}(.|\mathit{\boldsymbol{\psi}})$ is a multivariate polynomial of order $K$. Let $\mathit{\boldsymbol{\alpha}}$ be a multi-index, i.e., a $q \times 1$ vector with nonnegative integer elements: $\mathit{\boldsymbol{\alpha}} = (\alpha_{1}, ..., \alpha_{q})^{\top}$. Let $|\mathit{\boldsymbol{\alpha}}| = \sum_{k = 1}^{q}\alpha_{k}$ and $\mathit{\boldsymbol{z}}^{\mathit{\boldsymbol{\alpha}}} = \prod_{k = 1}^{q}\mathit{\boldsymbol{z}}_{k}^{\alpha_{k}}$. Then, the $K$ order multivariate polynomial can be written as follows:

Therefore, the densities of response and random effects can be expressed as follows:

where $z_{ij} = R^{-1} (y_{ij} - \mu_{ij}) + \eta$, and $\mathit{\boldsymbol{z}}_{b_i} = \mathit{\boldsymbol{R}}_b^{-1} \mathit{\boldsymbol{b}}_i + \mathit{\boldsymbol{\eta}}_b$; $\phi_q(.|\mathit{\boldsymbol{\gamma}}_b, \mathit{\boldsymbol{\Sigma}}_b)$ is q-variate normal density with mean $\mathit{\boldsymbol{\gamma}}_b$ and covariance matrix $\mathit{\boldsymbol{R}}_b\mathit{\boldsymbol{R}}_b^\top$ knowing $\mathit{\boldsymbol{\Sigma}}_b = \mathit{\boldsymbol{R}}_b\mathit{\boldsymbol{R}}_b^\top$. Thus, following ^[17], $E(Y_{ij}) = \mu_{ij}+\eta-\sigma E(Z_{ij})$; $var(Y_{ij}) = \sigma^2 var(Z_{ij})$; $E(b_{i}) = \eta_b+\Sigma_b^{1/2} E(Z_{i})$ and $var(b_{i}) = \Sigma_b var(Z_{i})$, where $E(Z_{ij}) = E(UP_{K}^{2}(\mathit{\boldsymbol{U}}; \mathit{\boldsymbol{\psi}}))$ and $var(Z_{ij}) = E(U^2P_{K}^{2}(\mathit{\boldsymbol{U}}; \mathit{\boldsymbol{\psi}}))-E^2(Z_{ij})$ and $U \sim \mathcal{N}(0, 1)$. The marginal likelihood for $\mathit{\boldsymbol{\theta}} = (\mathit{\boldsymbol{\beta}}^{\top}, R, \mathit{\boldsymbol{R}}_b, \mathit{\boldsymbol{\psi}}, \mathit{\boldsymbol{\psi}}_b)^{\top}$ for a fixed K is as folows:

where $\mathit{\boldsymbol{b}} = (\mathit{\boldsymbol{b}}_{1}^{\top}, \dots, \mathit{\boldsymbol{b}}_{m}^{\top})^{\top}$, $\mathit{\boldsymbol{y}} = (\mathit{\boldsymbol{y_{1}}}^{\top}, ..., \mathit{\boldsymbol{y}}_{m}^{\top})^{\top}$, and $\mathit{\boldsymbol{y}}_{i}$ = $(y_{i1}, ..., y_{im_i})^{\top}$, $i = 1, ..., n$. Moreover, the SNP-NLMM exploits the SNP approach to represent the density, where the degree of flexibility is controlled by the parameter K, which is chosen via inspection of the standard information criteria. Thus, the distribution density is parameterized in terms of a finite set of parameters $\theta$, where the dimension of $\theta$ (through $\mathit{\boldsymbol{\psi}}$ and $\mathit{\boldsymbol{\psi}}_b$) depends on K. When $K = 0$, the semi-nonparametric density reduces to a standard q-variate normal density, (3.13 and 3.14) are $N_{ni}({\eta}, R^2)$ and $N_{q}(\mathit{\boldsymbol{\eta_b}}, R_bR_b^\top)$ with ${\eta} = \mu_{ij}$ and $R^2 = \sigma^2$, so that the usual normal specification for NLMMs is a special case ^[9]. The approach is widely recommended for routine use in situations where a departure from the usual normal assumption for the random terms is either present or suspected. Often, K = 1 or 2 is sufficient to represent an adequate departure from normality ^[17]. Most studies considered the performance for varying the expansion order parameter $K$ and confirmed this underlined proposition ^[8,9,18,19]. A larger K gives more flexibility to represent the random terms distribution; however, choosing a K that is too large will result in an inefficient representation ^[17]. For facilitation and comparison purposes in this study, we just considered the SNP-NLMM with K = 1.

3.4.   Bayesian proposition

To complete the Bayesian specification, we need to consider prior distributions for all the unknown model parameters, $(\mathit{\boldsymbol{\beta}}^{\top}, \sigma^{2}, \mathit{\boldsymbol{\tau^\top}}, \mathit{\boldsymbol{\rho^\top}}, \delta_e, \mathit{\boldsymbol{\delta}}_{b}^{\top}, \nu, \nu_b, \psi, \mathit{\boldsymbol{\psi_{b}}})^{\top}$.

In this study, a Bayesian approach is adopted, thereby considering different levels of informative priors by fixing the mean hyperparameter for the prior distribution to represent relatively close, moderate, or distant values from the corresponding population value of the model parameters. These variations are specifically applied to the model's mixed effects parameters, particularly the fixed effects, while all other parameters (such as precision, scale, and structural model parameters) are completely assigned to diffuse priors, as described by ^[20]. Therefore, the primary factors manipulated in this study are the mean hyperparameter and the variance hyperparameter for each nonlinear parameter. A detailed description of these factors is provided below.

In Bayesian modeling, the elements of the parameter vector $\mathit{\boldsymbol{\theta}}$ are assumed to be independent. Consequently, the joint prior distribution ($\pi(\mathit{\boldsymbol{\theta}})$) for all unknown parameter densities can be expressed for the SMSN-NLMM and SNP-NLMM models, respectively, as follows:

By knowing the models for the observed data and the prior distributions for the unknown model parameters, a statistical inference for the parameters can be carried out using their posterior distributions under the Bayesian framework. Given the observed data, the joint posterior density of the parameter vector $\theta$ is expressed as follows:

In general, expression (3.18) does not have a closed form, and MCMC procedures can be used to sample based on (3.18) using the Hamiltonian Monte Carlo (HMC) algorithm. This latest is implemented in the Stan program along with the Bayesian Regression Models (brms) procedure ^[21,22]. This approach is particularly advantageous, as it allows the developed models to be implemented using the rstan package in R, which is a freely available software ^[23]. In this study, the tools and techniques are situated within a Bayesian framework and rely on the HMC algorithm implemented through rstan. This approach offers a practical alternative to the computationally expensive traditional MCMC methods, thus making it a more efficient option for Bayesian inference ^[4].

4.   Simulation experiment

4.1.   Data generating process

Considering the potentially asymmetric and skewed nature of infectious disease data ^[16], a natural extension of the modeling approach is to adopt a nonlinear mixed-effects model where the random components (residuals and random effects) are assumed to follow a flexible distribution. One such distribution is the skew Student-t, a specific case of the SMSN-NLMM framework ^[12,13], which incorporates the degrees of freedom ($\nu$) and the shape parameters ($\delta$) to capture skewness and heavy-tailed behavior:

where $\mu_{ij} = \mu_{ij}(\mathit{\boldsymbol{\phi}}_{i})$ is a non-linear function with the individual parameter vector

$\eta$ and $\mathit{\boldsymbol{\eta}}_b$ are location shifts used to zero means of residual errors and random effects, respectively. The intra-individual regression function considered is the derivative of the generalized logistic used to model COVID-19 death curves from American countries ^[13], as defined by the following:

In this equation, $\phi_{2i} = exp\{\beta_{2}+b_{2i}\}$, $\phi_{3i} = exp\{\beta_{3}+b_{3i}\}$, and $\phi_k = exp\{\beta_{k}\}$, for $k = 1; 4$, with the exponential transformation being used to ensure positiveness of the parameters; $x_{ij} = t_{ij}$ is the time of observation $j$ $(1 \leq j\leq n_i)$ on individual $i$ $(1 \leq i \leq n)$. The nonlinear mixed parameters is particularly express as follows: $\mathit{\boldsymbol{\phi}}_{i} = exp(\mathit{\boldsymbol{X}}_{i}\mathit{\boldsymbol{\beta}}+\mathit{\boldsymbol{b}}_{i})$.

Referring to ^[24], random effects are incorporated into (4.1) to facilitate a multivariate approach and enable information sharing across different time series. The parameter $\phi_3$ controls the infection rate, while $\phi_4$ acts as an asymmetry parameter. Additionally, $\phi_1$, $\phi_2$, and $\phi_4$ influence the curve's asymptote, with the total expected cases given by $\frac{\phi_1}{\phi_2^{\phi_4}}$. The peak of the curve occurs at the time point $t = -\frac{ln(\phi_2/\phi_4)}{\phi_3}$. The initial infection rate, $\phi_3$, of an epidemic is a critical measure of disease transmission, as it is often utilized to estimate the basic reproduction number $R_0$. This key epidemiological parameter, which represents the average number of secondary infections caused by a single infected individual in a fully susceptible population, is calculated as $R_0 = e^{\phi_3T}$ ^[25].

The parameters of the epidemic model are treated as random variables and assigned specific probability distributions. To generate data from the epidemic model, each parameter was initially simulated from its respective distribution (see population values in Table 1). The nonlinear parameters fixed effect population values are presented in the table with their standard errors in the parenthesis. Then, the resulting parameter values were used to simulate the dataset. These predefined distributions serve as benchmarks to evaluate the fitted posterior distributions derived from the simulated dataset.

We simulated $M$ datasets under the assumption that specific infectious disease measurements were reported daily over two timeframes: a maximum of 1 month for 6 clusters (representing a small dataset) and a maximum of 1 month for 12 clusters (representing a large dataset), thereby taking clusters for countries or areas of the country. A limit of quantification for the reported cases was set at 0.0 individuals, meaning any simulated measurements below this threshold were adjusted to 0.0. The parameter values used in the simulations were from estimates obtained through a nonlinear least squares (NLS) estimation for fixed effects and the maximum likelihood estimation proposed by ^[13] applied to some Latin American countries COVID-19 data (Table 1).

The primary focus of this paper has been to examine the effects of the varying degrees of (in) accurate prior distributions placed on parameters within the context of FMNLM. Meanwhile, it involves manipulating the mean (and variance) hyperparameter for the nonlinear parameter priors. Specifically, informative priors were only placed only on the nonlinear parameters ($\beta \sim \mathcal{N}_p(\mathit{\boldsymbol{\beta_0}}, \mathit{\boldsymbol{\Sigma_0}})$; $\mathit{\boldsymbol{\beta_0}}$, $\mathit{\boldsymbol{\Sigma_0}}$ are constituted of the defined hyperparameters) in this study (Table 1), while all other parameters were given completely diffuse priors (refering to ^[20], ^[12] and ^[19]).

Following the systematic evaluation framework of ^[3], the mean hyperparameter term was determined as a function of the fixed variance hyperparameter. Specifically, the mean hyperparameter was decreased by either 1 standard deviation (SD), 3SD, or 5SD based on the fixed variance hyperparameter term. However, the levels of variance hyperparameters were systematically set, taking either 10.00% (relatively low variance), 20.00% (relatively moderate variance), or 50.00% (relatively high variance) of the corresponding population value. Meanwhile, this setup resulted in nine distinct prior settings: 1SD10, 1SD20, 1SD50, 3SD10, 3SD20, 3SD50, 5SD10, 5SD20, and 5SD50, where $\alpha$SD$\gamma$ denotes an $\alpha$SD mean hyperparameter and a $\gamma$% variance hyperparameter level. This method was deemed the most straightforward and practical to define the critical conditions of interest in this investigation, thus giving the infeasibility of exploring all possible combinations of mean and variance hyperparameter conditions. We refer to ^[3] for further details.

Table 1 provides a summary of the population model values and their corresponding priors across all mean and variance hyperparameter levels.

As an example, all mean hyperparameter levels for the first parameter are illustrated in Figure 3. This figure demonstrates the prior distributions relative to the population value for all nine hyperparameter levels considered in this study.

For the population value $V$, the corresponding variance hyperparameters are computed as follows: $Var_{1} = \arrowvert V \arrowvert*10/100$, $Var_{2} = \arrowvert V \arrowvert*20/100$, $Var_{3} = \arrowvert V \arrowvert*50/100$, with the corresponding SDs: $SD_{1} = \sqrt{Var_{1}}$ for 10%, $SD_{2} = \sqrt{Var_{2}}$ for 20%, and $SD_{3} = \sqrt{Var_{3}}$ for 50%. For example, the mean hyperparameters for 10% variance are as follows: $V-(1*SD_{1})$, $V-(3*SD_{1})$, and $V-(5*SD_{1})$ for 1SD, 3SD, and 5SD, respectively.

4.2.   Analysis and performance criteria

Under each of the simulation scenarios, we fitted the Bayesian flexible semi-nonparametric multilevel nonlinear models and the Bayesian flexible parametric multilevel nonlinear models in addition to the actual normal NLMM to analyse the characteristics of each fitting model.

The Hamiltonian Monte-Carlo (HMC) algorithm has been proposed to be suitable for such a modeling approach ^[4]. Meanwhile, we took advantage of brms within the STAN interface in the R software ^[23] for the analysis. The population parameters were initialized at the values obtained through the nlme package setting and the random effects were initialized to zero. The scale parameters were initialized to 1, while the correlation, shape, and skewness parameters were initialized to zero. We ran 3 chains of 400 iterations for each simulated dataset, where we kept 200 iterations per chain for the posterior sample and the first 200 iterations were considered to ensure convergence (i.e., warm-up = 200 iterations). We obtained a posterior sample of L = 600 replications for each simulated dataset, where the trace plots of some simulated datasets were examined to check the convergence of the HMC algorithm based on the rstan package. In its initial conception, the latest package considered neither the SMSN family nor the SNP distribution of a random terms specification. However, we refer to the standard N-NLMM stan code to write the specific flexible stan codes by adding the required parameters following their stochastic representations (3.6 and 3.13).

We assessed how different priors affect the model convergence rate, parameter estimation, and computation time. We relied on the effective sample size $\hat{n}_{ESS}$ and its relative Gelman-Rubin statistics (classic split-Rhat) $\hat{R}$ to assess the impact of the prior information on the speed of convergence. Then, we computed the mean of the effective sample sizes and $\hat{R}$ over the M simulated datasets for each population parameter $\mathit{\boldsymbol{\theta}}$ under each prior scenario. In each simulation condition, the larger the mean of $\hat{n}_{ESS}$ quantities are, the better the posterior distribution is estimated, and an $\hat{R}$ indicator close to 1 shows a good convergence of the chains.

For the assessment of the impact of the prior information on the model estimation, we obtained the parameter estimate empirical standard error ($\widehat{SE}_k = [\hat{\theta}_k^{(r)} - E\theta_k]$), the average standard error ($\widehat{ASE}_k$), the mean squared error ($\widehat{MSE}_k$), and the recovery ratios ($RR_k$) for each parameter based on the converged simulation replications (see Table 2). The recovery ratios ($RR$) represent ratios of the estimate/population, where values closer to 1.0 indicate a proper logistic parameter ^[3], while the recovery quantity $CR_{(\widehat{HDI})_k}$ is the proportion of datasets for which the simulation mean value is included in the Bayesian highest posterior density interval ^[26,27]. For the $\widehat{MSE}_k$, the values far from 0.00, indicated that the population values were not properly recovered. Additionally, we recorded the estimation time (in minutes) for each replication, and then computed the average estimation time (AET) for all the converged replications.

4.3.   Simulation result

The sensitivity to the prior information was assessed, and we checked the conventional convergence diagnostics tools provided by Stan software. First, we validated the estimation of Bayesian flexible multilevel nonlinear models parameters using the HMC algorithm through a simulation study. We assessed the estimators accuracy in terms of the average standard errors (ASE) and the mean of squared errors (MSE), and handled incertitude using ratios of the estimate/population, where values closer to 1.0 indicate proper epidemic parameters recovery. The Bayesian flexible multilevel nonlinear models showed different sensitivities compared to the standard normal nonlinear mixed model.

Figure 5 shows the convergence appreciation based on the mean effective sample size for BFMNLM curve modeling with different nonlinear parameter priors. This figure clearly shows that the levels of non-informativeness of the priors harm the model convergence. The mean effective sample size over the 200 simulated datasets of parameters were varied (Figure 5). For parameters $\phi_{1}$, $\phi_{2}$, $\phi_{3}$, and $\phi_{4}$ fixed effects knowing the informativeness level, the more accurate the prior was, the faster the convergence was, translating into higher effective sample sizes under all models (for SMSN-NLMM: $mean(\hat{n}_{ESS}(\beta_1))$ = (136, 90, 8); $mean(\hat{n}_{ESS}(\beta_2))$ = (487,180, 45); $mean(\hat{n}_{ESS}(\beta_3))$ = (116, 28, 11); $mean(\hat{n}_{ESS}(\beta_4))$ = (292, 46, 21) for 1SD10, 1SD20 and 1SD50, respectively). On the other hand (knowing prior accuracy level), the more informative the prior was, the faster the convergence was, translating into higher effective sample sizes for SMSN-NLMM ($mean(\hat{n}_{ESS}(\beta_1, \beta_2, \beta_3, \beta_4))$ = (136,487,116,292)) under the highly informative prior (1SD10) against (146,196,170,179) under the weakly informative prior (3SD10) and (19, 83, 83, 56) for a highly weak informative prior while the contrast should be concluded on the standard multilevel nonlinear model (N-NLMM). The less informative the prior was, the faster the convergence was translating into higher effective sample sizes ($mean(\hat{n}_{ESS}(\beta_1, \beta_2, \beta_3, \beta_4))$ = (450,556,454,315)) under the non-informative prior (5SD10) against ($mean(\hat{n}_{ESS}(\beta_1, \beta_2, \beta_3, \beta_4))$ = (168,291,129, 87)), under the weakly informative prior (3SD10), and (123,321, 84, 33) under the highly informative prior (1SD10). The informative prior seemed to be the most favorable to good convergence for flexible multilevel nonlinear models, as also confirmed by Gelman Rubin statistic ($mean(\hat{R}(\beta_2, \beta_3))$ = (1.01, 1.09)) under the highly informative prior (1SD10) against ($mean(\hat{R}(\beta_2, \beta_3))$ = (1.20, 1.10)) and under the weakly informative prior (3SD10). However, the weakly (or non) informative prior seemed to be the most favorable to good convergence for the standard multilevel nonlinear model as $mean(\hat{R}(\beta_2, \beta_3))$ = (0.99, 1.02) under the non-informative prior (5SD10) against ($mean(\hat{R}(\beta_2, \beta_3))$ = (1.02, 1.61)), under the weakly informative prior (3SD10), and (1.03, 1.48) under the highly informative prior (1SD10).

For each model, the computation time remained stable from one prior scenario to another (Figure 4), thus suggesting a low impact of the prior on the time consumption. It's noteworthy that the SMSN-NLMM is about two time more time-consuming in all scenarios compared to SNP-NLMM and N-NLMM. The average estimation time (min-max) in minutes was 59.7 (55–72) of all scenarios combined for SMSN-NLMM, 36.3 (31–41) for SNP-NLMM, and 28.3 (20–34) for standard multilevel nonlinear model.

Tables 3–5 contain two pieces of information aimed at quantifying the difference between the parameter estimates and the population value. First, the differences between the estimated nonlinear parameters proportions and the population proportions are presented in terms of ratios, which were defined earlier. Next, the difference between the estimated model parameters and the population values was quantified in terms of the MSE and the ASE where the closer the value was to zero, the more appreciable the estimates.

Figure 6 shows the parameters estimates for different models along with the true parameters values, while Tables 3–5 present corresponding average standard error ($\widehat{ASE}_k$), the $\widehat{MSE}_k$, and the $RR$ for all models parameters. The first general result is that under the accurate informative prior, parametric flexible multilevel nonlinear models (SMSN-NLMM) recover better epidemic nonlinear parameters, as shown with values closer to the true simulation values (Figure 6). Tables 4 and 5 present results for the 3SD and 5SD mean hyperparameter levels, respectively, and show that the nonlinear parameters were overall poorly recovered with SMSN-NLMM.

However, the SNP-NLMM had facility in recovering epidemic parameters in case of such inaccurate informative priors (cases of 50% variance hyperparameters) or weakly informative priors (cases of 3SD and 5SD mean hyperparameters), where nonlinear parameters were consistently underestimated (Figure 6). Given the same mean hyperparameter setting, the uncertainty in estimates increased with the variance hyperparameters values ($RR_{\beta_1} = (0.797, 0.789, 0.452)$; $\widehat{ASE}_{\beta_1} = (0.144, 0.639, 1.538)$ for 10%, 20%, and 50% variance, respectively). However, the more the mean hyperparameter value is far from the true value, the more the estimates uncertainty is pronounced ($RR_{\beta_1} = (0.8, 0.3, -0.1)$; $\widehat{MSE}_{\beta_1} = (0.02, 0.30, 0.52)$ for 1SD, 3SD and 5SD mean respectively with 10% variance). These specific results are as indicated in the case of SMSN-NLMM for the SNP-NLMM and N-NLMM models (see Tables 3–5).

When considering an increased number of clusters, the results in Table 6 show the outperformance (estimates) of SMSN-NLMM on other models when more time consumption as known. Meanwhile, the results for varying the mean and variance hyperparameters are as the case of six clusters (SMSN-NLMM: $RR_{\beta_2} = (1.331, 1.505, 1.861)$; $RR_{\beta_3} = (0.865, 1.020, 0.969)$; SNP-NLMM: $RR_{\beta_2} = (1.671, 1.791, 2.067)$; $RR_{\beta_3} = (1.064, 1.192, 1.109)$; N-NLMM: $RR_{\beta_2} = (1.279, 1.496, 2.204)$; $RR_{\beta_3} = (0.723, 0.797, 0.478)$ for 10%, 20%, and 50% variance, respectively).

5.   Analysis of COVID-19 pandemic in Francophone West Africa

As described earlier, we considered the case of seven Francophone West African countries while estimating the average behavior of countries in the population and the variability among and within them during the first 7 months of the pandemic based on the studied approach.

Our modeling strategy was to fit the studied flexible multilevel nonlinear models and present epidemic interpretations of the results. We referred to our simulation result and Figure 2 to highlight that the flexible multilevel nonlinear models (SMSN-NLMM) with informative prior is more suitable model to the considered COVID-19 data. The code for the application of the proposed Bayesian flexible multilevel nonlinear models, with other necessary r codes, stan codes, and the considered COVID-19 data in the case of this study, is available through the GitHub link https://github.com/mathildeadeoti/BFMNLM-and-Prior. Figure 7 shows the fitting output, while the Table 7 presents keys epidemiological statistics that indicate the dynamics of COVID-19 in the seven west African countries and per country.

The figure confirmed how likely the model fits the data trajectory better by considering the inter-variability of countries. Table 7 includes the estimated total number of cases ($TEC:Total\_Est\_Cases$) at the end of the pandemic, the estimated peak time ($EPT:Est\_Peak\_Time$), and the basic reproduction numbers ($R_0$) based on flexible multilevel nonlinear models (SMSN-NLMM) with the informative prior. It showed a great heterogeneity in the region in terms of epidemic dynamic ($Total\_Est\_Cases$, $Est\_Peak\_Time$, etc.). This is clearly understandable as the population density and testing efforts were not homogeneous ^[28]. In most countries, the model estimated an average of one new case of infection caused by an infected individual during the infectious period ($R_0 \approx 1$). The estimated total number of cases ($Total\_Est\_Cases$) serves as a reference to determine whether one country is probably at the end of the pandemic or not, while $Est\_Peak\_Time$ indicates whether the worst phase has passed.

The results show that the final sizes of the epidemic were estimated to be relatively larger in Senegal, Ivory-Coast, and Guinea compared to all other countries with an estimated number of 16,607, 16,297, and 10,111 cases, respectively, before expecting the end of the pandemic. For the remaining countries, namely Burkina-Faso, Mali, Benin, and Niger, the epidemic reached its final size at 4798, 4255, 3846, and 2644 reported cases, respectively. The peak time was estimated at 74, 90,105,110,118, and 127 days from the first detected case date for Mali, Niger, Benin, Ivory-Coast, Senegal, and Guinea, respectively, whereas it was relatively larger for Burkina-Faso (around 153 days). This means that the epidemic progression is relatively slower in the latter country compared to the former (as confirmed in the fitted trajectories on Figure 7). In general, the estimated peak times range along with the observed true peak time, while their estimated reported peak sizes accounted for less than $1\%$ of their population size, on average (Table 7, Figure 7). This underlined that all the considered countries would reach their peak with less than $1\%$ of their population being affected.

The basic usefulness of the approach was confirmed from the estimates of the parameters ($\hat{\phi}_{k}$) of the logistic dynamic of the model. The result corroborate the assumption on the data, as the estimated $\hat{\phi}_{4} \gg 1$ indicates a right skewness for the West African countries pandemic data, thus suggesting that the increased phase occurs much faster than the decrease phase in many of the countries.

6.   Discussion

While compartmental models and other statistical tools are most commonly used in epidemic modeling ^[16,29,30], NLMMs have received particular attention due to their flexibility in handling the heterogeneous and unbalanced repeated measures data that can arise in infectious disease modeling. Although flexible multilevel nonlinear curve modeling has been proposed, especially in the Bayesian approach, the effect of the parameter prior distribution has not been fully studied.

A sensitivity analysis procedure relatively to the specification of priors can assess the sensitivity of the posterior densities to the choice of prior distributions ^[31]. Informative priors can have an impact on estimates, even if little methodological research has directly focused on this matter, as previously reported ^[3]. Basically, this process is a means to investigate the robustness of a model by systematically changing priors to assess the impact of those changes. That is, this paper simply reconsidered the proposed Bayesian multilevel nonlinear models by making reasonable modifications to the prior assumptions, recomputing the posterior quantities of interest, and observing whether they have changed in a way that significantly affects the resulting interpretations or conclusions.

Our results corroborate with the one of ^[32] for the fact that outliers influence the computation time. Therefore, the proposed model based on the scale mixture of skew normal distribution is more time consuming than the standard. Additionally, the results clearly showed that the levels of non-informativeness of the priors harmed the model convergence. For logistics parameters $\phi_{1}$, $\phi_{2}$, $\phi_{3}$, and $\phi_{4}$ fixed effects knowing informativeness level, the more accurate the prior, the faster the convergence translated into higher effective sample sizes for the studied flexible multilevel nonlinear model (SMSN-NLMM). Therefor, the informative prior seemed to be the most favorable to a good convergence for flexible multilevel nonlinear models as confirmed by the Gelman Rubin statistic, while the weakly (or non) informative prior would be the most favorable to a good convergence for either the semi-nonparametric nonlinear mixed model or the standard multilevel nonlinear model.

The analysis of such complex epidemic data requires more consideration due to some inherent features: countries at different stages of pandemic, outliers in the response, and the specific case of Francophone West African countries considered in this study, with several consecutive days with zero detected cases mainly due to their test capacity. As a flexible multilevel modeling, the fitted curve (Figure 7) showed an increased sensitivity to outliers and variations of the epidemic stages; this can either positively or negatively influence the analysis or interpretations. It's advisable to consider estimates and computable key epidemic parameters for rational analyses of the fitting trajectories.

One of the advantages of the Bayesian approach is that the uncertainty in all parameter values is taken into account through the use of prior distributions ^[33]. However, priors have also been pinpointed as being one of the main drawbacks to this framework as a result of the inherent subjectivity that is coupled with choosing prior distributions and the corresponding hyperparameters. Meanwhile, the actual impact of prior distributions is an important research topic to explore within Bayesian estimation. Although (non) informative priors can have a large impact on the estimates (see, e.g., ^[34]), little methodological research has directly focused on this matter. Most researches explore a feature of the sample size used for an estimation in the data analysis within any estimation framework. However, this issue of sample sizes is closely tied to the use of prior distributions within the Bayesian estimation. Research has indicated that even noninformative priors can impact the estimates using larger sample sizes. Then, it is important to never underestimate the impact of a prior distribution, regardless of the sample size, as previously underlined ^[35]. Meanwhile, this paper endorsed the use of the proposed methodology that jointly accommodated the different stages of the diseases while borrowing information of the different time series to provide a more robust and reliable fit and prediction. With such a joint accommodation, we strongly believe that this methodology of Bayesian flexible multilevel nonlinear models can provide more reliable and meaningful predictions that allow policymakers around the world to make effective decisions.

7.   Conclusions

Bayesian flexible multilevel nonlinear models are powerful tools to analyze infectious disease data with asymmetric and unbalanced structures, such as varying epidemic stages across countries. This study investigated how varying levels of prior informativeness can influence the model convergence, parameter estimation, and computation time in a Bayesian FMNLM. The key finding from this study is outlined as follows:

- Outliers influence the computation time, thus leading to more time consuming of flexible multilevel nonlinear model, though the levels of informativeness of the priors do not harm model computation time;

- The results clearly show that the levels of non-informativeness of the priors harm the model convergence;

- Knowing the informativeness level, the more accurate the prior, the faster the convergence for the studied flexible multilevel nonlinear model (SMSN-NLMM);

- The weakly (or non) informative prior is the most favorable to a good convergence of either the semi-nonparametric nonlinear mixed model or the standard multilevel nonlinear model;

- Under an accurate informative prior, the parametric flexible multilevel nonlinear models (SMSN-NLMM) recover the epidemic nonlinear parameters better;

- Giving the same mean hyperparameter setting, the uncertainty in estimates increases with the variance hyperparameters values;

- The cluster size or number of countries/areas have no impact on the estimates of Bayesian FMNLM.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

OMA acknowledges the support from Deutscher Akademischer Austauschdienst German Academic Exchange Service (DAAD) through the programme In-Country/In-Region Scholarship. OMA is also grateful to the Centre d'Excellence d'Afrique en Sciences Mathématiques, Informatique et Applications (CEA-SMIA) for supporting his works. RGK acknowledges the Sub-Saharan Africa Advanced Consortium for Biostatistics (SSACAB) Phase Ⅱ.

Author contributions

OMA and RGK conceptualized the problem. OMA and AD designed the Methodology. RGK and AD administrated the project. OMA did the resource and software handling. RGK and AD supervised the project. The original draft writing and editing involved all the authors.

Conflict of interest

The authors declare there is no conflict of interest.