Effects of vaccinations on seropositivity and symptoms prevalences
Summary
Here we first look at the relationship between vaccination and seroprevalence. Out of the 8 pathogens for which we have serological data, 5 have vaccination. None of these 5 pathogens showed a significant effect of vaccination on the seroprevalence. The analysis was done by correcting for the age of the chicken that tested for serology. Then, for the 8 pathogens, we characterized the trend in seroprevalence as a function of age. 6 did not show any significant trend as a function of age. Seroprevalence was significantly increasing with age for IBV and MG only.
In a second analysis we look at the relationship between vaccination and the prevalence of symptoms of respiratory infections, diarrhoea, CNS infections and sudden death. The prevalences of these 4 symptoms significantly decrease with age. The effect of the 6 vaccines on these four symptoms are somewhat puzzling in the sense that some vaccines do not have any effect on the symptoms that they are supposed to decrease. In some situations, they even increase them instead of decreasing them! Furthermore, other vaccines that are not supposed to have any effect on some symptoms do significantly increase and decrease them. This suggests either that the correspondance between the symptoms and the pathogens that are targeted by the vaccines is way too loose in this system or that there are other factors at play that we didn’t take into account into our analysis (e.g. AMU). The summary of the effects is detailed below:
respiratory infections: NDV and IBD significantly decrease the prevalence of symptoms of respiratory infections as we could expect (although not so much) but HPASI and IBV vaccine on the contrary significanlty increase the prevalence of symptols of respiratory infections. Furthermore, FP and PM which are not supposed to have any effect on the prevalence of respiratory infections significantly increase the prevalence of symptoms of respiratory infections.
sudden death: HPAI significantly decrease the prevalence of sudden death as we could explect but NDV and PM fail in doing so. Surprisingly, FP decreases the prevalence of sudden death whereas IBV, on the contrary, increases it.
diarrhoea: IBD decreases the prevalence of diarrhoea as expected and is the vaccine that does so the most. Surprisingly, IBV and PM also decrease the prevalence of diarrhoea but to a lesser extent than IBD. Still surprinsingly, HPAI increases the prevalence of diarrhoea.
CNS infections: NDV increases the prevalence of CNS a lot whereas it’s supposed to decrease it on the contrary. PM also increases the prevalence of CNS whereas it’s not supposed to have an effect. Surprisingly, IBV decreases the prevalence of CNS whereas it’s not supposed to have any effect on CNS.
Packages
The required packages:
> required_packages <- c("dplyr", "magrittr", "purrr", "readr", "readxl", "tibble")
> to_install <- setdiff(required_packages, rownames(installed.packages()))
> if (length(to_install)) install.packages(to_install)Laoding the packages for interactive use:
> library(magrittr)Loading and exploring the serological data
Let’s load and transform the data
> serodata <- readxl::read_excel("SeroData.xlsx") %>%
+ dplyr::mutate_at("Weight.Kg.", as.numeric) %>%
+ dplyr::mutate_at(dplyr::vars(dplyr::starts_with("Sero")), as.logical) %>%
+ dplyr::mutate_at(dplyr::vars(dplyr::ends_with("vaccinated")), as.logical) %>%
+ dplyr::rename(age = `Age(Week)`)
Warning: NAs introduits lors de la conversion automatiqueNumber of observations:
> nrow(serodata)
[1] 267Number of unique flocks:
> length(unique(serodata$FlockID))
[1] 267(i.e. one observation per flock). Number of farms:
> length(unique(serodata$`Farm ID`))
[1] 100Number of flocks per farm:
> serodata %>%
+ dplyr::group_by(`Farm ID`) %>%
+ dplyr::tally() %>%
+ purrr::pluck("n") %>%
+ table()
.
1 2 3 4 5 6 7 8 11
35 23 14 11 11 2 2 1 1 Ages:
> hist(serodata$age, col = "grey", main = NA, xlab = "age (weeks)", ylab = "number of chicken")
Weights:
> hist(serodata$Weight.Kg., col = "grey", main = NA, xlab = "weight (kg)", ylab = "number of chicken")
No correlation between age and weight:
> summary(lm(Weight.Kg. ~ age, serodata))
Call:
lm(formula = Weight.Kg. ~ age, data = serodata)
Residuals:
Min 1Q Median 3Q Max
-0.96443 -0.37720 -0.02499 0.33004 1.86397
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.12222 0.23855 8.896 <2e-16 ***
age -0.01105 0.01400 -0.789 0.431
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4943 on 262 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.002369, Adjusted R-squared: -0.001439
F-statistic: 0.6221 on 1 and 262 DF, p-value: 0.431Visually:
> with(serodata, plot(age, Weight.Kg., xlab = "age (weeks)", ylab = "weight (kg)"))
Flock sizes:
> hist(serodata$Farmsize, col = "grey", xlab = "farm size (nb of chickens)", ylab = "number of chicken")
No correlation between farm size and weight:
> summary(lm(Weight.Kg. ~ Farmsize, serodata))
Call:
lm(formula = Weight.Kg. ~ Farmsize, data = serodata)
Residuals:
Min 1Q Median 3Q Max
-0.94609 -0.39502 -0.03182 0.35437 1.86825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.960e+00 5.390e-02 36.365 <2e-16 ***
Farmsize -7.032e-05 1.274e-04 -0.552 0.582
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4946 on 262 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.001161, Adjusted R-squared: -0.002651
F-statistic: 0.3045 on 1 and 262 DF, p-value: 0.5815
> with(serodata, plot(Farmsize, Weight.Kg., xlab = "farm size (nb of chickens)", ylab = "weight (kg)"))
But age seems to increase with farm size:
> summary(lm(age ~ Farmsize, serodata))
Call:
lm(formula = age ~ Farmsize, data = serodata)
Residuals:
Min 1Q Median 3Q Max
-5.0276 -1.3600 -0.0966 1.3882 6.0929
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.572e+01 2.198e-01 71.500 < 2e-16 ***
Farmsize 3.324e-03 5.187e-04 6.408 6.72e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.024 on 265 degrees of freedom
Multiple R-squared: 0.1341, Adjusted R-squared: 0.1309
F-statistic: 41.06 on 1 and 265 DF, p-value: 6.72e-10Visually:
> mod <- loess(age ~ Farmsize, serodata)
> x_val <- with(serodata, seq(min(Farmsize), max(Farmsize), le = 512))
> pred <- predict(mod, data.frame(Farmsize = x_val), se = TRUE)
> y_val <- with(pred, c(fit + qt(.025, Inf) * se.fit, rev(fit + qt(.975, Inf) * se.fit)))
> with(serodata, plot(Farmsize, age, xlab = "farm size (nb of chickens)",
+ ylab = "age (weeks)",
+ ylim = c(min(c(age, y_val)), max(c(age, y_val)))))
> polygon(c(x_val, rev(x_val)), y_val, col = adjustcolor("blue", .1), border = NA)
> lines(x_val, pred$fit, col = "blue")
Effect of vaccinations on seroprevalence against pathogens
Making the data
The following function makes a data frame for a given disease var, with 3 columns: age, seropositive and vaccinated:
> make_data <- function(var) {
+ a <- dplyr::select(serodata, age, dplyr::contains(var))
+ sel <- grep("vaccinated", names(a))
+ if (length(sel) > 0) {
+ setNames(a, c("age", "seropositive", "vaccinated"))
+ }
+ dplyr::mutate(a, vaccinated = FALSE) %>%
+ setNames(c("age", "seropositive", "vaccinated"))
+ }The 8 pathogens:
> (vars <- sub("Sero.", "", grep("Sero", names(serodata), value = TRUE)))
[1] "AI" "CAV" "IBD" "ORT" "NDV" "IBV" "MG" "PM" Making the data sets for all the pathogens:
> datasets <- setNames(lapply(vars, make_data), vars)Checking which pathogens have vaccination:
> (vaccinated <- datasets %>%
+ sapply(function(x) length(unique(unlist(x[, "vaccinated"]))) > 1))
AI CAV IBD ORT NDV IBV MG PM
TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE Modeling the data
Methods: for the 5 pathogens for which we have both serological and vaccination data we use a logistic regression the look at the effect of age, vaccination status and interaction on the seroprevalence. Since none of these models show a significant effect of vaccination status, we then use logistic regressions to look at the effects of age on seroprevalence of the 8 pathogens. Only 2 (IBV and MG) show a significant increase of seroprevalence as a function of age with levels close to zero at hatching. Only for these two we perform extrapolation over the whole life time of the chicken. For the 6 other pathogen for which seroprevalence does not depend on age we don’t do interpolation over the age range for which we have data since we have no clue to know whether the oberved level of seroprevalence is reach after a quick increased from zero or was already as such at hatching.
The models:
> models <- lapply(datasets[vaccinated], function(x) glm(seropositive ~ age * vaccinated, binomial, x))The coefficients estimations:
> sapply(models, coef)
AI IBD NDV IBV PM
(Intercept) -1.13604430 6.8636373 1.78510672 -1.9528695 -10.691360
age 0.01330976 -0.2694324 -0.12184654 0.1580114 0.422604
vaccinatedTRUE 1.15331113 -6.1369275 -0.80875282 -1.2606290 8.624709
age:vaccinatedTRUE -0.06258336 0.4163490 0.09831679 0.1088592 -0.439410The significativities of age and vaccination on seropositivity against the 5 pathogens against which there is a vaccine:
> models %>%
+ sapply(function(x) anova(x, test = "LRT")$`Pr(>Chi)`) %>%
+ `[`(-1, ) %>%
+ `rownames<-`(c("age", "vaccinated", "age:vaccinated"))
AI IBD NDV IBV PM
age 0.5169695 0.6614363 0.6912190 0.004088135 0.04095817
vaccinated 0.7475972 0.1988441 0.2082602 0.101073655 0.19222618
age:vaccinated 0.7397933 0.2020851 0.7179156 0.507339495 0.07402805The only disease for which vaccination may have a slight effect is PM, but, not only the effect of vaccination would be small, it wouldn’t even be in the right directions…
> predict2 <- function(x) {
+ require(magrittr)
+ linkinv <- family(x)$linkinv
+ ages <- x$data$age
+ age_val <- seq(min(ages), max(ages), le = 512)
+ predict(x, data.frame(age = age_val), se.fit = TRUE) %>%
+ data.frame() %>%
+ dplyr::mutate(lower = linkinv(fit + qt(.025, Inf) * se.fit),
+ upper = linkinv(fit + qt(.975, Inf) * se.fit),
+ fit = linkinv(fit),
+ age = age_val) %>%
+ dplyr::select(-residual.scale, -se.fit)
+ }Let’s now focus on the effect of age on seroprevalence for all the pathogens:
> models <- lapply(datasets, function(x) glm(seropositive ~ age, binomial, x))Let’s compute additional data for predictions:
> linkinv <- family(models[[1]])$linkinv
> age_val1 <- seq(min(serodata$age), max(serodata$age), le = 512)
> age_val2 <- seq(0, min(serodata$age), le = 512)
> age_val3 <- range(serodata$age)
> newdata1 <- data.frame(age = age_val1)
> newdata2 <- data.frame(age = age_val2)
> age_pol1 <- c(age_val1, rev(age_val1))
> age_pol2 <- c(age_val2, rev(age_val2))
> age_pol3 <- c(age_val3, rev(age_val3))The following function uses a model as an input and computes a data frame with 512 rows and 4 columns: one for age, one for predicted seroprevalence and the last two for the lower and upper bounds of the confidence interval:
> predict2 <- function(x, newdata) {
+ require(magrittr)
+ predict(x, newdata, se.fit = TRUE) %>%
+ data.frame() %>%
+ dplyr::mutate(lower = linkinv(fit + qt(.025, Inf) * se.fit),
+ upper = linkinv(fit + qt(.975, Inf) * se.fit),
+ fit = linkinv(fit)) %>%
+ dplyr::select(-residual.scale, -se.fit)
+ }This function does the same as predict2() but in case there is no significant relationship with age:
> constant <- function(x) {
+ with(with(x,
+ binom.test(sum(seropositive), length(seropositive))),
+ c(estimate, conf.int))
+ }Let’s look at the significativities of the age effects:
> sapply(models, function(x) coef(summary(x))[2, "Pr(>|z|)"])
AI CAV IBD ORT NDV IBV MG PM
0.518343896 0.560203081 0.663876751 0.300271151 0.691107067 0.005384472 0.002566981 0.039601828 Age effect on PM after correction for multiple testing does not even remain significant:
> length(models) * sapply(models, function(x) coef(summary(x))[2, "Pr(>|z|)"])
AI CAV IBD ORT NDV IBV MG PM
4.14675117 4.48162465 5.31101401 2.40216921 5.52885654 0.04307578 0.02053585 0.31681462 Identifying the pathogens for which seroprevalence significantly increases with age:
> (ages_signif <- sapply(models, function(x) coef(summary(x))[2, "Pr(>|z|)"]) < .05 / length(models))
AI CAV IBD ORT NDV IBV MG PM
FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE Computing the models predictions, interpolation for the models where age effect is significant:
> predictions1 <- lapply(models, predict2, newdata1)Extrapolation for the models where age effect is significant:
> predictions2 <- lapply(models, predict2, newdata2)Interpolation for the models where the age effect is not significant:
> predictions3 <- lapply(datasets, constant)This function plots the model predictions in case the age effect is significant:
> plot_pred1 <- function(x, y) {
+ plot(age_val1, x$fit, xlim = c(0, 25), ylim = 0:1, type = "n",
+ xlab = "age (weeks)", ylab = "seroprevalence")
+ polygon(age_pol1, c(x$lower, rev(x$upper)), col = adjustcolor("blue", .2), border = NA)
+ lines(age_val1, x$fit, col = "blue")
+ lines(age_val1, x$lower, col = adjustcolor("blue", .2))
+ lines(age_val1, x$upper, col = adjustcolor("blue", .2))
+ polygon(age_pol2, c(y$lower, rev(y$upper)), 30, 135, NA, col = adjustcolor("blue", .2))
+ lines(age_val2, y$lower, col = adjustcolor("blue", .2))
+ lines(age_val2, y$upper, col = adjustcolor("blue", .2))
+ lines(age_val2, y$fit, col = "blue", lty = 2)
+ }where x and y are the interpolation and extrapolation predictions respectively. The following function plots the model interpolation predictions in case the age effect is not significant:
> plot_pred2 <- function(x) {
+ plot(age_val3, rep(x[1], 2), xlim = c(0, 25), ylim = 0:1,
+ type = "n", xlab = "age (weeks)", ylab = "seroprevalence")
+ polygon(age_pol3, x[c(2, 2, 3, 3)], col = adjustcolor("blue", .2), border = NA)
+ lines(age_val3, rep(x[1], 2), col = "blue")
+ lines(age_val3, rep(x[2], 2), col = adjustcolor("blue", .2))
+ lines(age_val3, rep(x[3], 2), col = adjustcolor("blue", .2))
+ }Plotting the models’ predictions:
> opar <- par(mfrow = c(2, 4), cex = .83, plt = c(.2, .85, .2, .85))
> for (i in seq_along(vars)) {
+ if (ages_signif[i]) plot_pred1(predictions1[[i]], predictions2[[i]])
+ else plot_pred2(predictions3[[i]])
+ mtext(vars[i], line = 0)
+ }
> par(opar)where the areas represent the 95% confidence interval (for extrapolation when hatched and for interpolation when full color). Extrapolations are shown only in cases where age effect is significant. Indeed, when age effect is not significant we have not much clue to know whether the level plateau corresponds to a level that has been reached after an icreased for which we don’t have data or to a level that was constant since hatching.
Effect of vaccinations on symptoms prevalences
Making the data
Loading the ViParc data for the week and the symptoms information:
> viparc <- readr::read_csv("https://raw.githubusercontent.com/viparc/clires_data/master/data/viparc.csv",
+ col_types = paste(c("ciillidddllllll", rep("d", 45)), collapse = "")) %>%
+ dplyr::mutate(FlockID = paste0(farm, flock))Laoding the data on the timings of vaccination:
> fig1 <- readr::read_csv("Figure1.csv", col_types = "cdddddddc") %>%
+ dplyr::select(-CODE) %>%
+ dplyr::mutate_if(is.numeric, as.integer) %>%
+ dplyr::arrange(WEEKNO)The following function splits the information on a given vaccine by flock:
> g <- function(vacc_name) {
+ fig1[, c("FlockID", "WEEKNO", vacc_name)] %>%
+ unique() %>%
+ split(.$FlockID)
+ }The following function recodes, for a data frame containing the vaccination information of a given flock, the vaccination information with FALSE for the weeks before vaccination and TRUE for the weeks after vaccination:
> f <- function(df) {
+ vacc_names <- grep("vaccinated", names(df))
+ out <- setNames(data.frame(df$FlockID[1], seq_len(max(df$WEEKNO)), 0L, stringsAsFactors = FALSE), names(df))
+ first <- which(df[, vacc_names] > 0)[1]
+ if (is.na(first)) return(out)
+ out[first:nrow(out), vacc_names] <- 1L
+ out
+ }Using the 2 datasets and the 2 functions defined above to generate the dataset for our analysis:
> (vacc_data <- names(fig1) %>%
+ setdiff(c("FlockID", "WEEKNO")) %>%
+ lapply(g) %>%
+ lapply(function(x) dplyr::bind_rows(lapply(x, f))) %>%
+ purrr::reduce(dplyr::left_join, c("FlockID", "WEEKNO")) %>%
+ dplyr::left_join(viparc, "FlockID") %>%
+ dplyr::select(farm, flock, week, respiratory, diarrhoea, cns, sudden_death, ends_with("vaccinated")) %>%
+ dplyr::mutate_at(dplyr::vars(dplyr::ends_with("vaccinated")), as.logical) %>%
+ tibble::as_tibble())
# A tibble: 83,062 x 13
farm flock week respiratory diarrhoea cns sudden_death Fowl.Pox.vaccin… HPAI.vaccinated NDV.vaccinated IBV.vaccinated
<chr> <int> <int> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl> <lgl>
1 75-0… 1 1 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
2 75-0… 1 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
3 75-0… 1 3 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
4 75-0… 1 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
5 75-0… 1 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
6 75-0… 1 6 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
7 75-0… 1 7 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
8 75-0… 1 8 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
9 75-0… 1 9 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
10 75-0… 1 10 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
# … with 83,052 more rows, and 2 more variables: IBD.vaccinated <lgl>, PM.vaccinated <lgl>Modeling the data
Methods: here we combine two data sets for the analysis: The first one gives the timings of vaccination and the second one give the presence of symptoms. The timing information is transformed into a boolean vector that says, for each week of each flock, whether we are before or after vaccination against a particular pathogen. Then, we use logistic regression to model, for each symptom of interest, the effect of being before or after vaccination for all the vaccines, as well age. The results show that, in absence of vaccination, the prevalences of symptoms significantly decrease with age. When accounting for this age effect, as well as correcting for multiple tests, there are a number of vaccination effects that are significant but not always in a way that is expected. This is detailed below and in the summary at the top of this document.
Estimating logistic regressions, one for each of the 4 symptoms:
> vaccine_names <- grep("vaccinated", names(vacc_data), value = TRUE)
> symptoms <- c("respiratory", "sudden_death", "diarrhoea", "cns")
> models <- setNames(lapply(
+ paste(symptoms, "~", paste(c("week", vaccine_names), collapse = "+")),
+ function(x) glm(as.formula(x), binomial, vacc_data)), symptoms)The coefficients of the 4 models:
> (coeff <- sapply(models, coef))
respiratory sudden_death diarrhoea cns
(Intercept) -3.26886509 -2.480021669 -1.59367353 -5.76842973
week -0.01034532 -0.135563230 -0.13279752 -0.05864383
Fowl.Pox.vaccinatedTRUE 0.71491124 -0.440324531 0.08814244 -0.03104626
HPAI.vaccinatedTRUE 0.17182484 -0.188275723 0.26538141 -0.18044706
NDV.vaccinatedTRUE -0.41561711 0.069989718 -0.07257318 1.58461512
IBV.vaccinatedTRUE 0.91847669 1.296963086 -0.16813184 -1.25034082
IBD.vaccinatedTRUE -0.17484649 -0.061754260 -0.34871743 0.20255575
PM.vaccinatedTRUE 0.23386012 0.009836652 -0.20486358 0.86785499The confidence intervals of the coefficients of the 4 models:
> coeff_ci <- lapply(models, confint)
Waiting for profiling to be done...
Waiting for profiling to be done...
Waiting for profiling to be done...
Waiting for profiling to be done...The significativities of the effects for the 4 models, corrected for potential confounding effects:
> models %>%
+ sapply(function(x) car::Anova(x, test = "LR")$`Pr(>Chisq)`) %>%
+ round(4) %>%
+ `rownames<-`(rownames(coeff)[-1])
respiratory sudden_death diarrhoea cns
week 7e-04 0.0000 0.0000 0.0000
Fowl.Pox.vaccinatedTRUE 0e+00 0.0000 0.0167 0.7144
HPAI.vaccinatedTRUE 0e+00 0.0000 0.0000 0.0490
NDV.vaccinatedTRUE 0e+00 0.4309 0.2668 0.0000
IBV.vaccinatedTRUE 0e+00 0.0000 0.0000 0.0000
IBD.vaccinatedTRUE 4e-04 0.2298 0.0000 0.0546
PM.vaccinatedTRUE 0e+00 0.8802 0.0000 0.0000Correcting for multiple tests:
> models %>%
+ sapply(function(x) car::Anova(x, test = "LR")$`Pr(>Chisq)`) %>%
+ round(4) %>%
+ `rownames<-`(rownames(coeff)[-1]) %>%
+ `*`(length(.))
respiratory sudden_death diarrhoea cns
week 0.0196 0.0000 0.0000 0.0000
Fowl.Pox.vaccinatedTRUE 0.0000 0.0000 0.4676 20.0032
HPAI.vaccinatedTRUE 0.0000 0.0000 0.0000 1.3720
NDV.vaccinatedTRUE 0.0000 12.0652 7.4704 0.0000
IBV.vaccinatedTRUE 0.0000 0.0000 0.0000 0.0000
IBD.vaccinatedTRUE 0.0112 6.4344 0.0000 1.5288
PM.vaccinatedTRUE 0.0000 24.6456 0.0000 0.0000The link inverse function of the models:
> linkinv <- family(models[[1]])$linkinvGenerating some new data for the models predictions:
> newdata <- data.frame(week = seq(min(vacc_data$week), max(vacc_data$week), le = 512),
+ Fowl.Pox.vaccinated = FALSE, HPAI.vaccinated = FALSE,
+ NDV.vaccinated = FALSE, IBV.vaccinated = FALSE,
+ IBD.vaccinated = FALSE, PM.vaccinated = FALSE)The following function takes a model as an input and returns a data frame with 3 columns: predicted value and the lower and upper bounds of the confidence interval of the predictions:
> predict2 <- function(x) {
+ require(magrittr)
+ predict(x, newdata, se.fit = TRUE) %>%
+ data.frame() %>%
+ tibble::as_tibble() %>%
+ dplyr::mutate(lower = linkinv(fit + qt(.025, Inf) * se.fit),
+ upper = linkinv(fit + qt(.975, Inf) * se.fit),
+ fit = linkinv(fit)) %>%
+ dplyr::select(-residual.scale, -se.fit)
+ }Generating the models’ predictions:
> predictions <- lapply(models, predict2)Reformatting the diseases names:
> leg <- names(predictions) %>%
+ sub("_", " ", .) %>%
+ sub("cns", "CNS infections", .) %>%
+ sub("y", "y infections", .)The prevalences of the diseases as functions of age:
> cols <- c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3")
> fills <- adjustcolor(cols, .2)
> x <- newdata$week
> x_val <- c(x, rev(x))
> plot(x, predictions[[1]]$fit, ylim = c(0,.17), type = "n", xlab = "age (weeks)",
+ ylab = "average prevalence")
> for (i in seq_along(predictions)) {
+ polygon(x_val, c(predictions[[i]]$lower, rev(predictions[[i]]$upper)),
+ col = fills[i], border = NA)
+ lines(x, predictions[[i]]$fit, col = cols[i])
+ }
> legend("topright", legend = leg, fill = fills, col = cols, lty = 1, bty = "n", border = NA)
The effects of the vaccines on the prevalences of the diseases:
> xs <- 1:6
> ys <- c(1 / seq(0, 10, 1), 1 + seq(0, 10, 1))
> f <- function(x) log2(exp(x))
>
> opar <- par(mfrow = c(1, 4), plt = c(.2, 1, .05, .95))
>
> for (i in seq_along(coeff_ci)) {
+ plot(xs, f(coeff[-(1:2), i]), pch = 19, xlim = c(0, 7), ylim = c(-3.1, 3.1),
+ axes = FALSE, xlab = NA, ylab = "odds ratio")
+ axis(2, log2(ys), ys)
+ text(3.5, 2.5, leg[i])
+ # text(xs, -2.5, c("NDV", "HPAI", "PM", "IBV", "IBD", "FP"))
+ text(xs, -2.5, c("FP", "HPAI", "NDV", "IBV", "IBD", "PM"))
+ arrows(xs, f(coeff_ci[[i]][-(1:2), 1]), xs, f(coeff_ci[[i]][-(1:2), 2]), .05, 90, 3)
+ abline(h = 0, lty = 2)
+ }
> par(opar)Note that all those for which the confidence interval touchs the 1 line are not significant after correction for multiple tests. All the other ones are significant.