Summary

This document contains data and code required to reproduce the analysis of immune response data presented in the manuscript:

Yeh, H.-Y., J. G. Frye, C. R. Jackson, Q. D. Read, J. E. Line, and A. Hinton. 2023. Use of automated capillary immunoassay for quantification of antibodies in chicken sera against recombinant Salmonella enterica serotype Heidelberg proteins. Journal of Microbiological Methods.

Also, additional supplementary results are presented here, not all of which appear in the manuscript.

In this document, we analyze data from the chicken immunoblot assay described in the manuscript, comparing how much IgG and IgM is produced by chicken serum exposed to the FliD and FimA proteins in an immunoblot assay. 10 non-immunized and 10 immunized chickens were each bled twice and the assay was done on each blood sample.

The statistical model is fit separately for the two response variables: IgG and IgM production. It takes into account the repeated measures within each individual chicken (2 blood samples, each of which is tested with each of the two proteins FliD and FimA, resulting in 4 samples per chicken per response variable). The response variables follow a gamma distribution with additional zero values where the immune response was below the limit of detection. Because of this, a Bayesian hurdle gamma model mixed was used. (We compared the hurdle gamma model with a hurdle lognormal model, and the hurdle gamma model was shown to be a better fit.)

We present results for each component of the model: the probability of observing a non-zero response, the magnitude of the response if there was one, and the final prediction averaging the zeros together with the non-zero values.

Abstracted results: Taking into account both the probability of a non-zero response and the size of the response, we have very strong evidence that the immunization provokes greater production of both IgG and IgM. Immunization causes high production of IgG after exposure to both FliD (~2000-fold increase averaged across bleedings) and FimA (18-fold increase), consistently across the first and second bleedings. In contrast, immunization only causes high production of IgM after exposure to FliD (1000-fold increase at first bleeding). IgM production in response to FimA is similar in immunized individuals versus non-immunized (no evidence for any increase at either bleeding). In addition, the effect of immunization on IgM production in response to FliD is weaker in the second bleeding compared to the first (120-fold increase at second bleeding, compared to ~1000-fold). Posterior parameter estimates and uncertainty intervals, as well as Bayes factors, are used to assess strength of evidence for these claims. See plots and tables for more detailed results.

Setup

Load any necessary analysis packages and read in the raw data. Use a manually created lookup table to assign treatments to columns.

library(tidyverse)

library(readxl)

library(brms)

library(emmeans)

library(tidybayes) 

library(glue)

library(gt)

library(ggnewscale)

library(bayestestR)

# Note: as of 2023-01-26 this requires you to update ggdist from github release with remotes::install_github('mjskay/ggdist')

# Otherwise the stat_interval plots do not display correctly



options(mc.cores = 4, brms.backend = 'cmdstanr', brms.file_refit = 'never')



raw <- read_xlsx('Chicken Sera FliD FimA.xlsx', skip = 4, sheet = 1)



trt_lookup <- data.frame(column = as.character(2:9), bleeding = rep(c('1st', '2nd'), each = 4), antibody = rep(c('IgG', 'IgM'), each = 2))

Clean up the data. Reshape the data to long form with a column for IgG and a column for IgM.

clean <- raw %>% 

  rename(id = ...1) %>%

  pivot_longer(-id) %>%

  separate(name, into = c('protein', 'column')) %>%

  mutate(treatment = if_else(grepl('32-', id), 'non-immunized', 'immunized') %>% factor(levels = c('non-immunized', 'immunized'))) %>%

  left_join(trt_lookup) %>%

  mutate(across(c(id, protein, bleeding, antibody), factor))



clean_wider <- clean %>%

  select(-column) %>%

  pivot_wider(names_from = antibody, values_from = value)

Exploratory plots

Plot distributions of the data by treatment, lumping the different bleedings together.

The plot shows that almost all the non-immunized show a very low immune response. Many of the immunized show an immune response orders of magnitude higher than the non-immunized individuals, but there are still a good number of immunized individuals with low or no response.

Plot the same histogram on a logarithmic scale with 1 added to all values for visualization purposes. This helps visualize differences at the lower end of the scale.

The logarithmic scale plot shows more clearly the separation between the treatments, especially for the IgG response, but still shows 4 or 5 cases from the immunized treatment where zero responses were observed for each antibody. These exploratory plots show us that we will need to fit a hurdle model, with either a gamma or a lognormal response distribution. We will try both and compare the models.

Statistical models

Fit two models for each response variable, one with hurdle-gamma response distribution and one with hurdle-lognormal. Uninformative priors are set on all parameters.

