§ INTRODUCTION

Verbal autopsy (VA) is a well established approach to ascertain
cause-of-death when medical certification and full autopsies are not
feasible or practical [Garenne, 2014, Taylor et al., 1983].
After a death is identified, a specially-trained fieldworker interviews
the caregivers (usually family members) of the decedent. A typical VA
interview includes a set of structured questions with categorical or
quantitative responses and a narrative section that records the `story'
of the death from the respondent's point of view
[World Health Organization and others, 2012]. Currently, there are multiple commonly used
questionnaires with overlapping, but not identical questions.

The process of inferring a cause from VA data consists of two
components. First, there must be some external information about the
relationship between causes and symptoms. This is similar to the role of
training data in supervised learning tasks. In the context of VA,
training dataset of cases with cause of death assigned by in person
autopsy is rare, and collecting such data is extremely time and
resource-intensive, and requires strong assumptions about the
generalizability of deaths in the training set to the population of
interest. A more common practice is to obtain training data using
clinically trained, experienced physicians read a fraction of the
interviews and determine causes. To address the fact that physicians
frequently do not agree on causes, VA interviews are often read by two
physicians, and sometimes three, and the final causes are determined
through a consensus mechanism <cit.>. Another
means of obtaining this information is directly through expert opinion.
A common practice is to ask groups of physicians to rate the likelihood
of a symptoms occurring given a particular cause of death, which can be
converted into a set of probabilities of observing a symptom given a
particular cause. Such expert knowledge can then be used for modeling.
The second component is the algorithmic or statistical method that
actually assigns the cause of death by combining the expert or learned
cause-symptom relationship with the symptoms observed in the VA
interviews from the target population.

In all, therefore, there are three key pieces required to analyze VA
data: (i) VA survey data itself, (ii) inputs that give information about
the association between symptoms and causes, and (iii) a statistical or
algorithmic method for assigning a likely cause. The current state of VA
literature usually does not distinguish between these three, in part
because existing software for algorithmic and statistical methods
require a specific set of inputs and survey format. This restriction
prevents robust comparison between methods and contexts. A health agency
in one region may, for example, want to analyze VA data using the same
VA algorithm as a neighboring region to ensure estimates are comparable.
Unless the agencies used the same survey format, however, this is not
possible with existing tools. The openVA package [Li et al., 2021]
addresses this issue through an open-source toolkit that (i) performs
data processing and conversion between existing survey formats, (ii)
implements multiple currently available algorithmic and statistical
methods, and (iii) provides visualization and tools for interpreting

§.§ The openVA package

The openVA comprises a suite of R [R Core Team, 2020]
packages for the analysis of verbal autopsy data. The goal of this
package is to provide researchers with an open-sourced tool that
supports flexible data input format and allows easier data manipulation
and model fitting. The openVA family consists of four core
required packages that are on CRAN, InterVA4, InterVA5,
InSilicoVA, and Tariff, and an optional package
nbc4va. Each of these packages implements one coding algorithm.
Three survey formats are supported currently: the WHO 2012 and WHO 2016
instrument after standard dichotomization procedures into binary
variables [World Health Organization and others, 2012, 26], the Institute for
Health Metrics and Evaluation (IHME) questionnaire in the format of the
Population Health Medical Research Consortium (PHMRC) dataset
[Murray et al., 2011], and customized binary data with training

The main focus of this paper is to provide a general introduction to the
implementation details of the included algorithms both in terms of the
underlying methodology, and through a series of case studies. For each
of the algorithms discussed, there is a standalone R package
available on CRAN, and three of them, InterVA-4, InterVA-5 and Tariff,
are also available in compiled software program distributed by the
original authors. The openVA package has four major contributions:

* It provides a standardized and user-friendly interface for analysts to fit and evaluate each method on different types of data input using standard syntax. Previously, most of the methods are designed to be used specifically with their own input formats and are usually incompatible with others. The openVA package closes this gap and allows easier and fair model comparison of multiple algorithms on the same data. This significantly facilitates further research on VA algorithms.

* It provides a series of functionalities to summarize and visualize results from multiple algorithms, which is helpful for analysts not familiar with data manipulation and visualization in R.

* It does not directly implement any algorithms for coding VA data[A special case is when we extend the InterVA-4 algorithm to the scenario where symptoms and causes are not in the pre-defined set provided by the original InterVA software, the extension is included in the openVA package instead of the InterVA4 and InterVA-5 packages for a faithful replication of the InterVA methods themselves.], so that it is possible for a research group to maintain their own algorithm implementations callable from the openVA package, while also make it available to the general users as a standalone piece of software. For example, the nbc4va was developed and maintained independently by researchers at the Center for Global Health Research in Toronto, but is designed so that it can be seamlessly integrated into the openVA package.

* It is fully open-sourced, and can be run on multiple platforms. The open-sourced nature of openVA significantly expands its potential for methodological research and its suitability for integration within a larger data analysis pipeline. Compared to the alternative implementations, the InterVA-4 and InterVA-5 software are distributed as compiled software that can only be run on Windows system. They provide the source codes as an additional code script, which are difficult to modify and re-compile.  Tariff, as implemented through the SmartVA-Analyze application [Serina et al., 2015], is also primarily distributed as a compiled application that can only be run on Windows system [Institute for Health Metrics and Evaluation, 2021]. Their source codes was recently made available under the open source MIT License on GitHub [Institute for Health Metrics and Evaluation, 2021].

The rest of this paper is organized as follows: In Section
<ref> we briefly introduce the main component packages and
the underlying algorithms. We demonstrate model fitting with the
openVA package for different input data formats in Section
<ref> to <ref>. We first show in Section
<ref> how to prepare the input data. We then demonstrate
fitting different VA models in Section <ref>, and
functionalities to summarize results in Section <ref>. We
then discuss how additional information can be incorporated into
modeling VA data in Section <ref>. Section <ref>
briefly surveys additional packages and software developments built
around the openVA package. We end in Section <ref>
with a discussion of remaining issues and limitations of the existing
automated VA algorithms and proposes new functionalities to be included
in openVA package.

§ STRUCTURE OF THE OPENVA SUITE

The openVA suite of packages currently consists of four standalone
packages that are available on the CRAN and one optional package hosted
on Github. In this section, we first provide a brief introduction to
these five packages, and we discuss the mathematical details behind each
algorithm in the next subsection.

* InterVA4 [Li et al., 2014, Li et al., 2019] is an R
package that implements the InterVA-4 model
[Byass et al., 2012]. It provides replication of both
InterVA software version 4.02 and the later release of version 4.03
update [Byass, 2015]. The standard input of InterVA4 is in
the form of a pre-defined set of indicators, based on the 2012 WHO VA
instrument (see
The default InterVA-4 algorithm cannot be applied to other data input
format because its internal built-in prior information is specific to
a fixed set of indicators and causes. The same restriction is also
maintained in InterVA4 package. However, the mathematical
formulation of InterVA-4 model is completely generalizable to other
binary input format. The generalized algorithm is described in Section
<ref>, and also implemented in the openVA

* InterVA5 [Thomas et al., 2021] is an R package that
implements the InterVA-5 model developed by Peter Byass
The InterVA-5 model updates the previous version in several ways.
First, the input data must adhere to the format of the 2016 WHO VA
instrument [26]. Second, changes have been made
to the data processing steps, which are described in Section
<ref>. It is also worth noting that the model outputs
have expanded by the inclusion of the most likely Circumstances Of
Mortality CATegory, or COMCAT, among the results – the categories
include: culture, emergency, health systems, inevitable, knowledge,
resources, or an indeterminant combination of multiple factors
<cit.>. Despite these
changes, the mathematical formulation of InterVA-5 is identical to
that of InterVA-4.

* InSilicoVA [Li et al., 2021] is an R package that
implements the InSilicoVA algorithm, a Bayesian hierarchical framework
for cause-of-death assignment and cause-specific mortality fraction
estimation proposed in McCormick et al., 2016. It is originally designed to
work with the WHO VA instrument, i.e., the same input data as in
InterVA software, but is also generalizable to other data input
format. It is a fully probabilistic algorithm and could incorporate
multiple sources of information, such as known sub-population in the
dataset, and physician coding when such data is available. The MCMC
sampler is implemented in Java for improved speed.

* Tariff [Li et al., 2018] is an R package that
implements the Tariff algorithm [James et al., 2011]. It most closely
reflects the description of Tariff 2.0 method
[Serina et al., 2015]. The Tariff algorithm is developed by the
Institute for Health Metrics and Evaluation (IHME) and officially
implemented in the SmartVA-Analyze software [Institute for Health Metrics and Evaluation, 2021].
However, as the developers of this R package are not
affiliated with the authors of the original algorithm, there are some
discrepancies in implementation. The source code of the two versions
of Tariff was not publicly available at the time when the Tariff
package was created, so the package was developed based solely on the
descriptions in the published work. Despite the difference in
implementation, Tariff is able to achieve comparable results as
the published work as demonstrated in McCormick et al., 2016. More detailed
descriptions of the Tariff implementations are also discussed in
the supplement of McCormick et al., 2016. The later released
Python source codes of SmartVA-Analyze has been
incorporated in the web application extension of the openVA
package, which we briefly discuss in Section <ref>.

* nbc4va [Miasnikof et al., 2015, Wen et al., 2018] is an R
package that implements the Naive Bayes Classifier for VA encoding. It
calculates the conditional probabilities of symptoms given causes of
death from training dataset instead of using physician provided
values. nbc4va is developed and maintained by researchers at the
Center for Global Health Research in Toronto, but is designed so that
it can be seamlessly integrated into openVA. Currently the
nbc4va package is hosted on Github and is an optional package
that users can choose to load separately.

The openVA package is hosted on CRAN and can be installed with the
following commands. For the analysis in this paper, we also install the
nbc4va package separately from Github. The versions of the
supporting packages can be checked in R using the
openVA_status() function.

R> library(openVA)
R> library(devtools)
R> install_github("rrwen/nbc4va")
R> library(nbc4va)

R> openVA_status()

§.§ Overview of VA cause-of-death assignment methods

The general modeling framework for model VA data consist of first
converting the surveys into a series of binary responses to questions
about each death, and then the two main goals of a typical VA analysis
are to estimate the population cause-specific mortality fractions
(CSMF), and the probability distribution or rankings of cause-of-death
(COD) for each individual death. In this section, we formally compare
the modeling approaches utilized by each algorithm. We adopt the
following notations. Consider \(N\) deaths, each with \(S\) binary
indicators of symptoms. Let \(s_{ij}\) denote the indicator for the
presence of \(j\)-th symptom in the \(i\)-th death, which can take
values from 0, 1, or NA (for missing data). We consider a pre-defined
set of causes is of size \(C\). For the \(i\)-th death, denote the COD
by \(y_i \in \{1, ..., C\}\) and the probability of dying from cause
\(k\) is denoted by \(P_{ik}\). For the population, the CSMF of cause
\(k\) is denoted as \(\pi_k\), with \(\sum_{k=1}^C \pi_k = 1\).

* InterVA4 [Byass et al., 2012] and InterVA5 [Byass, 2018] algorithm calculate the probability of each COD given the observed symptoms using the Bayes rule, so that

\[
P_{ik} = \frac{\pi_{k}^{(0)} \prod_{j=1}^S P(s_{ij}=1|y_{i}=k) \mathbf{1}_{s_{ij} = 1}}
{\sum_{k' = 1}^C \pi_{k'}^{(0)} \prod_{j=1}^S P(s_{ij}=1|y_{i}=k') \mathbf{1}_{s_{ij} = 1}}
\]

where both the prior distribution of each causes, $\pi_{k}^{(0)}$, and the conditional probability of each symptom given each cause, $P(s_{ij}=1|y_{i}=k)$, are fixed values provided in the algorithm. The conditional probabilities, $P(s_{ij}=1|y_{i}=k)$, used in InterVA algorithms are represented as rankings with letter grades instead of numerical value [Byass et al., 2012]. For example, $P(s_{ij}=1|y_{i}=k) = A+$ is translated into $P(s_{ij}=1|y_{i}=k) = 0.8$, etc. The standard InterVA software only supports the fixed set of symptoms and causes where such prior information is provided. For a different data input format, this formulation can be easily generalized if training data is available. We include in the openVA package an extension of the algorithm that calculates $\hat P(s_{ij}=1|y_{i}=k)$ from the empirical distribution in the training data and then maps to letter grades with different truncation rules. Details of the extended InterVA algorithm can be found in McCormick et al., 2016.

After the individual COD distributions are calculated, InterVA-4 utilizes a series of pre-defined rules to decide up to top three most likely COD assignments and truncates the probabilities for the rest of the CODs to 0 and adds an 'undetermined' category so that the probabilities sum up to 1 (See the user guide of Byass, 2015). Then the population-level CSMFs are calculated as the aggregation of individual COD distribution, such that

\[
\pi_k = \sum_{i=1}^N P^*_{ik}
\]

where $P^*_{ik}$ denotes the individual COD distribution after introducing the undetermined category.

* Naive Bayes Classifier [Miasnikof et al., 2015] is very similar to the InterVA algorithm with two major differences. First, instead of considering only symptoms that present, NBC algorithm also considers symptoms that are absent. Second, the conditional probabilities of symptoms given causes are calculated from training data instead of given by physicians, which is similar to our extension of InterVA discussed above. Similar to InterVA, the NBC method can be written as

\[
P_{ik} = \frac{\pi_{k}^{(0)} \prod_{j=1}^S (P(s_{ij}=1|y_{i}=k) \mathbf{1}_{s_{ij} = 1} + P(s_{ij} \neq 1|y_{i}=k) \mathbf{1}_{s_{ij} \neq 1})}
{\sum_{k' = 1}^C \pi_{k'}^{(0)} \prod_{j=1}^S (P(s_{ij}=1|y_{i}=k') \mathbf{1}_{s_{ij} = 1}+ P(s_{ij} \neq 1|y_{i}=k') \mathbf{1}_{s_{ij} \neq 1})}
\]

and the CSMFs are calculated by $\pi_k = \sum_{i=1}^N P_{ik}$.

* InSilicoVA algorithm [McCormick et al., 2016] assumes a generative model that characterizes both CSMF at the population level, and the COD distributions at the individual level. In short, the core generative process assumes

\begin{eqnarray} \nonumber
s_{ij} | y_i = k &\propto& \mbox{Bernoulli}(P(s_{ij} | y_i = k)) \\\nonumber
y_i | \pi_1, ..., \pi_C &\propto& \mbox{Multinomial}(\pi_1, ..., \pi_C) \\\nonumber
\pi_k &=& \exp \theta_k / \sum_{k=1}^C \exp \theta_k \\\nonumber
\theta_k &\propto& \mbox{Normal}(\mu, \sigma^2)
\end{eqnarray}

and hyperpriors are also placed on $P(s_{ij} | y_i = k)$, $\mu$, and $\sigma^2$. The priors for $P(s_{ij} | y_i = k)$ are set by the rankings used in InterVA-4 if the data is prepared into InterVA format, or learned from training data if otherwise. Parameter estimation is performed using Markov Chain Monte Carlo (MCMC), so that a sample of posterior distribution of $\pi_k$ can be obtained after the sampler converges. Additional heuristics of data processing steps introduced since the original publication of the method is discussed in Section <ref> of the Appendix.

* Tariff algorithm [James et al., 2011] differs from all other three methods in that it does not calculate an explicit probability distribution of COD for each death. Instead, for each death $i$, a Tariff score is calculated for each COD $k$ so that

\[
Score_{ik} = \sum_{j = 1}^{S} \mbox{Tariff}_{kj}\mathbf{1}_{s_{ij}=1}
\]

where the symptom-specific Tariff score $\mbox{Tariff}_{kj}$ is defined as

\[
\mbox{Tariff}_{kj} = \frac{n_{kj} - median(n_{1j}, n_{2j}, ..., n_{Cj})} {IQR(n_{1j}, n_{2j}, ..., n_{Cj})}
\]

where $n_{kj}$ is the count of how many deaths from cause $k$ contain symptom $j$ in the training data. The scores calculated are then turned into rankings by comparing to a reference distribution of scores calculated from re-sampling the training dataset to have a uniform COD distribution. It is worth noting that the Tariff algorithm produces the COD distribution for each death in terms of their rankings instead of the probability distributions. And thus the CSMF for each cause $k$ is calculated by the fraction of deaths with cause $k$ being the highest ranked cause, i.e.,

\[
\pi_k = \frac{\sum_{i=1}^N\mathbf{1}_{y_i = k}}{N}
\]

In addition to the different model specifications underlying each
algorithm, two major distinctions are most significant in understanding
the four methods for practitioners to interpret the model results.
First, the missing indicators are assumed to be equivalent to `absence'
in InterVA, NBC, and Tariff, but strictly as `missing' in InSilicoVA.
Second, the CSMFs are calculated in three different ways. Tariff
calculates CSMF as the proportion of the most likely COD assignments in
the dataset. InterVA-4 calculates CSMF as the aggregated distribution of
up to the top three most likely CODs. And InSilicoVA estimates CSMF
directly from the parametric model using MCMC. Some further feature
comparisons are summarized in Table <ref>.

Feature       InterVA     Tariff     NBC     InSilicoVA
Exact replication in the current openVA package     Yes     No    Yes     Yes

Implementable without training dataset      Yes     No     No     Yes

Can produce instantaneous results for single death      Yes     Yes     Yes     No

Accounts for absence of symptoms     No     No     Yes     Yes

Accounts for missing symptoms     No     No     No     Yes

Provides individual COD distribution     Yes     No     Yes     Yes

Direct estimation of CSMF and its uncertainty     No     No     No     Yes

Comparison the features of the four VA methods implemented in the openVA package.

§ DATA PREPARATION

In the openVA package, we consider two main forms of standardized
questionnaire: the WHO instrument and the IHME questionnaire. In this
section, we focus on describing these two data formats and tools to
clean and convert data into standardized forms for VA algorithms.
Pre-processing the raw data collected from the Open Data Toolkit is
usually performed with additional packages and software outside of the
analysis pipeline in R. We briefly mention software for data
preprocessing in Section <ref>.

§.§ Standardizing format: WHO

For users familiar with InterVA software and the associated data
processing steps, the standard input format from the WHO 2012 and 2016
instrument is usually well understood. For the 2012 instrument, the data
expected by the InterVA-4 software is organized into a data frame where
each row represents one death, and contains \(246\) fields, starting
from the first item being the ID of the death. The \(245\) items
following the ID each represent one binary variable of
symptom/indicator, where `presence' is coded by `Y', and `absence' is
coded by an empty cell. The details of this format, as well as the
translation from WHO 2012 instrument, could be found at

To accommodate updates for the WHO 2016 instrument
[26], InterVA-5 software accepts a data frame
with \(354\) columns that include \(353\) columns of symptom/indicators
followed by an additional column for the record ID. It should be noted
that the R package InterVA5 retains the format with the
record ID residing in the first column. Another important update with
InterVA-5 is that it acknowledge the difference between both “Yes” and
“No” (or “Y/y” and “N/n”, which is different from the coding
scheme in InterVA-4) are processed as relevant responses, while all
other responses are treated as missing values and ignored. With respect
to the list of causes of death, InterVA-5 utilizes the WHO 2016 COD
categories, which is nearly identical to the WHO 2012 COD categories
(used by InterVA-4) except that haemorrhagic fever and dengue fever are
two separate categories in the 2016 COD categories.

The same input format is inherited by the openVA package, except
for one modification. We further distinguish `missing' and `absent' in
the input data frame explicitly. This is conceptually easy to
understand: knowing that a symptom does not exist provides some
information to the possible cause assignment, while a missing symptom
does not. Missing data could arise from different stages of the data
collection process. Although in theory, most of the VA algorithms could
benefit from distinguishing `missing' from `absent', InSilicoVA is the
only algorithm that has implemented with missing data. We highly
recommend users to pre-process all the raw data so that a `missing'
value in the data spreadsheet is coded by a `.' (following the
stata practice), and an `absent' is coded by a empty cell, as
in the standard InterVA-4 software. For WHO 2016 data, both `.' and `-'
(following the default coding scheme of InterVA-5 software) are
interpreted as missing values. For methods other than InSilicoVA, the
`missing' and `absent' will be considered the same internally and thus
will not introduce compatibility problem.

