approach.qmd

---
execute: 
  echo: false
---
# APPROACH

## Data description

We will combine data from multiple influenza vaccination cohort studies for our
analyses. At the time of writing, we have data from two cohort studies conducted
by Ted Ross, one study conducted by Andrea Sant and Angela Branche (see
@gouma2020a), and a collaboration with Ben Cowling. Ben is the PI of several
similar studies and has offered to share data with us. The two studies conducted
by Ted Ross provide a unique opportunity to study heterologous antibody
responses to influenza vaccination, as a wide historical panel of assays was
conducted for each individual. Combining Ross' data with data from Sant and
Cowling will allow for us to compare the immunogenicity of multiple vaccines,
and increase the power of our analyses of host factors.

We will refer to both of Ross' cohort studies conjointly as the **UGAFluVac**
study. This data set consists of longitudinal HAI measurements taken at three
different study sites. From fall 2013 to spring 2016, participants were
recruited at the study site in either Pittsburgh, PA or Stuart, FL. Sample
collection is detailed in @nunez2017. Briefly, the study is a prospective open
cohort design. Participants were adults aged 18 and up who were allowed to
repeat each year, and the data includes a unique ID per participant that allows
for the identification of longitudinal measurements. Each participant received a
pre-vaccination blood draw, and was then administered a split-inactivated
standard dose (SD) Fluzone seasonal influenza vaccine (Sanofi Pasteur). Patients
aged 65 or older could opt to receive Fluzone High Dose (HD) instead. At the PA
study site in 2013, some patients were administered Fluzone intradermal rather
than the standard intramuscular Fluzone formulation. Followup whole blood draws
were targeted for 21 days post-vaccination.

Processed sera were used for HAI assays following standard protocols. HAI assays
were conducted using the homologous strain and a panel of heterologous strains.
The starting dilution was 1:10 and assays which did not agglutinate at the
starting dilution were coded as 1:5. For all further analyses, we will use the
reciprocal titer transformed as $\text{analysis titer} = \log_2 \left(\text{raw
reciprocal titer} / 5 \right).$ This transformation serves to set 0 as the limit
of detection on the log scale. Additionally, the following data were collected
from patients by a survey prior to vaccination: year of birth, age, gender (the
covariate is listed as gender but is coded as male/female only), and
race/ethnicity.

In January 2017, the study moved to the University of Georgia in Athens, GA. The
paper @wu2022 contains a description of the study, but at the time of writing,
no published paper contains a complete description of the study cohort. The
study design was similar, with additional covariate information collected:
complete date of birth, sex assigned at birth, race/ethnicity, BMI, height,
weight, and questions about smoking and comorbidities. Beginning in the
2017-2018 flu season, participants aged 10 and older were also
recruited for the study. The HAI assays were conducted in the same way. The
post-vaccination time point target was changed to Day 28 beginning in fall 2018.
Finally, additional subcohorts were administered other vaccines at certain
points during the study. See @fig-michael for details.

```{r}
#| label: fig-michael
#| fig-width: 6
#| fig-cap: "Study diagram showing the subcohorts and associated time points
#| collected during the ongoing UGA study. The diagram was created by Michael
#| Carlock, the CIVR-HRP research director."
knitr::include_graphics(here::here("figures", "study-diagram.png"))
```

We will refer to the data provided by Andrea Sant and Angela Branche as the
**RocFluVac** study. We have not yet finished cleaning and processing the data,
but it contains similar HAI measurements for homologous strains only,
participants were administered a variety of vaccines (including Fluzone HD
administered to participants aged 18 - 49), and also contains data on
T cell responses to the vaccine. See @gouma2020a for one brief
explanation of the data.

## Aim 1

For this aim, we will use the cohort data collected by Ted Ross.

### Preliminary results

