|
| 1 | +--- |
| 2 | +format: html |
| 3 | +toc: false |
| 4 | +--- |
| 5 | + |
| 6 | +```{r setup, include=FALSE} |
| 7 | +library(nimble) |
| 8 | +library(bayesTPC) |
| 9 | +library(dplyr) |
| 10 | +library(ggplot2) |
| 11 | +library(coda) |
| 12 | +library(HDInterval) |
| 13 | +library(patchwork) |
| 14 | +``` |
| 15 | + |
| 16 | +# Bluetongue Virus: Midge Fecundity Analysis |
| 17 | + |
| 18 | +**Source:** El Moustaid et al. 2021 Supplement | **Original experiment:** Lysyk & Danyk 2007 |
| 19 | + |
| 20 | +This analysis recreates the fecundity thermal performance curve from the Bluetongue virus transmission study, showing how midge reproduction varies with temperature. |
| 21 | + |
| 22 | +## Dataset |
| 23 | + |
| 24 | +We analyze 29 fecundity measurements (eggs per female per day) across 5 temperature groups: |
| 25 | + |
| 26 | +```{r fecundity-data, echo=FALSE, message=FALSE, warning=FALSE} |
| 27 | +fecundity <- tibble( |
| 28 | + T = c( |
| 29 | + rep(10, 6), # 10 C |
| 30 | + rep(15, 6), # 15 C |
| 31 | + rep(20, 6), # 20 C |
| 32 | + rep(25, 7), # 25 C |
| 33 | + rep(30, 7) # 30 C |
| 34 | + ), |
| 35 | + F = c( |
| 36 | + 5.528, 3.122, 13.11, 9.745, 6.206, 31.191, |
| 37 | + 19.034, 16.0, 1.361, 1.242, 11.08, 13.961, |
| 38 | + 17.93, 41.531, 60.535, 39.856, 51.724, 69.951, |
| 39 | + 12.731, 1.154, 36.417, 0.465, 5.844, 7.048, 19.469, |
| 40 | + 31.938, 21.195, 12.255, 22.332, 0.365, 1.703, 1.536 |
| 41 | + ) |
| 42 | +) |
| 43 | +
|
| 44 | +# Display data summary |
| 45 | +cat("Temperature groups:\n") |
| 46 | +print(fecundity %>% count(T)) |
| 47 | +cat("\nFecundity range:", round(min(fecundity$F), 1), "to", round(max(fecundity$F), 1), "eggs per female per day\n") |
| 48 | +``` |
| 49 | + |
| 50 | +## Thermal Performance Curve |
| 51 | + |
| 52 | +We fit a Brière model to capture the temperature-dependent fecundity pattern: |
| 53 | + |
| 54 | +**Model:** F(T) = k × T × (T - T_min) × √(T_max - T) |
| 55 | + |
| 56 | +```{r fit-fecundity, echo=FALSE, message=FALSE, warning=FALSE} |
| 57 | +# Prepare data for bayesTPC |
| 58 | +data_list <- list(Temp = fecundity$T, Trait = fecundity$F) |
| 59 | +
|
| 60 | +# Fit Brière model with biologically realistic priors |
| 61 | +set.seed(123) |
| 62 | +fit_fec <- b_TPC( |
| 63 | + data = data_list, |
| 64 | + model = "briere", |
| 65 | + priors = list( |
| 66 | + T_min = "dunif(5, 15)", # Lower thermal limit |
| 67 | + T_max = "dunif(32, 36)", # Upper thermal limit |
| 68 | + q = "dunif(0, 200)" # Scaling parameter |
| 69 | + ), |
| 70 | + nchains = 4, |
| 71 | + burn = 6000, |
| 72 | + niter = 18000 |
| 73 | +) |
| 74 | +``` |
| 75 | + |
| 76 | +## Results |
| 77 | + |
| 78 | +### Fecundity Thermal Performance Curve |
| 79 | + |
| 80 | +```{r plot-fecundity, echo=FALSE, message=FALSE, warning=FALSE, fig.width=8, fig.height=5} |
| 81 | +# Generate temperature grid for predictions |
| 82 | +temp_grid <- seq(0, 50, by = 0.1) |
| 83 | +
|
| 84 | +# Get posterior samples |
| 85 | +samples <- as.matrix(fit_fec$samples[[1]]) |
| 86 | +
|
| 87 | +# Generate predictions for each sample |
| 88 | +predictions <- apply(samples, 1, function(params) { |
| 89 | + T_min <- params["T_min"] |
| 90 | + T_max <- params["T_max"] |
| 91 | + q <- params["q"] |
| 92 | + |
| 93 | + # Brière function |
| 94 | + ifelse(temp_grid >= T_min & temp_grid <= T_max, |
| 95 | + q * temp_grid * (temp_grid - T_min) * sqrt(pmax(T_max - temp_grid, 0)), |
| 96 | + 0) |
| 97 | +}) |
| 98 | +
|
| 99 | +# Calculate median and HPD intervals |
| 100 | +median_pred <- apply(predictions, 1, median) |
| 101 | +hpd_lower <- apply(predictions, 1, function(x) HPDinterval(as.mcmc(x), prob = 0.95)[1]) |
| 102 | +hpd_upper <- apply(predictions, 1, function(x) HPDinterval(as.mcmc(x), prob = 0.95)[2]) |
| 103 | +
|
| 104 | +# Create single plot matching screenshot style |
| 105 | +par(mfrow = c(1, 1), mar = c(4, 4, 2, 1)) |
| 106 | +plot(temp_grid, median_pred, |
| 107 | + type = "l", |
| 108 | + lwd = 2, |
| 109 | + col = "black", |
| 110 | + xlab = "T (°C)", |
| 111 | + ylab = "Eggs per Female per Day", |
| 112 | + main = "", |
| 113 | + ylim = c(0, 80), |
| 114 | + xlim = c(0, 50), |
| 115 | + cex.lab = 1.1, |
| 116 | + cex.axis = 1.0, |
| 117 | + bty = "l") |
| 118 | +
|
| 119 | +# Add HPD interval as dashed lines |
| 120 | +lines(temp_grid, hpd_lower, lty = 2, lwd = 1.5, col = "black") |
| 121 | +lines(temp_grid, hpd_upper, lty = 2, lwd = 1.5, col = "black") |
| 122 | +
|
| 123 | +# Add data points as solid black circles |
| 124 | +points(fecundity$T, fecundity$F, pch = 16, col = "black", cex = 1.0) |
| 125 | +``` |
| 126 | + |
| 127 | +### Key Parameters |
| 128 | + |
| 129 | +```{r param-summary, echo=FALSE, message=FALSE, warning=FALSE} |
| 130 | +# MAP Estimates (Maximum A Posteriori) |
| 131 | +map_params <- MAP_estimate(fit_fec) |
| 132 | +cat("**MAP Estimates:**\n") |
| 133 | +cat("• T_min =", round(map_params["T_min"], 2), "°C\n") |
| 134 | +cat("• T_max =", round(map_params["T_max"], 2), "°C\n") |
| 135 | +cat("• q =", round(map_params["q"], 3), "\n") |
| 136 | +cat("• sigma.sq =", round(map_params["sigma.sq"], 4), "\n\n") |
| 137 | +
|
| 138 | +# Extract key parameter estimates |
| 139 | +param_summary <- summary(fit_fec) |
| 140 | +key_params <- param_summary$parameters[c("T_min","T_max","q"), ] |
| 141 | +
|
| 142 | +cat("**Posterior Medians:**\n") |
| 143 | +cat("• Lower limit (T_min):", round(as.numeric(key_params["T_min", "50%"]), 1), "°C\n") |
| 144 | +cat("• Upper limit (T_max):", round(as.numeric(key_params["T_max", "50%"]), 1), "°C\n") |
| 145 | +cat("• Scaling factor (q):", round(as.numeric(key_params["q", "50%"]), 1), "\n\n") |
| 146 | +
|
| 147 | +cat("**Optimal Temperature:**\n") |
| 148 | +# Calculate optimal temperature (where fecundity peaks) |
| 149 | +T_opt <- (as.numeric(key_params["T_min", "50%"]) + as.numeric(key_params["T_max", "50%"])) / 2 |
| 150 | +cat("• Peak fecundity at:", round(T_opt, 1), "°C\n") |
| 151 | +``` |
| 152 | + |
| 153 | +## Biological Interpretation |
| 154 | + |
| 155 | +The thermal performance curve reveals critical insights about midge reproduction: |
| 156 | + |
| 157 | +### **Thermal Limits** |
| 158 | +- **Lower limit (T_min):** Below this temperature, midges cannot reproduce |
| 159 | +- **Upper limit (T_max):** Above this temperature, reproduction ceases due to heat stress |
| 160 | + |
| 161 | +### **Optimal Range** |
| 162 | +- **Peak performance:** Maximum fecundity occurs around the midpoint of the thermal range |
| 163 | +- **Temperature sensitivity:** The curve shows how reproduction drops off rapidly at thermal extremes |
| 164 | + |
| 165 | +### **Transmission Implications** |
| 166 | +These results inform disease transmission models by quantifying how temperature affects vector population growth and disease spread potential. |
| 167 | + |
0 commit comments