§.§ Standardizing format: PHMRC

For studies that need to compare with existing work using PHMRC gold
standard dataset [Murray et al., 2011], a practical challenge is
that although the data is publicly accessible, the pre-processing steps
described from the relevant publications are not clear enough nor easy
to implement. The openVA package provides functions to clean up
the PHMRC gold standard data in the best way we could replicate from

First, we allow users to download part or all of the PHMRC gold standard
data directly from its on-line repository:

R> PHMRC_first1000 <- read.csv(getPHMRC_url("adult"), nrows = 1000)

The data can then be cleaned up into a set of dichotomized symptoms, as
described in Section 2 of the supplement material of McCormick et al., 2016.

This dichotomization process was originally described in
Murray et al., 2011. However, following the procedures
described in the original paper, we are not able to reproduce the
thresholds provided with the on-line supplement materials. Therefore we
have included both the thresholds used in Murray et al., 2011
and the thresholds calculated as in their description based on the data
that user inputs. See the following codes for the summary of the two
different transformation results.

R> convert.default <- ConvertData.phmrc(PHMRC_first1000, phmrc.type = "adult",
+                                      cutoff = "default", cause = "va34")

The first column is site, assign IDs to each row by default
1000 deaths in input. Format: adult
168 binary symptoms generated

Number of Yes        21023
Number of No         124644
Number of Not known  22333