We have calculated several metrics of antigenic distance, and will compare the
year-based method [@auladell2022], the dominant *p*-epitope sequence based
method [@gupta2006], and a distance based on antigenic cartography [@cai2012;
@fonville2014]. To compute the dominant *p*-epitope distance, the pairwise
Hamming distances are computed between the amino acid sequences of each of the
five hemagglutinin head epitopes, and these are divided by the length of their
respective sequences. The dominant *p*-epitope measurement is the maximum of
each of these probabilities, representing the different between the two strains
at their most different epitope.

Antigenic cartography analyses were conducted by Amanda Skarlupka, PhD, who
continues to work with us on this project. In short, antigenic cartography uses
a matrix where each row represents an individual in the study and each column
represents an influenza strain. The cells of the matrix are populated with
individual $i$'s titer to strain $j$. Multidimensional scaling (MDS) is used to
reduce the matrix to a specified column dimension, while minimizing the change
in the Euclidean distances between measurements. After performing iterative MDS
on a variety of target column dimensions, we found that two dimensions was
satisfactory, and our maps were similar to those in @fonville2014. After
calculating the MDS maps, we obtained antigenic distances as the Euclidean
distance between map coordinates of the vaccine strain and the other strains in
the panel following the method of @cai2012. For the purposes of this analysis,
we ignored longitudinal measurements between individuals, and treated each
observation of an individual as a unique measurement. Finally, all of the
antigenic distance were normalized *for each vaccine strain*. Therefore, each
vaccine strain had a distance of 0 with itself, and the most distant historical
strain had a distance of 1 with that strain. (We normalized all three of the
distance measurements in this way.)

After computing the normalized distance measurements, we fit simple linear
regression models with either post-vaccination titer (measuring the absolute
immune response post-vaccination) or fold-change in titer (measuring the
relative boost post-vaccination) against antigenic distance. @fig-amandafits
shows sample linear regression models following this protocol for two
vaccine strains.

Our proposed metrics for evaluating the **magnitude**, **breadth**, and
**overall strength** are the intercept of the regression line, slope of the
regression line, and area under the regression line respectively. The intercept
measures how strong the response is to the homologous strain, the slope
describes how the response weakens linearly as antigenic distance increases, and
the AUC provides a framework for combining the magnitude and breadth measurements
into one measurement of overall strength of the vaccine-induced response.

```{r}
#| label: fig-amandafits
#| out-width: 6.25in
#| fig-cap: "Fitted models with titer outcomes and antigenic distance as the
#| only predictor. The metrics reported are for titer increase for simplicity."

knitr::include_graphics(here::here("Figures", "dist-linear-fits.png"))
```

As a case study to motivate our framework, we plan to analyze differences
in these three metrics between the Fluzone SD and Fluzone HD vaccines, to determine
if the high dose vaccine elicits a stronger or broader response in our cohort, as shown in @fig-distdose. Our preliminary results suggests that this effect
differs qualitatively across vaccines, and we plan to analyze which vaccines show
the most notable differences between strains.

```{r}
#| label: fig-distdose
#| out-width: 6.25in
#| fig-cap: "Linear models for titer increase vs. antigenic distance, stratified
#| by vaccine dose. The same metrics are reported for both doses."

knitr::include_graphics(here::here("Figures", "dist-dose-fits.png"))
```

### Proposed studies

Going forward, we plan to consider more flexible models than the simple linear
regression models we used in this step of the analysis. We will consider
restricted spline GAMs and LOESS to fit a potentially nonlinear effect of
antigenic distance, which we would expect in the presence of original antigenic
sin if our panel contains strains which are distant enough to confer no
cross-protection. We will also implement Bayesian multilevel linear and spline
models, in order to incorporate between-subject variability into the overall
model fit.