fit_igg_gamma <- brm(

  bf(IgG ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id), 

     hu ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id)),

  data = clean_wider, family = hurdle_gamma(link = 'log', link_shape = 'log', link_hu = 'logit'),

  prior = c(

    prior(normal(0, 3), class = 'b'),

    prior(normal(0, 3), class = 'b', dpar = 'hu')

  ),

  chains = 4, iter = 2000, warmup = 1000, seed = 10171,

  file = 'fit_igg_gamma'

)



fit_igm_gamma <- brm(

  bf(IgM ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id), 

     hu ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id)),

  data = clean_wider, family = hurdle_gamma(link = 'log', link_shape = 'log', link_hu = 'logit'),

  prior = c(

    prior(normal(0, 3), class = 'b'),

    prior(normal(0, 3), class = 'b', dpar = 'hu')

  ),

  chains = 4, iter = 2000, warmup = 1000, seed = 10172,

  file = 'fit_igm_gamma'

)



fit_igg_lnorm <- brm(

  bf(IgG ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id), 

     hu ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id)),

  data = clean_wider, family = hurdle_lognormal(link = 'identity', link_sigma = 'log', link_hu = 'logit'),

  prior = c(

    prior(normal(0, 3), class = 'b'),

    prior(normal(0, 3), class = 'b', dpar = 'hu')

  ),

  chains = 4, iter = 2000, warmup = 1000, seed = 10173,

  file = 'fit_igg_lnorm'

)



fit_igm_lnorm <- brm(

  bf(IgM ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id), 

     hu ~ treatment + protein + bleeding + treatment:protein + treatment:bleeding + protein:bleeding + treatment:protein:bleeding + (1|id)),

  data = clean_wider, family = hurdle_lognormal(link = 'identity', link_sigma = 'log', link_hu = 'logit'),

  prior = c(

    prior(normal(0, 3), class = 'b'),

    prior(normal(0, 3), class = 'b', dpar = 'hu')

  ),

  chains = 4, iter = 2000, warmup = 1000, seed = 10174,

  file = 'fit_igm_lnorm'

)

Ensure all models converged by looking at the maximum R-hat statistic for each model. If all are < 1.05, the models converged. All are < 1.01 indicating convergence.

max(rhat(fit_igg_gamma))
## [1] 1.001283
max(rhat(fit_igm_gamma))
## [1] 1.000847
max(rhat(fit_igg_lnorm))
## [1] 1.003332
max(rhat(fit_igm_lnorm))
## [1] 1.009373

Comparing models

To determine whether the hurdle-gamma or hurdle-lognormal is more appropriate, we first look at the posterior predictive check plot for each model. The black line is the observed data and the thin blue lines are draws from the posterior predictive distribution. We want a model where the blue lines approximate the black line (the model is a good fit to the data).

The plots show no obvious difference in the model fit quality. Both response distributions look like reasonably good fits. We can also calculate information criteria for each model using leave-one-out (LOO) cross validation. We will use the distribution that gives us the best LOO value.

fit_igg_gamma <- add_criterion(fit_igg_gamma, 'loo')

fit_igg_lnorm <- add_criterion(fit_igg_lnorm, 'loo')



fit_igm_gamma <- add_criterion(fit_igm_gamma, 'loo')

fit_igm_lnorm <- add_criterion(fit_igm_lnorm, 'loo')



loo_compare(fit_igg_gamma, fit_igg_lnorm)
##               elpd_diff se_diff

## fit_igg_gamma  0.0       0.0   

## fit_igg_lnorm -6.4       3.0
loo_compare(fit_igm_gamma, fit_igm_lnorm)
##               elpd_diff se_diff

## fit_igm_gamma  0.0       0.0   

## fit_igm_lnorm -5.1       2.2

Both comparisons show that the hurdle gamma model is a better fit, by 5 or 6 information criteria units (on the log scale). Therefore we will use the hurdle gamma model to generate predictions and effect sizes.

Use model to test effects

We now estimate the effect of the immunization treatment on provoking an immune response for both proteins, and compare the effect across proteins and between the first and second bleedings. Because a hurdle model was used, we will look at both the effect on the size of the immune response (if one was detected) and the effect on the probability of observing a zero (no immune response detected). In addition, we will look at the effect averaging the zeros together with the responses to assess the overall difference in response accounting for both probability of any response and size of the response. In all cases, we will take these predictions from the posterior distributions of the parameters, resulting in a median estimate and quantile-based credible intervals for the uncertainty.

Estimated marginal means

comparison_list = list(

  treatment = ~ treatment,

  treatment_by_protein = ~ treatment | protein,

  treatment_by_bleeding = ~ treatment | bleeding,

  three_way = ~ treatment | protein + bleeding

)



