Change log

  • 17 Aug 2023: add raw data to figures where possible
  • 14 Aug 2023: changes to model as a result of statisticians’ review
  • 15 Mar 2023: version for archiving on Ag Data Commons
  • 28 Feb 2023: improve figure formatting for manuscript
  • 27 Oct 2022: cleaned up version, includes tables at the end
  • 05 Oct 2022: make nicer looking figures and compute odds ratio for the interaction model
  • 27 April 2022: update with new data
  • 26 April 2022: model includes interactions with variety
  • 19 April 2022: completely updated version with new data and logistic regression instead of null model
  • 23 February 2022: first version

Summary

This notebook includes the R code needed to reproduce the main statistical model fitting and post-processing of model output from the manuscript:

Penn, H. J. and Q. D. Read. 2023. Stem borer damage dependent on interactions of plant variety, age, and prior conspecific herbivory. Functional Ecology, in review. update this when MS is accepted.

A logistic regression was used to look at whether the probability of damage to an internode is a function of the number of internodes previously damaged on that same stalk, as well as whether that probability changes over time, with variety, or with crop cycle. Interactions are also included. We only used varieties that had data from >= 5 years and excluded treated plots. We generated posterior estimated marginal means and trends that are shown as figures and tables.

A follow-up analysis exploring the influence of variety-level traits on probability of internode damage is included in a second notebook 02_trait_covariate_analysis.Rmd.

Setup

Load packages, set options for the brms modeling package, and set plotting theme to a simpler one. Load data, including lookup table of variety resistance classes. Need to specify column types because the automatic column type checker fails with so many empty cells.

Note: The Rutilitybelt package is one I wrote for myself. You can install it yourself by calling remotes::install_github('qdread/Rutilitybelt').

library(tidyverse)

library(brms)

library(emmeans)

library(Rutilitybelt)

library(gt)

library(brmsmargins)



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



theme_set(

  theme_bw() + theme(panel.grid = element_blank(), strip.background = element_blank())

)



dat <- read_csv('project/data/BoredInternodes_26April2022_no format.csv', col_types = paste0('icccccffi', strrep('i', 30), 'cc'))

lookup <- read_csv('project/data/variety_lookup.csv') %>%

  mutate(check = grepl('\\*', `SCB resistance`),

         resistance = gsub('\\*', '', `SCB resistance`))



dat_long <- dat %>%

  select(-c(`# Internodes`, Notes)) %>%

  pivot_longer(cols = `1`:`30`, names_to = 'internode', values_to = 'damage') %>%

  filter(!is.na(damage)) %>%

  mutate(internode = as.integer(internode))

Correct typos and formatting issues in data.

dat_long <- dat_long %>% mutate(

  damage = if_else(damage == 22, 2L, damage),

  Stubble = if_else(Stubble == 'Plant cane', 'Plant Cane', Stubble),

  Block = factor(substr(as.character(Plot), 1, 1)),

  Variety = gsub(' ', '', gsub('-', '_', Variety))

)



lookup <- mutate(lookup, Variety_nospace = gsub(' ', '', gsub('-', '_', Variety)))

Exploring data

Are there trends over time in average percent of damaged internodes by stalk, and do these trends vary by variety? Remove data from experiments with different treatments. First get proportion damage by stalk then average the proportion over plots.

damage_over_time_by_stalk <- dat_long %>%

  filter(grepl('YRS', Experiment)) %>%

  group_by(Year, Variety, Plot, Stalk) %>%

  summarize(prop_damaged = sum(damage > 0)/length(damage))



damage_over_time_by_plot <- damage_over_time_by_stalk %>%

  group_by(Year, Variety, Plot) %>%

  summarize(prop_damaged = mean(prop_damaged))

Plot trends over time by variety. Show only those varieties (16 of them) that were planted in at least five different years. Put some regression lines on there just to guide the eye to see if a trend may exist. It looks like damage is decreasing over time in many varieties.

varieties_5yrs <- damage_over_time_by_plot %>%

  group_by(Variety) %>%

  filter(length(unique(Year)) >= 5) %>%

  pull(Variety) %>%

  unique