We will then compare our metrics to traditional metrics, notably the mean titer
increase, HAI composite score [@abreu2020], and proportion of seroconverted
strains, all calculated for each individual. In order to determine whether
our proposed metrics are more robust than traditional metrics, we will
subsample measurements from our cohort in order to mimic the use of
different virus panels across labs. For $k$ labs, we will subsample $p$ strains
from our panel, plus the homologous strain. For each of these panels, we will
calculate the metrics for each individual, and then analyze the variability of
mean metrics across labs.

Finally, we will explore weighting schema for our metrics. The unweighted
AUC assigns equal weight to strains of all antigenic distances, but perhaps
we would prefer to weight the response to distant strains higher or
lower than similar strains--for example, if we are primarily considered with
boosting the response to a specific pandemic strain, we could assign less
weight to distant strains. But if we are considering candidates for a
broadly reactive vaccine, we could weight distant strains higher in order to
favor vaccine candidates that induce responses to the most distant strains.
Examples of weighting schema are shown in @fig-weights.

```{r}
#| label: fig-weights
#| out-width: 6.25in
#| fig-cap: "Example weighting schema for AUC measurements. This figure shows
#| the responses to the H1N1-California-2009 vaccine in the UGAFluVac study.
#| In the linear weighting
#| scheme, the titer increase is multipled by the reciprocal antigenic distance
#| in order to construct weights that decrease linearly as strains become
#| more distance. In the 2
#| antigenic unit scheme, all strains which are within two antigenic units
#| on the cartography map (or equivalent with the other distances) are
#| weighted equally, while strains that are further away are not included
#| in calculating the AUC."

knitr::include_graphics(here::here("Figures", "Weighting-schema.png"))
```

### Expected outcomes

Our primary outcome for this goal will be a set of metrics that we propose
for the evaluation of broadly reactive vaccine candidates. These metrics
will be calculated on the UGAFluVac data and will be supported by our
subsampling analysis to estimate the robustness of our metrics. Developing these
metrics will improve our understanding of the functional relationship between
antigenic distance and immune response. Finally, our robustness analysis
will provide insight into the amount of error induced into metrics for
post-vaccination immune response by variability within virus panels used
across different research groups.

## Aim 2

For this aim, we will use the UGAFluVac data, specifically the portion
collected at the University of Georgia. The PA/FL data does not have
prior vaccination history and thus we can only use the UGA data for this aim.

### Preliminary results

We first explored the homologous case in order to identify simple first-order
effects of each of the covariates of interest. @fig-mada-eda shows plots of fold
change in titer (titer increase) plotted against available covariates. The only
covariate which independently explained a substantial proportion of the variance
in the outcome was prevaccination titer. Note that these results do not consider
interaction effects (or subgroup specific effects that can be modeled as
interactions): for example, we would not expect to see an effect of dose in the
entire population since otherwise healthy adults under the age of 65 were not
offered high dose vaccine. The $R^2$ metrics also only capture linear trend--if
a predictor acts through a primarily nonlinear (especially true for nonmonotonic
effects) the $R^2$ will drastically underestimate the strength of association
between the two variables.

```{r}
#| label: fig-mada-eda
#| out-width: 6.25in
#| fig-cap: "Homologous titer increase to the
#| influenza A(H1N1) Michigan-2015 strain in the UGAFluVac data, plotted
#| against a variety of covariates. The size of the point corresponds to the
#| number of overlapping measurements at that location. The purple vertical
#| lines indicate medians, and the horizontal purple lines mark the 2.5th and
#| 97.5th quantiles. While many of the univariate associations with titer
#| increase are weak or difficult to see, many of these effects become
#| much stronger when interactions are taken into effect. These trends
#| also vary widely across other strains which are not shown."

knitr::include_graphics(here::here("Figures", "Single_var_titer_panel.png"))
```

From our exploratory analyses, we also know that the effect of prevaccination
titer on the homologous immune response varies across vaccine strain and dose
(@fig-pvt-lm). We fit stratified linear models with both parallel slopes (A and
C) and varying slopes (B and D) with these two factors. We further expect to see
a third-order interaction between vaccine strain and dose. Note that these
preliminary figures show all individuals, and the effect of dose will be further
modified by age, since high dose vaccine was only offered to elderly
participants.