R> convert.adapt <- ConvertData.phmrc(PHMRC_first1000, phmrc.type = "adult",
+                                    cutoff = "adapt", cause = "va34")

The first column is site, assign IDs to each row by default
1000 deaths in input. Format: adult
168 binary symptoms generated

Number of Yes        21711
Number of No         123956
Number of Not known  22333

Notice that the original PHMRC data are useful for comparing and
validating new methods, as well as for using as training data, but the
cleaning functions require only the columns to be exactly as the PHMRC
gold standard dataset on-line, so they could also be used for new data
that are pre-processed into the same
format[However, since calculating the thresholds adaptively requires the knowledge of underlying cause-of-deaths, the latter approach could not be applied to unlabeled datasets.].

§.§ Other formats

In addition to the standard formats discussed previously, researchers
might also be interested in using customized dichotomous symptoms data
in their analysis. The openVA package also supports customized
input as long as they are dichotomous. In such case, neither the
built-in conditional probability matrix of InterVA nor the PHMRC gold
standard dataset could be used to learn the relationship between
training and testing data, thus some training data with known causes of
death is necessary for all three algorithm. The ConvertData()
function can be used to convert data in a different coding scheme into
the format used by the openVA package. This is demonstrated in the
following small example.

R> toyid <- c("d1", "d2")
R> toycause <- c("A", "B")
R> toyData <- matrix(c("Yes", "No", "Don't know",
+                 "Yes", "Refused to answer", "No"),
+                 byrow = TRUE, nrow = 2, ncol = 3)
R> toyData <- cbind(toyid, toycause, toyData)
R> colnames(toyData) <- c("ID", "Cause", "S1", "S2", "S3")

R> toyData

ID   Cause S1    S2                  S3
[1,] "d1" "A"   "Yes" "No"                "Don't know"
[2,] "d2" "B"   "Yes" "Refused to answer" "No"

R> toyDataNew <- ConvertData(toyData, yesLabel = "Yes", noLabel = "No",
+                   missLabel = c("Don't know", "Refused to answer"))
R> toyDataNew

ID Cause S1 S2 S3
1 d1     A  Y     .
2 d2     B  Y  .

§.§ Convert data between standardized formats

The openVA package supports customized set of symptoms to be used
as input for all methods, thus it is usually less important to convert
data from one format to another, since the conversion inevitably creates
loss of information. We recognize the exact mapping of symptoms between
different format can be useful for some applications. For example, a
full mapping of PHMRC dataset into the InterVA format enables the use of
physician provided conditional probabilities included in the InterVA
software. This remains as an important feature to be added to the
package in the future.

§ FITTING VA CAUSE-OF-DEATH ASSIGNMENT MODELS

In this section, we demonstrate model fitting process in the
openVA package using two types of datasets: (1) a random sample of
\(1,000\) deaths from ALPHA network without gold standard labels
collected with the WHO 2012 instrument, and (2) the adult VA records in
the PHMRC gold standard data, with the \(1,554\) records from Andhra
Pradesh, India used as testing set and the rest as training dataset. In
the first case without gold standard training data, only InterVA and
InSilicoVA can be fitted. All four methods can be fitted on the second

§.§ Modeling data collected with WHO 2012 questionnaire

The randomly sampled VA records from ALPHA network sites are already
included in the openVA package and can be loaded directly.

R> data(RandomVA1)
R> dim(RandomVA1)

[1] 1000  246

R> head(RandomVA1[, 1:10])

ID elder midage adult child under5 infant neonate male female
1 d1     Y                                             Y
2 d2     Y                                                    Y
3 d3            Y                                      Y
4 d4                  Y                                       Y
5 d5                  Y                                Y
6 d6                  Y                                       Y