emm_igg <- emmeans(fit_igg_gamma, specs = comparison_list, type = 'response', epred = FALSE) 

emm_igg_combined <- emmeans(fit_igg_gamma, specs = comparison_list, type = 'response', epred = TRUE)

emm_igg_probzero <- emmeans(fit_igg_gamma, specs = comparison_list, type = 'response', dpar = 'hu')



emm_igm <- emmeans(fit_igm_gamma, specs = comparison_list, type = 'response', epred = FALSE)

emm_igm_combined <- emmeans(fit_igm_gamma, specs = comparison_list, type = 'response', epred = TRUE)

emm_igm_probzero <- emmeans(fit_igm_gamma, specs = comparison_list, type = 'response', dpar = 'hu')

Contrast ratios

In addition, take pairwise contrasts of all the treatment effects: overall, by protein, by bleeding, and by protein/bleeding combination. This will show us evidence for the treatment differences and interactions. Again we will get a median posterior estimate of the differences, plus quantile-based credible intervals.

Because all the analyses are done on a log scale, the contrasts are ratios. The ratios are response of immunized / response of non-immunized, so the value is the fold change in response resulting from immunization treatment. A ratio of 1 between the treatments would indicate no difference. A ratio of 10, for example, would indicate the immunization treatment results in 10 times greater immune response.

con_igg <- map(emm_igg, contrast, method = 'revpairwise')

con_igg_combined <- map(emm_igg_combined, ~ contrast(regrid(., transform = 'log'), 'revpairwise'))

con_igg_probzero <- map(emm_igg_probzero, contrast, method = 'revpairwise')



con_igm <- map(emm_igm, contrast, method = 'revpairwise')

con_igm_combined <- map(emm_igm_combined, ~ contrast(regrid(., transform = 'log'), 'revpairwise'))

con_igm_probzero <- map(emm_igm_probzero, contrast, method = 'revpairwise')

Calculate the medians and quantile-based credible intervals of the posterior mean and contrast estimates.

qi_widths <- c(.66, .90, .95)