```{r}
#| label: fig-pvt-lm
#| out-width: 6.25in
#| fig-cap: "Titer increase plotted against pre-vaccination titer. Models
#| are stratified by vaccine strain type (A and B) or by vaccine dose (C and D).
#| Figures A and C show a parallel slopes model where only the intercepts
#| are allowed to vary (corresponding to the model with no interaction
#| term, only main effects for both variables), while Figures B and D
#| show the varying slopes (interaction and main effects) models. The
#| varying slopes models improve R-squared marginally, but there are
#| several limitations to these analyses that may hide larger
#| differences."

knitr::include_graphics(here::here("Figures", "prevac_titer_lm_panel.png"))
```

We further conducted an analysis at the strain-specific level, fitting
separate models for each combination of vaccine strain and assay strain
represented in the data (@fig-ns-panel).

```{r}
#| label: fig-ns-panel
#| out-width: 6.5in
#| fig-cap: "Relationship between titer increase and prevaccination titer
#| for each of the assay strains which were used as part of a historical panel.
#| Only assays for individuals who received H1N1-Michigan-2015 containing
#| vaccine are shown for simplicity. We can see that the intercept and
#| slope have the largest magnitude for the homologous strain and for
#| the previous vaccine strain (H1N1-California-2009), and decreases with
#| other assay strains."

knitr::include_graphics(here::here("Figures", "michigan-ns-panel.png"))
```

We can further stratify these strain specific analyses by other factors like
dose (@fig-dose-panel). These analyses reveal a non-constant modification of
the effect of prevaccination titer by dose, which differs across vaccines
and across assay strains within vaccines.

```{r}
#| label: fig-dose-panel
#| out-width: 6.5in
#| fig-cap: "The same analyses as shown in @fig-ns-panel, but stratified by
#| vaccine dose. Additionally, only individuals aged 65 and older are included
#| in these models, since younger individuals were not eligible to
#| receive high dose vaccines. Some relationships are distorted by small
#| sample effects, but in general the dose appeared to make little difference
#| for this particular vaccine. However, we see varying, potentially real
#| effects in the MI/15 and CA/09 effects. The effect of dose also
#| varies across vaccine strains (data not shown)."

knitr::include_graphics(here::here("Figures", "michigan-dose-panel.png"))
```

Differences in strain-specific models for the same vaccine demonstrate
the interaction between antigenic distance, prevaccination titer, and
other stratifying factors (like dose and prior vaccination history).
We conducted a preliminary analysis of the strength of this interaction by
compiling the slopes of each model (as shown in @fig-ns-panel, but incorporating
data from all vaccine strains that were used in the study), and plotting these
slopes against the antigenic distance between the vaccine strain and the
assay strain for that model, shown in @fig-slope-plots. Note that in this framework,
vaccines were only compared against assays of the same type. HAI assays for
influenza B strains which predated the divergence of the Victoria and Yamagata
lineages were compared against both B-Yamagata and B-Victoria vaccines.

For the type A influenza assays, we see a strong positive correlation between
antigenic distance and slope of the strain-specific model. As the vaccine
strain and the assay strain become more distant, pre-vaccination titer has
less of an effect on the amount of boosting produced by the vaccine.

```{r}
#| label: fig-slope-plots
#| out-width: 6.25in
#| fig-cap: "Slopes from each of the individual regression analyses (similar to
#| those depicted
#| in @fig-ns-panel) plotted against the dominant p-Epitope antigenic distance
#| between the vaccine strain and the assay strain. We performed these individual
#| regression analyses for every combination of vaccine strain and assay strain
#| in the UGAFluVac data. Point colors in this plot represent different
#| assay strains, while point shapes represent different vaccine strains."

knitr::include_graphics(here::here("Figures", "slopes-plot.png"))
```