The codeVA() function provides a standardized syntax to fit
different VA models. Internally, the codeVA() function organizes
the input data into standardized formats, checks for incompatibility of
the data and specified model, and calls the specified model fitting
functions. It returns a classed object of the specified model class. In
this example, we use the version 4.03 of InterVA software, which is the
latest release of the original software compatible with the WHO 2012
instrument. Any additional model-specific parameters can be passed
through the arguments of codeVA(). Here we specify the HIV and
malaria prevalence levels required by the InterVA model to be `high'.

R> fit_inter_who <- codeVA(data = RandomVA1, data.type = "WHO2012",
+                         model = "InterVA", version = "4.03",
+                         HIV = "h", Malaria = "h")

R> summary(fit_inter_who)

InterVA-4 fitted on 1000 deaths
CSMF calculated using reported causes by InterVA-4 only
The remaining probabilities are assigned to 'Undetermined'

Top 5 CSMFs:
cause                     likelihood
Undetermined              0.156
HIV/AIDS related death    0.122
Stroke                    0.072
Reproductive neoplasms MF 0.058
Pulmonary tuberculosis    0.055

We can implement InSilicoVA method with the similar syntax. We use the
default parameters and run the MCMC for \(10,000\) iterations. Setting
the auto.length argument to FALSE specifies that the algorithm
does not automatically increase the length of the chain when convergence
failed. In addition to this simple interface, the algorithm
implementation in the InSilicoVA packages allows many more types
of model customization and extensions that cater to specific practical
modeling needs. We will briefly discuss some unique features in Section

R> fit_ins_who <- codeVA(RandomVA1, data.type = "WHO2012", model = "InSilicoVA",
+                     Nsim = 10000, auto.length = FALSE)

R> summary(fit_ins_who)

InSilicoVA Call:
1000 death processed
10000 iterations performed, with first 5000 iterations discarded
250 iterations saved after thinning
Fitted with re-estimated conditional probability level table
Data consistency check performed as in InterVA4

Top 10 CSMFs:
Mean Std.Error Lower Median Upper
Other and unspecified infect dis  0.266    0.0168 0.235  0.265 0.301
HIV/AIDS related death            0.102    0.0091 0.085  0.102 0.119
Renal failure                     0.101    0.0108 0.084  0.101 0.123
Other and unspecified neoplasms   0.062    0.0089 0.046  0.061 0.080
Other and unspecified cardiac dis 0.058    0.0076 0.044  0.058 0.075
Digestive neoplasms               0.050    0.0077 0.033  0.050 0.065
Acute resp infect incl pneumonia  0.048    0.0073 0.034  0.049 0.063
Pulmonary tuberculosis            0.039    0.0068 0.025  0.039 0.054
Stroke                            0.038    0.0061 0.027  0.038 0.052
Other and unspecified NCD         0.034    0.0089 0.018  0.034 0.052

§.§ Modeling the PHMRC data

For the PHMRC gold standard dataset, we first load the complete dataset
from its on-line repository, and then organize them into training and
testing splits.

R> PHMRC_all <- read.csv(getPHMRC_url("adult"))
R> is.test <- which(PHMRC_all$site == "AP")
R> test <- PHMRC_all[is.test, ]
R> train <- PHMRC_all[-is.test, ]
R> dim(test)
\end{CodeInput}
\begin{CodeOutput}
[1] 1554  946
\end{CodeOutput}
\begin{CodeInput}
R> dim(train)
\end{CodeInput}
\begin{CodeOutput}
[1] 6287  946
\end{CodeOutput}
\end{CodeChunk}

In order to fit the models on the PHMRC data, we specify
\code{data.type = "PHMRC"}, and the column of the causes of death label
in the training data. The rest of the syntax is similar to the previous
example. When the input consist of both training and testing data, the
InterVA and InSilicoVA algorithm requires a small modification that map
the conditional probabilities of symptoms given causes in the training
dataset to a letter grade system assumed by the model. In such case, the
\code{version} argument for the InterVA algorithm is no longer needed.
There are several ways to do this conversion, specified by the
\code{convert.type} argument. The \code{convert.type = "quantile"}
performs the mapping so that the percentile of each rank stays the same
as the original \(P_{s|c}\) matrix in InterVA software. Alternatively we
can also use the original fixed values of translation, and assign letter
grades closest to each entry in \(\hat{P}_{s|c}\). This conversion is
specified by \code{convert.type = "fixed"}, and is more closely aligned
to the original InterVA and InSilicoVA setting. Finally, we can also
directly use the values in the \(\hat{P}_{s|c}\) without converting them
to ranks and re-estimating the values associated with each rank. This
can be specified by \code{convert.type = "empirical"}. In this
demonstration, we assume the fixed value conversion.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_inter <- codeVA(data = test, data.type = "PHMRC", model = "InterVA",
+                      data.train = train, causes.train = "gs_text34",
+                      phmrc.type = "adult", convert.type = "fixed")
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> fit_ins <- codeVA(data = test, data.type = "PHMRC", model = "InSilicoVA",
+                     data.train = train, causes.train = "gs_text34",
+                     phmrc.type = "adult", convert.type = "fixed",
+                     Nsim=10000, auto.length = FALSE)
\end{CodeInput}
\end{CodeChunk}

The NBC and Tariff method, on the other hand, does not need to perform
such conversion, and only need to change the \code{model} argument.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_nbc <- codeVA(data = test, data.type = "PHMRC", model = "NBC",
+                    data.train = train, causes.train = "gs_text34",
+                    phmrc.type = "adult")
\end{CodeInput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> fit_tariff <- codeVA(data = test, data.type = "PHMRC", model = "Tariff",
+                      data.train = train, causes.train = "gs_text34",
+                      phmrc.type = "adult")
\end{CodeInput}
\end{CodeChunk}

Notice that we do not need to transform the PHMRC data manually as
described in the previous section. Data transformations are performed
automatically within the \code{codeVA()} function. The arguments in
\code{ConvertData.phmrc()} can also be passed into \code{codeVA()}.

\section{Summarizing results}
\label{sec:results}

In this section we demonstrate how to summarize results, extract output,
and visualize and compare fitted results. All the fitted object returned
by \code{codeVA()} are \code{S3} objects, where readable summary of
model results can be obtained with the \code{summary()} function as
shown in the previous section. In addition, several other metrics are
commonly used to evaluate and compare VA algorithms at either the
population or individual levels. In the rest of this section, we show
how to easily calculate and visualize some of these metrics with the
\pkg{openVA} package.

\subsection{CSMF Accuracy}

We can extract the CSMFs directly using the \code{getCSMF()} function.
The function returns a vector of the point estimates of the CSMFs, or a
matrix of posterior summaries of the CSMF for the InSilicoVA algorithm.

\begin{CodeChunk}
\begin{CodeInput}
R> csmf_inter <- getCSMF(fit_inter)
R> csmf_ins <- getCSMF(fit_ins)
R> csmf_nbc <- getCSMF(fit_nbc)
R> csmf_tariff <- getCSMF(fit_tariff)
\end{CodeInput}
\end{CodeChunk}

One commonly used metric to evaluate the CSMF estimates is the so-called
CSMF accuracy, defined as \[
CSMF_{acc} = 1 - \frac{\sum_i^C CSMF_i - CSMF_i^{(true)}}{2(1 - \min CSMF^{(true)})}
\] We use the empirical distribution in the test data to calculate the
true CSMF distribution and evaluate the CSMF accuracy using the
\code{getCSMF_accuracy()} function. As discussed previously, the default
CSMF calculation is slightly different for different methods. For
example, the InterVA algorithm creates the additional category of
\texttt{Undetermined} by default, which is not in the true CSMF
categories and needs to be specified. The creation of the undetermined
category can also be suppressed by \code{interVA.rule = FALSE} in the
\code{getCSMF()} function call. For the InSilicoVA algorithm, we will
use the posterior mean to calculate the point estimates of the CSMF

\begin{CodeChunk}
\begin{CodeInput}
R> csmf_true <- table(c(test$gs_text34, unique(PHMRC_all$gs_text34))) - 1
R> csmf_true <- csmf_true / sum(csmf_true)
R> c(getCSMF_accuracy(csmf_inter, csmf_true, undet = "Undetermined"),
+   getCSMF_accuracy(csmf_ins[, "Mean"], csmf_true),
+   getCSMF_accuracy(csmf_nbc, csmf_true),
+   getCSMF_accuracy(csmf_tariff, csmf_true))
\end{CodeInput}
\begin{CodeOutput}
[1] 0.58 0.74 0.77 0.68
\end{CodeOutput}
\end{CodeChunk}

\subsection{Individual COD summary}

At each individual level, we can extract the most likely cause-of-death
assignment from the fitted object using the \code{getTopCOD()} function.

\begin{CodeChunk}
\begin{CodeInput}
R> cod_inter <- getTopCOD(fit_inter)
R> cod_ins <- getTopCOD(fit_ins)
R> cod_nbc <- getTopCOD(fit_nbc)
R> cod_tariff <- getTopCOD(fit_tariff)
\end{CodeInput}
\end{CodeChunk}

With the most likely COD assignment, other types of metrics based on
individual COD assignment accuracy can be similarly constructed by
users. The summary methods can also be called for each death ID. For
example, using Tariff method, we can extract the fitted rankings of
causes for the death with ID 6288 by

\begin{CodeChunk}
\begin{CodeInput}
R> summary(fit_inter, id = "6288")
\end{CodeInput}
\begin{CodeOutput}
InterVA-4 fitted top 5 causes for death ID: 6288

Cause                     Likelihood
Stroke                    0.509
Pneumonia                 0.318
COPD                      0.081
Other Infectious Diseases 0.064
Renal Failure             0.013
\end{CodeOutput}
\end{CodeChunk}

Notice that the direct call of summary for InSilcoVA does not provide
uncertainty estimates for individual COD assignments. This is because
the calculation of individual posterior probabilities of COD
distribution is relatively time-consuming and memory-intense. To obtain
individual-level uncertainty measurement, we can either run the MCMC
chain with the additional argument \code{indiv.CI = 0.95} when calling
\code{codeVA()}, or update the fitted object directly with the saved
posterior draws.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_ins <- updateIndiv(fit_ins, CI = 0.95)
R> summary(fit_ins, id = "6288")
\end{CodeInput}
\begin{CodeOutput}
InSilicoVA fitted top  causes for death ID: 6288
Credible intervals shown: 95%
Mean  Lower Median  Upper
Stroke                      0.5043 0.3485 0.5083 0.6361
Pneumonia                   0.4116 0.2615 0.4083 0.5834
Other Infectious Diseases   0.0660 0.0411 0.0642 0.0966
Epilepsy                    0.0099 0.0064 0.0097 0.0142
COPD                        0.0053 0.0031 0.0052 0.0079
Malaria                     0.0007 0.0005 0.0007 0.0011
Diabetes                    0.0005 0.0003 0.0005 0.0009
Acute Myocardial Infarction 0.0004 0.0003 0.0004 0.0006
Falls                       0.0004 0.0001 0.0004 0.0013
Renal Failure               0.0004 0.0002 0.0003 0.0005
\end{CodeOutput}
\end{CodeChunk}

For \(N\) deaths, \(C\) causes, the posterior mean of individual COD
distributions returned by the InSilicoVA model, along with median and
with credible intervals can be represented by a
\((N \times C \times 4)\)-dimensional array. The function
\code{getIndivProb()} extracts this summary in the form of a list of
\(4\) matrices of dimension \(N\) by \(C\), which can then be saved to
other formats to facilitate further analysis. For other methods, the
point estimates of individual COD distribution are returned as the \(N\)
by \(C\) matrix.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_prob <- getIndivProb(fit_inter)
R> dim(fit_prob)
\end{CodeInput}
\begin{CodeOutput}
[1] 1554   34
\end{CodeOutput}
\end{CodeChunk}

\subsection{Visualization}

The previous sections discuss how results could be extracted and
examined in \proglang{R}. In this subsection, we show some visualization
tools provided in the \pkg{openVA} package for presenting these results.
The fitted CSMFs for the top causes can be easily visualized by the
\code{plotVA()} function. The default graph components are specific to
each algorithm and individual package implementations, with options for
further customization.

\begin{CodeChunk}
\begin{CodeInput}
R> plotVA(fit_inter, title = "InterVA")
R> plotVA(fit_ins, title = "InSilicoVA", bw = TRUE)
R> plotVA(fit_nbc, title = "NBC")
R> plotVA(fit_tariff, title = "Tariff")
\end{CodeInput}
\begin{figure}[ht]

{\centering \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/vis-1-1} \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/vis-1-2} \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/vis-1-3} \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/vis-1-4}

\caption[Plot of top 10 CSMFs estimated by InterVA, InSilicoVA, NBC, and Tariff]{Plot of top 10 CSMFs estimated by InterVA, InSilicoVA, NBC, and Tariff.}\label{fig:vis-1}
\end{figure}
\end{CodeChunk}

The CSMFs can also be aggregated for easier visualization of groups of
causes. For InterVA-4 cause list, we included an example grouping built
into the package, so the aggregated CSMFs can be compared directly. In
practice, the grouping of causes of deaths often needs to be determined
according to context and the research question of interest. Changing the
grouping can be easily achieved by modifying the \code{grouping}
argument in \code{stackplotVA()} function. For example, to facilitate
the new category of \texttt{Undetermined} returned by InterVA, we first
modify the grouping matrix to include it as a new cause:

\begin{CodeChunk}
\begin{CodeInput}
R> data(SampleCategory)
R> grouping <- SampleCategory
R> grouping[,1] <- as.character(grouping[,1])
R> grouping <- rbind(grouping, c("Undetermined", "Undetermined"))
R> compare <- list(InterVA4 = fit_inter_who,
+                 InSilicoVA = fit_ins_who)
R> stackplotVA(compare, xlab = "", angle = 0,  grouping = grouping)
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=.7\linewidth,]{openVA-JSS_files/figure-latex/vis-2-1}

\caption[Comparing aggregated CSMF for InterVA-4 and InSilicoVA, adding undetermined category]{Comparing aggregated CSMF for InterVA-4 and InSilicoVA, adding undetermined category.}\label{fig:vis-2}
\end{figure}
\end{CodeChunk}

The ordering of the stacked bars can also be changed to reflect the
structures within the aggregated causes.

\begin{CodeChunk}
\begin{CodeInput}
R> order.group <- c("TB/AIDS",  "Communicable", "NCD", "External", "Maternal",
+             "causes specific to infancy", "Undetermined")
R> stackplotVA(compare, xlab = "", angle = 0, grouping = grouping,
+             order.group = order.group)
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=.7\linewidth,]{openVA-JSS_files/figure-latex/vis-3-1}

\caption[Comparing aggregated CSMF for InterVA-4 and InSilicoVA, with the causes reordered]{Comparing aggregated CSMF for InterVA-4 and InSilicoVA, with the causes reordered.}\label{fig:vis-3}
\end{figure}
\end{CodeChunk}

\section{Incorporating additional information}
\label{sec:insilico}

Among the VA methods discussed in this paper, the InSilicoVA algorithm
[McCormick \emph{et~al.}, 2016] allows more flexible modifications to the Bayesian
hierarchical model structure when additional information is available.
In this section, we illustrate two features unique to the InSilicoVA
method: jointly estimating CSMFs from multiple population, and
incorporating partial and potentially noisy physician codings in the

\subsection{Sub-population specific CSMFs}

In practice researchers may want to estimate and compare CSMFs for
different regions, time periods, or demographic groups in the
population. Running separate models on subsets of data can be
inefficient and does not allow each parameter estimation to borrow
information across different groups. The generative framework adopted by
InSilicoVA allows the specification of sub-population in analyzing VA
data. Consider an input dataset with \(G\) different sub-population, we
can estimate different CSMFs \(\pi^{(g)}\) for \(g = 1, ..., G\) for
each sub-population, while assuming the same conditional probability
matrix, \(P_{s|c}\) and other hyperpriors. As an example, we show how to
estimate different CSMFs for sub-populations specified by sex and age
groups, using a randomly sampled ALPHA dataset with additional columns
specifying the sub-population each death belongs to.

\begin{CodeChunk}
\begin{CodeInput}
R> data(RandomVA2)
R> head(RandomVA2[, 244:248])
\end{CodeInput}
\begin{CodeOutput}
stradm smobph scosts   sex age
1      .      .      .   Men 60+
2      .      .      . Women 60-
3      .      .      . Women 60-
4      .      .      . Women 60+
5      .      .      . Women 60-
6      .      .      . Women 60-
\end{CodeOutput}
\end{CodeChunk}

Then we can fit the model with one or multiple additional columns
specifying sub-population membership for each observation.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_sub <- codeVA(RandomVA2, model = "InSilicoVA",
+               subpop = list("sex", "age"),  indiv.CI = 0.95,
+               Nsim = 10000, auto.length = FALSE)
\end{CodeInput}
\end{CodeChunk}

Functions discussed in the previous sections works in the same way for
the fitted object with multiple sub-populations. Additional
visualization tools are also available. By specify
\code{type = "compare"}, we can plot the CSMFs for two sub-populations
on the same plot.

\begin{CodeChunk}
\begin{CodeInput}
R> plotVA(fit_sub, type = "compare", title = "Comparing CSMFs", top = 3)
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=0.9\linewidth,]{openVA-JSS_files/figure-latex/ins3-1}

\caption[Comparing CSMF for different sub-population]{Comparing CSMF for different sub-population.}\label{fig:ins3}
\end{figure}
\end{CodeChunk}

By default, the comparison plots will select all the CODs that are
included in the top \(K\) for each of the sub-populations. We can also
plot only subsets of them by specifying with causes are of interest.

\begin{CodeChunk}
\begin{CodeInput}
R> plotVA(fit_sub, type = "compare", title = "Comparing CSMFs",
+                 causelist = c("HIV/AIDS related death",
+                               "Pulmonary tuberculosis",
+                               "Other and unspecified infect dis",
+                               "Other and unspecified NCD"))
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=0.8\linewidth,]{openVA-JSS_files/figure-latex/ins4-1}

\caption[ Comparing the fraction of deaths due to selected CODs for different sub-populations]{ Comparing the fraction of deaths due to selected CODs for different sub-populations.}\label{fig:ins4}
\end{figure}
\end{CodeChunk}

If a single sub-population is selected by the \code{which.sub} argument,
the top CSMFs of the chosen sub-population will be plotted in the same
way as before.

\begin{CodeChunk}
\begin{CodeInput}
R> plotVA(fit_sub, which.sub = "Women 60-", title = "Women 60-")
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=0.8\linewidth,]{openVA-JSS_files/figure-latex/ins5-1}

\caption[Top 10 CSMF for a specified sub-population]{Top 10 CSMF for a specified sub-population.}\label{fig:ins5}
\end{figure}
\end{CodeChunk}

The \code{stackplot()} function can also be used to compare different
sub-populations in aggregated cause groups.

\begin{CodeChunk}
\begin{CodeInput}
R> stackplot(fit_sub)
R> stackplot(fit_sub, type = "dodge")
\end{CodeInput}
\begin{figure}[H]

{\centering \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/ins6-1} \includegraphics[width=0.49\linewidth,]{openVA-JSS_files/figure-latex/ins6-2}

\caption[Comparing aggregated CSMF for four different sub-population]{Comparing aggregated CSMF for four different sub-population.}\label{fig:ins6}
\end{figure}
\end{CodeChunk}

\subsection{Physician coding}

When physician coded cause of death is available for all or a subset of
the death, we can incorporate such information in the InSilicoVA model.
The physician coded causes can be either the same as the CODs used for
the final algorithm, or a higher level aggregation of them. When there
are more than one physician codes for each death and the physician
identity is known, we can first de-bias the multiple codes provided from
different physicians using the process described in McCormick \emph{et~al.}, 2016.
For the purpose of implementation, we only need to specify which columns
are physician IDs, and which are their coded causes respectively.

\begin{CodeChunk}
\begin{CodeInput}
R> data(SampleCategory)
R> data(RandomPhysician)
R> head(RandomPhysician[, 245:250])
\end{CodeInput}
\begin{CodeOutput}
smobph scosts        code1 rev1        code2 rev2
1      .      .          NCD doc9          NCD doc6
2      .      .          NCD doc4          NCD doc3
3      .      .          NCD doc1          NCD doc5
4      .      .      TB/AIDS doc4      TB/AIDS doc7
5      .      .      TB/AIDS doc5      TB/AIDS doc9
6      .      . Communicable doc9 Communicable <NA>
\end{CodeOutput}
\end{CodeChunk}

\begin{CodeChunk}
\begin{CodeInput}
R> doctors <- paste0("doc", c(1:15))
R> causelist <- c("Communicable", "TB/AIDS", "Maternal",
+               "NCD", "External", "Unknown")
R> phydebias <- physician_debias(RandomPhysician,
+           phy.id = c("rev1", "rev2"), phy.code = c("code1", "code2"),
+           phylist = doctors, causelist = causelist, tol = 0.0001, max.itr = 100)
\end{CodeInput}
\end{CodeChunk}

The de-biased step essentially creates a prior probability distribution
for each death over the broader categories of causes. Then to run
InSilicoVA with the de-biased physician coding, we can simply pass the
fitted object from the previous step to the model. Additional arguments
are needed to specify the external cause category, since they are
handled by separated heuristics, and the unknown category, which is
equivalent to an uniform probability distribution over all other
categories, i.e., the same as the case where no physician coding exist.

\begin{CodeChunk}
\begin{CodeInput}
R> fit_ins_phy <- codeVA(RandomVA1, model = "InSilicoVA",
+              phy.debias = phydebias, phy.cat = SampleCategory,
+              phy.external = "External", phy.unknown = "Unknown",
+              Nsim = 10000, auto.length = FALSE)
\end{CodeInput}
\end{CodeChunk}

This can be compared with the previous results without including
physicians codes.

\begin{CodeChunk}
\begin{CodeInput}
R> plotVA(fit_ins_who, title = "Without physician coding")
R> plotVA(fit_ins_phy, title = "With physician coding")
\end{CodeInput}
\begin{figure}[!h]

{\centering \includegraphics[width=0.7\linewidth,]{openVA-JSS_files/figure-latex/phy3-1} \includegraphics[width=0.7\linewidth,]{openVA-JSS_files/figure-latex/phy3-2}

\caption[Comparing fitted CSMF with and without physicians]{Comparing fitted CSMF with and without physicians}\label{fig:phy3}
\end{figure}
\end{CodeChunk}

\section{Other related software packages}
\label{sec:other}

Since the release of the \pkg{openVA} package on CRAN, there have been
many new developments in both methodology and software implementations
that build on the \pkg{openVA} suite of packages and further extended
its functionalities. Here we briefly survey some of the related methods,
software packages, and ongoing work in the VA community that are related
to \pkg{openVA}.

First, the ability to easily fit and compare existing methods using the
\pkg{openVA} package facilitated the development of several new VA
methods in the last several years. Most of the development focuses on
building statistical models to relax the conditional independence
assumption of symptoms given a cause of death
\citep[e.g.,][]{li2017mix, tsuyoshi2017, moran2021bayesian}. These
methods tend to be computationally more demanding compared to the
algorithms currently included in the \pkg{openVA} package, but usually
provide improved inference. It is future work to include some of these
latest developments to the \pkg{openVA} package for routine use. Most of
these methods have publicly available implementations, such as the
\pkg{farva} package [Moran, 2020] on Github. In another direction of
research, a class of transfer learning methods focuses on building
models to correct for bias in existing methods when predicting out of
domain. These methods take the predicted cause of death assignments and
distributions obtained with the \pkg{openVA} package and learn an
ensemble of the predictions calibrated to the target population
[Datta \emph{et~al.}, 2020, Fiksel \emph{et~al.}, 2021]. The \pkg{calibratedVA} package
[Fiksel and Datta, 2018] is available to implement these models.

Outside of research community that develops new VA algorithms,
\pkg{openVA} has also been used extensively by governments and health
care organizations, particularly in locations that lack a strong vital
registration system and use VA data to identify the leading causes of
death. To facilitate the work of these users, \pkg{openVA} has been
wrapped into a web application, the openVA App [Thomas, 2021], using
the \pkg{shiny} package. The open source openVA App is available on
GitHub and provides an intuitive interface to \pkg{openVA} that does not
require one to learn \proglang{R}, but provides visual and tabular
output produced by the different VA algorithms. It also runs the
official Tariff algorithm by calling the \proglang{Python} source code
of SmartVA-Analyze [Institute for Health Metrics and Evaluation}, 2021] and processing the output to be
consistent with the other algorithms. The \pkg{openVA} \proglang{R}
package has also been a key component in a larger data analysis pipeline
that pulls VA data from an Open Data Kit (ODK) Collect server and
deposits the assigned causes of death to another server running the
District Health Information Software 2 (DHIS2), which is a common
information management system used in low and middle-income countries.
This open source software is implemented as a \proglang{Python} package,
\pkg{openva-pipeline}, and is available on GitHub and the Python Package
Index [Thomas \emph{et~al.}, 2021]. Finally, the \proglang{R} package
\pkg{CrossVA} [Thomas \emph{et~al.}, 2020] and the \proglang{Python} package
\pkg{pyCrossVA} [Choi and Thomas, 2021] provide additional toolkit to convert
raw VA data from its original format from the ODK server to the
standardized formats discussed before. Both packages are open source and
available on GitHub. The \pkg{pyCrossVA} package is also available on
the Python Package Index.

\section{Conclusion}
\label{sec:conclusion}

In this paper, we introduce the \pkg{openVA} package. This is the first
open-sourced software that implements and compares the major VA methods.
The \pkg{openVA} package allows researchers to easily access the most
recent tools that are previously difficult to work with or unavailable.
It also enables the compatibility of multiple data input formats and
significantly reduces the tedious work to pre-process different data
formats specific to each algorithm. The software framework of the
\pkg{openVA} package allows the integration of new methods in the
future. The \pkg{openVA} package makes all the steps involved in
analyzing VA data, i.e., data processing, model tunning and fitting,
result summarization, and evaluation metrics, transparent and
reproducible. This contributes significantly to the public health
community using VA.

Finally, we note that several features that will be helpful for future
development. First, many users of the \pkg{openVA} package may not be
familiar with the command line tools or does not have access to
\proglang{R} on their local machines. A well designed graphical user
interface can be very useful in such settings. The work on the shiny web
application, openVA app, is a first step towards making the package more
accessible. The authors intend to extend it to a better front end hosted
on secure centralized servers. Second, although in this paper, we aim to
provide users with all the available methods for assigning causes of
death, without comparing the accuracy and robustness between them, much
future work is needed to systematically assess, compare, and combine
these methods in a better analysis framework. Finally, the development
of VA algorithms is still an active area of research and it would be
possible to extend the \pkg{openVA} suite to incorporate better VA
algorithms and new types of data such as free-text narratives.

\section*{Acknowledgement}

This work was supported by grants K01HD078452, R01HD086227, and
R21HD095451 from the Eunice Kennedy Shriver National Institute of Child
Health and Human Development (NICHD).

\begin{thebibliography}{35}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Byass, 2015]
Byass P (2015).
\newblock \enquote{{InterVA-4 Software [Windows Executable]}.}
\newblock \emph{www.interva.net}.

[Byass, 2018]
Byass P (2018).
\newblock \enquote{{InterVA-5 Software [Windows Executable]}.}
\newblock \emph{www.interva.net}.

[Byass \emph{et~al.}, 2012]
Byass P, Chandramohan D, Clark SJ, D'Ambruoso L, Fottrell E, Graham WJ, Herbst
AJ, Hodgson A, Hounton S, Kahn K, \emph{et~al.} (2012).
\newblock \enquote{Strengthening Standardised Interpretation of Verbal Autopsy
Data: The New {InterVA-4} Tool.}
\newblock \emph{Global Health Action}, \textbf{5}.

[Choi and Thomas, 2021]
Choi E, Thomas J (2021).
\newblock \enquote{Prepare data from WHO and PHRMC instruments for verbal
autopsy algorithms.}
\newblock \url{https://pypi.org/project/pycrossva/}; Accessed July 21, 2021.

[5]
D'Ambruoso L, Boerma T, Byass P, Fottrell E, Herbst K, Kallander K, Mullan Z
\newblock \enquote{{The case for verbal autopsy in health systems
\newblock \emph{LANCET GLOBAL HEALTH}, \textbf{{5}}({1}), {E20--E21}.
\newblock ISSN {2214-109X}.

[Datta \emph{et~al.}, 2020]
Datta A, Fiksel J, Amouzou A, Zeger SL (2020).
\newblock \enquote{Regularized Bayesian transfer learning for population-level
etiological distributions.}
\newblock \emph{Biostatistics}, p. ISSN 1465–4644.

[Fiksel and Datta, 2018]
Fiksel J, Datta A (2018).
\newblock \emph{calibratedVA: Locally calibrated cause specific mortality
fractions using verbal autopsy data}.
\newblock \urlprefix\url{https://github.com/jfiksel/CalibratedVA}.

[Fiksel \emph{et~al.}, 2021]
Fiksel J, Datta A, Amouzou A, Zeger S (2021).
\newblock \enquote{Generalized {Bayesian} Quantification Learning.}
\newblock \emph{Journal of the American Statistical Association}, pp. 1--19.

[Garenne, 2014]
Garenne M (2014).
\newblock \enquote{Prospects for automated diagnosis of verbal autopsies.}
\newblock \emph{BMC Medicine}, \textbf{12}(1), 18.

[Institute for Health Metrics and Evaluation}, 2021]
{Institute for Health Metrics and Evaluation} (2021{\natexlab{a}}).
\newblock \enquote{{SmartVA-Analyze} Desktop Application.}
\newblock \url{https://github.com/ihmeuw/SmartVA-Analyze}; Accessed July 21,

[Institute for Health Metrics and Evaluation}, 2021]
{Institute for Health Metrics and Evaluation} (2021{\natexlab{b}}).
\newblock \enquote{Verbal Autopsy tools.}
\newblock \url{http://www.healthdata.org/verbal-autopsy/tools}; Accessed July
21, 2021.

[James \emph{et~al.}, 2011]
James SL, Flaxman AD, Murray CJ, {Consortium Population Health Metrics
Research} (2011).
\newblock \enquote{Performance of the Tariff Method: validation of a simple
additive algorithm for analysis of verbal autopsies.}
\newblock \emph{Population Health Metrics}, \textbf{9}(31).

[Kahn \emph{et~al.}, 2012]
Kahn K, Collinson MA, Gómez-Olivé FX, Mokoena O, Twine R, Mee P, Afolabi SA,
Clark BD, Kabudula CW, Khosa A, Khoza S, Shabangu MG, Silaule B, Tibane JB,
Wagner RG, Garenne ML, Clark SJ, Tollman SM (2012).
\newblock \enquote{Profile: {Agincourt} Health and Socio-Demographic
Surveillance System.}
\newblock \emph{International Journal of Epidemiology}, \textbf{41}(4),

[Kunihama \emph{et~al.}, 2020]
Kunihama T, Li ZR, Clark SJ, McCormick TH (2020).
\newblock \enquote{{Bayesian} factor models for probabilistic cause of death
assessment with verbal autopsies.}
\newblock \emph{The Annals of Applied Statistics}, \textbf{14}(1), 241--256.

[Li \emph{et~al.}, 2021]
Li ZR, McCormick T, Clark S (2021{\natexlab{a}}).
\newblock \emph{\pkg{openVA}: Automated Method for Verbal Autopsy}.
\newblock R package version 1.0.12,
\urlprefix\url{https://cran.r-project.org/web/packages/openVA/index.html}.

[Li \emph{et~al.}, 2014]
Li ZR, McCormick TH, Clark SJ (2014).
\newblock \enquote{InterVA4: An {R} package to analyze verbal autopsy data.}
\newblock \emph{Center for Statistics and the Social Sciences Working Paper,

[Li \emph{et~al.}, 2018]
Li ZR, McCormick TH, Clark SJ (2018).
\newblock \emph{\pkg{Tariff}: Replicate Tariff Method for Verbal Autopsy}.
\newblock \proglang{R} package version 1.0.5,
\urlprefix\url{http://CRAN.R-project.org/package=Tariff}.

[Li \emph{et~al.}, 2020]
Li ZR, McCormick TH, Clark SJ (2020).
\newblock \enquote{Using {Bayesian} latent {Gaussian} graphical models to infer
symptom associations in verbal autopsies.}
\newblock \emph{Bayesian Analysis}, \textbf{15}(3), 781.

[Li \emph{et~al.}, 2021]
Li ZR, McCormick TH, Clark SJ (2021{\natexlab{b}}).
\newblock \emph{\pkg{InSilicoVA}: Probabilistic Verbal Autopsy Coding with
'InSilicoVA' Algorithm}.
\newblock \proglang{R} package version 1.3.4,
\urlprefix\url{http://CRAN.R-project.org/package=InSilicoVA}.

[Li \emph{et~al.}, 2019]
Li ZR, McCormick TH, Clark SJ, Byass P (2019).
\newblock \emph{\pkg{InterVA4}: Replicate and Analyse 'InterVA4'}.
\newblock \proglang{R} package version 1.7.6,
\urlprefix\url{http://CRAN.R-project.org/package=InterVA4}.

[McCormick \emph{et~al.}, 2016]
McCormick TH, Li ZR, Calvert C, Crampin AC, Kahn K, Clark SJ (2016).
\newblock \enquote{Probabilistic Cause-of-Death Assignment Using Verbal
\newblock \emph{Journal of the American Statistical Association},
\textbf{111}(515), 1036--1049.

[Miasnikof \emph{et~al.}, 2015]
Miasnikof P, Giannakeas V, Gomes M, Aleksandrowicz L, Shestopaloff AY, Alam D,
Tollman S, Samarikhalaj A, Jha P (2015).
\newblock \enquote{Naive Bayes classifiers for verbal autopsies: comparison to
physician-based classification for 21,000 child and adult deaths.}
\newblock \emph{BMC Medicine}, \textbf{13}(1), 1.

[Moran, 2020]
Moran KR (2020).
\newblock \enquote{FActor Regression for Verbal Autopsy.}
\newblock \url{https://github.com/kelrenmor/farva}; Accessed July 21, 2021.

[Moran \emph{et~al.}, 2021]
Moran KR, Turner EL, Dunson D, Herring AH (2021).
\newblock \enquote{Bayesian hierarchical factor regression models to infer
cause of death from verbal autopsy data.}
\newblock \emph{Journal of the Royal Statistical Society: Series C (Applied

[Murray \emph{et~al.}, 2011]
Murray CJ, Lopez AD, Black R, Ahuja R, Ali SM, Baqui A, Dandona L, Dantzer E,
Das V, Dhingra U, \emph{et~al.} (2011).
\newblock \enquote{Population Health Metrics Research Consortium gold standard
verbal autopsy validation study: design, implementation, and development of
analysis datasets.}
\newblock \emph{Population Health Metrics}, \textbf{9}(1), 27.

[26]
Nichols EK, Byass P, Chandramohan D, Clark SJ, Flaxman AD, Jakob R, Leitao J,
Maire N, Rao C, Riley I, Setel PW, Grp WVAW ({2018}).
\newblock \enquote{{The WHO 2016 verbal autopsy instrument: An international
standard suitable for automated analysis by InterVA, InSilicoVA, and Tariff
\newblock \emph{PLOS Medicine}, \textbf{{15}}({1}).
\newblock ISSN {1549-1676}.

[R Core Team}, 2020]
{R Core Team} (2020).
\newblock \emph{R: A Language and Environment for Statistical Computing}.
\newblock R Foundation for Statistical Computing, Vienna, Austria.
\newblock \urlprefix\url{https://www.R-project.org/}.

[Serina \emph{et~al.}, 2015]
Serina P, Riley I, Stewart A, James SL, Flaxman AD, Lozano R, Hernandez B,
Mooney MD, Luning R, Black R, \emph{et~al.} (2015).
\newblock \enquote{Improving performance of the Tariff Method for assigning
causes of death to verbal autopsies.}
\newblock \emph{BMC Medicine}, \textbf{13}(1), 1.

[Taylor \emph{et~al.}, 1983]
Taylor CE, Parker R, Reinke W, Faruquee R, \emph{et~al.} (1983).
\newblock \emph{Child and maternal health services in rural India. The
Narangwal experiment. 2. Integrated family planning and health care.}
\newblock Johns Hopkins University Press.

[Thomas, 2021]
Thomas J (2021).
\newblock \enquote{Software for automating the processing of verbal autopsy
\newblock \url{https://github.com/verbal-autopsy-software/openVA_App}; Accessed
July 21, 2021.

[Thomas \emph{et~al.}, 2020]
Thomas J, Choi E, Li Z, Maire N, McCormick T, Byass P, Clark S (2020).
\newblock \emph{\pkg{CrossVA}: Verbal Autopsy Data Transformation for
{InSilicoVA} and {InterVA5} Algorithms}.
\newblock R package version 1.0.0,
\urlprefix\url{https://CRAN.R-project.org/package=CrossVA}.

[Thomas \emph{et~al.}, 2021]
Thomas J, Clark SJ, Bratschi MW (2021{\natexlab{a}}).
\newblock \enquote{Software for automating the processing of verbal autopsy
\newblock \url{http://openva-pipeline.readthedocs.io/}; Accessed July 21, 2021.

[Thomas \emph{et~al.}, 2021]
Thomas J, Li Z, Byass P, McCormick T, Boyas M, Clark S (2021{\natexlab{b}}).
\newblock \emph{{InterVA5}: Replicate and Analyse {'InterVA5'}}.
\newblock R package version 1.1.3,
\urlprefix\url{https://CRAN.R-project.org/package=InterVA5}.

[Wen \emph{et~al.}, 2018]
Wen R, Miasnikof P, Giannakeas V, Gomes M (2018).
\newblock \emph{\pkg{nbc4va}: Bayes Classifier for Verbal Autopsy Data}.
\newblock \proglang{R} package version 1.1,
\urlprefix\url{https://github.com/rrwen/nbc4va}.

[World Health Organization and others}, 2012]
{World Health Organization and others} (2012).
\newblock \enquote{Verbal Autopsy Standards: The 2012 WHO Verbal Autopsy
Instrument Release Candidate 1.}
\newblock \emph{Technical report}, World Health Organization (WHO), Geneva.

\end{thebibliography}

\appendix

\section{New changes to the InSilicoVA algorithm}
\label{app:insilico}

Since the original publication of McCormick \emph{et~al.}, 2016, several major
changes have been introduced in the InSilicoVA algorithm implementation
to improve cause-of-death assignments. Two major changes are summarized

\subsection{Removal of physically impossible causes}\label{sec:impossible}

The originally proposed InSilicoVA assumes all causes of death are
possible for each observation. The impact form such assumption is mild
when data are abundant, but could be problematic when either sample size
is small or the proportion of missing data is high. In both cases,
physically impossible causes might get assigned with non-ignorable
posterior mass. Since the version 1.1.5 of \pkg{InSilicoVA}, the
algorithm automatically checks and removes impossible causes before
fitting the model. The \(k\)-th cause is defined as physically
impossible for the \(i\)-th death if \(P(s_{ij}=1|y_{i}=k) = 0\) for any
\(j\) representing either sex or age group that presents. We then
consider a cause to be physically impossible to the underlying
population if it is impossible for all the observations of the input
data. For example, with the new implementation, CSMF for
pregnancy-related causes will not be estimated if the input data consist
of only male deaths.

\subsection{Structured symptom dependence}\label{sec:dependence}

Another change in InSilicoVA algorithm is the implementation of the data
consistency check. Inherited from the InterVA-4 algorithm, InSilicoVA
performs a data consistency check before fitting the algorithm as
described in Algorithm \ref{alg:check}. This consistency check consists
of two parts. First, it removes indicators whose presence contradicts
the existence of a lower-level indicator. For example, if the question
\textit{pregnant at the time of death} is answered with an `Yes', then
the algorithm makes sure that the item \textit{child} is set to `No'.
The second part of the algorithm checks the existence of each
higher-level indicator when lower-level indicators that implied by it
exists. Again, if the question \textit{pregnant at the time of death} is
answered with an `Yes', then the algorithm makes sure that the item
\textit{female} is set to `Yes' as well.

\begin{algorithm}
\caption{Data Check}
\label{alg:check}
\begin{algorithmic}[1]

\Statex

\Require
\Statex $\mathbf{s}_j$ = vector of $j=1 \mbox{ to } 245$ indicators for one death \vspace{3pt}
\Statex $\mbox{notask}(s)$: higher-level indicator(s) whose presence ensures indicator $s$ does not exist \vspace{3pt}
\Statex $\mbox{anc}(s)$: higher-level indicator(s) that must exist if indicator $s$ exists \vspace{3pt}

\Statex

\Ensure $\mathbf{s}^*$ = modified vector of $j=1 \mbox{ to } 245$ indicators for one death \vspace{3pt}

\Statex

\For{each indicator $j = 1 \mbox{ to } 245$}
\State Set $ \mathbf{s}^*_j \leftarrow \mathbf{s}_j$ \Comment{Initialize}
\EndFor

\For{each indicator $j = 1 \mbox{ to } 245$} \label{line:start}
\If{there exists $\mathbf{s}^*_{\ell}$ in $\mbox{notask}(\mathbf{s}^*_j)$ and $\mathbf{s}^*_{\ell}$ is `yes'}\label{line:firstif}
\State $\mathbf{s}^*_j \leftarrow$ `no' \Comment{This indicator is nonsensical, set to `No'} \label{line:changetomissing}
\EndIf
\If{there exists $\mathbf{s}^*_{\ell}$ in $\mbox{anc}(\mathbf{s}_j)$ and $\mathbf{s}^*_j$ is `yes'}\label{line:secondif}
\State $\mathbf{s}^*_\ell \leftarrow$ `yes' \Comment{Set more general version of this indicator to `yes'}
\EndIf
\EndFor \label{line:end}
\State Repeat loop in lines \ref{line:start} to line \ref{line:end} again \Comment{Processes 2-level hierarchies}