ggplot(damage_over_time_by_plot %>% filter(Variety %in% varieties_5yrs), aes(x = Year, y = prop_damaged)) +

  geom_point(alpha = 0.5) +

  facet_wrap(~ Variety) +

  labs(y = 'proportion damaged stalks in each plot') +

  stat_smooth(method = 'lm', se = FALSE)

Statistical model

We will model damage as a binary variable (yes or no) and ask whether the probability of damage on a stalk is a function of variety, stubble type, year, and damage on lower-numbered internodes on that same stalk. This is a logistic regression with multiple covariates.

The following interaction terms are also included:

  • interaction between variety and crop cycle (plant cane or stubble)
  • interaction between variety and previous damage

We have random intercepts for each plot and each stalk (not sure whether the random intercept for each stalk is needed but it is probably helpful). Reasonable priors are set on the fixed effects.

We need to create an additional covariate, previous_damage, which is the number of damaged internodes (binary yes/no) lower on the same stalk. That will be one of the fixed effects in the model. Also calculate this as a proportion of damaged internodes lower on the stalk.

dat_long_withcov <- dat_long %>%

  arrange(Year, Date, Location, Experiment, Stubble, Variety, Block, Plot, Stalk, Treatment, internode) %>%

  group_by(Year, Date, Location, Experiment, Stubble, Variety, Block, Plot, Stalk, Treatment) %>%

  mutate(previous_damage = lag(cumsum(damage > 0), default = 0),

         previous_damage_proportion = c(0, (previous_damage / 0:(length(damage)-1))[-1])) %>%

  ungroup

Now fit the model as described above.

dat_fit <- dat_long_withcov %>% 

  filter(grepl('YRS', Experiment)) %>% 

  dplyr::select(-c(Location, Experiment, Treatment)) %>%

  mutate(damage_binary = as.integer(damage > 0),

         Year = Year - min(Year))



dat_fit_16vars <- filter(dat_fit, Variety %in% varieties_5yrs)

The following model was fit on the USDA SciNet computing cluster.

binom_fit_interactions <- brm(

  damage_binary ~ internode + previous_damage + Year + Stubble + Variety + Variety:internode + Variety:previous_damage + Variety:Year + Variety:Stubble + (1|Year) + (1|Plot:Year) + (1|Stalk:Plot:Year), 

  data = dat_fit_16vars, family = bernoulli(link = 'logit'),

  prior = c(

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

  ),

  chains = 4, iter = 4000, warmup = 3000, seed = 1338,

  file = 'project/binom_fit_intxns_updated_only5yrs'

)

Here is the R-hat convergence diagnostic showing that the model converged:

max(rhat(binom_fit_interactions))
## [1] 1.011187

Results

Posterior estimate plots

Now let’s take a look at some of the posterior estimates. Because the response variable is probability of damage of an internode, we used the logit link function in the model. We back-transform the estimates with plogis() (inverse of the logit) to get it to an interpretable probability (ranging from 0 to 1 where 1 is the highest probability of being damaged).

In all the figures, the point estimate is the median of the posterior samples, the thick error bar is the 66% quantile credible interval (roughly \(\pm 1 SD\)) and the thin error bar is the 95% credible interval (roughly \(\pm 2 SD\)).

Here are the estimated marginal means for each variety. For those that don’t have any first stubble data, first stubble estimates are not accounted for.

emm_variety <- emmeans(binom_fit_interactions, ~ Variety, weights = 'flat')

emm_variety_quant <- emm_quantile(emm_variety) %>%

  mutate(across(starts_with('q'), plogis)) %>%

  mutate(Variety = fct_reorder(Variety, q0.5))

Varieties are colored by their resistance level in different shades of blue and red. Data for each variety are plotted as semitransparent points (each point is the mean proportion damaged for a single plot in a single year).

lookup_short <- lookup %>% select(Variety = Variety_nospace, Variety_label = Variety, check, resistance)



plotdat_variety <- left_join(emm_variety_quant, lookup_short) %>%

  mutate(Variety = fct_reorder(Variety_label, q0.5),

         resistance = factor(resistance, levels = c('Susceptible', 'Moderately Susceptible', 'Moderately Resistant', 'Resistant'))) 



p_year_var_plot <- left_join(damage_over_time_by_plot, lookup_short) %>%

  filter(Variety %in% varieties_5yrs) %>%

  mutate(Variety = factor(Variety_label, levels = levels(plotdat_variety$Variety)))



