text:
\section{Introduction}

When it comes to the global burden of disease, the weight falls most heavily on the shoulders of low-income countries. As measured by Disability Adjusted Life Years (DALYs) lost, estimated global rates range from 40,000 to 70,000 DALYs per 100,000 individuals across low-income countries. On the other hand, rates in most developed countries tend to fall between 10,000 and 30,000 DALYs per 100,000 individuals \citep{GBDCN2017, owidburdenofdisease}. Furthermore, many deaths in low-income countries occur without registration, recording, or notice by the health system \citep{nichols20182016}. Hospitals, community health workers, and public health planners are hindered in their ability to treat new patients, allocate resources, and plan for the future when the landscape of cause of death (COD) is poorly mapped.

Medical certification of cause of death without an autopsy is difficult in most cases. Unfortunately, performing an autopsy is often infeasible or impossible. Decedents often pass away in the home. If deaths occur in the hospital, the next of kin may not agree to a full or even minimally invasive autopsy. If the next of kin does agree, the hospital may lack the resources to perform the autopsy, or the findings may be inconclusive. In cases where no medical certification occurs, verbal autopsy (VA) offers a practical alternative approach for assessing cause of death. VA involves exploring the signs and symptoms (hereafter referred to only by symptoms) a decedent experienced before death by structured interview with a relative or caregiver of the deceased. Questions can include details such as ``For how long was [decedent name] ill before they died?'' and ``Did [decedent name] have a fever in the three weeks leading up to death?'' The 2016 World Health Organization (WHO) VA instrument \citep{va2016}, available online at \url{https://www.who.int/healthinfo/statistics/verbalautopsystandards/en/}, offers a standardized form by which VA data may be recorded.

One option for analyzing these VA records is to have physicians go through and decide on a COD for each decedent. However, this process is time consuming and takes resources away from existing patients in already resource-limited settings. In addition, such expert labeling may differ when multiple physicians are presented with the same data. An alternative approach is computer-coded VA, in which the COD is assigned via an algorithm or probabilistic model. Existing computer-coded VA algorithms include those for which the relationship between symptoms and cause of death is encoded by experts (InterVA \citep{byass2012strengthening} and InSilicoVA \citep{mccormick2016probabilistic}) and those for which it is learned by relying on a labeled subset of the data having known COD (the King and Lu method \citep{king2008verbal}, the Tariff method \citep{james2011performance}, the Simplified Symptom Pattern method \citep{murray2011simplified}, the naive Bayes classifier \citep{miasnikof2015naive}, the Bayesian factor model \citep{kunihama2018bayesian}, and latent Gaussian graphical model \citep{li2018using}).

The majority of the work on VA algorithms relies on a conditional independence assumption \citep{mccormick2016probabilistic, byass2012strengthening, james2011performance, miasnikof2015naive}. That is, assume that once the cause of death is known, the knowledge that an individual had symptom $A$ should not impact belief that an individual had symptom $B$ for ANY combination of symptoms $A$ and $B$. Recent work \citep{li2018using, kunihama2018bayesian} has shown that this assumption is not valid in general. Logically, this finding is unsurprising; for example, the knowledge that someone had difficulty breathing would increase belief that they had a cough even given that their COD is already known.

\citet{li2018using} and \citet{kunihama2018bayesian} have relaxed the conditional independence assumption by probabilistically addressing conditional associations of symptoms given causes. The Bayesian latent Gaussian graphical model of \citet{li2018using} assumes the conditional dependence structure is common across causes, i.e., the symptom-level correlation does not vary with cause. On the other hand, the latent factor model of \citet{kunihama2018bayesian} models the symptom-level association independently for each cause. Neither allows for covariates (e.g. age, season, malaria endemicity of region) to be included in the model unless they are treated as ``symptoms'' themselves. The capability to allow covariates to affect cause assignment is potentially quite useful. For example, the information that an individual had balance and memory issues may be highly informative for certain CODs given an individual is young, but be uninformative among elder decedents for whom such symptoms are quite common.

The Population Health Metrics Research Consortium (PHMRC) created an open source ``Gold Standard'' VA database for training and testing VA models \citep{murray2011population}. This database uses a standardized VA questionnaire developed by the PHMRC based on WHO standards. These data include 7,840 ``adults'' (defined by both the WHO and PHMRC questionnaires to mean decedents aged 12 and older), for whom potential causes for analysis number 34 for adults (see Figure \ref{fig:phmrc_causes}). \begin{figure}[!htb]
	\centering
	\includegraphics[width=1.0\textwidth,keepaspectratio]{figs/extra/phmrc_causes.pdf}
	\caption{Number of deaths due to each of the 34 broad cause groupings in the adult PHMRC data set ($n=7840$). See the online supplementary materials for an analogous plot broken down by age.}
	\label{fig:phmrc_causes}
\end{figure} Locations include six sites across four countries, and data were collected from 2007-2010. The ensuing examination of the PHMRC data illustrates the importance of allowing both the mean and the covariance to vary with both COD and covariate information, highlighting the importance of the methodological development in this article.

Figures \ref{fig:symptom_mean} and \ref{fig:symptom_cov} offer illustrative data-driven examples from the adult PHMRC data set using the five symptoms of fever, cough, chest pain, weight loss, and headaches and the three broad COD clusters AIDS/TB, cardiovascular conditions (including heart attack and other cardiovascular diseases), and injuries (including road traffic, drowning, falls, fires, homicides, suicides, and other injuries). These COD clusters encompass 11 of the 34 PHMRC-defined COD categories. Figure \ref{fig:symptom_mean} shows that the prevalence of symptoms varies significantly with age. \begin{figure}[!htb]
	\centering
	\includegraphics[width=1.0\textwidth,keepaspectratio]{figs/1_motivate_corr_meanplot_multiple_6by7p.pdf}
	\caption{The average proportion of symptoms present varies across both age and cause of death (COD) in the adult PHMRC data. Each subplot corresponds to the proportion of ``Yes'' responses given by interviewees when asked about each symptom with proportion point estimates (points) and 95\% confidence intervals (bands) shown. Responses of don't know or refuse to answer are omitted, and individuals having unknown age (123) are omitted from age-specific subplots. A more detailed figure showing the number of observations contributing to each data point is provided in the online supplementary materials.}
	\label{fig:symptom_mean}
\end{figure} For example, across all causes (seen in the bottom row of the figure) the prevalence of the interviewee reporting the decedent experienced weight loss increases with decedent age. However, for AIDS/TB deaths the prevalence of reported weight loss remains relatively stable across ages, and is much higher than in the full population. Estimates are more precise for CODs having more observations.

Figure \ref{fig:symptom_cov} illustrates that the association between the set of example symptoms varies with both age group and COD. \begin{figure}[hbpt]
	\centering
	\includegraphics[width=1.0\textwidth,keepaspectratio]{figs/1_motivate_corr_assocplot_multiple_6by7l.pdf}
	\caption{The association between symptoms varies across both age and cause of death (COD). Each subplot corresponds to the association between pairs of select symptoms in the adult PHMRC data as measured by Yule's Q (possible range $-1$ to $1$). Rows correspond to COD groups and columns correspond to age groups. ``Don't know'' or ``refused to answer'' responses were omitted for the Yule's Q calculations, and individuals having unknown age (123) are omitted from age-specific subplots. A more detailed figure showing the number of observations contributing to the calculations in each subplot symptom-by-symptom square is provided in the online supplementary materials.
	\label{fig:symptom_cov}
\end{figure} Across all CODs (seen in the bottom row of the figure), the association between symptoms weakens with age. Using an extension of McNemar's test \citep{zhao2014testing}, the null hypothesis of homogeneous age group effects is rejected at the $0.05$ level for all $10$ of these illustrative symptom pairs (and for $72\%$ of the 14432 testable symptom pairs in the data set, i.e. those pairs of symptoms for which at least two of the category-level $2\times2$ tables have all nonzero margins-- results not shown, but code available in the online supplementary materials). This finding could be due to the fact that older people tend to have more ailments in general (imagine asking a grandmother about her everyday symptoms), while younger people likely experience clusters of symptoms relating to a specific illness. Across all ages (seen in the far right column of the figure), there are differences in associations between symptoms by cause. Again using the test of \citet{zhao2014testing}, the null hypothesis of homogeneous COD effects is rejected at the $0.05$ level for all $10$ of these illustrative symptom pairs (and for $88\%$ of all 14065 testable pairs in the data set). In summary, the data presented in Figures \ref{fig:symptom_mean} and \ref{fig:symptom_cov} show that both symptom prevalence and pairwise association vary by decedent age.

To the authors' knowledge, no existing methods consider potential covariate impacts on the prevalence or association of symptoms. Furthermore, few move beyond the conditional independence assumption (that is, that the probability of observing symptom are independent given COD). Only two current works explicitly model the conditional covariance between symptoms: \citet{li2018using} and \citet{kunihama2018bayesian}.The former models one common covariance structure across causes, while the latter models cause-specific symptom-level association via a latent factor model. The FActor Regression for Verbal Autopsy (FARVA) model described in the current paper has the most in common with \citet{kunihama2018bayesian}, which also models symptom-level association via a latent factor model. FARVA extends \citet{kunihama2018bayesian} in several key ways. It defines both the mean and the covariance of the latent variable associated with the symptom vector hierarchically so as to share information while still allowing flexibility across disparate causes. Furthermore, it allows the mean occurrence of symptoms and covariance between symptoms to depend on individual-level covariates (e.g., age). Section \ref{previous_work} discusses existing methods and distinct advantages of FARVA relative to the model of \citet{kunihama2018bayesian}. Sections \ref{prop_model} and \ref{sampler} describe the FARVA model and sampler. Section \ref{sim} describes a simulation experiment designed to separate aspects of performance of existing VA algorithms, and the model testing performed using real data. Section \ref{results} reports performance of select methods using simulated and PHMRC data with regard to CSMF accuracy and top cause assignment accuracy. Section \ref{disc} discusses future directions for this model and for the field of probabilistic COD assignment using VA data.

\section{Model}\label{model}

Let $y_i \in \{1,\ldots,C\}$ denote the underlying cause of death of person $i, i=1,\ldots,N$ and $\bm{s}_i = (s_{i1},\ldots,s_{iP})'$ be a vector of symptoms as measured by a verbal autopsy instrument. Consider the model:
\begin{equation}
\displaystyle \pi(y_i=c| \bm{s}_i) = \frac{\pi(\bm{s}_i | y_i=c)\pi(y_i=c)}{\Sigma_{c^{*}=1}^C \pi(\bm{s}_i | y_i=c^{*})\pi(y_i=c^{*})}.
\end{equation}
Here the probability that the COD of person $i$ equals $c$ is modeled conditional on the reported symptoms. The predictive goal of the model is to infer COD in cases where only symptoms are observed and the true cause is unknown. Using the above model, one can also calculate the posterior cause specific mortality fraction (CSMF, the population proportion of deaths attributed to each cause; see equation (\ref{popdist}) in the Appendix) for the population of interest. The modeling choices remaining to be made are in how to treat the prior distribution over causes $\pi(y_i)$ and how to model the likelihood of symptoms given causes $\pi(\bm{s}_i | y_i).$

\subsection{Previous work}\label{previous_work}

Often in regions in which VA methods are adopted, there is little knowledge about the population CSMF. Even in cases where clinical cases in the living are studied, causes of disease are difficult to pin down \citep{crump2013etiology}. Furthermore, it is not necessarily known which diseases are most contributory to deaths, especially outside of the hospital setting. Thus the major distinguishing feature in the existing body of work on VA algorithms is in the treatment of $\pi(\bm{s}_i | y_i),$ rather than that of $\pi(y_i)$.

Three commonly cited non-probabilistic methods are the Tariff method \citep{james2011performance}, the King-Lu (KL) method \citep{king2008verbal}, and the Simplified Symptom Pattern (SSP) method \citep{murray2011simplified}. The Tariff method is a score-based system that uses a heuristic (that tends to do well) rather than a probabilistic formulation. Tariff gives a score  to each combination of cause and symptom. These scores are used to assign a value to each possible death-cause combination, and this value is then used to rank CODs for each individual. An improved version of the Tariff method \citep{serina2015improving} is the foundation for the SmartVA-Analyze Application, available at \url{http://www.healthdata.org/verbal-autopsy/tools}. The KL method is designed to only estimate the CSMF and uses the assumption that the probability of symptoms given causes is the same in both the set of data for which COD is known (usually referred to as the training set) and that for which COD is unknown; it cannot infer individual cause of deaths (CODs). The closer this assumption is to the truth, and the more training data are available, the better this method will do. The SSP method uses the CSMF calculated from the KL method and calculates individual CODs using averages for the probability of symptoms given causes calculated across multiple random draws of symptoms. In each non-probabilistic method, some set of decedents having known COD are required for the model to learn a mapping between patterns of symptoms and COD.

InterVA was the first VA model to frame the relationship between symptoms and causes in a more probabilistic light \citep{byass2012strengthening}. Rather than relying on some training data for which COD is known, InterVA relies on a matrix of physician-generated scores assigned to the probability of observing each symptom given each cause (denoted by the $\text{P}(s|c)$ matrix, referring to the probability of symptom s given cause c). For example, physicians were asked to rate the probability that someone would answer ``Yes'' to ``Did the patient have a cough'' given that they died of TB, stroke, a motorcycle accident, etc. This matrix is useful because it allows models to be run in settings in which no training data are available. It is problematic because the values in this $\text{P}(s|c)$ matrix are likely not generalizable to all settings, are difficult to elicit (much of physician experience is likely in clinical cases and not deaths, and often the cause is not known), or may be internally inconsistent (because there are so many entries to fill out, it is very easy to, e.g., rank the probability of a patient having a cough for at least a month as higher than the probability of the patient having a cough for at least a week, even though the former case is a subset of the latter). InterVA models the probability of symptoms given causes as conditionally independent for all symptoms and causes.

There are three main issues with InterVA. First, the score used by InterVA to assign causes is based only on the presence of symptoms; it disregards symptoms when they are absent. Second, it is not possible to quantify uncertainty because the model does not contain features that are allowed to vary probabilistically. Third, InterVA is unable to incorporate physician-coded VA cases (i.e., gold standard data) into its algorithm. InSilicoVA (McCormick, 2016) is an extension of InterVA in both name and spirit. It addresses the three main issues with InterVA and provides a valid probabilistic framework in which the main ideas of InterVA can live. That is, InSilicoVA returns valid uncertainty estimates for both the COD and CSMF values. However, it still relies on the conditional independence assumption of InterVA. An even simpler framing for VA modeling directly uses a naive Bayes classifier (NBC) to assign the probability of a death being due to a given cause for each individual and each cause \citep{miasnikof2015naive}; the COD with the highest probability is then treated as the cause. Again, this framing relies on the conditional independence assumption.

Two more recent models that have relaxed the conditional independence assumption are the latent Gaussian graphical model of \citet{li2018using} and the factor model of \citet{kunihama2018bayesian}. Both models allow the prevalence of symptoms to vary by cause. The model of \citet{li2018using} assumes a shared covariance matrix between symptoms across all causes. On the other hand, the model of \citet{kunihama2018bayesian} defines separate covariance matrices for each cause. The former has the advantage of allowing for mixed data (i.e., not just binary indicators), but neither model allows for demographic, spatial, or temporal information to inform the prevalence of symptoms nor their covariance.

\subsection{Proposed model}\label{prop_model}

FARVA is a probabilistic model for VA data in which some decedents have known COD (i.e., there exists some labeled training data) and there is interest in learning about individual COD, population CSMF, and/or the mean and covariance structure of responses in the symptom questionnaire. The explicit goals of the FARVA model proposed here are to capture dependence of symptoms given a cause; share information across causes via hierarchical modeling to improve estimates associated with causes having few observed deaths; allow both the conditional prevalence and the conditional association between symptoms to vary with covariates (e.g., age of patient, time of year, geographic region); probabilistically predict cause of death for a new individual given their symptoms; and improve on COD and CSMF estimation relative to current state-of-the-art VA algorithms. A graphical representation of the model is shown in Figure \ref{fig:graph_model}, the parameters are summarized in Tables \ref{tab:model_terms1} and \ref{tab:model_terms2}, and a toy example to explain how the model works in lower dimensions is described in the online supplementary materials.

\subsubsection{Prior over causes}\label{prior}

Using a Bayesian approach to characterize uncertainty in the proportions of deaths in each COD category, let:
\begin{equation}
\{\Pr(y_i=1), \ldots, \Pr(y_i=C)\} \sim \text{Dirichlet}(a_1,...,a_C).
\end{equation}
Under the assumption that little is known about the CSMF in the region of interest and the number of possible causes is high relative to the number of common causes, set $a_1 = \ldots = a_C < 1.$ This choice leads to equal probability of each cause a priori (i.e., $\pi(y_i=c) \propto 1 \ \forall \ c \in \{1,\ldots,C\}$ where $\pi(y_i)$ is a probability distribution over causes of death) and encourages concentration near sparse subvectors. This assumption is likely valid in real applications; e.g., the gold standard PHMRC data set exhibits many causes having few observed deaths (see Figure \ref{fig:phmrc_causes}). However, this prior can easily be modified to incorporate prior knowledge about the CSMF or to encourage uniformity in the CSMF. For example, if a pilot study performed in the region of interest found counts $n_c$ attributable to COD $c, c=1\ldots C$, then an option for an informed prior is $a_1 = n_1, \ldots = a_C = n_C.$ If uniformity was desired, one could set $a_1 = \ldots = a_C = 1.$

\subsubsection{Likelihood of symptoms given a cause}\label{likelihood}

Recall that the likelihood takes the form $\pi(\bm{s}_i | y_i),$ i.e. the distribution of the symptom vector $\bm{s}_i $ is conditioned on COD $y_i$. Throughout the likelihood description, the motivation behind each modeling decision and a short description of said decision will be stated in bold preceding each mathematical description.

\textbf{Motivation: Symptom data includes questions on duration of symptoms in addition to their presence/absence. Therefore, allow the model to encompass data of mixed type.} In order to allow the FARVA model to encompass data of binary, continuous, count, and categorical type, in the specification of $\pi(\bm{s}_i | y_i)$ define:
\begin{align}
s_{ij} &= f_j(z_{ij}), \quad  j=1,\ldots,P,
\end{align}
where $\bm{z}_i=(z_{i,1},\ldots,z_{iP})'$ is a vector of latent continuous symptoms. The particular link function $f_j$ depends on the symptom specification, allowing for mixed scale data via selection of an appropriate link by scale and type. Let $f_j(z_{ij}) = z_{ij}$ or $f_j(z_{ij}) = \text{log}(z_{ij})$ for continuous $s_{ij}$, with the latter chosen for strictly positive and positively skewed cases. Scale continuous symptoms (after a log transformation if applicable) to have unit variance. Let $f_j(z_{ij}) = 1(z_{ij}>0)$ for binary $s_{ij},$ where $1(\cdot)$ is an indicator function taking the value of 1 when the argument is true and 0 when the argument is false. Finally, for (potentially zero-inflated) count $s_{ij}$ let $f_j(z_{ij})$ be a rounding operator such that $f_j(z_{ij}) = 0$ if $z_{ij}<0$ and $f_j(z_{ij}) = k$ if $k-1 \leq z_{ij} < k$ \citep{canale2013nonparametric}. Categorical variables may be handled by choosing one category as the baseline and transforming the $T$ categories into $T-1$ binary variables indicating whether or not the categorical value for that individual took that non-baseline category value (so either none or one of the $T-1$ variables will take on a value of 1).  Other specifications of $f$ would be needed for ordered categorical data to fit within this framework, if desired.  Note that of existing models, only \citet{li2018using} allows for non-binary data, specifically continuous data.

\textbf{Motivation: It is likely that the full set of symptoms are observations of some lower dimensional syndromic state. Thus, model symptom mean and covariance structure via a smaller set of underlying (unobserved) factors.} The number of measured symptoms, $P$, may be quite large for verbal autopsies (e.g., the 2016 WHO VA instrument \citep{va2016} asks up to 253 questions for adult and child deaths). Furthermore, many of the causes of interest may have very few observed deaths (e.g., see deaths by cause in the gold standard PHMRC data set in Figure \ref{fig:phmrc_causes}). Rather than estimating $P(P+1)/2$ parameters for each of the $C$ conditional covariance matrices, FARVA introduces a novel hierarchical latent factor model for the latent continuous vector $\bm{z}_i$ underlying $\bm{s}_i.$ To begin, define:
\begin{align}\label{eq:z_i_1}
\bm{z}_i &= \Lambda_{c[i]}(\bm{x}_i)\bm{\eta}_i + \bm{\epsilon}_i, \quad \bm{\eta}_i \sim \text{N}(\bm{\psi}_{c[i]}(\bm{x}_i),I_K),  \quad \bm{\epsilon}_i \sim \text{N}(\bm{0}_P, \Sigma_0),
\end{align}
where $c[i]$ indexes the COD for person $i$, $\bm{x}_i$ are predictors (i.e. covariates) associated with person $i$, and $K$ is the number of latent factors; in practice $K$ is allowed to be unknown. Here, $\Lambda_{c[i]}(\bm{x}_i)$ is a predictor-dependent $P \times K$ factor loadings matrix, $\bm{\eta}_i$ is a predictor-dependent $K \times 1$ latent factor vector, and $\bm{\epsilon}_i$ is the $P \times 1$ independent noise vector. The noise covariance $\Sigma_0$ is a diagonal matrix having entries $\sigma_j^2, j=1\ldots P,$ with $\sigma_j^2$ fixed to equal 1 when $s_{ij}$ is binary or categorical and estimated otherwise. In the above model, $\bm{\eta}_i$ captures a set of latent syndromes, $\bm{\psi}_{c[i]}(\bm{x}_i)$ the covariate-dependent mean of these latent symptoms for COD $c[i]$, and $\Lambda_{c[i]}(\bm{x}_i)$ is a low-rank description of the symptom correlations for COD $c[i]$ with covariates $\bm{x}_i.$ Note that throughout the paper the convention $\text{N}(\text{mean}, \text{variance/covariance})$ is used in parameterizing the normally distributed terms.

Since $\bm{\eta}_i$ is not actually known, dependence is induced in the latent $z_i$ via the marginalization over $\bm{\eta}_i.$ This induced dependence in elements of $z_i$ creates the possibility for dependence to be captured in the likelihood of symptoms given a cause. Specifically, the prior induced on the latent $z_i$ by integrating out $\bm{\eta}_i$ is:
\begin{equation}
\bm{z}_i|y_i \sim \text{N}(\Lambda_{c[i]}(\bm{x}_i)\bm{\psi}_{c[i]}(\bm{x}_i), \Lambda_{c[i]}(\bm{x}_i) \Lambda_{c[i]}(\bm{x}_i)' + \Sigma_0),
\end{equation}
leading to a parsimonious representation of both the mean and the dependence between the $P$ symptoms. An alternative is to define the model as $\bm{z}_i = \bm{\mu}_{c[i]}(\bm{x}_i) + \Lambda_{c[i]}(\bm{x}_i)\bm{\eta}_i + \bm{\epsilon}_i$, with $\bm{\eta}_i \sim \text{N}(\bm{0}_K,I_K)$, i.e. with a $P$-dimensional mean regression. This direct parameterization of the mean has the advantage of allowing a more direct/interpretable way to encode prior information about conditional symptom prevalence, but in the case of large $P$ it is more computationally intensive.

\textbf{Motivation: The association between symptoms likely depends on factors such as season and decedent age, but modeling a number of covariate-dependent parameters whose size depends on the number of symptoms is computationally intensive. Define a computationally efficient structure by which to model covariate-dependent covariance structure.} In order to estimate covariate-dependent factor loadings for large $P$, decompose $\Lambda_{c[i]}(\bm{x}_i)$ to be a weighted combination of a smaller set of basis elements, as in \citet{fox2015bayesian}. Namely, let
\begin{align}\label{eq:z_i_2}
\begin{split}
\Lambda_{c[i]}&(\bm{x}_i) = \Theta_{c[i]} \bm{\xi}_{c[i]}(\bm{x}_i), \\
&\Theta_{c[i]} \in \mathbb{R}^{p \times L}, \quad \bm{\xi}_{c[i]}(\bm{x}_i) = \{\xi_{c[i],lk}(\bm{x}_i), l=1,\ldots,L, \ k=1,\ldots,K\},
\end{split}
\end{align}
where $\Theta_{c[i]}$ is the coefficient matrix with weights mapping the smaller set of predictor-dependent basis functions $\bm{\xi}_{c[i]}(\bm{x}_i)$ to the higher-dimensional loadings matrix $\Lambda_{c[i]}(\bm{x}_i).$ The number of covariate-dependent parameters to be estimated is $P(P+1)/2$ in direct modeling of the conditional covariance matrix, $KP$ in (\ref{eq:z_i_1}), and $KL$ in (\ref{eq:z_i_2}). Here, $K$ is the dimension of the subspace that is assumed to capture the statistical variability in $\bm{z}$ and $L$ is the maximum size of the basis for any given value of $K$. Typically $K\ll P$ and $L \ll P,$ yielding significant computational savings by using the formulation in (\ref{eq:z_i_2}). If the dependence structure underlying the relationship between symptoms were actually highly complex it would be expected that $K$ and $L$ would need to approach $P$ in size, negating computational advantages of the factor model framework.

\textbf{Motivation: The size of the underlying basis defining relationships between symptoms isn't actually known. Furthermore, these relationships may be difficult to estimate when the number of observed deaths of a given cause is small. Thus, learn (rather than fix) the basis size while sharing information across causes hierarchically.} To share information across causes define the entries of each coefficient matrix $\Theta_{c}, c=1,\ldots, C,$ to share a common population level mean, $\Delta,$ across causes. This hierarchical structure serves the purpose of linking factor loadings and increasing robustness when estimating the factor loadings for uncommon causes. Thus similar causes should have similar factor loadings, and loading estimates for rare causes will tend toward the overall population-level mean loadings. Sparsity is induced on the population mean parameter for each entry in the coefficient matrix via the adaptive shrinkage prior of \citet{bhattacharya2011sparse}. The hierarchical shrinkage prior on each element of $\Theta_{c[i]}$ is shown below:
\begin{align}
\begin{split}
\theta_{c[i], jl} \sim & \text{N}(\Delta_{jl}, \phi_{jl}^{-1} \tau_{l}^{-1}), \\
\Delta_{jl} \sim & \text{N}(0, \phi_{jl}^{-1} \tau_{l}^{-1}), \quad  \phi_{jl} \sim \text{Ga}(\gamma/2,\gamma/2), \quad  \tau_{l}=\prod_{h=1}^l \delta_{h}, \\
&  j=1,\ldots,P, l=1,\ldots,L,
\end{split}
\end{align}
where $L$ is a conservative upper bound on the basis dimension. The shrinkage prior of \citet{bhattacharya2011sparse} with specification suggested by \citet{durante2017note} is used on the entries of both the within-cause coefficient matrix $\Theta$ and across-causes mean matrix $\Delta.$ This prior specifies $\phi_{jl}$ to be a local precision specific to element $j,l$ and $\tau_{l}$ to be a column-specific multiplier. Specifically, this shrinkage prior pulls entries in later columns more strongly toward the mean by letting $\delta_{1} \sim \text{Ga}(d_{1},1)$ and $\delta_{h} \sim \text{Ga}(d_{2},1), h \geq 2,$ with $d_{2} > 1.$  The aim of this prior is to increase the degree of shrinkage as the column index of the matrix grows. In the case of $\Theta,$ cause-specific entries are shrunk towards the population mean $\Delta$. In the case of $\Delta,$ the population mean is shrunk towards zero. Because the precision terms $\phi_{jl}$ and $\tau_{l}$ are shared across $\theta_{c[i], jl}$ and $\Delta_{jl}$, the overall effect is shrinkage of columns of $\Theta$ toward columns of $\Delta$, which approach 0. The result of this shrinkage is an effective truncation of the number of latent factors because the corresponding rows of the basis $\bm{\xi}$ will then have insignificant effect in the definition of the resulting mean and covariance of $\bm{z}$.  A hallmark of $L$ not being large enough would be columns of $\Delta$ close to the $L$th column remaining large in spite of the shrinkage prior on its entries.

\textbf{Motivation: The association between symptoms likely depends on covariates, but this covariate dependent relationship may be difficult to learn for causes with a small number of observed deaths. Therefore, model covariate-dependent covariance structure while sharing information across causes hierarchically.} The choice of how to model each $\bm{\xi}_{c[i]}(\bm{x}_i),$ the cause-specific $L \times K$ matrix of basis functions of equation (\ref{eq:z_i_2}), to include covariates is flexible. A non-hierarchical nonparametric model can be introduced by choosing independent Gaussian process priors on entries of $\bm{\xi},$ i.e. $\xi_{c[i],lk} \sim \text{GP}(0,\sigma_{\xi}(\cdot))$ with $\sigma_{\xi}(\bm{x},\bm{x}')$ an appropriate kernel function. If the covariate included is spatial location, for example $\bm{x}=(lat, lon)'$, a squared exponential covariance kernel with unit variance could be used (the restriction that the kernel has unit variance arises due to identifiability issues with the multiplication with $\Theta$). A separable covariance structure across space and time could also be chosen. The advantage of such a formulation is the ability to model symptom dependence that varies smoothly across space and time, and to infer dependence while predicting COD at locations for which no VAs were recorded. Parametric models for each entry could also be considered by letting $\xi_{c[i],lk}(\bm{x}_i) = \bm{\beta}_{c[i],lk}^T \bm{x}_i.$ Here $\bm{\beta}$ is a $B \times 1$ vector, where $B$ is the number of covariates included in the model. This form would be suitable for, e.g., including an indicator variable for whether the $i$th person died in the hospital, whether the death occured during the malaria season, or some other relevant information. Aspects of the interviewer, e.g., level of experience, could also be included as predictors. Smoothing splines, piecewise splines, or combinations of parametric and nonparametric forms for predictors can be included, but the level of complexity it is possible and useful to include will be limited by the expressiveness and dimension of the data.

In order to borrow information in the case of small sample sizes within a given cause, introduce a hierarchical prior on elements of the predictor-dependent basis functions $\bm{\xi}$ and parametric regression parameters $\bm{\beta}$. Rather than modeling $\bm{\beta}_{c[i],lk}$ independently for each $c$, model $\bm{\beta}_{c[i],lk}$ as coming from a shared population mean parameter $\bm{\mu}_{\beta_{lk}}.$ For convenience of notation, focus on the case when the covariates of interest are the same for the mean and covariance components of the model. Define the following prior on elements $\xi_{c[i],lk}$ of $\bm{\xi}_{c[i]}:$
\begin{align}\label{eq:xi_regression}
\begin{split}
\xi_{c[i],lk}(\bm{x}_i) = &\bm{\beta}_{c[i],lk}^T \bm{x}_i,  \\% + \epsilon_{\xi_{lk}i}, \quad \epsilon_{\xi_{lk}i} \sim \text{N}(0,1) \\
& \bm{\beta}_{c[i],lk} \sim \text{N}_B(\bm{\mu}_{\beta_{lk}}, \Sigma_{\beta_{lk}}) ,  \\
& \bm{\mu}_{\beta_{lk}} \sim \text{N}_B(\mu_0, \Lambda_0), \quad \Sigma_{\beta_{lk}} \sim \text{IW}(\nu_0, S_0), \\
&  l=1,\ldots,L, k=1,\ldots,K.
\end{split}
\end{align}

\textbf{Motivation: Average symptom observations likely depend on covariates, but this covariate dependent relationship may be difficult to learn for causes with a small number of observed deaths. Model covariate-dependent mean structure while sharing information across causes hierarchically.} The choice of how to model $\bm{\eta}_{i}$, the cause-specific $K \times 1$ vector of latent syndromes of equation (\ref{eq:z_i_1}), to include covariate information is similarly flexible to that of how to model $\bm{\xi}$. As with the model for $\bm{\xi}_{c[i]}(\bm{x}_i),$ a non-hierarchical nonparametric model can be introduced by choosing independent Gaussian process priors on entries of $\bm{\psi}_{c[i]}.$ A parametric model is defined by fixing $\psi_{c[i],k} = \bm{\alpha}_{c[i],k}^T \bm{x}_i.$  This choice induces the following hierarchical formulation for elements $\eta_{i,k}$ of $\bm{\eta}_{i}$ :
\begin{align}\label{eq:eta_regression}
\begin{split}
\eta_{i,k}(\bm{x}_i) = &\psi_{c[i],k} + \epsilon_{\eta_{k}i} \\
= &\bm{\alpha}_{c[i],k}^T \bm{x}_i + \epsilon_{\eta_{k}i}, \quad \epsilon_{\eta_{k}i} \sim \text{N}(0,1) \\
& \bm{\alpha}_{c[i],k} \sim \text{N}_B(\bm{\mu}_{\alpha_{k}}, \Sigma_{\alpha_{k}}) ,  \\
& \bm{\mu}_{\alpha_{k}} \sim \text{N}_B(A_0, L_0), \quad \Sigma_{\alpha_{k}} \sim \text{IW}(v_0, D_0), \\
&  k=1,\ldots,K.
\end{split}
\end{align}

When there are no covariates, model elements of $\bm{\eta}_{i}$ and $\bm{\xi}_{c[i]}$ using a mean-only model. Finally, it is possible to include covariates in the model for $\bm{\xi}_{c[i]}(\cdot)$ that are not included in that for $\bm{\eta}_{i}(\cdot)$ or vice versa.

\begin{table
	\footnotesize
	\bgroup
	\def\arraystretch{1.2
	\begin{center}
		\begin{tabular}{||l l||}
			\hline
			Variable & Description \\ [0.5ex]
			\hline\hline
			$s_{ij}$ & Observed symptom $j$ of decedent $i$ \\
			$\bm{x}_{i}$ & Observed covariates for decedent $i$ \\
			$z_{ij}$ & Latent continuous symptom $j$ of decedent $i$ \\
			$\sigma_j^2$ & Variance for latent continuous symptom $j$ \\
			$\Lambda_{ijk}$ & $jk$-th entry of decedent $i$'s factor loadings matrix \\
			$\eta_{ik}$ & $k$-th entry of decedent $i$'s latent factor vector \\
			$\psi_{ik}$ & $k$-th entry of decedent $i$'s latent factor mean \\
			$\bm{\alpha}_{ck}$ & Regression coefficients for $k$-th entry of $\bm{\psi}$ for decedents dying of cause $c$ \\
			$\bm{\mu}_{\alpha,k}$ & Global mean of cause-specific regression coefficients for $k$-th entry of $\bm{\psi}$ \\
			$\Sigma_{\alpha,k}$ & Global covariance of cause-specific regression coefficients for $k$-th entry of $\bm{\psi}$ \\
			$\xi_{ilk}$ & $lk$-th entry of decedent $i$'s latent predictor-dependent basis functions \\
			$\bm{\beta}_{clk}$ & Regression coefficients for $lk$-th entry of $\xi$ for decedents dying of cause $c$ \\
			$\bm{\mu}_{\beta,lk}$ & Global mean of cause-specific regression coefficients for $lk$-th entry of $\xi$ \\
			$\Sigma_{\beta,lk}$ & Global covariance of cause-specific regression coefficients for $lk$-th entry of $\xi$ \\
			$\Theta_{cjl}$ & $jl$-th entry of coefficient matrix mapping $\bm{\xi}$ to $\Lambda$ for decedents dying of cause $c$ \\
			$\Delta_{jl}$ & $jl$-th entry of global mean of cause-specific coefficient matrices $\Theta_{c}$ \\
			$\phi_{jl}$ & Local precision specific to element $jl$ of the global and cause-specific coefficient matrices \\
			$\tau_{l}$ & Column-specific precision for $l$-th column of the global and cause-specific coefficient matrices \\
			$\delta_{l}$ & $l$-th component of the multiplicative gamma process forming column-specific precision $\tau_{l}$ \\ [0.5ex]
			\hline
		\end{tabular}
	\end{center}
	\caption{Summary of the parameters of the FARVA model.}
	\label{tab:model_terms1}
	\egroup
\end{table}

\begin{table
	\footnotesize
	\bgroup
	\def\arraystretch{1.2
	\begin{center}
		\begin{tabular}{||l l||}
			\hline
			Plate & Description \\ [0.5ex]
			\hline\hline
			$N$ & Decedents $i=1,\ldots,N$ \\
			$P$ & Symptoms $j=1,\ldots,P$ \\
			$K$ & Latent factors $k=1,\ldots,K$ \\
			$C$ & Causes of death $c=1,\ldots,C$ \\
			$L$ & Maximum basis size $l=1,\ldots,L$ \\
			$B$ & Covariates $b=1,\ldots,B$ \\ [0.5ex]
			\hline
		\end{tabular}
	\end{center}
	\caption{Summary of the indices (plates in Figure \ref{fig:graph_model}) of the FARVA model.}
	\label{tab:model_terms2}
	\egroup
\end{table}

\begin{figure
	\centering
	 \tikz{ %
		\node[obs] (s) {$\bm{s}_{i}$} ; %
		\node[latent, above=of s] (z) {$\bm{z}_{i}$} ; %
		\node[latent, left=of s, xshift=-1.4cm] (sig) {$\sigma^2$} ; %
		\node[latent, above=of z, yshift=0.0cm] (lambda) {$\Lambda$} ; %
		\node[latent, above right=of z, yshift=0.5cm] (eta) {$\bm{\eta}_i$} ; %
		\node[latent, above=of lambda, xshift=-1.9cm] (xi) {$\xi$} ; %
		\node[latent, above=of eta] (psi) {$\bm{\psi}$} ; %
		\node[obs, below right=of psi] (x) {$\bm{x}_{i}$} ; %
		\node[latent, left=of xi, xshift=0.3cm] (beta) {$\bm{\beta}$} ; %
		\node[latent, above=of beta, yshift=0.7cm] (alpha) {$\bm{\alpha}$} ; %
		\node[latent, above=of alpha, yshift=-0.3cm] (mualph) {$\bm{\mu}_{\alpha}$} ; %
		\node[latent, right=of mualph, xshift=-0.8cm] (sigalph) {$\bm{\Sigma}_{\alpha}$} ; %
		\node[latent, above=of xi, xshift=-0.2cm, yshift=-0.5cm] (mubet) {$\bm{\mu}_{\beta}$} ; %
		\node[latent, right=of mubet, xshift=-0.8cm] (sigbet) {$\bm{\Sigma}_{\beta}$} ; %
		\node[latent, left=of z, xshift=-1.6cm] (theta) {$\Theta$} ; %
		\node[latent, left=of theta, xshift=0.2cm] (delta) {$\Delta$} ; %
		\node[latent, above=of delta, yshift=-0.0cm] (phithet) {$\phi$} ; %
		\node[latent, above=of phithet, yshift=-0.5cm] (tauthet) {$\tau$} ; %
		\node[latent, above=of tauthet, yshift=-0.0cm] (delthet) {$\delta$} ; %

		\plate[inner sep=0.25cm, xshift=0.0cm, yshift=0.1cm] {plateC} {(beta) (alpha) (theta)} {$C$}; %
		\plate[inner sep=0.25cm, xshift=0.1cm, yshift=-0.1cm] {plateK} {(beta) (alpha) (xi) (psi) (lambda) (eta) (mualph) (sigalph)} {$K$}; %
		\plate[inner sep=0.25cm, xshift=0.0cm, yshift=0.1cm] {plateN} {(z) (s) (lambda) (eta) (xi) (psi) (x)} {$N$}; %
		\plate[inner sep=0.25cm, xshift=0.1cm, yshift=0.0cm] {plateP} {(theta) (z) (lambda) (s) (sig) (phithet) (delta)} {$P$}; %
		\plate[inner sep=0.25cm, xshift=-0.0cm, yshift=0.0cm] {plateL} {(delta) (beta) (xi) (phithet) (tauthet) (delta) (delthet) (mubet) (sigbet)} {$L$}; %

		\edge {mubet} {beta} ; %
		\edge {sigbet} {beta} ; %
		\edge {mualph} {alpha} ; %
		\edge {sigalph} {alpha} ; %
		\edge {delthet} {tauthet} ; %
		\edge {phithet} {delta} ; %
		\path (tauthet) edge[bend left, ->]  (delta) ;

		\edge {phithet} {theta} ; %
		\edge {tauthet} {theta} ; %
		\edge {delta} {theta} ; %
		\edge {theta} {lambda} ; %
		\path (alpha) edge[bend left, ->]  (psi) ;
		\edge {beta} {xi} ; %
		\edge {x} {xi, psi} ; %
		\edge {xi} {lambda} ; %
		\edge {psi} {eta} ; %
		\edge {lambda, eta} {z} ; %
		\path (sig) edge[bend right, ->]  (z) ;
		\edge {z} {s} ; %

	}
	\caption{Graphical representation of the FARVA model. Nodes shaded in gray are data (observed), unshaded nodes are parameters (learned by the model). Note $B,$ the number of covariates, is not shown but is the dimension of vectors $\bm{x}_i, \bm{\alpha}_{ck}, \bm{\mu}_{\alpha,k}, \bm{\beta}_{clk}, \bm{\mu}_{\beta,lk}$ and square matrices $\Sigma_{\alpha,k}, \Sigma_{\beta,lk}$.}
	\label{fig:graph_model}
\end{figure}

\subsubsection{Model assumptions and limitations}

The FARVA model assumes the symptom means can be specified via the linear projection of latent factors into the symptom space ($\Lambda_{c[i]}(\bm{x}_i)\bm{\psi}_{c[i]}(\bm{x}_i)$). In a practical sense, this is analogous to assuming that the full set of symptoms arise from some smaller set of underlying lower-dimensional syndromes. If symptom means do not arise from  such a process, or if the number of latent factors required to capture symptom means well is large, it is preferable to simply use a separate mean vector as discussed in Section \ref{likelihood}.

The model also assumes the covariance between symptoms can be specified by a structured low-rank matrix ($\Lambda_{c[i]}(\bm{x}_i) \Lambda_{c[i]}(\bm{x}_i)'$) and a diagonal noise matrix ($\Sigma_0$). If the rank of the structured component needs to be large, i.e. as $K$ approaches $P$, the computational benefits of the factor model specification disappear and the amount of data required to learn all of the model parameters becomes impractically high for the VA setting.

Recall that $K$ need only be specified as a conservative upper bound on the number of latent factors believed to be necessary for modeling. If the chosen $K$ is too small (diagnosable by observing that later columns of $\Delta$ do not converge to 0) model performance and convergence may suffer. If the chosen $K$ is too large (diagnosable by observing that many columns of $\Delta$, not only the final few, have converged to 0) computation time will be higher than necessary but performance and convergence should not be impacted.

Modelers face decisions as to the inclusion/specification of (1) symptoms and (2) covariates. Regarding (1) one could, for example, either exclude or include information on coughing as a symptom. If the latter, one could decide to use a binary specification (e.g., did or didn't cough) or a numeric time specification (e.g., coughed for $x$ days before death). Regarding (2) one could, for example, either exclude or include age as a covariate. If the latter, one could use a binary specification (e.g., elder vs. non-elder) or a numeric year specification (e.g., 68 vs. 22 years old) in a linear model. As noted in Section \label{likelihood} it would also be possible to include certain covariates via nonlinear models (e.g., Gaussian processes).

A possible Bayesian solution to model selection would be to compute posterior model probabilities and Bayes factors for each symptom/covariate specification of interest. These quantities rely on the calculation of the marginal likelihood, an analytically intractable integral. Approximating the marginal likelihood for FARVA, e.g. via bridge sampling, could aid in understanding and formally deciding which symptoms and/or decedent characteristics are important when modeling CODs. A simpler alternative for exploring model predictive performance under different possible sets of symptoms/covariates is to perform cross validation, in which the model is trained/tested using each input set.

\subsection{Posterior computation}\label{sampler}

The posterior for the FARVA model is not available in closed form, but it may be approximated via samples obtained from a Markov chain Monte Carlo (MCMC) algorithm. Closed form full conditional distributions of the parameters associated with the model allow the use of a straightforward Gibbs sampler for these draws. Full details on each Gibbs step, posterior computations, hyperparameter settings, and the practical considerations on how the steps are modified for real data are provided in the online supplementary materials.

\subsubsection{Performance metrics}\label{perf_met}

The performance metrics used to assess model accuracy are top cause accuracy $\text{ACC}_1$ (\ref{eq:acc}), a measure for how well CODs for individuals are predicted, and CSMF accuracy $\text{ACC}_{\text{CSMF}}$ (\ref{eq:acccsmf}), a measure for how well the population CSMF is predicted. The top cause accuracy and CSMF accuracy take possible values from 0 to 1, with 1 being best.
\begin{equation}\label{eq:acc}
	\text{ACC}_1 = \frac{\text{\# of correct COD being the top cause assignment}}{N}
\end{equation}
\begin{equation}\label{eq:acccsmf}
\text{ACC}_{\text{CSMF}} = 1 - \frac{\sum_{c=1}^C |\text{CSMF}_c^{true} - \text{CSMF}_c^{pred}| }{2(1-\text{min}(\text{CSMF}^{true}))}
\end{equation}
The CSMF accuracy formula was defined in \citet{murray2011robust}, with the idea being that the worst case scenario for CSMF prediction is to put all predicted CSMF weight on the least common cause. This scenario corresponds to a total absolute error of $2(1-\text{min}(\text{CSMF}^{true}))$. Note that these metrics all rely on there being some hold-out data set for which information on COD is known but can be hidden from the model during training. Results for chance-corrected concordance, a metric with a 0 value corresponding to the performance expected of random uniform guessing, are included in the online supplementary materials.

\subsection{Simulations}\label{sim}

Commonly, work on VA based COD algorithms jump immediately to training/testing using samples from the PHMRC data set. This approach has the advantage of allowing questions such as ``How do algorithms trained in one context and tested in another perform?'' and ``What causes are easier/harder to predict?''; it places the emphasis on performance of algorithms with real data. However, the underlying mean/covariance of symptoms in these resampled PHMRC data sets are unknown. This limitation makes it difficult to manipulate and explore how differing symptom level mean and covariance structures impact model performance. To address this question, the current work includes simulated scaled down data sets having just 928 observations, 21 ``symptoms'', and 4 ``causes''. Via different generative processes, it is possible to explore various structures of cause-specificity and covariate dependence for the simulated datasets' mean and covariance.

Table \ref{sim_table} shows detailed information about the various simulation configurations. For each simulation configuration sample 1000 data sets, a number chosen to balance the considerations of Monte Carlo error and computation time, and run FARVA, Kunihama's Bayesian factor (BF) model, open source Tariff, InsilicoVA, and the naive Bayes classifier (NBC). Broadly speaking, there are two classes of configuration:
\begin{enumerate}
	\item No covariate information included: Simulations \textit{a} through \textit{d}, with various combinations of cause-specificity in mean/covariance. FARVA is likely to only do as well as existing methods.
	\item Covariate information included: Simulations \textit{e} through \textit{g}, with various combinations of cause-specificity in mean/covariance, with covariate effects. FARVA is expected to outperform existing methods.
\end{enumerate}
Note that although simulation \textit{g} includes mixed and continuous versions of the data set for testing with FARVA, all other methods compared can use only binary data. More details of how data are simulated and visuals of the impact of various simulation settings are provided in the online supplementary materials.

\begin{table}
	\footnotesize
	\centering
	\begin{tabular}{ l c||c|c|c  }
		\toprule
		\multicolumn{2}{c||}{} &\multicolumn{3}{c}{\textbf{Artificial data structure}} \\
		\multicolumn{1}{c}{} & \multicolumn{1}{c||}{Simulation} &
		\multicolumn{1}{c}{Mean structure} &
		\multicolumn{1}{c}{Covariance structure} &
		\multicolumn{1}{c}{Data type} \\

		\midrule
		\multirow{2}{*}{\begin{minipage}{2cm}No covariate information included\end{minipage}} & \textit{a}  & S & C, I &  B \\
		& \textit{b} &  S & C, D   & B \\
		& \textit{c} & C & S, D &  B \\
		& \textit{d} & S & S, D & B \\
		\midrule
		\multirow{2}{*}{\begin{minipage}{1.8cm}Covariate information included\end{minipage}} & \textit{e} & S, V & C, I & B \\
		& \textit{f} & C & S, D, V & B  \\
		& \textit{g}$_1$ & S, V & S, D, V & B \\
		& \textit{g}$_2$ & S, V & S, D, V & M \\
		& \textit{g}$_3$ & S, V & S, D, V & T \\
		\bottomrule
	\end{tabular}
	\caption{\label{sim_table} In the table above S denotes cause-specific, C denotes common across causes; I denotes independent, D denoted non-independent; V denotes covariate-dependent; B denotes binary, M denotes mixed, and T denotes continuous. For example, simulation $d$ has cause-specific (S) mean structure and cause-specific (S) dependent (D) covariance structure, and the data type is binary (B). The simulation letter denotes that 1000 data sets are sampled for that letter. For example, simulations $g_1,$ $g_2,$ and $g_3$ are based on the same underlying data set so the same inputs are provided for each to the BF, InSilicoVA, NBC, and Tariff methods, but different inputs are given to FARVA ($g_1$ uses binary inputs, $g_2$ uses mixed inputs, and $g_3$ uses continuous inputs).}

\end{table}

Simulations using the PHMRC data are also performed. To begin, the data are read into \textbf{R} using the command \texttt{read.csv(openVA::getPHMRC\_url("adult"))} (yielding 7841 observations) and the one observation for which a decedent has age $<12$ is removed (leaving 7840 adult deaths). The data are then divided by site (Andhra Pradesh, India ; Bohol, Philippines; Dar es Salaam, Tanzania; Mexico City, Mexico; Pemba Island, Tanzania; and Uttar Pradesh, India). For each site the following is repeated 100 times: First, site-specific data are split into 75\% training, 25\% test. Then, data cleaning steps used in the \textbf{OpenVA} software were performed, i.e. all variables converted to dichotomous symptoms matching those used in InterVA algorithm. Specifically, the adult PHMRC data was processed using the \texttt{ConvertData.phmrc()} command with \texttt{cutoff = "default"} and \texttt{cause = "va34"}. Variables were then cleaned such that symptoms having a  high missingness rate (specifically, over 95\% missing) or zero variability (i.e., all respondents had the same answer) were removed from analysis. Only the subset of causes included in the training data set were included as possible causes in each analysis. Finally, each model is run, with FARVA including whether or not each decedent was an elder ($\geq 65$) as a covariate. The competitor methods included are the same as those used in the simulated data runs. The InterVA and King Lu methods were not included as competitors, as both were found by \citet{kunihama2018bayesian} to have consistently poor performance relative to state of the art.

 To explore cause-specific and potentially covariate-dependent prevalence and patterns of association between symptoms, the FARVA model is run using the full PHMRC data set (i.e., including data from across all sites and with no held-out test data) with (1) no covariates included and (2) a covariate for whether or not each decedent was an elder included. This exploration is a novel contribution to the VA literature. \citet{li2018using} described overall patterns of association between symptoms, but their model does not allow for covariance structure to differ across causes. No previous paper has explored the difference in patterns of association between symptoms across CODs (although the BF model would also allow such an analysis, none was performed in \citet{kunihama2018bayesian}). Furthermore, as no previous models have allowed for covariate inclusion, the question of how the mean or covariance structure differs within (or across) COD(s) for a given covariate has yet to be addressed for any covariate. Note that results here are meant to be illustrative rather than exhaustive, highlighting the utility of the FARVA model for exploring both the symptom level mean structure and patterns of association. Further simulation run details are provided in the online supplementary materials.

\section{Results and discussion}\label{results}

Performance of both the simulated and the PHMRC data runs are discussed below. In summary, explicitly accounting for covariate dependence in the model does improve predictive performance in both simulated and real data. The degree to which performance improves depends on the strength of the relationship between a covariate and symptom prevalence and/or associations.

\subsection{Simulation performance}

Two sets of simulations are considered. In the first, no covariate information is included or modeled (simulations $a$ through $d$). In the second, a single binary covariate is introduced as a modifier of the mean and/or covariance structure of symptoms (simulations $e$ through $g_3$). See the online supplementary materials for visual examples of the symptom mean and covariance structure from each simulation setting, and for performance as a function of degree of cause-specificity, and for plots of performance as a function of covariate dependence and proportion of data allowed to be continuous.

\subsubsection{Simulations with no covariate dependence}

Simulation \textit{a} fixes cause-specific symptom mean structure for each cause and shared \textbf{independent} covariance structure across causes. This simulation is designed to mimic a scenario in which the conditionally independent assumption is valid, and the only information about cause of death comes from symptom prevalence. The models allowing the specification of a non-independent conditional covariance matrix (FARVA and BF) perform comparably to models assuming conditional independence. Simulation \textit{b} is similar to simulation \textit{a} in that the only information on COD is found in the differing mean structure across causes. The difference is that symptoms are assumed to have a shared \textbf{dependent} covariance structure across causes. Both FARVA and BF are able to learn this shared covariance structure, and both perform slightly better than models assuming conditional independence.

\begin{table}
	\footnotesize
	\centering
	\begin{tabular}{ cc||ccccc  }
		\toprule
		\multicolumn{2}{c||}{} &\multicolumn{5}{c}{\textbf{VA model}} \\
		\multicolumn{1}{c}{} & \multicolumn{1}{c||}{Simulation} &
		\multicolumn{1}{c}{FARVA} &
		\multicolumn{1}{c}{BF} &
		\multicolumn{1}{c}{InSilico} &
		\multicolumn{1}{c}{NBC} &
		\multicolumn{1}{c}{Tariff} \\
		\midrule
		\multirow{2}{*}{\begin{minipage}{1cm}$\text{ACC}_1$\end{minipage}} & \textit{a}   & 0.524 (0.047) & 0.518 (0.047) & 0.493 (0.069) & \textbf{0.528} (0.046) & 0.50 (0.045) \\
		& {b} &  \textbf{0.532} (0.045) & 0.527 (0.045) & 0.470 (0.068) & 0.522 (0.048) & 0.498 (0.045) \\
		& \textit{c} & \textbf{0.635} (0.039) & 0.634 (0.040) & 0.268 (0.065) & 0.275 (0.035) & 0.262 (0.040) \\
		& \textit{d} & 0.711 (0.036) & \textbf{0.712} (0.037) & 0.420 (0.070) & 0.468 (0.042) & 0.444 (0.043) \\
		\midrule
		\multirow{2}{*}{\begin{minipage}{1cm}$\text{ACC}_\text{CSMF}$\end{minipage}} & \textit{a}   & 0.894 (0.013) & 0.895 (0.014) & 0.841 (0.104) & \textbf{0.908} (0.024) & 0.903 (0.027) \\
		& \textit{b} & 0.895 (0.014) & 0.898 (0.016) & 0.831 (0.104) & 0.901 (0.034) & \textbf{0.902 (0.030)} \\
		& \textit{c} & 0.913 (0.014) & \textbf{0.915} (0.015) & 0.650 (0.148) & 0.858 (0.040) & 0.808 (0.087) \\
		& \textit{d} & 0.926 (0.015) & \textbf{0.929 (0.016)} & 0.784 (0.139) & 0.896 (0.030) & 0.892 (0.031) \\
		\bottomrule
	\end{tabular}
	\caption{\label{sim_table_res_1} Above are the mean (SD) of the top cause accuracy ($\text{ACC}_1$) and CSMF accuracy ($\text{ACC}_\text{CSMF}$) for the 1000 data sets sampled for each simulation having no covariate dependence. The highest performing model is indicated by a bolded mean. Note that FARVA and BF perform quite similarly with regards to both $\text{ACC}_1$ and $\text{ACC}_\text{CSMF}$, as expected. See Table \ref{sim_table} forc data structure for each simulation setting.}
\end{table}

Simulation \textit{c} fixes a common symptom mean structure across causes and symptom covariance structure specific to each cause. In other words, all of the information differentiating CODs is found in the covariance structure. In this setting, models making the conditional independence assumption have no information by which to make a COD determination, and as expected FARVA and BF are the only models to perform better than simple uniform guessing. Simulation \textit{d} fixes a symptom mean and covariance structure specific to each cause. In other words, information about the COD may be found in both the mean and the covariance structure of symptoms. As expected, FARVA and BF outperformed the models that only utilize information about symptom means. Also, FARVA and BF outperformed their own metrics in previous simulation runs, illustrating that both the mean and the covariance have the potential to contain valuable information.

\subsubsection{Simulations with covariate dependence}

Simulation \textit{e} fixes cause-specific symptom mean structure that varies with covariates for each cause and shared independent covariance structure across causes. This simulation is designed to mimic a scenario in which the conditionally independent assumption is valid, and COD-specific symptom prevalence is modified by a single binary covariate. FARVA can utilize this covariate information to more accurately capture cause-specific symptom mean structure for all individuals, and as expected it has the highest performance. Simulation \textit{f} fixes a common symptom mean structure across causes and symptom covariance structure that varies with covariates specific to each cause. Again, only FARVA and BF are expected to do better than chance, and of the two only FARVA can utilize this covariate information to more accurately capture cause-specific symptom association structure for all individuals.  Simulation \textit{g} sets both the mean and covariance of symptoms to depend on cause and vary with a binary covariate. This simulation illustrates that FARVA does best when both the mean and the covariance depend on a covariate, and its performance improves when truly continuous data are included in the model as such.

\begin{table}
	\footnotesize
	\centering
	\begin{tabular}{ cc||ccccc  }
		\toprule
		\multicolumn{2}{c||}{} &\multicolumn{5}{c}{\textbf{VA model}} \\
		\multicolumn{1}{c}{} & \multicolumn{1}{c||}{Simulation} &
		\multicolumn{1}{c}{FARVA} &
		\multicolumn{1}{c}{BF} &
		\multicolumn{1}{c}{InSilico} &
		\multicolumn{1}{c}{NBC} &
		\multicolumn{1}{c}{Tariff} \\
		\midrule
		\multirow{2}{*}{\begin{minipage}{1cm}$\text{ACC}_1$\end{minipage}} & \textit{e} & \textbf{0.442} (0.039) & 0.397 (0.036) & 0.332 (0.079) & 0.389 (0.037) & 0.372 (0.037) \\
		& \textit{f} &  \textbf{0.578} (0.038) & 0.496 (0.038) & 0.255 (0.071) & 0.264 (0.032) & 0.258 (0.038) \\
		& \textit{g}$_1$ & \textbf{0.653} (0.035) & 0.552 (0.038) & 0.333 (0.075) & 0.392 (0.036) & 0.375 (0.038) \\
		& \textit{g}$_2$ & \textbf{0.674} (0.035) & | & |& |&  |\\
		& \textit{g}$_3$ & \textbf{0.749} (0.033)  & | & |& |&  |\\
		\midrule
		\multirow{2}{*}{\begin{minipage}{1cm}$\text{ACC}_\text{CSMF}$\end{minipage}} & \textit{e}   & 0.883 (0.010) & 0.881 (0.011) & 0.647 (0.184) & \textbf{0.884} (0.030) & 0.883 (0.031) \\
		& \textit{f} &  \textbf{0.902} (0.013) & 0.894 (0.014) & 0.647 (0.161) & 0.862 (0.037) & 0.809 (0.086) \\
		& \textit{g}$_1$ & \textbf{0.913} (0.013) & 0.902 (0.014) & 0.649 (0.180) & 0.882 (0.030) & 0.882 (0.031) \\
		& \textit{g}$_2$ & \textbf{0.915} (0.013)  & | & |& |&  |\\
		& \textit{g}$_3$ & \textbf{0.924} (0.013) & | & |& |&  |\\
		\bottomrule
	\end{tabular}
	\caption{\label{sim_table_res_2} Above are the mean (SD) of the top cause accuracy ($\text{ACC}_1$) and CSMF accuracy ($\text{ACC}_\text{CSMF}$)  for the 1000 data sets sampled for each simulation having covariate dependence in the mean, covariance, or both. Recall that different inputs are given to FARVA for the $g$ simulations ($g_1$ uses binary inputs, $g_2$ uses mixed inputs, and $g_3$ uses continuous inputs); the results for the non-FARVA models are shared across each $g$. See Table \ref{sim_table} for other information about the data structure for each simulation setting. The highest performing model is indicated by a bolded mean. The proportion of datasets for which FARVA outperformed BF for simulations $e$ through $g_3$ is 0.91, 1.00, 1.00, 1.00, and 1.00 for $\text{ACC}_1$ and 0.61, 0.82, 0.87, 0.88, and 0.94 for $\text{ACC}_\text{CSMF}$, respectively.}
\end{table}

\subsection{PHMRC predictive performance}

Top cause accuracy results for the 100 test/train splits of each site-specific PHMRC data set are shown in Figure \ref{fig:res_phmrc_ACC1}. FARVA outperforms its competitors with regards to top cause accuracy at each site. Notably, FARVA shows improvement relative the BF model in each setting.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1.0\textwidth,keepaspectratio]{figs/res_phmrc_ACC1.png}
	\caption{Top cause accuracy for PHMRC data over 100 test/train splits at each site. The centered solid point (horizontal bar) shows the mean (median) across test/train splits, for which individual performance results are shown as semi translucent horizontally jittered points. Each method shares data processing steps, but decedent age ($\geq 65$) is included in FARVA as a binary covariate. For each site (read left to right, top to bottom), the FARVA model outperforms the BF model for 88, 75, 81, 97, 61, and 99\% of the test/train splits. }
	\label{fig:res_phmrc_ACC1}
\end{figure}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1.0\textwidth,keepaspectratio]{figs/res_phmrc_CSMF.png}
	\caption{CSMF accuracy for PHMRC data over 100 test/train splits at each site. The centered solid  point (horizontal bar) shows the mean (median) across test/train splits, for which individual performance results are shown as semi translucent horizontally jittered points. Each method shares data processing steps, but decedent age ($\geq 65$) is included in FARVA as a binary covariate. For each site (read left to right, top to bottom), the FARVA model outperforms the BF model for 72, 68, 87, 86, 64, and 97\% of the test/train splits.}
	\label{fig:res_phmrc_CSMF}
\end{figure}

CSMF accuracy results are shown in Figure \ref{fig:res_phmrc_CSMF}. FARVA outperforms its competitors with regards to CSMF accuracy in each site. As before, FARVA shows an improvement in performance over the BF model in each setting. Unsurprisingly, the regions of highest individual performance gain for FARVA are also those with highest CSMF performance gain (e.g., Mexico City and Uttar Pradesh).  However, in some contexts the BF model is able to allocate probability across true causes in test data nearly as well as FARVA (e.g., Pemba Island and to a lesser extent Bohol). See the online supplementary materials for histograms showing the difference in top cause and CSMF accuracy between FARVA and BF in each train/test split

\subsection{Inference on VA questionnaire data}

Of interest in addition to prediction is model-based inference on the structure of responses within VA questionnaires. The FARVA model allows for both the mean and the covariance structure of the latent symptom vector $\bm{z}_i$ to be interrogated. In the following examples the same set of binary symptoms as those used to train the model in the PHMRC predictive ability assessment are considered. The PHMRC variable names and questions are included in full in the online supplementary materials.

\subsubsection{Symptom-level mean}

First consider an intercept-only FARVA model, i.e. one in which $x_i$ is set to $1$ in equations (\ref{eq:xi_regression}) and (\ref{eq:eta_regression}), run using all 7840 adult deaths in the PHMRC data set. Figure \ref{fig:inference_full_example} illustrates model inference on $\mathbb{E}[\bm{z}_{c[i]}],$ the mean of the latent symptom vector $\bm{z}_{c[i]}$, for some select symptoms and causes. The most important thing to notice about the example is that the latent symptom posterior means differ by cause. For example, samples of $\mathbb{E}[z_{c[i],\text{Cough}}]$ tend to be negative for CODs cirrhosis, and homicide, and positive for COD pneumonia, implying that coughing occurs more frequently for pneumonia deaths. These learned cause-specific patterns of latent symptom means, which differ across CODs, are a part of what drives the model's prediction ability.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=0.8\textwidth,keepaspectratio]{figs/inference/all/symp_mean_example.pdf}
	\caption{Violin plot visualizing samples of $\mathbb{E}[z_{c[i],j}]$ for CODs cirrhosis ($n=313$), homicide ($n=167$), and pneumonia ($n= 540$) for select symptoms. An expected latent symptom mean being positive (negative) corresponds to more (less) expected observed ``Yes'' responses for those symptoms.}
\label{fig:inference_full_example}
\end{figure}

Next consider a FARVA model run using all training data with the binary covariate of  decedent age ($\geq 65$) included, i.e. one in which $\bm{x_i}$ is set to $(\bm{1}_{\text{age}[i]\geq 65},\bm{1}_{\text{age}[i]< 65})'$ in equations (\ref{eq:xi_regression}) and (\ref{eq:eta_regression}). Figure \ref{fig:inference_age_example} provides an example of FARVA's potential to elucidate differences in individual response patterns within COD by covariate, showing $\mathbb{E}[z_{c[i],j}]$ for select symptoms and CODs. The included non-injury-related latent symptom means for homicide deaths are smeared across zero for older decedents, likely due to the small number of observed elder deaths due to homicide combined with the general trend across CODs for senior decedents to experience more constellations of symptoms. Pneumonia appears to present with small differences across age groups for the select symptoms shown, with fever (weight loss) being more (less) common among non-elder decedents. Cirrhosis shows the least age-specific differentiation in these example symptoms.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=0.95\textwidth,keepaspectratio]{figs/inference/age/symp_mean_example.pdf}
	\caption{Violin plot of posterior draws of $\mathbb{E}[z_{c[i],j}]$ for CODs cirrhosis ($n_{\geq 65}=49$, $n_{< 65}=260$), homicide ($n_{\geq 65}=5$, $n_{< 65}=161$), and pneumonia ($n_{\geq 65}=193$, $n_{< 65}=341$) for select symptoms.}
	\label{fig:inference_age_example}
\end{figure}

When age is included as a covariate, it acts as a moderator of the latent mean pattern for a given cause and symptom. The most notable feature of the pattern of symptoms for elders is that they appear less pronounced than symptom patterns of younger decedents (note, e.g., how the 95\% credible interval for cough and fever spans zero for elders even in cirrhosis and homicide CODs). This is consistent with the hypothesis that elders tend to have more diffuse symptoms from a myriad of ailments, not just a pointed set of symptoms leading to a specific diagnosis. The latent mean $\mathbb{E}[z_{c[i],j}]$ is significantly greater than (less than) 0 for 5\% (71\%) of the 4658 possible symptom-cause combinations (137 included symptoms $\times$ 34 causes)  for non-elders, but only 4\% (67\%) of the possible symptom-cause combinations for elders. This phenomenon, although slight, illustrates that FARVA is able to be more focused for younger individuals and spread probability more diffusely for elders, who tend to exhibit a wider array of symptoms across causes, contributing to the improvement in top cause accuracy and CSMF accuracy seen in the previous section.

\subsubsection{Symptom-level association}

\citet{li2018bayesian} consider inference on overall symptom-level covariance in verbal autopsy data for continuous symptoms, but the results are not cause-specific. \citet{kunihama2018bayesian} discuss a model-based version of Cramer's V, but this metric is only applicable to binary symptoms. To the authors' knowledge, no work has the ability to perform model-based inference on cause-specific covariance in VA data of mixed type.

Again, first consider an intercept-only FARVA model run using all training data. Figure \ref{fig:inference_full_example_cov} provides an example of model inference on $\text{Cov}(z_{c[i],j},z_{c[i],k})$ for pairs of symptoms $j$ and $k$, which can be thought of as the latent symptom covariance. These learned patterns of latent symptom covariances capture additional information beyond just symptom prevalence that the model can use to differentiate between CODs.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=0.85\textwidth,keepaspectratio]{figs/inference/all/symp_cov_example.pdf}
	\caption{Posterior mean and lower/upper 95\% credible interval of $\text{Cov}(z_{c[i],j},z_{c[i],k})$ for CODs cirrhosis ($n=313$), homicide ($n=167$), and pneumonia ($n= 540$) for select symptoms. An expected latent covariance mean being positive (negative) corresponds to observed ``Yes'' or ``No'' responses for those symptoms to (not) co-occur. For example, samples of $\text{Cov}(z_{c[i],\text{Cough}},z_{c[i],\text{Fever}})$ tend to be positive across all shown CODs, implying that the symptoms cough and fever tend to co-occur. On the other hand, violent injury tends not to be associated with the other symptoms shown, and is negatively associated with these symptoms in the case of homicide deaths.}
	\label{fig:inference_full_example_cov}
\end{figure}

When age is included as a covariate, it acts as a weak moderator of the latent covariance pattern for given pair of symptoms within a cause (figures included in  online supplementary materials). Of note is that the patterns of association within a COD are similar across age groups, but the strength of association tends to be weaker for younger decedents. Furthermore, the credible intervals tend to be wider for older decedents. This finding differs from the exploratory finding that, across all causes of death, age tends to weaken symptom associations. The difference is due to the conditional nature of these covariance matrices; rather than marginalizing over COD, these symptom associations are conditional on a COD. It is likely that the combination of slightly stronger average covariance along with a wider credible interval is due to the tendency for elderly individuals to be ``hit harder'' by an illness, i.e. to suffer from a large number of the common symptoms for that illness. At the same time, they likely also tend to have many other symptoms unrelated to their COD, and make up a smaller proportion of deaths for each COD, hence the increased uncertainty.

\section{Discussion}\label{disc}

The PHMRC data set, while useful for model testing and validation, has some drawbacks. Deaths clear enough to be considered ``gold standard'' by the PHMRC may not be representative of the population of community deaths, relationships between symptoms/causes learned in one context may not transfer to new contexts (location/time), and, relevant to this work, few potential covariates are recorded. More generally, information in VA data is only as reliable as the interviewee's knowledge. Finally, the question of how to merge different questionnaires (e.g., the adult/child PHMRC form) and select/prioritize symptoms for inquiry remains an open area of research.

For comparability to recent competing methods, the real data analysis of this paper considered only the binary symptoms initially chosen by \citet{byass2012strengthening}. However, the selection of which symptoms to include and whether to include them as binary or continuous variables is an area of open and important research. This issue is not unique to FARVA, and a rigorous consideration of the importance of various symptoms in assigning COD across models would be useful for both model performance and to lower the burden placed on interviewees (a reduced interview set that takes less time would be preferred). An area of potential improvement to the FARVA model itself would be incorporating a prior that could encourage row-level sparsity in the cause-specific loadings matrix $\Lambda_{c[i]}$. This addition could help the model better capture the likely realistic idea that many symptoms may be unimportant to assigning a cause, and could also help assess the relevance of symptoms for inclusion on future questionnaires.

Training on data from one (or multiple) setting(s) and testing in another is common practice for VA algorithms (e.g., using data from Dar es Salaam to train a model, then using test data from Pemba Island). Recent work by \citet{clark2018quantifying} has quantified the effect on predictive performance of this practice. They show that the impact of the information about the joint distribution between symptoms and causes is at least as high as that of the underlying algorithm's logic. In other words, existing algorithms cannot correct for discrepancy between the joint distribution between symptoms and causes in a new site relative to the body of training data. Work by \citet{datta2018local} has focused on accounting for local discrepancy when learning the CSMF at a site with limited training data. It is likely that some (although almost certainly not all) of the variation between sites is systematic and related to measurable features such as seasonal difference, endemicity of various pathogens, urbanicity, etc. If so, some of this decline in predictive performance could be corrected for explicitly in a model like FARVA via the inclusion of region-specific or individual-level covariates in a model. Exploring which covariates help model transferability is a potential future area of research.

The work of McCormick, Clark, Li, and the rest of their team has been pivotal in universalizing algorithm-assigned COD using VA data. In large part, this is due to his team's excellent documentation and code base. The FARVA team strives to maintain similarly high standards of documentation and usability of the code under hopes that it will be friendly to the broader VA community. Algorithm code and a user manual for the FARVA model are available at \url{https://github.com/kelrenmor/farva}.

\section{Acknowledgments}\label{ack}

The authors are grateful to John A. Crump, Manuela Carugati, Matthew P. Rubach, and Michael J. Maze for introducing us to the problem and inspiring us to think about methodological developments, and for helpful discussions with them and the rest of the Investigating Febrile Deaths in Tanzania (INDITe) team. This research was partly supported by the US National Institutes of Health (R01AI121378) for the INDITe project and by the Department of Energy Computational Science Graduate Fellowship (DE-FG02-97ER25308). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

\section{Introduction}

Figures 1, 2, and 3 are created in the \path{intro_vis.R} file. These exploratory plots are based on the publicly available PHMRC data, which is loaded into \textbf{R} in the manner described in the manuscript.

\section{Model}\label{model}

No code in the files is explicitly used in the manuscript in this section. The Gibbs sampler described in Section 2.3 may be viewed in code form in the FARVA package at \url{https://github.com/kelrenmor/farva}.

\section{Results and discussion}\label{results}

Each of the following subsections has its own specific run folder.

\subsection{Simulation performance}

The simulation lives in the \path{simulation/} folder. Functions for simulated dataset generation and for running the Gibbs sampler are included in the FARVA \textbf{R} package, which is called from the \path{_run_simulation.R} run file. If the user has the capability to run code in parallel, they can set \texttt{num\_cores} to some value greater than 1 in the preamble section. If the user doesn't have the capability to run Ox code, they should set \texttt{ox\_path=NULL}. It is recommended that users interested in replicating the results of the paper on a much smaller scale set \texttt{num\_datasets=10} and examine the results for this smaller sample, and then scale back up if desired.

Simulated data sets will go into the \path{simulation/sim_data/} folder (one data set will yield one \texttt{.Rdata} file), FARVA model run output will go into the \path{simulation/sim_farva_runs/} folder (one model run will yield one \texttt{.Rdata} file), and results will go into the \path{simulation/sim_results/} file (one simulation setting will yield one \texttt{.Rdata} file, so there will be 9 in total). Note that the results from the 1000 runs included in the manuscript are maintained in the \path{simulation/sim_results/} folder so users can explore the comparative performance more carefully without needing to re-run the entire simulation.

\subsection{PHMRC predictive performance}

The PHMRC resampled performance assessment lives in the \path{phmrc_performance/} folder. The main function for running the Gibbs sampler is included in the FARVA \textbf{R} package, which is called from the \path{_run_phmrc.R} run file. If the user doesn't have the capability to run Ox code, they should set \texttt{ox\_path=NULL}. It is recommended that users interested in replicating the results of the paper on a much smaller scale set \texttt{nsims=5} and examine the results for this smaller sample, and then scale back up if desired.

Indices for each test/train split will go into the \path{phmrc_performance/phmrc_inds/} folder (one test/train split will yield one \texttt{.Rdata} file per location), FARVA model run output will go into the \path{phmrc_performance/phmrc_farva_runs/} folder (one model run will yield one \texttt{.Rdata} file), and results will go into the \path{phmrc_performance/phmrc_results/} file (one location will yield one \texttt{.Rdata} file, and there will be one complete file with all locations so there will be 7 in total). Note that the results from the 100 runs included in the manuscript are maintained in the \path{phmrc_performance/phmrc_results/} folder so users can explore the comparative performance more carefully without needing to re-run the entire simulation.

\subsection{Inference on VA questionnaire data}

The PHMRC inference code lives in the \path{phmrc_inference/} folder. The main function for running the Gibbs sampler is included in the FARVA \textbf{R} package, which is called from the \path{run_inference.R} file. Note that running FARVA on the full data set is fairly time-consuming and may take 8-12 hours.

\section{Discussion}\label{disc}

No code is used in this section of the manuscript.

\section{Acknowledgments}\label{ack}

No code is used in this section of the manuscript.

redpajama_set_name: RedPajamaArXiv


text:
STAGE REVIEW: Les Miserables @ Theatre Royal, Newcastle (20.08.19)
Mark Grainger reviews the epic theatrical spectacle
By Mark Grainger on Wednesday, August 21st, 2019
Image: The Barricade by Helen Maybanks
Les Miserables is at Theatre Royal, Newcastle until Saturday 5th October
It's more than fair to say that Les Miserables is a cross-cultural phenomenon. A giant of the theatrical landscape with a name that provokes a sense of awe in anyone who's seen it.
For those few who have managed to avoid either the novel, this long-running musical, or Tom Hooper's lavish, Oscar-winning film adaptation, Les Mis is an epic piece of storytelling. Initially it's the tale of Jean Valjean, a man imprisoned for the best part of twenty years for stealing a loaf of bread, as he tries to live apart from the relentless pursuit of obsessive and self-righteous police officer, Javert. This is a work that spans decades though, taking in the cruelty and poverty that befall the innocent Fantine, Valjean's shot at redemption with the child, Cosette, and her eventual romance with Marius, whose desperate stand at the barricades with his revolutionary friends forms the thrilling bulk of act two.
The thing with Les Mis is there's not much point in going further into specifics. Most people know the plot and the key songs, and there's no way that the performances and technical aspects of such a monolithic show are going to be anything less than top class. So it is at the Theatre Royal, but let's be clear, that is the highest compliment. If I Dreamed A Dream doesn't make you well up, and Do You Hear The People Sing doesn't stir your soul, then you should question where you misplaced it.
In truth, there's no show quite like Les Mis in terms of scope, scale and staging. The streets of revolutionary France, its sewers and taverns, are all brought vividly to life thanks to complex machinations that seemingly appear out of the air and sparingly used, but effective, projections. The characters are what Les Mis does best, and will always hit home. Whether in the zealous, black and white world-view of Javert, the unrequited yearning of Eponine, the irascible youth of Gavroche, or in any of the supporting roles, maximum skill and commitment simply pours from every actor on stage to intoxicating and emotional effect.
But then, that's the magic of Les Mis. A true theatrical spectacle for sure, but so much more. It's an experience that you owe it to yourself to see at least once in your life.
News: Newcastle Theatres – Photo Prints Fundraiser
INTERVIEW: David Baddiel
LIVE REVIEW: The Wildhearts, Backyard Babies, CKY @ O2 Academy, Newcastle (28.01.20)
LIVE REVIEW: Holy Moly & The Crackers, GGAllan Partridge @ Sage Gateshead (02.10.19)
ALBUM REVIEWl: The Hold Steady – Thrashing Thru The Passion
ALBUM REVIEW: The National – I Am Easy To Find
FEATURE: Five Reasons To Watch balletLORENT's After Dark @ Theatre Royal
STAGE REVIEW: A Christmas Carol @ Northern Stage, Newcastle (06.12.18)
STAGE REVIEW: The Kite Runner @ Theatre Royal, Newcastle (30.04.18)
ALBUM REVIEW: Eels – The Deconstruction
NEWS: Hedda Gabler @ Theatre Royal

redpajama_set_name: RedPajamaCommonCrawl


text:
Will fascism triumph in the Brazilian elections?
Federico Finchelstein
Former President Donald Trump officially endorsed Jair Bolsonaro in his re-election bid. In a context of economic crisis and general disgust with the Brazilian president's extremism, the rival candidate, former president Lula da Silva, leads in all polls by a difference of 15 points with a support of 45%. In this context, to what extent can Trump's praise serve Bolsonaro?
Trump explained that the "tropical Trump," as he is affectionately called, has done a great job for the wonderful people of Brazil." "President Bolsonaro loves Brazil above all else", Trump wrote. "He is a wonderful man and has my full and complete endorsement!"
On a personal level, love is reciprocal. Bolsonaro went so far as to say he "loves Trump" but it is important to think about the centrality of the cult of the leader and the expectation of total obedience from followers. This type of relationship is not uncommon among fascists and would-be fascists. In 1923, Hitler said, "If a German Mussolini were delivered to Germany, the people would kneel down and worship him more than Mussolini ever did."
Bolsonaro also agrees in making vulgar and exacerbated machismo a central campaign theme. For example, just as Trump had done in the 2016 campaign, Bolsonaro recently praised his penis and sexual potency to present his suitability as a candidate.
It makes sense that a president under investigation by the judiciary, who failed to deal with the pandemic and who was rejected by a majority of voters and then lied and continues to lie about the outcome of the election, would support a colleague who presents exactly his characteristics.
Worse, Trump then supported the January 6, 2021, takeover of parliament on his behalf on the basis of the big lie that the election had been stolen from him. The purpose of such mobilization was to keep Trump permanently in power regardless of the outcome of the election. It is at this point and in the repeated threats that, if he loses, Bolsonaro will not recognize the election results where a fascist danger for Brazil is possible.
Nowadays, the world is undergoing a profound transformation: an attempted return from populism to fascism. During the last century, fascism has evolved, leaders have reformulated its appearance. Although explicit fascism disappeared from power after World War II, its anti-democratic ideas survived, often intertwined with various currents of populism.
Despite his lack of originality, Bolsonaro stands out among contemporary authoritarian leaders as being, like Donald Trump, closer to fascism than other populists.
There are four elements that differentiate fascism from populism. These are fundamental pillars that fascism has and populism does not. The first is violence and the militarization of politics. This occurred in Spain in the Civil War and also in Fascist Italy and Hitler Germany. All catastrophic examples of political problems being solved through violence.
The second is the lie, the totalitarian, extreme lie, which not only deforms reality and creates alternative realities in people's minds, but attempts to create that new reality. The third is xenophobia, racism, the total demonization of the imagined external enemy. Any problem that exists (real or invented) is explained and solved through hatred of what is different. And the fourth is dictatorship. There is no fascism without dictatorship, although there can be dictatorship without fascism.
Historically, populists have been elected, for example, Juan Perón, Getulio Vargas, Hugo Chávez, Cristina Kirchner or Silvio Berlusconi. All of them, despite their authoritarianism, recognized the legitimacy of elections even when they did badly in them.
Any observer of Brazilian politics can note that Bolsonaro, a deep admirer of dictatorships and dictators, lacks in practice the fourth element (dictatorship) to become a full-fledged fascist.
In his Independence Day speech, and this is not irrelevant, Bolsonaro classified the other side as 'evil'. The president's words were as follows: "The evil that lasted 14 years in our country, that almost broke our homeland and wants to return to the scene of the crime. They will not, the people are on the side of good". Would this apocalyptic argument allow him to justify in the name of "good" a coup d'état?
Like Donald Trump, Bolsonaro represents a new type of global autocratic ruler who is legally elected, but also adopts elements that populist figures like Perón considered too controversial or even toxic: totalitarian lies, racism and illegal means as coups to destroy democracy from within.
Trump might best be considered a "wannabe fascist." Trump is a populist who aspires to return to a form of fascism. His administration was not full-blown fascism because it did not descend into dictatorship. But it could have been if his attempts to retain power after the 2020 elections had succeeded.
This is the danger facing Brazilian democracy. Bolsonaro is a sort of mini-Trump. The "tropical Trump" not only shares Trump's failed fascist desires but also follows the Trumpist model to destroy democracy from within.
Translated from Spanish by Janaína Ruviaro da Silva
Professor of History at New School for Social Research (New York). Taught also at Brown University. PhD at Cornell Univ. He is the author of several books on fascism, populism, dictatorships and the Holocaust. His new book is "A Brief History of Fascist Lies" (2020).


text:
Tippi Hedren tells Michael Jackson's ...
Bangladeshi man arrested for beating ...
Tippi Hedren tells Michael Jackson's former pet tigers of his death
The Birds star says she informed the two tigers she adopted from Jackson's Neverland ranch of the singer's death, telepathically
Wednesday 1 July 2009 17.10 BST
Tippi Hedren has called upon mourning fans of Michael Jackson to honour his memory by coming to the aid of the exotic animals rescued from his Neverland ranch.
In 2005, Hedren adopted Jackson's two tigers, Thriller and Sabu, after the singer gave up the private zoo at his Santa Barbara estate. She cares for them at her animal sanctuary, the Los Angeles Shambala Preserve, home to 68 tigers, lions and other big cats who were formerly employed in the entertainment industry, or who were unwanted pets.
The star of Alfred Hitchcock's The Birds says she informed Thriller and Sabu about the singer's death at the end of last week. "I went up and sat with them for a while and let them know that Michael was gone," she explained. "You don't know what mental telepathy exists from the human to the animal. But I hope they understood."
Hedren has called upon grieving members of the public to donate money to her organisation and help her call for a law that bans the keeping of big cats as pets.
Other former inhabitants of Jackson's private zoo were taken on by various sanctuaries and individuals throughout America. His most famed pet, the chimpanzee Bubbles, lives at the Centre for Great Apes in Florida, whose owners are also urging fans for support.
http://www.guardian.co.uk/film/2009/jul/01/tippi-hedren-michael-jackson-tigers
https://bigcatrescue.org/
Bangladeshi man arrested for beating tiger to death
Tripartite Mou To Ensure Effective Tiger Conservation
Ocelots move in to Louisville Zoo
12:24 PM EDT on Thursday, June 7, 2007 Louisville Zoo Press Release LOUISVILLE, Ky. ...
BCR June 8, 2007
'Female presence' to lure runaway tiger back to Panna
'Female presence' to lure runaway tiger back to Panna Neha SinhaPosted: Monday , Dec ...
Sportsmanship in Hunting Act
Putting a Stop to Pay-Per-View and Pay-to-Kill Hunting Michael Markarian Posted: June 4, 2009 ...


This recipe is a great success, because it does not require any effort and time for baking. But there is only one nuance: the cake for cheesecake should be sand. Therefore, if you can not find one in the store, it's better to do it yourself, because Biscuit cake for cheesecake is not suitable.
3. Pour the remaining milk into a saucepan and boil.
4. Pour the egg mixture into the boiling milk with a thin trickle, constantly stirring the milk with a whisk. Brew the mixture for 30-45 seconds.
9. put the cake on the bottom of the mold, pour the cream and put the cheesecake in the freezer for 3 hours.
10. Get cheesecake from the form, decorate with slices, berries and mint.


text:
The Dead Zone; Monk
By Bruce Fretts
Updated August 09, 2002 at 04:00 AM EDT
Pinterest Mail Email iphone Send Text Message Print Comment
When the USA Network lost its rights to telecast the WWF two years ago, it was widely perceived as a body slam to the basic cabler. Now, with the WWF — excuse me, WWE — in a ratings free fall, it seems like this stroke of good luck forced the net to develop better programming. In fact, USA has bounced back off the Nielsen canvas with two of the summer's highest-rated new cable series, The Dead Zone and Monk.
Ironically, these shows are the antithesis of Vince McMahon's aesthetic: They're about average guys with superhuman mental strength. Traumatic auto-related events trigger their psychological powers: The Dead Zone's Johnny Smith (Anthony Michael Hall) awakens from a car-crash-induced six-year coma to discover he can see visions of people's past, present, and future. And in the aftermath of his wife's death in a car bombing, detective Adrian Monk (Tony Shalhoub) develops an extreme case of obsessive-compulsive disorder that enables him to analyze evidence in mind-boggling detail.
Both dramas are built around can't-take-your-eyes-off-'em lead performances. With his naturally battle-scarred complexion and spooked-out orbs, Hall looks like he's been in a long-term coma. Since his days as John Hughes' resident geek god, he's gained heft as an actor, literally and figuratively. His Johnny isn't as creepy as Christopher Walken's was in David Cronenberg's 1983 film version (inspired, as is the series, by Stephen King's novel), nor should he be; would you want to invite a Walken zombie into your home every week? Hall even tosses off a few one-liners with an old pro's comic timing. Temporarily rehired at the Maine high school where he used to teach science, Johnny disarms unruly students with this quip: "I may just be your worst nightmare — a substitute teacher with extra-sensory perception."
Shalhoub, familiar to USA viewers from years of endless Wings reruns, has finally found a TV vehicle to display the dazzling talents he's shown in big-screen roles from slithery alien Jack Jeebs in the MIB movies to silver-tongued lawyer Freddy Riedenschneider in The Man Who Wasn't There. He doesn't overplay Monk's OCD symptoms for cheap laughs; instead, he subtly conveys the quiet misery of a man who's trapped in vicious cycles of irrational behavior. Deftly shifting between scenes of pathos and physical comedy, he renders Monk as an endearingly original hero, a germaphobe who's fighting to clean up the gritty streets of San Francisco without getting his hands dirty.
Alas, the series' supporting casts aren't so uniformly spotless. A League of Their Own's Bitty Schram, who's like a less strident Debi Mazar, plays off Shalhoub amusingly as Sharona Fleming, Monk's tough-tawkin' nurse/gal Friday. But Ted Levine (The Silence of the Lambs) deserves better than the half-note role of the hotheaded police captain who's eternally skeptical of Monk's unorthodox methods.
The Dead Zone's most affecting subplot involves the sprightly Nicole deBoer (Star Trek: Deep Space Nine) as Johnny's ex-fiancee, Sarah, who gave up hope and married the town's foursquare sheriff (Chris Bruno) — but not before bearing Johnny's namesake son. Such a triangle may be a soapy cliche, but it's played with real conviction by Hall, deBoer, and All My Children alum Bruno. If only John L. Adams (as Johnny's physical therapist/man Friday) and M*A*S*H vet David Ogden Stiers (as a shady reverend who became the coma victim's legal guardian) had equally well-defined parts.
Like all shows, USA's duo will rise or fall based on the quality of their scripts, which have been uneven so far. Created by father-and-son Star Trek scribes Michael and Shawn Piller, The Dead Zone has used Hall's powers to spin a clever variation on 12 Angry Men — Johnny must convince 11 other jurors of a defendant's innocence based solely on his psychic insight. Yet another episode relied too heavily on corny 1940s flashbacks as Johnny helped a WWII vet find his long-lost love (the dialogue was straight out of a Jon Lovitz-Phil Hartman SNL old-movie parody: "Man oh man, a dolly just walked in to beat the band!"). And while Monk's pilot, penned by Rat Race's Andy Breckman, was fresh and funny, the second episode felt like a rerun, right down to the shadow-of-the-Golden Gate Bridge epilogue.
Still, compared to previous USA Network dreck like Weird Science (based on Hall's 1985 flick), Silk Stalkings, and Duckman, The Dead Zone and Monk pack more creative muscle than The Rock and Triple H combined. And if you don't like them, you can always watch the WWE on TNN. The Dead Zone: B Monk: B+


Running a business can be quite the challenge, especially a restaurant. A group's opinion on food can greatly vary, making it difficult to exactly discern their likes or distastes for the dish itself or the way it was prepared. Nevertheless, a good chef uses quality ingredients and products and ingredients that will make it rich in flavor and texture, all of which can be found at our store.
No matter what restaurant supplies you're looking for, we're sure to have it. From condiment holders, to plates and utensils, we'll make sure your restaurant is supplied with everything needed to strengthen the flavor of your food and ambiance of the setting.