emm_igg_draws <- map(emm_igg, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

emm_igg_combined_draws <- map(emm_igg_combined, ~ gather_emmeans_draws(.)) 

emm_igg_probzero_draws <- map(emm_igg_probzero, ~ gather_emmeans_draws(.) %>% mutate(.value = plogis(.value)))

emm_igm_draws <- map(emm_igm, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

emm_igm_combined_draws <- map(emm_igm_combined, ~ gather_emmeans_draws(.))

emm_igm_probzero_draws <- map(emm_igm_probzero, ~ gather_emmeans_draws(.) %>% mutate(.value = plogis(.value)))



emm_igg_quants <- map(emm_igg_draws, median_qi, .width = qi_widths)

emm_igg_combined_quants <- map(emm_igg_combined_draws, median_qi, .width = qi_widths)

emm_igg_probzero_quants <- map(emm_igg_probzero_draws, median_qi, .width = qi_widths)

emm_igm_quants <- map(emm_igm_draws, median_qi, .width = qi_widths)

emm_igm_combined_quants <- map(emm_igm_combined_draws, median_qi, .width = qi_widths)

emm_igm_probzero_quants <- map(emm_igm_probzero_draws, median_qi, .width = qi_widths)



con_igg_draws <- map(con_igg, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

con_igg_combined_draws <- map(con_igg_combined, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

con_igg_probzero_draws <- map(con_igg_probzero, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

con_igm_draws <- map(con_igm, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

con_igm_combined_draws <- map(con_igm_combined, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))

con_igm_probzero_draws <- map(con_igm_probzero, ~ gather_emmeans_draws(.) %>% mutate(.value = exp(.value)))



con_igg_quants <- map(con_igg_draws, median_qi, .width = qi_widths)

con_igg_combined_quants <- map(con_igg_combined_draws, median_qi, .width = qi_widths)

con_igg_probzero_quants <- map(con_igg_probzero_draws, median_qi, .width = qi_widths)

con_igm_quants <- map(con_igm_draws, median_qi, .width = qi_widths)

con_igm_combined_quants <- map(con_igm_combined_draws, median_qi, .width = qi_widths)

con_igm_probzero_quants <- map(con_igm_probzero_draws, median_qi, .width = qi_widths)

Bayes factors

We may also be interested in Bayes Factors (BF) which provide a measure of the strength of the evidence that the posterior distribution of each model parameter is different than our initial assumed “null” prior distribution. The BF is the ratio of evidence for the posterior / evidence for the prior. The greater the BF, the more evidence we have for the posterior distribution of the parameter compared to the prior null distribution. Here we compute the BF for all parameters in each of the two models.

bf_igg <- bayesfactor_parameters(fit_igg_gamma)

bf_igm <- bayesfactor_parameters(fit_igm_gamma)

Plots and tables of results

All the plots show the median posterior estimate as a black dot, with a set of shaded orange (non-immunized) and blue (immunized) bars representing the 66%, 90%, and 95% credible intervals. The tables contain the exact same information as the plots.

Plot: estimated marginal means of responses

The plot above illustrates that the immune response is almost always higher in the immunized individuals. However there are some interactions. The IgG response is elevated in immunized individuals across both proteins (FliD and FimA) and at both the first and second bleedings. However, the IgM response is most elevated in FliD in the first bleeding (elevated response is detected for FimA but is weak). By the second bleeding, the effect of immunization on IgM production is much less. The response is still present for FliD but weaker, and has completely disappeared for FimA.

The ratios (contrasts between immunized and non-immunized) are presented below in a table, along with their credible intervals.

Plot: estimated marginal means of probability of zero immune response

The previous plot compared the size of the immune response between the treatments, if there was one. The next plot looks at the probability of observing a zero response.

Plot: estimated marginal means combining size and probability of response

This plot puts together the previous two, averaging together the sizes of the responses and the 0 values to look at the overall differences.

The plot above is not that different from the first plot, except that we now see a larger effect of immunization on IgM response triggered by FliD in both bleedings 1 and 2. FimA still is not shown to produce very much response of IgM.

Plot for manuscript

This is similar to the above plot but with a different orientation, and eliminating the overall for the protein column, and with different labels.

Alternative plot for manuscript

This combines the jitter plot of raw data with the modeled predictions.

Figure legend: Observations and modeled estimates of immune response, measured as chemiluminescence in arbitrary units, in immunized and non-immunized birds. Observed data are plotted as small semi-transparent points; modeled estimates are shown as large black points (median expected values of posterior predictive distributions), with 66%, 90%, and 95% quantile credible intervals shown as shaded green bars. Separate panels are shown for all combinations of IgG and IgM response (rows), FliD and FimA protein, and time points (two and three weeks post-vaccination, as well as overall modeled estimates averaged across both time points).

Table: estimated marginal means and contrasts of size of response

protein bleeding treatment median estimate 66% QI 90% QI 95% QI
IgG
FimA 1st immunized 32,400 (23400, 45300) (18500, 57600) (16700, 65600)
FimA 1st non-immunized 3,850 (2970, 5100) (2500, 6430) (2300, 7070)
FimA 1st ratio immunized/non-immunized 8.43 (5.48, 12.7) (3.92, 17.6) (3.37, 20)
FimA 2nd immunized 134,000 (94000, 192000) (74700, 254000) (67700, 286000)
FimA 2nd non-immunized 4,180 (3120, 5650) (2580, 7120) (2370, 8050)
FimA 2nd ratio immunized/non-immunized 32.2 (20, 50.9) (13.9, 73.5) (12, 83.8)
FimA overall immunized 66,600 (50200, 88100) (40600, 107000) (37100, 118000)
FimA overall non-immunized 4,020 (3240, 5050) (2770, 6130) (2580, 6750)
FimA overall ratio immunized/non-immunized 16.5 (11.3, 23.6) (8.45, 30.3) (7.29, 34.2)
FliD 1st immunized 748,000 (566000, 989000) (471000, 1220000) (432000, 1350000)
FliD 1st non-immunized 811 (595, 1130) (478, 1460) (435, 1650)
FliD 1st ratio immunized/non-immunized 909 (599, 1410) (435, 1910) (371, 2130)
FliD 2nd immunized 657,000 (503000, 877000) (423000, 1100000) (387000, 1220000)
FliD 2nd non-immunized 386 (47.5, 3210) (10.7, 14700) (5.45, 28400)
FliD 2nd ratio immunized/non-immunized 1,680 (200, 13800) (46.4, 62900) (24.8, 118000)
FliD overall immunized 706,000 (568000, 877000) (489000, 1030000) (451000, 1100000)
FliD overall non-immunized 575 (194, 1700) (88.4, 3610) (63.5, 5120)
FliD overall ratio immunized/non-immunized 1,230 (403, 3770) (185, 8460) (134, 11700)
overall 1st immunized 156,000 (124000, 198000) (103000, 235000) (94000, 254000)
overall 1st non-immunized 1,790 (1430, 2250) (1220, 2680) (1130, 2900)
overall 1st ratio immunized/non-immunized 87.5 (63.3, 122) (48.7, 153) (42.3, 170)
overall 2nd immunized 298,000 (236000, 382000) (198000, 462000) (181000, 505000)
overall 2nd non-immunized 1,300 (442, 3890) (200, 8470) (143, 11500)
overall 2nd ratio immunized/non-immunized 231 (74.2, 708) (35.2, 1570) (24.2, 2190)
overall overall immunized 218,000 (177000, 264000) (152000, 301000) (139000, 321000)
overall overall non-immunized 1,530 (862, 2670) (571, 4010) (472, 4900)
overall overall ratio immunized/non-immunized 143 (77.7, 260) (51, 402) (40.8, 483)
IgM
FimA 1st immunized 15,700 (12100, 20500) (10000, 24900) (8990, 27100)
FimA 1st non-immunized 8,410 (6550, 10900) (5460, 13200) (4970, 14300)
FimA 1st ratio immunized/non-immunized 1.86 (1.31, 2.66) (1.02, 3.54) (0.89, 4.08)
FimA 2nd immunized 17,400 (13200, 23200) (10800, 28800) (9820, 31800)
FimA 2nd non-immunized 16,500 (12500, 22100) (10100, 27200) (9200, 30400)
FimA 2nd ratio immunized/non-immunized 1.06 (0.709, 1.56) (0.533, 2.12) (0.459, 2.48)
FimA overall immunized 16,600 (13200, 21000) (11000, 25200) (10100, 27100)
FimA overall non-immunized 11,800 (9360, 15100) (7850, 17800) (7280, 19400)
FimA overall ratio immunized/non-immunized 1.40 (1.01, 1.94) (0.792, 2.51) (0.702, 2.81)
FliD 1st immunized 155,000 (122000, 2e+05) (101000, 240000) (92100, 258000)
FliD 1st non-immunized 176 (135, 235) (113, 294) (102, 331)
FliD 1st ratio immunized/non-immunized 877 (603, 1250) (457, 1650) (393, 1910)
FliD 2nd immunized 96,400 (74200, 126000) (60300, 153000) (55400, 167000)
FliD 2nd non-immunized 14,200 (8190, 24100) (5600, 36600) (4620, 46100)
FliD 2nd ratio immunized/non-immunized 6.82 (3.72, 12.4) (2.37, 19.2) (1.9, 22.6)
FliD overall immunized 123,000 (98500, 153000) (82800, 178000) (75300, 194000)
FliD overall non-immunized 1,580 (1140, 2210) (908, 2850) (804, 3210)
FliD overall ratio immunized/non-immunized 77.2 (51.9, 114) (37.9, 152) (32.5, 170)
overall 1st immunized 49,500 (39600, 61900) (33500, 72600) (30500, 78500)
overall 1st non-immunized 1,220 (977, 1540) (824, 1860) (756, 2020)
overall 1st ratio immunized/non-immunized 40.5 (29.4, 55.3) (23.4, 69.9) (20.5, 80)
overall 2nd immunized 40,900 (32700, 51800) (27200, 62000) (24700, 67500)
overall 2nd non-immunized 15,400 (10900, 21600) (8570, 28100) (7630, 31200)
overall 2nd ratio immunized/non-immunized 2.68 (1.76, 4.07) (1.28, 5.41) (1.12, 6.3)
overall overall immunized 45,100 (36900, 55600) (31200, 64400) (28500, 69800)
overall overall non-immunized 4,320 (3420, 5530) (2850, 6610) (2610, 7270)
overall overall ratio immunized/non-immunized 10.4 (7.57, 14.4) (5.92, 18.2) (5.34, 20.3)

Table: estimated means and contrasts of probability of zero immune response

protein bleeding treatment median estimate 66% QI 90% QI 95% QI
IgG
FimA 1st immunized 0.0570 (0.0173, 0.142) (0.00611, 0.233) (0.00379, 0.293)
FimA 1st non-immunized 0.0224 (0.00676, 0.0653) (0.00247, 0.127) (0.00159, 0.16)
FimA 1st ratio immunized/non-immunized 2.44 (0.677, 9.27) (0.238, 23.1) (0.153, 37.3)
FimA 2nd immunized 0.241 (0.123, 0.399) (0.0675, 0.527) (0.0485, 0.576)
FimA 2nd non-immunized 0.133 (0.0568, 0.251) (0.0289, 0.358) (0.0195, 0.411)
FimA 2nd ratio immunized/non-immunized 2.08 (0.708, 6.51) (0.329, 15.3) (0.222, 22.6)
FimA overall immunized 0.118 (0.0583, 0.208) (0.0312, 0.3) (0.023, 0.34)
FimA overall non-immunized 0.0549 (0.024, 0.112) (0.0116, 0.178) (0.00783, 0.214)
FimA overall ratio immunized/non-immunized 2.26 (0.873, 6.2) (0.4, 12.9) (0.293, 19.1)
FliD 1st immunized 0.00558 (0.00063, 0.0342) (0.00012, 0.0788) (5.69e-05, 0.114)
FliD 1st non-immunized 0.178 (0.0807, 0.312) (0.0427, 0.428) (0.0302, 0.486)
FliD 1st ratio immunized/non-immunized 0.0267 (0.00317, 0.173) (0.000656, 0.507) (0.000305, 0.76)
FliD 2nd immunized 0.0209 (0.00342, 0.0794) (0.000547, 0.157) (0.000194, 0.203)
FliD 2nd non-immunized 0.966 (0.894, 0.992) (0.812, 0.997) (0.761, 0.999)
FliD 2nd ratio immunized/non-immunized 0.000619 (6.83e-05, 0.0037) (1.13e-05, 0.0117) (4.61e-06, 0.0203)
FliD overall immunized 0.0104 (0.00216, 0.0363) (0.000556, 0.0696) (0.000279, 0.0935)
FliD overall non-immunized 0.708 (0.528, 0.857) (0.399, 0.92) (0.346, 0.943)
FliD overall ratio immunized/non-immunized 0.00394 (0.000714, 0.0173) (0.000167, 0.0413) (9.06e-05, 0.0608)
overall 1st immunized 0.0172 (0.00462, 0.0522) (0.00152, 0.0984) (0.000844, 0.122)
overall 1st non-immunized 0.0642 (0.0294, 0.127) (0.0148, 0.191) (0.0106, 0.226)
overall 1st ratio immunized/non-immunized 0.260 (0.0596, 0.941) (0.0193, 2.3) (0.0106, 3.73)
overall 2nd immunized 0.0732 (0.027, 0.158) (0.0103, 0.237) (0.00656, 0.278)
overall 2nd non-immunized 0.674 (0.484, 0.833) (0.358, 0.905) (0.309, 0.931)
overall 2nd ratio immunized/non-immunized 0.0359 (0.00968, 0.113) (0.00321, 0.237) (0.00199, 0.355)
overall overall immunized 0.0347 (0.0143, 0.0739) (0.00637, 0.111) (0.00413, 0.134)
overall overall non-immunized 0.276 (0.168, 0.412) (0.111, 0.522) (0.0911, 0.576)
overall overall ratio immunized/non-immunized 0.0945 (0.0323, 0.244) (0.0135, 0.461) (0.00939, 0.651)
IgM
FimA 1st immunized 0.0234 (0.00491, 0.0817) (0.00143, 0.165) (0.000732, 0.217)
FimA 1st non-immunized 0.0233 (0.00635, 0.0675) (0.0022, 0.136) (0.00133, 0.178)
FimA 1st ratio immunized/non-immunized 1.04 (0.211, 4.49) (0.0606, 12.7) (0.0329, 20.4)
FimA 2nd immunized 0.226 (0.0879, 0.432) (0.038, 0.602) (0.0229, 0.681)
FimA 2nd non-immunized 0.273 (0.125, 0.472) (0.0645, 0.63) (0.0471, 0.704)
FimA 2nd ratio immunized/non-immunized 0.776 (0.199, 2.84) (0.065, 7.52) (0.0379, 11.7)
FimA overall immunized 0.0776 (0.0276, 0.171) (0.0105, 0.274) (0.00629, 0.326)
FimA overall non-immunized 0.0852 (0.0349, 0.178) (0.0166, 0.285) (0.0107, 0.34)
FimA overall ratio immunized/non-immunized 0.905 (0.252, 2.83) (0.0921, 6.65) (0.054, 10.6)
FliD 1st immunized 0.00314 (0.000344, 0.0207) (5.8e-05, 0.0649) (2.53e-05, 0.0898)
FliD 1st non-immunized 0.120 (0.0444, 0.267) (0.0197, 0.424) (0.0133, 0.505)
FliD 1st ratio immunized/non-immunized 0.0227 (0.00259, 0.155) (0.000441, 0.62) (0.000178, 1.09)
FliD 2nd immunized 0.0511 (0.0121, 0.151) (0.00318, 0.278) (0.00165, 0.355)
FliD 2nd non-immunized 0.946 (0.849, 0.986) (0.735, 0.996) (0.659, 0.997)
FliD 2nd ratio immunized/non-immunized 0.00276 (0.000339, 0.0158) (5.8e-05, 0.0472) (2.45e-05, 0.0763)
FliD overall immunized 0.0127 (0.00273, 0.0441) (0.000683, 0.0913) (0.00035, 0.131)
FliD overall non-immunized 0.610 (0.401, 0.796) (0.266, 0.893) (0.214, 0.924)
FliD overall ratio immunized/non-immunized 0.00809 (0.00127, 0.0379) (0.000254, 0.106) (0.000121, 0.165)
overall 1st immunized 0.00864 (0.00167, 0.0323) (0.00043, 0.0736) (0.000225, 0.0962)
overall 1st non-immunized 0.0543 (0.0198, 0.12) (0.00936, 0.204) (0.00583, 0.262)
overall 1st ratio immunized/non-immunized 0.155 (0.0294, 0.691) (0.00726, 1.93) (0.00346, 2.92)
overall 2nd immunized 0.111 (0.0409, 0.231) (0.0158, 0.357) (0.00965, 0.414)
overall 2nd non-immunized 0.720 (0.532, 0.866) (0.389, 0.934) (0.325, 0.952)
overall 2nd ratio immunized/non-immunized 0.0471 (0.0106, 0.166) (0.00303, 0.366) (0.00161, 0.545)
overall overall immunized 0.0317 (0.0102, 0.0768) (0.00336, 0.13) (0.00204, 0.165)
overall overall non-immunized 0.277 (0.152, 0.435) (0.093, 0.587) (0.0715, 0.654)
overall overall ratio immunized/non-immunized 0.0862 (0.0222, 0.262) (0.00665, 0.568) (0.00351, 0.813)

Table: estimated marginal means and contrasts of combined response

As discussed above, this combines the two components above (size of response and probability of a zero response) to a single measure of response.

protein bleeding treatment median estimate 66% QI 90% QI 95% QI
IgG
FimA 1st immunized 29,700 (21100, 42100) (16600, 54100) (15100, 60600)
FimA 1st non-immunized 3,700 (2840, 4920) (2390, 6210) (2180, 6920)
FimA 1st ratio immunized/non-immunized 8.07 (5.19, 12.2) (3.65, 16.9) (3.12, 19.4)
FimA 2nd immunized 97,600 (64700, 147000) (48300, 199000) (42100, 223000)
FimA 2nd non-immunized 3,480 (2560, 4890) (2050, 6170) (1830, 6950)
FimA 2nd ratio immunized/non-immunized 27.6 (16.4, 46.6) (10.9, 68.9) (9.2, 81)
FimA overall immunized 64,600 (46000, 91200) (36400, 119000) (32700, 131000)
FimA overall non-immunized 3,700 (2930, 4700) (2480, 5710) (2260, 6260)
FimA overall ratio immunized/non-immunized 17.5 (11.5, 26.6) (8.32, 36.4) (7.16, 42.2)
FliD 1st immunized 734,000 (554000, 967000) (461000, 1210000) (422000, 1340000)
FliD 1st non-immunized 648 (458, 931) (350, 1230) (304, 1380)
FliD 1st ratio immunized/non-immunized 1,110 (721, 1780) (512, 2490) (439, 2900)
FliD 2nd immunized 631,000 (480000, 847000) (395000, 1040000) (362000, 1190000)
FliD 2nd non-immunized 12.5 (1.02, 128) (0.153, 670) (0.0734, 1440)
FliD 2nd ratio immunized/non-immunized 51,400 (4890, 605000) (932, 4150000) (444, 8610000)
FliD overall immunized 698,000 (557000, 868000) (475000, 1020000) (438000, 1110000)
FliD overall non-immunized 354 (246, 543) (184, 839) (163, 1150)
FliD overall ratio immunized/non-immunized 1,960 (1200, 3050) (749, 4140) (549, 4920)
overall 1st immunized 383,000 (292000, 5e+05) (242000, 621000) (223000, 685000)
overall 1st non-immunized 2,210 (1740, 2860) (1460, 3500) (1370, 3850)
overall 1st ratio immunized/non-immunized 172 (121, 245) (91.1, 321) (80.7, 368)
overall 2nd immunized 369,000 (287000, 481000) (238000, 588000) (222000, 652000)
overall 2nd non-immunized 1,800 (1310, 2550) (1040, 3320) (941, 3860)
overall 2nd ratio immunized/non-immunized 203 (132, 312) (96.2, 417) (81.1, 493)
overall overall immunized 383,000 (308000, 472000) (267000, 551000) (248000, 595000)
overall overall non-immunized 2,050 (1630, 2590) (1410, 3170) (1290, 3520)
overall overall ratio immunized/non-immunized 187 (135, 255) (106, 318) (95.2, 354)
IgM
FimA 1st immunized 14,900 (11500, 19600) (9340, 24000) (8400, 26200)
FimA 1st non-immunized 8,100 (6260, 10400) (5140, 12700) (4670, 13900)
FimA 1st ratio immunized/non-immunized 1.84 (1.29, 2.65) (0.99, 3.54) (0.869, 4.08)
FimA 2nd immunized 12,800 (8530, 18200) (5990, 23200) (4960, 25800)
FimA 2nd non-immunized 11,600 (7550, 16400) (5220, 21300) (4380, 23900)
FimA 2nd ratio immunized/non-immunized 1.12 (0.648, 1.91) (0.428, 2.94) (0.358, 3.69)
FimA overall immunized 14,100 (10700, 18300) (8740, 22100) (7970, 24200)
FimA overall non-immunized 9,830 (7360, 13000) (5870, 16000) (5290, 17700)
FimA overall ratio immunized/non-immunized 1.42 (0.971, 2.11) (0.725, 2.89) (0.627, 3.29)
FliD 1st immunized 153,000 (120000, 198000) (99400, 237000) (89900, 256000)
FliD 1st non-immunized 149 (108, 204) (82.8, 255) (71, 287)
FliD 1st ratio immunized/non-immunized 1,030 (696, 1560) (515, 2170) (449, 2540)
FliD 2nd immunized 88,700 (66200, 116000) (52400, 143000) (46800, 158000)
FliD 2nd non-immunized 720 (179, 2360) (56.4, 4820) (29.9, 6750)
FliD 2nd ratio immunized/non-immunized 121 (36.4, 506) (17.5, 1750) (12.6, 3220)
FliD overall immunized 122,000 (97300, 153000) (81200, 181000) (73700, 195000)
FliD overall non-immunized 440 (165, 1260) (91.2, 2490) (74.1, 3450)
FliD overall ratio immunized/non-immunized 271 (94.5, 760) (47.5, 1400) (33.9, 1780)
overall 1st immunized 84,100 (66600, 108000) (55800, 129000) (50500, 139000)
overall 1st non-immunized 4,130 (3200, 5320) (2630, 6460) (2400, 7050)
overall 1st ratio immunized/non-immunized 20.5 (14.4, 29) (11.2, 38.1) (9.85, 43)
overall 2nd immunized 51,000 (38400, 66000) (31200, 80500) (28000, 88500)
overall 2nd non-immunized 6,410 (4180, 9200) (2850, 11900) (2290, 13700)
overall 2nd ratio immunized/non-immunized 8.00 (4.99, 13.2) (3.53, 20.6) (3.01, 25.3)
overall overall immunized 68,200 (54700, 85400) (45900, 99900) (41600, 107000)
overall overall non-immunized 5,300 (3900, 7080) (3060, 8660) (2760, 9640)
overall overall ratio immunized/non-immunized 12.9 (8.85, 18.8) (6.78, 25.4) (6.01, 29)

Table: Bayes Factors

The BFs for each parameter in the model for each antibody are presented below in a table.

response component parameter IgG IgM
size of response treatment 159 0.526
size of response protein 38.2 1.82 × 1012
size of response bleeding 0.112 1.9
size of response treatment:protein 1.51 × 106 2.35 × 1011
size of response treatment:bleeding 3.8 0.474
size of response protein:bleeding 0.841 6.13 × 104
size of response treatment:protein:bleeding 0.74 4.02 × 104
probability of zero response treatment 0.546 0.493
probability of zero response protein 2.23 1.29
probability of zero response bleeding 1.58 6.91
probability of zero response treatment:protein 17.9 7.43
probability of zero response treatment:bleeding 0.488 0.582
probability of zero response protein:bleeding 3.51 1.36
probability of zero response treatment:protein:bleeding 2.6 1.09

We can take the following key results from this table:

  • There is strong evidence for an overall effect of immunization on IgG but due to the interaction effects for IgM there is not evidence for an overall effect.
  • There is somewhat moderate evidence for an overall difference between the proteins for IgG, but extremely strong evidence for a difference between the proteins for IgM. There is also extremely strong evidence for an immunization by protein interaction for both Ig responses.
  • There is no or very weak evidence for any interactions with bleeding for IgG, but there is for IgM. We see a protein by bleeding interaction, as well as a three-way interaction (because the response was biggest for immunized, FliD, first bleeding).
  • There is much less evidence overall for effects on the probability of zero component of the response. We do see moderately strong evidence that there is an immunization effect that differs by protein for both IgG and IgM.