dx <- 0.15



(p <- ggplot(plotdat_variety, aes(x = Variety, y = q0.5)) +

  geom_point(aes(x = Variety, y = prop_damaged), data = p_year_var_plot, alpha = .25, position = position_nudge(x = dx), size = 0.8) +

  geom_errorbar(aes(ymin = q0.17, ymax = q0.83, color = resistance), size = 2.5, width = 0, position = position_nudge(x = -dx)) +

  geom_errorbar(aes(ymin = q0.025, ymax = q0.975, color = resistance), size = 0.8, width = 0, position = position_nudge(x = -dx)) +

  geom_point(size = 1.5, position = position_nudge(x = -dx)) +

  labs(y = 'Estimated proportion of bored internodes', x = 'Variety') +

  coord_flip() +

  theme(legend.position = c(0.72, 0.2), legend.title = element_blank()) +

  scale_color_brewer(palette = 'RdBu') + scale_fill_brewer(palette = 'RdBu')

)

#ggsave('project/MS/figs_revision/2b.png', p, height = 4, width = 5, dpi = 400)

Are there any interactions between variety and crop cycle? We take the estimated marginal means for each combination of variety and crop cycle and plot them. Not all the varieties are represented in the first stubble category (6 of them were not part of studies where first stubble was tracked) so varieties that have only a single crop year are excluded. We see very little effect of crop cycle.

varieties_1ststubble <- dat_fit %>%

  filter(Stubble %in% '1st Stubble') %>%

  pull(Variety) %>%

  unique



emm_variety_cropcycle <- emmeans(binom_fit_interactions, ~ Variety + Stubble)

emm_variety_cropcycle_quant <- emm_quantile(emm_variety_cropcycle) %>%

  mutate(across(starts_with('q'), plogis)) %>%

  filter(Variety %in% varieties_1ststubble) %>%

  mutate(Variety = factor(Variety, levels = levels(emm_variety_quant$Variety)))

In the figure, raw data are shown as semitransparent points (the points are proportion damaged internodes for a single plot in a single year, colored by crop cycle)

pd <- position_dodge(width = 0.2)

pdwide <- position_dodge(width = 1)



plotdat_1ststubble <- left_join(emm_variety_cropcycle_quant, lookup_short)%>% 

  mutate(Variety = factor(Variety_label, levels = levels(plotdat_variety$Variety)),

         Stubble = factor(Stubble, levels = c('Plant Cane', '1st Stubble'), labels = c('plant cane', '1<sup>st</sup> ratoon')))



p_year_var_ccycle <- dat_fit_16vars %>%

  group_by(Year, Variety, Plot, Stubble) %>%

  summarize(p = sum(damage_binary)/n()) %>%

  mutate(Stubble = factor(Stubble, levels = c('Plant Cane', '1st Stubble'), labels = c('plant cane', '1<sup>st</sup> ratoon'))) %>%

  left_join(lookup_short) %>%

  filter(Variety %in% varieties_5yrs) %>%

  mutate(Variety = factor(Variety_label, levels = levels(plotdat_variety$Variety))) %>%

  filter(Variety %in% plotdat_1ststubble$Variety)



emm_stubble_withint <- emmeans(binom_fit_interactions, ~ Stubble)

emm_stubble_quant_withint <- emm_quantile(emm_stubble_withint) %>%

  mutate(across(starts_with('q'), plogis),

         Stubble = factor(Stubble, levels = c('Plant Cane', '1st Stubble'), labels = c('plant cane', '1<sup>st</sup> ratoon')))



(p <- ggplot(plotdat_1ststubble, aes(x = Variety, y = q0.5, color = Stubble, group = interaction(Variety, Stubble))) +

  geom_point(aes(y = p), data = p_year_var_ccycle, alpha = .25, position = pdwide, size = 1) +

  geom_errorbar(aes(ymin = q0.17, ymax = q0.83), size = 2.5, width = 0, position = pd) +

  geom_errorbar(aes(ymin = q0.025, ymax = q0.975), size = 0.8, width = 0, position = pd) +

  geom_point(size = 1.5, color = 'black', position = pd) +

  labs(y = 'Proportion of bored internodes', x = 'Variety') +

  coord_flip() +

  scale_color_brewer(palette = 'Dark2') +

  theme(legend.position = c(0.8, 0.12), legend.title = element_blank(), legend.text = ggtext::element_markdown())

)