\Statex

\end{algorithmic}
\end{algorithm}

This data check procedure ensures the collected data are internal
consistent, yet it might introduce additional complications under the
common assumption that symptoms are conditionally independent given any
cause of death. This can be illustrated with a simple example below. In
calculating the conditional probability of observing two symptoms, e.g.,
`infant', and `recent abortion', we can decompose
\begin{eqnarray}\nonumber
\mbox{Pr}(\mbox{infant} = Y \;\&\; \mbox{recent abortion} = N \;|\; \mbox{some cause})
\mbox{Pr}(\mbox{infant} = Y \;|\; \mbox{some cause})\\\nonumber
\mbox{Pr}(\mbox{recent abortion} = N \;|\; \mbox{some cause}).
\end{eqnarray} However, such independence does not hold since the two
symptoms are mutually exclusive, i.e., the presence of symptom `infant'
implies no recent abortion. In fact, the structure of the WHO instrument
assures that the lower-level question will not be asked at all if the
interview is about an infant death. Thus in practice, the above
probability calculation yields an underestimate of the target value and
it should be calculated simply with \[
\mbox{Pr}(\mbox{infant} = Y \;\&\; \mbox{recent abortion} = N \;|\; \mbox{some cause})
\mbox{Pr}(\mbox{infant} = Y \;|\; \mbox{some cause})
\] To reflect the different joint probability decompositions when the
known symptom hierarchy is available, it turns out we can easily modify
the data check algorithm by imputing missing instead of absence to the
lower level symptoms in line \ref{line:changetomissing}. We have found
that such a change usually yield more reasonable CSMF estimates,
especially for child and neonate deaths, in practice. Also it should be
noticed that this problem stems from the conditional independence
assumption that has been adopted in all other methods besides
InSilicoVA. The computational trick here only eliminates the bias from
mutually exclusive symptoms known from the hierarchical structure of VA
questionnaires. More systematic ways to deal with symptom correlations
are needed to further improve the estimation.

With the update of the WHO VA questionnaire to the 2016 format, three
major changes were made to the data check procedure implemented by
InterVA-5. First, similar to InSilicoVA, InterVA-5 also imputes missing
instead of absence to the lower level symptoms in line
\ref{line:changetomissing}. Second, the hierarchical structure of the
2016 VA questionnaire makes use of both ``Yes'' and ``No'' (or ``Y'' and
``N'') as substantive values that trigger the asking of lower-level
questions. For example, if the indicator/question ``Did the baby ever
cry?'' has the response ``No'', then the following indicator/question,
``Did the baby cry immediately after birth, even if only a little bit?''
should not be asked -- as implemented with the notask(\(s\)) function in
Algorithm \ref{alg:check}. Thus, additional steps are needed in the
InterVA-5 data check which test the value of each indicator with the
appropriate substantive response that implies subsequent action (i.e.,
asking or not asking a different indicator/question) in line
\ref{line:firstif} and \ref{line:secondif}. The last major change is
that an additional check is made for indicators/questions that only
apply to neonates. If a record for a non-neonate contains a value for an
indicator/question that only applies to neonates, then that value is
cleared in the working copy of the data file. The same check rules are
used by InSilicoVA for WHO 2016 data as well. The implementation is
slightly different between InSilicoVA and InterVA-5. InterVA-5 updates
\(s_j^*\) to be missing in line \ref{line:changetomissing} only when it
is recorded to be the corresponding substantive value, since the
non-substantive value is considered the same as missing values in the
calculation. InSilicoVA, on the other hand, make the update regardless
of the value of \(s_j^*\), since the non-substantive value and missing
are treated differently in the algorithm.

\end{document}


# Servicifying zk-SNARKs Execution for
Verifiable Off-chain Computations ††thanks:

Alvaro Alonso Domenech, Jonathan Heiss, Stefan Tai  Information System
Engineering
TU Berlin
Berlin, Germany
<EMAIL_ADDRESS>
###### Abstract

Zk-SNARKs help scale blockchains with Verifiable Off-chain Computations (VOC).
zk-SNARK DSL toolkits are key when designing arithmetic circuits but fall
short of automating the subsequent proof-generation step in an automated
manner. We emphasize the need for portability, interoperability, and
manageability in VOC-based solutions and introduce a Proving Service that is
designed to provide a scalable and reusable solution for generating zk-SNARK
proofs leveraging clouds.

###### Index Terms:

Blockchain, Zero-knowledge Proofs, zk-SNARKs, ZoKrates, Proving, Service
Architecture, Cloud Computing

##  I Introduction

The succinctness property and the short verification times of zk-SNARK s come
at the cost of large computational complexity, static execution models, and
memory overheads during the proof generation process [1]. These limitations
represent a problem in VOC applications which require handling large and
varying workloads, like rollups. This type of application would benefit from
scalable, interoperable, and manageable system environments like clouds.
However, DSL toolkits like ZoKrates [2] or Circom [3] concentrate on circuit
development and put little emphasis on how to integrate these circuits into
production systems. Bridging the gap between advanced cryptography and modern
systems engineering also helps to create standardized benchmarks for proving
systems and tools.

##  II Proving Service

We introduce a service-oriented approach for VOC that facilitates the use of
the cryptographic procedures of zk-SNARK s within cloud system architectures.
Our system allows for executing arithmetic circuits as encapsulated
application logic in containers that are deployable to different machines
realizing scalability, provide interoperability with other services, and are
better manageable. For that, we give an overview of the service-oriented
system architecture and describe the internals of the proving service.

###  II-A Service Oriented Verifiable Off-chain Computation

The problems above can be fulfilled through a service-oriented approach as
depicted in Figure 1. Starting from a higher-level system’s perspective, we
treat the proving service as a black box which, upon a request, returns the
proof together with the computation’s output. Following the VOC model [2], we
distinguish between the blockchain infrastructure hosting the verifier
contract and a cloud-based off-chain infrastructure that runs outside the
consensus protocol and hosts the consumer and the proving service.

Figure 1: Proving Service Model

The consumer service is responsible for managing the proving service’s inputs
and outputs and interacting with the relying verifier contract. It receives
data from external sources and translates them into a proof request. Upon a
request through the prover client, the proving service executes the VOC and
returns the ZKP attesting to the computational integrity of the VOC. The
consumer service then submits the ZKP to the verifier contract through its
verifier client. On submission, the verifier contract verifies the proof
computed by the proving service.

Figure 2: Proving Service Architecture Figure 3: Proving Time and Memory
Consumption

###  II-B Service Architecture

The Proving Service depicted in Figure 2 is an application service that runs
stand-alone on the prover’s off-chain infrastructure. It serves proof requests
by exposing the proving-related operations through an Application Program
Interface (API) as minimal functionality to the consumer service.

The procedure can be summarized in four simple steps: First, the proving
service receives the Proof Request containing the public ($x$) and private
($x^{\prime}$) proof arguments and the identifier of the addressed Executable
Constraint System ($ecs$)($id$). Second, using the $id$ the API fetches the
corresponding $ecs$ and private key ($pk$) from the Proof Registry which is
the persistent storage component containing these large, recurrently requested
files needed for proving. The proof registry manages different reusable pairs
of $ecs$ and $pk$, each addressable through a unique $id$.

Third, the Proving Instance is executed in two stages: the witness
computations and the proof generation. Fourth, the ZKP is returned to the
consumer service through the API upon successful execution.

##  III ZoKrates-API

The previous service-oriented architecture serves as a technology-agnostic
blueprint for building proving services for DSL circuits. For evaluation, we
technically instantiate the proving service for ZoKrates [2] and present the
ZoKrates-API111https://github.com/ZK-Plus/zokrates-api as a ready-to-use open-
source software. We servicify ZoKrates by wrapping an API around the ZoKrates
interpreter, the central component of the ZoKrates software that previously
has only been addressable through a Command Line Interface and a Javascript
library. The ZoKrates-API exposes the methods of the ZoKrates interpreter
through HTTP endpoints. We containerized the ZoKrates-API using Docker, making
the services easily deployable among a wide range of machines and allowing us
to further leverage cloud-native tools like Kubernetes for horizontal
scalability, manageability, and observability. Furthermore, the ZoKrates-API
supports multi-threading so a single instance can compute multiple proofs in
parallel.

##  IV Evaluation

To test the presented implementation, we deployed the containerized ZoKrates-
API on a Kubernetes cluster. As a workload for our experiments, we generated a
large number of EdDSA signatures[4] which amounts for more than $2^{29}$ of
circuit constraints similar to [5].

We conducted three experiments (see Figure 3) to measure the average proving
time per signature in [sec] and the memory consumption in [gb] using various
cluster configurations. Figure 3A) shows a 33% improvement in proving time for
a single machine when choosing an appropriate machine size. An argument in
favour of vertical escalability. For a single machine as well, Figure 3B)
shows that the proving time can be brought down drastically when few parallel
threads are enabled, though the gains plateaued rapidly due to the increasing
resources needed. Running the same experiments in parallel VMs instead of
threads, Figure 3C) demonstrate a better approach to scaling proving as the
computational burden is distributed over several machines, proportional
increasing the processing time.

##  V conclusion

As the experiments demonstrate, significant performance improvements can be
gained from horizontal (more nodes) and vertical (larger nodes) scalability of
an arbitrary zk-SNARK proof. We facilitate this process by leveraging modern
cloud-native architectures such as Docker and Kubernetes.

## References

  * [1] E. Ben-Sasson, A. Chiesa, E. Tromer, and M. Virza, “Scalable zero knowledge via cycles of elliptic curves,” _Algorithmica_ , vol. 79, pp. 1102–1160, 2017.
  * [2] J. Eberhardt and S. Tai, “Zokrates - scalable privacy-preserving off-chain computations,” in _IEEE International Conference on Blockchain_. IEEE, 2018.
  * [3] M. Bellés-Muñoz, M. Isabel, J. L. Muñoz-Tapia, A. Rubio, and J. Baylina, “Circom: A circuit description language for building zero-knowledge applications,” _IEEE Transactions on Dependable and Secure Computing_ , 2022.
  * [4] J. Heiss, A. Busse, and S. Tai, “Trustworthy pre-processing of sensor data in data on-chaining workflows for blockchain-based iot applications,” in _Service-Oriented Computing: 19th International Conference, ICSOC 2021, Virtual Event, November 22–25, 2021, Proceedings 19_. Springer, 2021, pp. 133–149.
  * [5] A. Chiesa, R. Lehmkuhl, P. Mishra, and Y. Zhang, “Eos: Efficient private delegation of zksnark provers,” in _USENIX Security Symposium. USENIX Association_ , 2023.