Fall armyworm analysis
Quentin D. Read
November 13, 2024
This notebook contains all statistical analysis presented in the manuscript:
Woolfolk, S., G. Matthews, and Q. D. Read. 2024. Comparison of infestation rates of fall armyworm (Lepidoptera: Noctuidae) neonates for maize resistance screening. Journal of Insect Science.
Summary
This is an analysis of the data from the FAW leaf feeding damage experiments, 2023 and 2024. We correct the issues with the statistics identified by the editor from the original submission of the manuscript. First, we now correctly treat the response variable (leaf damage rating) as an ordered categorical response. We include a random effect to account for the split plot design, and include three-way interactions between fixed effects of genotype, treatment, and year. We present means for each genotype and treatment based on the underlying modeled probabilities of each category (weighted averages of the modeled probabilities).
This document contains all code to produce quantitative results from the tables where means comparisons for treatment, genotype, and treatment by genotype interaction are presented, both within each year and averaged across years. Code to produce the figure is also included.
Setup
Load packages and import data from CSV.
library(tidyverse)
library(ordinal)
library(emmeans)
library(multcomp)
genotype_IDs <- c('Ab24E', 'Tx601', 'T173', 'Va35', 'Mp496', 'Mp704', 'Mp705', 'Mp706', 'Mp707', 'Mp708', 'Mp714')
faw <- read_csv('FAW_combinedyears.csv') %>%
mutate(across(c(Year, Row, Rep, Entry, Treatment, Plant), factor),
Pedigree = factor(Pedigree, levels = genotype_IDs),
Rating = ordered(Rating))Fit statistical model
We need to fit a model with an ordered multinomial response variable.
First we will do a quick check to diagnose the proportional odds assumption by genotype and by treatment. Results are not shown. The assumption is not grossly violated so we can proceed.
logodds_successive <- function(x) {
logits <- sapply(levels(x), function(i) qlogis(mean(x >= i)))
diffs <- diff(logits)
tibble(level = names(diffs), `diff(level-1)` = diffs)
}
faw %>%
group_by(Year, Treatment) %>%
reframe(logodds_successive(Rating)) %>%
pivot_wider(names_from=Treatment, values_from=`diff(level-1)`)
faw %>%
group_by(Year, Pedigree) %>%
reframe(logodds_successive(Rating)) %>%
pivot_wider(names_from=Pedigree, values_from=`diff(level-1)`)Next, let’s try to fit the model as a cumulative logistic mixed model. We will treat rep nested within year as a random effect. Fixed effects will be treatment, pedigree, year, and their interactions.
clmm_fit <- clmm(Rating ~ Pedigree * Treatment * Year + (1|Rep:Year), link = 'logit', threshold = 'flexible', data = faw)Results
Analysis of deviance tables
First, we do joint chi-squared tests for each model term to get an analysis of deviance table. Because this is a GLMM, we are now doing chi-squared tests instead of F-tests. In the original manuscript, this is done overall, broken down only by year (Tables 1 and 2), then by genotype and by treatment within each year (Tables 3 and 4). It is now also done by genotype and by treatment averaged across years.
Note: In some of these tables we see a message that the degrees of freedom of the Chi-squared distribution have been reduced due to non-estimability, and there is a row for confounded variance. This is related to not every genotype being present in every year.
Overall analysis of deviance table
This table has the Chi-squared values that correspond to the “combined-years” F-tests for genotype and treatment.
joint_tests(clmm_fit)## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 356.458 2851.664 <.0001 e
## Pedigree:Treatment 32 Inf 5.470 175.040 <.0001 e
## Pedigree:Year 8 Inf 29.267 234.136 <.0001 e
## Pedigree:Treatment:Year 32 Inf 7.625 244.000 <.0001 e
## (confounded) 19 Inf 14.202 3244.750 <.0001
##
## e: df1 reduced due to non-estimabilityAnalysis of deviance tables by year
This table has the Chi-squared tests that correspond to the 2023 and 2024 F-tests for genotype and treatment.
joint_tests(clmm_fit, by = 'Year')## Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 261.503 2353.527 <.0001 e
## Pedigree:Treatment 36 Inf 5.376 193.536 <.0001 e
## (confounded) 4 Inf 32.301 2546.856 <.0001
##
## Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 181.210 1630.890 <.0001 e
## Pedigree:Treatment 36 Inf 6.839 246.204 <.0001 e
## (confounded) 4 Inf 50.123 1943.318 <.0001
##
## e: df1 reduced due to non-estimabilityAnalysis of deviance for treatment by genotype combined across years
Note that neither Mp706 nor Mp707 are included here because they were not present in the data for both years.
joint_tests(clmm_fit, by = 'Pedigree', at = list(Pedigree = genotype_IDs[!genotype_IDs %in% c('Mp706', 'Mp707')]))## Pedigree = Ab24E:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 7.357 29.428 <.0001
## Year 1 Inf 33.528 33.528 <.0001
## Treatment:Year 4 Inf 1.979 7.916 0.0947
##
## Pedigree = Tx601:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 23.545 94.180 <.0001
## Year 1 Inf 59.706 59.706 <.0001
## Treatment:Year 4 Inf 5.627 22.508 0.0002
##
## Pedigree = T173:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 36.268 145.072 <.0001
## Year 1 Inf 10.909 10.909 0.0010
## Treatment:Year 4 Inf 4.894 19.576 0.0006
##
## Pedigree = Va35:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 41.186 164.744 <.0001
## Year 1 Inf 29.113 29.113 <.0001
## Treatment:Year 4 Inf 2.421 9.684 0.0461
##
## Pedigree = Mp496:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 25.332 101.328 <.0001
## Year 1 Inf 1.917 1.917 0.1662
## Treatment:Year 4 Inf 31.862 127.448 <.0001
##
## Pedigree = Mp704:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 29.387 117.548 <.0001
## Year 1 Inf 1.458 1.458 0.2273
## Treatment:Year 4 Inf 11.661 46.644 <.0001
##
## Pedigree = Mp705:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 24.054 96.216 <.0001
## Year 1 Inf 6.519 6.519 0.0107
## Treatment:Year 4 Inf 2.150 8.600 0.0719
##
## Pedigree = Mp708:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 24.987 99.948 <.0001
## Year 1 Inf 0.014 0.014 0.9046
## Treatment:Year 4 Inf 2.760 11.040 0.0261
##
## Pedigree = Mp714:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 56.011 224.044 <.0001
## Year 1 Inf 0.012 0.012 0.9112
## Treatment:Year 4 Inf 3.991 15.964 0.0031Analysis of deviance for treatment by genotype within each year
This table has the Chi-squared tests that correspond to the F-tests for treatment within each genotype within each year. The genotype that was not present in each of the two years is not included in these tests.
2023
joint_tests(clmm_fit, by = 'Pedigree', at = list(Pedigree = genotype_IDs[!genotype_IDs %in% 'Mp706'], Year = '2023'))## Pedigree = Ab24E:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 1.478 5.912 0.2059
##
## Pedigree = Tx601:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 20.231 80.924 <.0001
##
## Pedigree = T173:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 23.717 94.868 <.0001
##
## Pedigree = Va35:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 20.605 82.420 <.0001
##
## Pedigree = Mp496:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 28.653 114.612 <.0001
##
## Pedigree = Mp704:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 19.423 77.692 <.0001
##
## Pedigree = Mp705:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 14.129 56.516 <.0001
##
## Pedigree = Mp707:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 7.963 31.852 <.0001
##
## Pedigree = Mp708:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 9.169 36.676 <.0001
##
## Pedigree = Mp714:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 20.611 82.444 <.00012024
joint_tests(clmm_fit, by = 'Pedigree', at = list(Pedigree = genotype_IDs[!genotype_IDs %in% 'Mp707'], Year = '2024'))## Pedigree = Ab24E:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 9.305 37.220 <.0001
##
## Pedigree = Tx601:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 4.862 19.448 0.0006
##
## Pedigree = T173:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 16.694 66.776 <.0001
##
## Pedigree = Va35:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 23.293 93.172 <.0001
##
## Pedigree = Mp496:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 30.408 121.632 <.0001
##
## Pedigree = Mp704:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 25.455 101.820 <.0001
##
## Pedigree = Mp705:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 12.586 50.344 <.0001
##
## Pedigree = Mp706:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 9.692 38.768 <.0001
##
## Pedigree = Mp708:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 20.334 81.336 <.0001
##
## Pedigree = Mp714:
## model term df1 df2 F.ratio Chisq p.value
## Treatment 4 Inf 41.659 166.636 <.0001Analysis of deviance for genotype by treatment combined across years
joint_tests(clmm_fit, by = 'Treatment')## Treatment = 10:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 144.838 1158.704 <.0001 e
## Pedigree:Year 8 Inf 6.389 51.112 <.0001 e
## (confounded) 3 Inf 17.954 1224.846 <.0001
##
## Treatment = 20:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 117.662 941.296 <.0001 e
## Pedigree:Year 8 Inf 13.637 109.096 <.0001 e
## (confounded) 3 Inf 14.518 1073.935 <.0001
##
## Treatment = 30:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 126.320 1010.560 <.0001 e
## Pedigree:Year 8 Inf 19.253 154.024 <.0001 e
## (confounded) 3 Inf 21.979 1131.904 <.0001
##
## Treatment = 40:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 129.746 1037.968 <.0001 e
## Pedigree:Year 8 Inf 8.041 64.328 <.0001 e
## (confounded) 3 Inf 17.244 1106.093 <.0001
##
## Treatment = 50:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 8 Inf 116.588 932.704 <.0001 e
## Pedigree:Year 8 Inf 16.924 135.392 <.0001 e
## (confounded) 3 Inf 31.346 1061.299 <.0001
##
## e: df1 reduced due to non-estimabilityAnalysis of deviance for genotype by treatment within each year
This table has the Chi-squared tests that correspond to the F-tests for genotype within each combination of treatment and year.
joint_tests(clmm_fit, by = c('Treatment', 'Year'))## Treatment = 10, Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 89.734 807.606 <.0001 e
##
## Treatment = 20, Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 77.513 697.617 <.0001 e
##
## Treatment = 30, Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 88.952 800.568 <.0001 e
##
## Treatment = 40, Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 84.223 758.007 <.0001 e
##
## Treatment = 50, Year = 2023:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 87.164 784.476 <.0001 e
##
## Treatment = 10, Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 64.344 579.096 <.0001 e
##
## Treatment = 20, Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 54.937 494.433 <.0001 e
##
## Treatment = 30, Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 50.954 458.586 <.0001 e
##
## Treatment = 40, Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 52.679 474.111 <.0001 e
##
## Treatment = 50, Year = 2024:
## model term df1 df2 F.ratio Chisq p.value note
## Pedigree 9 Inf 44.096 396.864 <.0001 e
##
## e: df1 reduced due to non-estimabilityComparison of means
We do the following sets of multiple comparisons:
- Compare each level of genotype, averaged over treatment levels, within 2023, within 2024, and combined across years
- Compare each level of treatment, averaged over genotype levels, within 2023, within 2024, and combined across years
- Compare each level of genotype within each treatment level within 2023
- Compare each level of treatment within each genotype level within 2023
- Compare each level of genotype within each treatment level within 2024
- Compare each level of treatment within each genotype level within 2024
- Compare each level of genotype within each treatment level combined across years
- Compare each level of treatment within each genotype level combined across years
The means presented here are actually weighted averages of the modeled probabilities. Of course, the results should be interpreted with the caution that these means are derived from a categorical rating scale. We are using the Tukey test to adjust the p-values accounting for multiple hypothesis testing. The means are presented with the 95% confidence intervals adjusted for multiple comparisons (Sidak adjustment), in brackets after each mean.
Note that the genotype Mp707 is only included in the 2023 comparison, Mp706 is only included in the 2024 comparison, and neither of them are included in the combined years comparison. This also applies to the overall treatment means separation: the genotypes that aren’t present in both years can’t be included in the overall treatment mean that averages over both years.
emm_g_by_year <- emmeans(clmm_fit, ~ Pedigree | Year, mode = 'mean.class')
emm_g_overall <- emmeans(clmm_fit, ~ Pedigree, mode = 'mean.class')
cld_g_by_year <- cld(emm_g_by_year, adjust = 'tukey', Letters = letters, decreasing = TRUE)
cld_g_overall <- cld(emm_g_overall, adjust = 'tukey', Letters = letters, decreasing = TRUE)
emm_t_2023 <- emmeans(clmm_fit, ~ Treatment, mode = 'mean.class', at = list(Year = '2023', Pedigree = genotype_IDs[!genotype_IDs %in% c('Mp706')]))
emm_t_2024 <- emmeans(clmm_fit, ~ Treatment, mode = 'mean.class', at = list(Year = '2024', Pedigree = genotype_IDs[!genotype_IDs %in% c('Mp707')]))
emm_t_overall <- emmeans(clmm_fit, ~ Treatment, mode = 'mean.class', at = list(Pedigree = genotype_IDs[!genotype_IDs %in% c('Mp706', 'Mp707')]))
cld_t_2023 <- cld(emm_t_2023, adjust = 'tukey', Letters = letters, decreasing = TRUE)
cld_t_2024 <- cld(emm_t_2024, adjust = 'tukey', Letters = letters, decreasing = TRUE)
cld_t_overall <- cld(emm_t_overall, adjust = 'tukey', Letters = letters, decreasing = TRUE)
emm_t_by_gxyear <- emmeans(clmm_fit, ~ Treatment | Pedigree + Year, mode = 'mean.class')
emm_g_by_txyear <- emmeans(clmm_fit, ~ Pedigree | Treatment + Year, mode = 'mean.class')
cld_t_by_gxyear <- cld(emm_t_by_gxyear, adjust = 'tukey', Letters = letters, decreasing = TRUE)
cld_g_by_txyear <- cld(emm_g_by_txyear, adjust = 'tukey', Letters = LETTERS, decreasing = TRUE)
emm_t_by_g_overall <- emmeans(clmm_fit, ~ Treatment | Pedigree, mode = 'mean.class')
emm_g_by_t_overall <- emmeans(clmm_fit, ~ Pedigree | Treatment, mode = 'mean.class')
cld_t_by_g_overall <- cld(emm_t_by_g_overall, adjust = 'tukey', Letters = letters, decreasing = TRUE)
cld_g_by_t_overall <- cld(emm_g_by_t_overall, adjust = 'tukey', Letters = LETTERS, decreasing = TRUE)Means by genotype
means_g_table <- bind_rows(as_tibble(cld_g_by_year), as_tibble(cld_g_overall)) %>%
filter(!is.na(mean.class)) %>%
mutate(Year = fct_na_value_to_level(Year, 'Combined years'), .group = trimws(.group),
across(c(mean.class, asymp.LCL, asymp.UCL), ~ round(., 2)),
display = glue::glue('{mean.class}{.group} [{asymp.LCL}, {asymp.UCL}]')) %>%
pivot_wider(id_cols = Pedigree, names_from = Year, values_from = display) %>%
arrange(Pedigree)
knitr::kable(means_g_table)| Pedigree | 2023 | 2024 | Combined years |
|---|---|---|---|
| Ab24E | 6.75a [6.43, 7.07] | 5.72ab [5.34, 6.09] | 6.23a [5.99, 6.48] |
| Tx601 | 6.85a [6.54, 7.17] | 5.48b [5.08, 5.87] | 6.17a [5.92, 6.42] |
| T173 | 6.35b [6.01, 6.69] | 5.77a [5.4, 6.14] | 6.06a [5.81, 6.31] |
| Va35 | 6.14b [5.79, 6.5] | 5.11c [4.72, 5.5] | 5.63b [5.36, 5.89] |
| Mp496 | 4.43c [4.02, 4.84] | 4.14d [3.74, 4.54] | 4.28c [4, 4.57] |
| Mp704 | 3.64de [3.23, 4.06] | 3.39f [3.01, 3.76] | 3.52d [3.24, 3.79] |
| Mp705 | 3.83d [3.41, 4.25] | 3.31f [2.94, 3.69] | 3.57d [3.29, 3.85] |
| Mp706 | NA | 3.7e [3.29, 4.1] | NA |
| Mp707 | 3.29ef [2.91, 3.67] | NA | NA |
| Mp708 | 3.17f [2.81, 3.53] | 3.22f [2.87, 3.57] | 3.2e [2.95, 3.44] |
| Mp714 | 3.03f [2.72, 3.35] | 3.12f [2.79, 3.46] | 3.08e [2.85, 3.3] |
Means by treatment
Note that the means separation that spans both years is the average of the nine genotypes that were present in both years, not including Mp706 and Mp707.
means_t_table <- bind_rows(as_tibble(cld_t_2023), as_tibble(cld_t_2024), as_tibble(cld_t_overall)) %>%
mutate(Year = rep(c('2023', '2024', 'Combined years'), each = 5), .group = trimws(.group),
across(c(mean.class, asymp.LCL, asymp.UCL), ~ round(., 2)),
display = glue::glue('{mean.class}{.group} [{asymp.LCL}, {asymp.UCL}]')) %>%
pivot_wider(id_cols = Treatment, names_from = Year, values_from = display) %>%
arrange(rev(Treatment))
knitr::kable(means_t_table)| Treatment | 2023 | 2024 | Combined years |
|---|---|---|---|
| 50 | 5.22a [4.89, 5.54] | 4.87a [4.53, 5.22] | 5.18a [4.95, 5.42] |
| 40 | 5.03b [4.71, 5.35] | 4.49b [4.15, 4.82] | 4.89b [4.66, 5.12] |
| 30 | 4.93b [4.62, 5.24] | 4.35bc [4.02, 4.68] | 4.74c [4.51, 4.96] |
| 20 | 4.62c [4.3, 4.95] | 4.3c [3.97, 4.62] | 4.56d [4.33, 4.79] |
| 10 | 3.95d [3.67, 4.22] | 3.48d [3.19, 3.77] | 3.82e [3.61, 4.02] |
Means by treatment x genotype in 2023
# Remove combinations where all letters are a.
cld_t_by_gxyear_table <- as_tibble(cld_t_by_gxyear) %>%
group_by(Pedigree, Year) %>%
mutate(ngrps = length(unique(.group)),
group_t = if_else(ngrps == 1, '', trimws(.group))) %>%
dplyr::select(-c(SE, df, .group, ngrps)) %>%
filter(!is.na(mean.class))
cld_g_by_txyear_table <- as_tibble(cld_g_by_txyear) %>%
group_by(Treatment, Year) %>%
mutate(ngrps = length(unique(.group)),
group_g = if_else(ngrps == 1, '', trimws(.group))) %>%
dplyr::select(-c(SE, df, .group, ngrps)) %>%
filter(!is.na(mean.class))
cld_gt_table <- left_join(cld_t_by_gxyear_table, cld_g_by_txyear_table, by = join_by(Treatment, Pedigree, Year)) %>%
mutate(across(c(mean.class.x, asymp.LCL.x, asymp.UCL.x), ~ round(., 2)),
display = glue::glue('{mean.class.x}{group_t}{group_g} [{asymp.LCL.x}, {asymp.UCL.x}]'))
means_gt_2023_table <- cld_gt_table %>%
filter(Year == '2023') %>%
pivot_wider(id_cols = Treatment, names_from = Pedigree, values_from = display) %>%
arrange(Treatment)
knitr::kable(means_gt_2023_table)| Treatment | Ab24E | Tx601 | T173 | Va35 | Mp496 | Mp704 | Mp705 | Mp707 | Mp708 | Mp714 |
|---|---|---|---|---|---|---|---|---|---|---|
| 10 | 6.53A [6.08, 6.98] | 6.15dAB [5.66, 6.64] | 5.48bBC [4.99, 5.97] | 5.04cC [4.5, 5.58] | 3.41cD [2.93, 3.88] | 2.39cF [2.02, 2.76] | 2.97cDE [2.57, 3.36] | 2.65bEF [2.27, 3.04] | 2.62cEF [2.22, 3.02] | 2.22cF [1.96, 2.48] |
| 20 | 6.65A [6.24, 7.06] | 6.39cdAB [5.95, 6.83] | 5.9bB [5.44, 6.37] | 6.49abAB [6.07, 6.9] | 3.65cCD [3.15, 4.15] | 3.87abCD [3.16, 4.58] | 3.71bC [3.2, 4.22] | 3.37aCD [2.82, 3.93] | 3.01bcD [2.52, 3.5] | 3.2abCD [2.8, 3.61] |
| 30 | 6.91A [6.48, 7.33] | 6.72bcA [6.29, 7.14] | 6.54aA [6.14, 6.94] | 6.57aA [6.17, 6.97] | 5.48aB [4.98, 5.98] | 3.66bCD [3.11, 4.21] | 4.15abC [3.49, 4.81] | 3.35aCDE [2.83, 3.86] | 3.04bcDE [2.56, 3.53] | 2.9bE [2.51, 3.28] |
| 40 | 6.76A [6.35, 7.18] | 7.23bA [6.8, 7.67] | 7aA [6.61, 7.4] | 6.08bB [5.66, 6.51] | 4.95abC [4.41, 5.5] | 4.53aCD [3.91, 5.15] | 3.8bDE [3.24, 4.36] | 3.36aE [2.78, 3.95] | 3.34abE [2.87, 3.81] | 3.2abE [2.8, 3.61] |
| 50 | 6.9B [6.49, 7.32] | 7.78aA [7.44, 8.11] | 6.82aB [6.46, 7.18] | 6.53abB [6.1, 6.96] | 4.65bC [4.07, 5.23] | 3.76bDE [3.16, 4.37] | 4.53aCD [3.95, 5.1] | 3.73aE [3.23, 4.23] | 3.84aDE [3.29, 4.4] | 3.64aE [3.15, 4.12] |
Means by treatment x genotype in 2024
means_gt_2024_table <- cld_gt_table %>%
filter(Year == '2024') %>%
pivot_wider(id_cols = Treatment, names_from = Pedigree, values_from = display) %>%
arrange(Treatment)
knitr::kable(means_gt_2024_table)| Treatment | Ab24E | Tx601 | T173 | Va35 | Mp496 | Mp704 | Mp705 | Mp706 | Mp708 | Mp714 |
|---|---|---|---|---|---|---|---|---|---|---|
| 10 | 5.16cA [4.67, 5.66] | 5.13bA [4.58, 5.68] | 4.88cA [4.35, 5.41] | 3.75bB [3.22, 4.28] | 3.4cBC [2.9, 3.9] | 2.74cDE [2.36, 3.12] | 2.54bDE [2.21, 2.87] | 2.93bCD [2.51, 3.35] | 2.36cE [2.03, 2.68] | 1.9cF [1.69, 2.11] |
| 20 | 5.49bcA [5, 5.97] | 5.38bA [4.9, 5.87] | 5.76bA [5.3, 6.21] | 5.49aA [5.04, 5.94] | 4.55bB [4.02, 5.08] | 2.58cE [2.2, 2.96] | 3.72aC [3.17, 4.28] | 3.7aC [3.21, 4.2] | 2.87bDE [2.5, 3.23] | 3.41abCD [2.95, 3.88] |
| 30 | 5.81abAB [5.35, 6.27] | 5.26bBC [4.75, 5.78] | 6.07abA [5.66, 6.49] | 5.2aC [4.73, 5.66] | 3.17cE [2.73, 3.6] | 3.61bDE [3.13, 4.08] | 3.31aE [2.79, 3.83] | 4.19aD [3.68, 4.71] | 3.64aDE [3.14, 4.14] | 3.23bE [2.81, 3.66] |
| 40 | 6.2aA [5.79, 6.62] | 5.67abAB [5.19, 6.15] | 5.8bAB [5.35, 6.24] | 5.47aB [4.99, 5.95] | 4.34bC [3.81, 4.87] | 3.53bCD [2.87, 4.2] | 3.34aD [2.85, 3.82] | 3.8aCD [3.27, 4.33] | 3.49aD [3.01, 3.98] | 3.21bD [2.8, 3.62] |
| 50 | 5.93abAB [5.47, 6.38] | 5.94aAB [5.49, 6.4] | 6.35aA [5.95, 6.75] | 5.64aBC [5.18, 6.1] | 5.25aC [4.77, 5.73] | 4.48aD [3.92, 5.05] | 3.67aE [3.08, 4.25] | 3.87aDE [3.32, 4.42] | 3.76aE [3.27, 4.25] | 3.87aDE [3.34, 4.39] |
Means by treatment x genotype combined across years
Note that neither Mp706 nor Mp707 are included here because they do not have means that span both years.
# Remove combinations where all letters are a.
cld_t_by_g_overall_table <- as_tibble(cld_t_by_g_overall) %>%
group_by(Pedigree) %>%
mutate(ngrps = length(unique(.group)),
group_t = if_else(ngrps == 1, '', trimws(.group))) %>%
dplyr::select(-c(SE, df, .group, ngrps)) %>%
filter(!is.na(mean.class))
cld_g_by_t_overall_table <- as_tibble(cld_g_by_t_overall) %>%
group_by(Treatment) %>%
mutate(ngrps = length(unique(.group)),
group_g = if_else(ngrps == 1, '', trimws(.group))) %>%
dplyr::select(-c(SE, df, .group, ngrps)) %>%
filter(!is.na(mean.class))
cld_gt_overall_table <- left_join(cld_t_by_g_overall_table, cld_g_by_t_overall_table, by = join_by(Treatment, Pedigree)) %>%
mutate(across(c(mean.class.x, asymp.LCL.x, asymp.UCL.x), ~ round(., 2)),
display = glue::glue('{mean.class.x}{group_t}{group_g} [{asymp.LCL.x}, {asymp.UCL.x}]'))
means_gt_overall_table <- cld_gt_overall_table %>%
pivot_wider(id_cols = Treatment, names_from = Pedigree, values_from = display) %>%
arrange(Treatment)
knitr::kable(means_gt_overall_table)| Treatment | Ab24E | Tx601 | T173 | Va35 | Mp496 | Mp704 | Mp705 | Mp708 | Mp714 |
|---|---|---|---|---|---|---|---|---|---|
| 10 | 5.85cA [5.51, 6.18] | 5.64cAB [5.27, 6.01] | 5.18cB [4.82, 5.54] | 4.4bC [4.02, 4.77] | 3.4dD [3.06, 3.75] | 2.56dE [2.3, 2.83] | 2.75cE [2.5, 3.01] | 2.49dE [2.23, 2.75] | 2.06cF [1.89, 2.23] |
| 20 | 6.07bcA [5.75, 6.39] | 5.89cA [5.56, 6.21] | 5.83bA [5.5, 6.16] | 5.99aA [5.68, 6.3] | 4.1cB [3.74, 4.46] | 3.23cCD [2.82, 3.63] | 3.72abBC [3.34, 4.09] | 2.94cD [2.63, 3.24] | 3.31bCD [3, 3.62] |
| 30 | 6.36abA [6.05, 6.67] | 5.99cAB [5.66, 6.33] | 6.31aA [6.02, 6.6] | 5.88aB [5.57, 6.19] | 4.32bcC [3.99, 4.65] | 3.63bcD [3.27, 4] | 3.73abD [3.31, 4.15] | 3.34bDE [2.99, 3.69] | 3.06bE [2.78, 3.35] |
| 40 | 6.48aA [6.19, 6.78] | 6.45bA [6.13, 6.78] | 6.4aA [6.1, 6.7] | 5.78aB [5.46, 6.1] | 4.65abC [4.27, 5.03] | 4.03abD [3.58, 4.49] | 3.57bDE [3.2, 3.94] | 3.42abE [3.08, 3.75] | 3.2bE [2.92, 3.49] |
| 50 | 6.41abBC [6.11, 6.72] | 6.86aA [6.58, 7.14] | 6.59aAB [6.31, 6.86] | 6.08aC [5.77, 6.4] | 4.95aD [4.58, 5.32] | 4.12aE [3.71, 4.54] | 4.1aE [3.68, 4.51] | 3.8aE [3.43, 4.17] | 3.75aE [3.39, 4.11] |
Figure
The multiple-comparison letters for 2023 are in lowercase and the letters for 2024 are in uppercase, to make it clear that the comparisons are done within each year, not between.
fig_dat <- as_tibble(cld_g_by_txyear) %>%
mutate(label = if_else(Year == 2023, tolower(trimws(.group)), trimws(.group)))
pd <- position_dodge(width = .5)
ggplot(fig_dat, aes(x = Pedigree, y = mean.class, ymin = asymp.LCL, ymax = asymp.UCL, label = label, group = interaction(Pedigree, Year))) +
geom_errorbar(position = pd, width = .25) +
geom_point(aes(color = Year), position = pd, size = 2) +
geom_text(aes(y = asymp.UCL), position = pd, vjust = -.4, size = 3) +
facet_wrap(~ Treatment, labeller = as_labeller(function(x) paste(x, 'Neonates infestation rate')), ncol = 2) +
scale_color_grey() +
theme(legend.position = 'inside',
legend.position.inside = c(0.75, 0.1),
legend.title = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = 'Genotype', y = 'Leaf-feeding damage score')