#ggsave('project/MS/figs_revision/3b.png', p, height = 4, width = 5, dpi = 400)

Looking at the main effect of crop cycle in the interaction model, we see essentially no difference.

(p <- ggplot(emm_stubble_quant_withint, aes(x = Stubble, y = q0.5, color = Stubble)) +

  geom_jitter(aes(y = p), data = p_year_var_ccycle, alpha = .25, size = 0.8, width = 0.1, height = 0) +

  geom_errorbar(aes(ymin = q0.17, ymax = q0.83), size = 2.5, width = 0) +

  geom_errorbar(aes(ymin = q0.025, ymax = q0.975), size = 0.8, width = 0) +

  geom_point(size = 1.5, color = 'black') +

  scale_color_brewer(palette = 'Dark2') +

  coord_flip() +

  scale_x_discrete(limits = rev) +

  labs(y = 'Proportion of bored internodes') +

  theme(axis.title.y = element_blank(), legend.position = 'none', axis.text.y = ggtext::element_markdown())

)

#ggsave('project/MS/figs_revision/3a.png', p, height = 4, width = 4, dpi = 400)

The odds ratio is almost exactly 1 indicating no difference.

contrast(emm_stubble_withint, 'pairwise', type = 'response')
##  contrast                 odds.ratio lower.HPD upper.HPD

##  1st Stubble / Plant Cane      0.947     0.555      1.51

## 

## Results are averaged over the levels of: Variety 

## Point estimate displayed: median 

## Results are back-transformed from the log odds ratio scale 

## HPD interval probability: 0.95

In the interaction model, we see the main effects of previous damage and the time trend over year. This is now done with marginal effects instead of conditional effects.

marginal_draws_prevdamage <- brmsmargins(binom_fit_interactions, at = data.frame(previous_damage = -2:16))

marginal_summ_prevdamage <- Rutilitybelt::pred_quantile(x_pred = data.frame(previous_damage = -2:16), y_pred = marginal_draws_prevdamage$Posterior)



marginal_draws_year <- brmsmargins(binom_fit_interactions, at = data.frame(Year = -2:30)) # 1991 to 2023 (year 0 is 1993)

marginal_summ_year <- Rutilitybelt::pred_quantile(x_pred = data.frame(Year = -2:30), y_pred = marginal_draws_year$Posterior)



save(marginal_draws_prevdamage, marginal_summ_prevdamage, marginal_draws_year, marginal_summ_year, file = 'project/marginal_trends.RData')

In the following plots, distributions of the proportion damaged internodes are shown at the year x variety x plot level for the time trend plot and at the year x variety x previous damage class for the previous damage trend plot.

# Bin up the previous damage values by year and variety

p_year_var_prevdmg <- dat_fit_16vars %>%

  group_by(Year, Variety, previous_damage) %>%

  summarize(p = sum(damage_binary)/n())



(ptrend_prevdamage <- ggplot(marginal_summ_prevdamage, aes(x = previous_damage, y = q0.5)) + 

   geom_boxplot(aes(y = p, x = previous_damage, group = previous_damage), data = p_year_var_prevdmg %>% filter(previous_damage %in% 0:13), width = 0.5) +  

   geom_ribbon(aes(ymin = q0.025, ymax = q0.975), alpha = 0.3, color = NA) + 

   geom_ribbon(aes(ymin = q0.17, ymax = q0.83), alpha = 0.3, color = NA) + 

   geom_line(linewidth = 0.8) +

   coord_cartesian(xlim = c(-0.5, 13.5)) +

   scale_x_continuous(name = 'Number of previously damaged internodes on stalk', expand = expansion(mult = c(0, 0)), breaks = 0:13) +

   labs(y = 'Proportion of bored internodes') +

   theme(plot.margin = unit(c(2, 4, 2, 2), 'mm')))