We also used machine learning models to predict response to the homologous
vaccine strain based on all of the covariate data that was available to us.
We used permutation variable importance to score the importance of the
included predictors (@fig-vip). We found that prevaccination titer was, by far,
the most important covariate in predicting vaccine boosting. However, we
also observed differences between influenza seasons, individuals of different
ages, and the vaccine strain type. Future models will include antigenic distance
measurements instead of indicators for strain types.

```{r}
#| label: fig-vip
#| out-width: 6.25in
#| fig-cap: "Permutation variable importance (VI) for each feature included in the
#| final tuned random forest model. Permutation VI is calculated by randomly
#| permuting the variable of interest (which should destroy any signal), and
#| subtracting the estimated RMSE with the permuted variable from the
#| estimated RMSE when the variable is not permuted. The model was trained
#| using a subset of 70% of the UGAFluVac data, and these RMSE for
#| both the permutated and non-permuted data were calculated using the
#| holdout set of the remaining 30% of data. Despite using a random forest
#| model, which can detect interactions and nonlinear effects, any effects
#| that interact with age or dose may be artificially dampened in this
#| analysis due to class imbalance."

knitr::include_graphics(here::here("Figures", "rf-vip-edit.png"))
```

### Proposed studies

First, we will expand and refine our preliminary machine learning approach to include
all predictors of interest. We will apply similar models
to quantify the predictive power of each variable using permutation importance.
For the final models, we will perform more in-depth feature engineering to ensure
that predictors have optimal predictive power while retaining biological meaning. These models can learn nonlinear and interaction effects without prior specification and
thus will allow us to understand how
well all we can expect inferential models to capture trends in the outcome. We
will also fit models using both the post-vaccination titer and the titer
increase as outcomes, to determine if the drivers of the raw immune response
are less influenced by pre-vaccination titer. We will also fit separate
models for elderly participants in order to gain a better understanding of the
effect of vaccine dose in this subgroup.

After we understand how much of the post-vaccination titer can be explained
by the data we have, we will build Bayesian hierarchical models to take
advantage of the clustered structure of the data. These hierarchical models
will be allow us to make inferences about the relationships between
variables, rather than only quantify the predictive power. Furthermore, we
can implement random effects in these models to partially absorb unmeasured
confounding, which is not possible within a predictive machine learning
framework. However, we can compare the overall predictive power of our
inferential models to the best predictive models.

Finally, we will modify previously developed mechanistic models [@zarnitsyna2015; @zarnitsyna2016] which model steric hindrance between multiple epitopes of the same antigen. The model for steric hindrance is mathematically similar to how we would model differences in binding avidity between antigenically distance strains, so by incorporating a notion of antigenic similarity, we can use these models to predict how the relationship between
dose, prior immunity, and vaccine response varies with antigenic distance. We will compare the model predictions with the predictions made by our machine learning and inferential models. A notable limitation of the mechanistic models in this case is that we will not be able to make predictions about the effect of serial repeated vaccinations as we hope to do for the hierarchical models.

### Expected outcomes

Our work on this aim will produce fitted inferential models, with estimates
of the strength of the effect for pre-vaccination titer, prior vaccination
history, and other drivers of the immune response. Our models will also
provide estimates for the strength of interactions between these covariates.
We will also produce a set of tuned predictive models. Estimates of the
vaccine response for each individual will be computed from the inferential
and the predictive models, and by comparing the two sets of predictions,
we will estimate the performance of our inferential models. Finally, we will
update the mechanistic models and compare patterns in mechanistic model
results to the patterns observed in our data.

## Aim 3

For this aim, we will combine the UGAFluVac data, RocFluVac data, and
any applicable data provided by Ben Cowling.

### Preliminary results