(ptrend_year <- ggplot(marginal_summ_year, aes(x = Year + 1993, y = q0.5)) + 

    geom_boxplot(aes(y = prop_damaged, x = Year, group = Year), data = p_year_var_plot, width = 0.5) +

    geom_ribbon(aes(ymin = q0.025, ymax = q0.975), alpha = 0.3, color = NA) + 

    geom_ribbon(aes(ymin = q0.17, ymax = q0.83), alpha = 0.3, color = NA) + 

    geom_line(linewidth = 0.8) + 

    coord_cartesian(xlim = c(1992, 2022)) +

    scale_x_continuous(name = 'Year', expand = expansion(mult = c(0, 0)), breaks = seq(1995,2020,5)) +

    labs(y = 'Proportion of bored internodes') +

    theme(plot.margin = unit(c(2, 4, 2, 2), 'mm'))) 

#ggsave('project/MS/figs_revision/2a.png', ptrend_year, height = 4, width = 4, dpi = 400)

#ggsave('project/MS/figs_revision/4a.png', ptrend_prevdamage, height = 4, width = 4, dpi = 400)

What about the interaction between variety and amount of previous damage? We can plot the fitted slope of the previous damage trend for each variety. After revision, we see more uncertainty in this relationship. One variety even has a negative relationship, though most have positive or near zero relationship.

emt_var_damage <- emtrends(binom_fit_interactions, ~ Variety, var = 'previous_damage')

emt_var_damage_quant <- emm_quantile(emt_var_damage) %>%

  mutate(Variety = fct_reorder(Variety, q0.5))
plotdat_damage <- left_join(emt_var_damage_quant, lookup_short) %>%

  mutate(Variety = fct_reorder(Variety_label, q0.5),

         resistance = factor(resistance, levels = c('Susceptible', 'Moderately Susceptible', 'Moderately Resistant', 'Resistant'))) 



(p <- ggplot(plotdat_damage, aes(x = Variety, y = q0.5, color = resistance)) +

  geom_errorbar(aes(ymin = q0.17, ymax = q0.83), size = 2.5, width = 0) +

  geom_errorbar(aes(ymin = q0.025, ymax = q0.975), size = 0.8, width = 0) +

  geom_point(size = 1.5, color = 'black') +

  labs(y = 'Strength of previous damage effect', x = 'Variety') +

  geom_hline(yintercept = 0, linetype = 'dotted', size = .8, color = 'slateblue') +

  coord_flip() +

  theme(legend.position = c(0.72, 0.2), legend.title = element_blank()) +

  scale_color_brewer(palette = 'RdBu')

)

#ggsave('project/MS/figs_revision/4b.png', p, height = 4, width = 5, dpi = 400)

Now, we can also look at the trend over time in total amount of damage, separately by variety (interaction between variety and year). Almost all show a negative trend though the magnitude varies quite a lot.

emt_var_time <- emtrends(binom_fit_interactions, ~ Variety, var = 'Year')

emt_var_time_quant <- emm_quantile(emt_var_time) %>%

  mutate(Variety = fct_reorder(Variety, q0.5))
plotdat_time <- left_join(emt_var_time_quant, lookup_short) %>%

  mutate(Variety = fct_reorder(Variety_label, q0.5),

         resistance = factor(resistance, levels = c('Susceptible', 'Moderately Susceptible', 'Moderately Resistant', 'Resistant'))) 



(p <- ggplot(plotdat_time, aes(x = Variety, y = q0.5, color = resistance)) +

  geom_errorbar(aes(ymin = q0.17, ymax = q0.83), size = 2.5, width = 0) +

  geom_errorbar(aes(ymin = q0.025, ymax = q0.975), size = 0.8, width = 0) +

  geom_point(size = 1.5, color = 'black') +

  labs(y = 'Change in damage probability over time', x = 'Variety') +

  geom_hline(yintercept = 0, linetype = 'dotted', size = .8, color = 'slateblue') +

  coord_flip() +

  theme(legend.position = c(0.22, 0.87), legend.title = element_blank(), legend.background = element_blank()) +

  scale_color_brewer(palette = 'RdBu')

)

#ggsave('project/MS/figs_revision/newfigure.png', p, height = 4, width = 5, dpi = 400)

Posterior estimate tables

These are the same values as in the plots above, except they are in table form. 95% quantile credible interval bounds are provided. In addition, these tables provide the odds ratios for first stubble damage versus plant cane damage for each variety separately, as well as the uncertainty intervals, which is not presented in any figure.

Marginal means: variety

plotdat_variety %>%

  select(Variety_label, resistance, q0.5, q0.025, q0.975) %>%

  gt() %>%

  cols_label(Variety_label = 'variety', q0.5 = 'median estimate', q0.025 = '95% QI lower', q0.975 = '95% QI upper') %>%

  fmt_number(columns = 3:5, n_sigfig = 3) %>%

  tab_options(column_labels.font.weight = 'bold')
variety resistance median estimate 95% QI lower 95% QI upper
CP 70-321 Resistant 0.0320 0.00947 0.103
Ho 07-613 Moderately Resistant 0.0565 0.0253 0.123
Ho 95-988 Susceptible 0.155 0.108 0.213
HoCP 00-950 Susceptible 0.172 0.113 0.253
HoCP 04-838 Resistant 0.0241 0.0124 0.0479
HoCP 09-804 Moderately Susceptible 0.0658 0.0287 0.147
HoCP 85-845 Resistant 0.0617 0.0420 0.0893
HoCP 91-555 Susceptible 0.103 0.0601 0.170
HoCP 96-540 Susceptible 0.0920 0.0599 0.134
L 01-283 Moderately Resistant 0.0718 0.0458 0.108
L 01-299 Moderately Resistant 0.0505 0.0300 0.0828
L 97-128 Susceptible 0.105 0.0722 0.145
L 99-226 Moderately Resistant 0.0534 0.0365 0.0753
L 99-233 Susceptible 0.0983 0.0685 0.135
LCP 85-384 Susceptible 0.0798 0.0412 0.144
LHo 83-153 Resistant 0.0418 0.0108 0.143

Marginal means: variety by crop year

plotdat_1ststubble %>%

  select(Variety_label, resistance, Stubble, q0.5, q0.025, q0.975) %>%

  arrange(Variety_label, resistance, Stubble) %>%

  mutate(Stubble = map(as.character(Stubble), html)) %>%

  gt() %>%

  cols_label(Variety_label = 'variety', Stubble = 'crop year', q0.5 = 'median estimate', q0.025 = '95% QI lower', q0.975 = '95% QI upper') %>%

  fmt_number(columns = 4:6, n_sigfig = 3) %>%

  tab_options(column_labels.font.weight = 'bold')
variety resistance crop year median estimate 95% QI lower 95% QI upper
CP 70-321 Resistant plant cane 0.0295 0.00936 0.0887
CP 70-321 Resistant 1st ratoon 0.0347 0.00933 0.130
Ho 07-613 Moderately Resistant plant cane 0.0681 0.0319 0.137
Ho 07-613 Moderately Resistant 1st ratoon 0.0465 0.0166 0.127
HoCP 00-950 Susceptible plant cane 0.161 0.111 0.222
HoCP 00-950 Susceptible 1st ratoon 0.183 0.102 0.315
HoCP 04-838 Resistant plant cane 0.0298 0.0168 0.0512
HoCP 04-838 Resistant 1st ratoon 0.0196 0.00732 0.0544
HoCP 09-804 Moderately Susceptible plant cane 0.0943 0.0399 0.210
HoCP 09-804 Moderately Susceptible 1st ratoon 0.0459 0.0163 0.123
HoCP 85-845 Resistant plant cane 0.0563 0.0380 0.0784
HoCP 85-845 Resistant 1st ratoon 0.0681 0.0381 0.121
HoCP 96-540 Susceptible plant cane 0.0792 0.0548 0.109
HoCP 96-540 Susceptible 1st ratoon 0.107 0.0577 0.182
L 01-299 Moderately Resistant plant cane 0.0446 0.0280 0.0694
L 01-299 Moderately Resistant 1st ratoon 0.0573 0.0294 0.106
LCP 85-384 Susceptible plant cane 0.0838 0.0502 0.134
LCP 85-384 Susceptible 1st ratoon 0.0762 0.0288 0.176
LHo 83-153 Resistant plant cane 0.0495 0.0149 0.151
LHo 83-153 Resistant 1st ratoon 0.0349 0.00715 0.146

Marginal means: crop year