While working with the UGAFluVac data, we have constructed a small DAG exploring
only the necessary adjustment variables that affect the causal pathway
between dose and post-vaccination titer (@fig-dag). In this reduced framework,
the total and direct causal effects of dose would be the same, and the only
factor we would need to adjust for to obtain the true causal estimate of dose
is age. Although age may be a common cause of several other factors that affect
post-vaccination titer, vaccine dose is a quality of the intervention (vaccination)
and thus the only factors that affect vaccine dose are factors that determine
participation in the study.

```{r}
#| label: fig-dag
#| out-width: 6.25in
#| fig-cap: "A directed acyclic graph (DAG) showing a graphical causal model
#| for the effect of dose (the exposure, in green) on post-vaccination titer
#| (the outcome, in blue with an I symbol). Age is a confounder and is colored
#| pink in the DAG."

knitr::include_graphics(here::here("Figures", "dag-aim3.png"))
```

We have also conducted preliminary comparisons of the UGAFluVac and RocFluVac
data, primarily to compare the patterns in antibody responses observed in
the UGAFluVac data with the T cell responses observed in the RocFluVac data
(@fig-rocted). We also have HAI data for the RocFluVac data, and thus we can
compare the effect of dose on the immune response in the 18 - 49 year olds
in the RocFluVac study with the effect in the 65+ year olds in the UGAFluVac study.

```{r}
#| label: fig-rocted
#| out-width: 6.5in
#| fig-cap: "Antibody immune response measured via HAI in the UGAFluVac cohort
#| compared to T cell immune response measured via Elispot in the RocFluVac
#| cohort."

knitr::include_graphics(here::here("Figures", "ted-roc-h1.png"))
```

### Proposed studies

First, we will ensure that our theoretical causal model is accurate, and we will
update our model to include other common causes which are observed in both
data sets. Our causal model will be formulated as a DAG, which we will use
to generate a minimal sufficient adjustment set for the causal effect of
vaccine dose on post-vaccination titer and on titer increase. We will
analyze the conditional independencies of the DAG by comparing the observed
standardized covariance matrix to the covariance matrix predicted by the DAG
in order to determine how consistent our data are with the identified DAG.

We will apply traditional regression methods using either cluster-robust
standard error estimation or hierarchical modeling to obtain a causal estimate
of the effect of dose, using a minimally sufficient adjustment set identified by
our DAG. We will use quantitative bias analysis methodology to estimate the
effect of unmeasured confounding on our estimate. Finally, we will obtain the
same estimate using a targeted maximum likelihood estimation (TMLE) framework,
which is doubly robust to model misspecification, and we can compare the TMLE
estimate to the regression estimate.

We will explore the interaction of age and dose by estimating the causal effect
of dose in both the UGAFluVac and RocFluVac studies separately and comparing
these estimates to the combined study estimate. We note that differences in
these estimates could be due to sampling variation, systematic variation between
the two source populations, or due to the interaction effect with age. We can
furthermore estimate interaction effects in all three situations.

### Expected outcomes

The outcomes of this aim will be estimates of the effect
of vaccine dose on post-vaccination titer, controlling for age. We will
estimate the amount of confounding bias in our causal effects, and compare
the estimates from the two individual studies to the overall estimate. The
individual study and confounding analyses will allow us to understand
the limitations of our causal estimates.

## Timeline

@fig-timeline shows the expected timeline for our project. The expected
completion date for our project is Friday, March 14, 2025, and a dissertation
draft will be sent to the committee no later than the following Monday. My
dissertation defense will be planned for early April to accommodate
Graduate School deadlines.

```{r}
#| label: fig-timeline
#| out-width: 6.5in
#| fig-cap: "Proposed project timeline. We expect the majority of research
#| to be conducted during the 2023-2024 academic year, extending into 2025
#| as necessary. All research will be complete by March 1, 2025 in order to
#| finalize results for the written dissertation."

knitr::include_graphics(here::here("Figures", "timeline.png"))
```


{{< pagebreak >}}
<!-- END -->