emm_stubble_quant_withint %>%

  select(Stubble, q0.5, q0.025, q0.975) %>%

  arrange(Stubble) %>%

  mutate(Stubble = map(as.character(Stubble), html)) %>%

  gt() %>%

  cols_label(Stubble = 'crop year', q0.5 = 'median estimate', q0.025 = '95% QI lower', q0.975 = '95% QI upper') %>%

  fmt_number(columns = 2:4, n_sigfig = 3) %>%

  tab_options(column_labels.font.weight = 'bold')
crop year median estimate 95% QI lower 95% QI upper
plant cane 0.0726 0.0506 0.0989
1st ratoon 0.0685 0.0395 0.120

Odds ratios between first stubble and plant cane

Only the varieties that have data for both first stubble and plant cane are presented, as well as the overall which is also presented above. All odds ratios indicate little to no difference.

emm_cycle_byvariety <- emmeans(binom_fit_interactions, ~ Stubble | Variety, type = 'response')

contrast_cycle_byvariety <- contrast(emm_cycle_byvariety, 'pairwise')

contrast_cycle_overall <- contrast(emm_stubble_withint, 'pairwise', type = 'response')
contrast_cycle_byvariety %>%

  as.data.frame %>%

  filter(Variety %in% varieties_1ststubble) %>%

  bind_rows(as.data.frame(contrast_cycle_overall)) %>%

  left_join(lookup_short) %>%

  select(variety = Variety_label, resistance, odds.ratio, lower.HPD, upper.HPD) %>%

  mutate(variety = if_else(is.na(variety), 'Overall', variety)) %>%

  gt() %>%

  cols_label(odds.ratio = 'odds ratio', lower.HPD = '95% QI lower', upper.HPD = '95% QI upper') %>%

  fmt_number(columns = 3:5, n_sigfig = 3) %>%

  tab_options(column_labels.font.weight = 'bold') %>%

  fmt_missing(columns = 2)
variety resistance odds ratio 95% QI lower 95% QI upper
CP 70-321 Resistant 1.19 0.609 2.06
Ho 07-613 Moderately Resistant 0.663 0.192 1.39
HoCP 00-950 Susceptible 1.18 0.556 2.05
HoCP 04-838 Resistant 0.662 0.193 1.46
HoCP 09-804 Moderately Susceptible 0.452 0.144 0.979
HoCP 85-845 Resistant 1.24 0.599 2.09
HoCP 96-540 Susceptible 1.39 0.705 2.32
L 01-299 Moderately Resistant 1.30 0.737 2.04
LCP 85-384 Susceptible 0.897 0.301 1.75
LHo 83-153 Resistant 0.685 0.242 1.41
Overall — 0.947 0.555 1.51

Previous damage effects: variety

plotdat_damage %>%

  select(Variety_label, resistance, q0.5, q0.025, q0.975) %>%

  gt() %>%

  cols_label(Variety_label = 'variety', q0.5 = 'median estimate', q0.025 = '95% QI lower', q0.975 = '95% QI upper') %>%

  fmt_number(columns = 3:5, n_sigfig = 3) %>%

  tab_options(column_labels.font.weight = 'bold')
variety resistance median estimate 95% QI lower 95% QI upper
CP 70-321 Resistant −0.0704 −0.132 −0.0123
Ho 07-613 Moderately Resistant 0.239 0.0179 0.467
Ho 95-988 Susceptible 0.0283 −0.0132 0.0692
HoCP 00-950 Susceptible 0.0329 −0.00349 0.0696
HoCP 04-838 Resistant 0.210 0.0104 0.391
HoCP 09-804 Moderately Susceptible 0.444 0.309 0.583
HoCP 85-845 Resistant 0.0297 −0.0204 0.0831
HoCP 91-555 Susceptible 0.0357 −0.0109 0.0832
HoCP 96-540 Susceptible 0.0924 0.0408 0.145
L 01-283 Moderately Resistant 0.219 0.134 0.299
L 01-299 Moderately Resistant 0.218 0.116 0.318
L 97-128 Susceptible 0.0796 0.0295 0.131
L 99-226 Moderately Resistant 0.112 0.0509 0.169
L 99-233 Susceptible 0.0490 −0.00241 0.103
LCP 85-384 Susceptible 0.00229 −0.0386 0.0426
LHo 83-153 Resistant −0.00611 −0.0809 0.0671