Analysis of highbush blueberry cover crop experimental data
Quentin D. Read
January 03, 2025
This document contains all code necessary to reproduce the statistical analyses presented in the paper:
Smith, B. J., E. T. Stafne, and Q. D. Read. 2024. Establishment of southern highbush blueberry cultivars and suppression of Phytophthora root rot using cover crop and soil amendment treatments. PhytoFrontiers.
Summary
This is an analysis of all experimental data presented in the manuscript. In all cases, we fit linear mixed models (LMMs) aligned with the experimental design. We used generalized linear mixed models (GLMMs) in cases when we could not assume normally distributed residuals. Bayesian methods are used where necessary, and in some cases prior distributions are assumed within frequentist analyses. Tukey’s HSD is used as the p-value adjustment in all cases. We provide ANOVA tables or an ANOVA-like table where possible, but this was not possible for all models.
For Table 1, response variables are log-transformed where appropriate and a Gaussian LMM is fit to each one. In the Table 2 analysis there are repeated measures with a non-linear effect of time. Therefore we fit a generalized additive mixed model (GAMM), which includes a spline term that captures the non-linear time trend within each cover crop treatment. The Table 3 analysis includes several variables measured on a categorical scale. For those variables, a cumulative logistic mixed model (CLMM) is fit with an ordered multinomial response distribution. The binary variables in Table 3 are analyzed with a binomial GLMM. The CLMMs were fit with a Bayesian approach, using prior distributions on the fixed effect parameters to deal with the issue of complete separation (some treatment combinations that had no variation). Prior distributions were also used for some of the binomial GLMMs, without using a fully Bayesian approach, again to deal with the issue of complete separation. In Table 4, appropriate GLMMs are used (Poisson, Beta, and Gamma). For the variables that had only a few zeros, the zeros are removed, but for the live length variable, there are many zeros so a hurdle model was used to explicitly account for the probability of observing zero. In Table 5, LMMs with log-transformed response and a binomial GLMM were fit. In Table 6, a Poisson GLMM and linear mixed models with log-transformed response were fit. In Table 7, Poisson GLMM, binomial GLMM, and linear mixed model with log(x+0.1)-transformed response were fit.
Means separation tables for soil amendment treatment, cultivar (where present), and soil amendment treatment within each level of cultivar are presented.
Setup
Load packages and import data from each of the seven tables. Order treatment factor levels so that control is the reference level.
library(tidyverse)
library(readxl)
library(glmmTMB)
library(emmeans)
library(multcomp)
library(DHARMa)
library(ordinal)
library(mgcv)
library(gt)
library(brms)
library(tidybayes)
library(bayestestR)
library(patchwork)
library(rstan)
options(brms.backend = 'cmdstanr', brms.file_refit = 'on_change', mc.cores = 4)
excel_file <- 'Smith BB cover crop 7 data sets 10.9.24 .xlsx'
table1 <- read_excel(excel_file, sheet = 'Table 1', skip = 5) %>%
rename(Rep = rep, Date = date) %>%
mutate(across(c(Rep, Date), factor),
Cover = relevel(factor(Cover), ref = 'Control'))
table2 <- read_excel(excel_file, sheet = 'Table2', skip = 6) %>%
mutate(across(c(Rep, Date, Plot), factor),
Cover = relevel(factor(Cover), ref = 'Control'))
table3 <- read_excel(excel_file, sheet = 'Table 3', skip = 8, na = '.') %>%
rename(Rep = Row) %>%
mutate(across(c(Rep, Plant, Cultivar, Live921, Live622), factor),
across(c(Leaf, Flower, Crop, Ripe, Vig921, Vig622), ordered),
Cover = relevel(factor(Cover), ref = 'Control'))
table4 <- read_excel(excel_file, sheet = 'Table 4', skip = 7) %>%
rename(Cover = Crop) %>%
mutate(across(c(Plant, Cultivar, Rep), factor),
Cover = relevel(factor(Cover), ref = 'Control'))
table5 <- read_excel(excel_file, sheet = 'Table5', skip = 6) %>%
rename(Cover = Crop, Root_bound = `Root bound`) %>%
mutate(across(c(Cultivar, Rep, Root_bound), factor),
Cover = relevel(factor(Cover), ref = 'Control'))
table6 <- read_excel(excel_file, sheet = 'Table6', skip = 9) %>%
rename(weight_g = `Weight(g)`, wt_per_seedling = `Wt/seedling`) %>%
mutate(across(c(Assay, Rep), factor),
Cover = relevel(factor(Cover), ref = 'Control'))
table7 <- read_excel(excel_file, sheet = 'Table7', skip = 6) %>%
mutate(Rep = factor(Rep),
Cover = factor(Cover, levels = c('Control', 'Bark', 'Brassica', 'Grass', 'Legume', 'Ridomil', 'Ckpos', 'Ckneg', 'Ckneut')))Table 1 analysis
For these analyses, a two-way ANOVA is fit (cover by date including the interaction), with a random intercept for each replicate. The pH model assumes normally distributed residuals. For the models for elemental concentration and SS%, the response variable is log-transformed before fitting, and the predictions are back-transformed.
t1_pH_fit <- glmmTMB(pH ~ Cover * Date + (1 | Rep), data = table1)
t1_P_fit <- glmmTMB(log(P) ~ Cover * Date + (1 | Rep), data = table1)
t1_K_fit <- glmmTMB(log(K) ~ Cover * Date + (1 | Rep), data = table1)
t1_Mg_fit <- glmmTMB(log(Mg) ~ Cover * Date + (1 | Rep), data = table1)
t1_Zn_fit <- glmmTMB(log(Zn) ~ Cover * Date + (1 | Rep), data = table1)
t1_Ca_fit <- glmmTMB(log(Ca) ~ Cover * Date + (1 | Rep), data = table1)
t1_SS_fit <- glmmTMB(log(`SS%`) ~ Cover * Date + (1 | Rep), data = table1)Model residual diagnostics are acceptable (not shown).
plot(simulateResiduals(t1_pH_fit, n = 1000))
plot(simulateResiduals(t1_P_fit, n = 1000))
plot(simulateResiduals(t1_K_fit, n = 1000))
plot(simulateResiduals(t1_Mg_fit, n = 1000))
plot(simulateResiduals(t1_Zn_fit, n = 1000))
plot(simulateResiduals(t1_SS_fit, n = 1000))Table 1 results: ANOVA tables
Generate the F-statistics (ANOVA tables) using joint F-tests for each model term. Some of these have significant interactions between cover and date. Note here and elsewhere the p-values below 0.00001 are given as “< 0.00001” in the tables.
t1_fits <- tibble(variable = c('pH', 'P concentration', 'K concentration', 'Mg concentration', 'Zn concentration', 'Ca concentration', 'SS%'),
fit = list(t1_pH_fit, t1_P_fit, t1_K_fit, t1_Mg_fit, t1_Zn_fit, t1_Ca_fit, t1_SS_fit)) %>%
mutate(ANOVA = map(fit, joint_tests))
t1_ANOVAtables <- with(t1_fits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y))))
t1_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(F.ratio, decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts| model term | df1 | df2 | F.ratio | p.value |
|---|---|---|---|---|
| pH | ||||
| Cover | 5 | 46 | 25.95 | < 0.00001 |
| Date | 1 | 46 | 70.45 | < 0.00001 |
| Cover × Date | 5 | 46 | 3.32 | 0.012 |
| P concentration | ||||
| Cover | 5 | 46 | 21.26 | < 0.00001 |
| Date | 1 | 46 | 4.35 | 0.043 |
| Cover × Date | 5 | 46 | 0.66 | 0.66 |
| K concentration | ||||
| Cover | 5 | 46 | 19.37 | < 0.00001 |
| Date | 1 | 46 | 63.83 | < 0.00001 |
| Cover × Date | 5 | 46 | 0.62 | 0.69 |
| Mg concentration | ||||
| Cover | 5 | 46 | 13.06 | < 0.00001 |
| Date | 1 | 46 | 83.36 | < 0.00001 |
| Cover × Date | 5 | 46 | 3.78 | 0.0060 |
| Zn concentration | ||||
| Cover | 5 | 46 | 1.78 | 0.14 |
| Date | 1 | 46 | 17.78 | 0.00011 |
| Cover × Date | 5 | 46 | 1.34 | 0.26 |
| Ca concentration | ||||
| Cover | 5 | 46 | 5.03 | 0.00093 |
| Date | 1 | 46 | 43.40 | < 0.00001 |
| Cover × Date | 5 | 46 | 2.99 | 0.020 |
| SS% | ||||
| Cover | 5 | 46 | 60.18 | < 0.00001 |
| Date | 1 | 46 | 30.55 | < 0.00001 |
| Cover × Date | 5 | 46 | 19.39 | < 0.00001 |
Table 1 results: Means separation
For consistency we present means separation for all variables within each of the two dates. Here and elsewhere the means and 95% confidence interval endpoints are provided to three significant digits. The 95% confidence intervals are displayed in brackets after the mean and have been adjusted with the Sidak multiple comparison adjustment. The p-values on which the multiple comparison letters are based have been adjusted with the Tukey adjustment.
t1_fits <- t1_fits %>%
mutate(EMM = map(fit, ~ emmeans(., pairwise ~ Cover | Date, type = 'response', inv.lbl = 'emmean')),
CLD = map(EMM, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)))
t1_meantables <- with(t1_fits, map2_dfr(variable, CLD, ~ tibble(variable = .x, as_tibble(.y))))
t1_contrasttables <- with(t1_fits, map2_dfr(variable, EMM, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
rename(difference = estimate) %>%
relocate(ratio, .after = 'difference') %>%
mutate(null = if_else(is.na(null), 0, 1))
t1_meantables %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover, Date), names_from = variable, values_from = display) %>%
group_by(Date) %>%
arrange(Cover) %>%
gt() %>%
tabopts()| Cover | pH | P concentration | K concentration | Mg concentration | Zn concentration | Ca concentration | SS% |
|---|---|---|---|---|---|---|---|
| 417 | |||||||
| Control | 4.26 [4.08, 4.44] b | 117 [75.8, 180] b | 143 [114, 180] bc | 127 [100, 162] ab | 6.14 [5.05, 7.46] a | 769 [566, 1040] ab | 0.2 [0.168, 0.239] a |
| Bark | 4.64 [4.46, 4.82] a | 85.6 [55.6, 132] b | 175 [139, 221] bc | 166 [131, 211] a | 7.21 [5.93, 8.76] a | 953 [702, 1290] a | 0.1 [0.0838, 0.119] b |
| Brassica | 3.94 [3.76, 4.12] c | 238 [155, 367] a | 185 [147, 233] bc | 80 [62.9, 102] c | 6.13 [5.04, 7.45] a | 596 [439, 809] bc | 0.2 [0.168, 0.239] a |
| Grass | 3.94 [3.76, 4.12] c | 201 [130, 310] a | 272 [216, 343] a | 97.4 [76.6, 124] bc | 6.4 [5.27, 7.78] a | 652 [480, 885] bc | 0.235 [0.197, 0.281] a |
| Legume | 3.98 [3.8, 4.16] c | 224 [146, 346] a | 196 [156, 247] ab | 84.1 [66.1, 107] c | 5.67 [4.66, 6.89] a | 536 [395, 728] c | 0.217 [0.182, 0.259] a |
| Ridomil | 4.26 [4.08, 4.44] b | 115 [74.4, 177] b | 133 [106, 168] c | 125 [98, 158] ab | 5.93 [4.88, 7.21] a | 748 [551, 1020] abc | 0.2 [0.168, 0.239] a |
| 916 | |||||||
| Control | 4.6 [4.42, 4.78] ab | 125 [81.3, 193] c | 192 [153, 242] c | 162 [127, 206] ab | 6.77 [5.57, 8.23] a | 906 [667, 1230] a | 0.1 [0.0838, 0.119] c |
| Bark | 4.66 [4.48, 4.84] a | 126 [81.6, 194] bc | 249 [197, 313] bc | 185 [145, 235] a | 7.8 [6.42, 9.48] a | 1070 [791, 1460] a | 0.1 [0.0838, 0.119] c |
| Brassica | 4.3 [4.12, 4.48] c | 245 [159, 377] a | 267 [212, 336] abc | 127 [100, 162] b | 7.73 [6.36, 9.39] a | 830 [612, 1130] a | 0.174 [0.146, 0.208] b |
| Grass | 4.34 [4.16, 4.52] bc | 206 [134, 318] ab | 380 [302, 478] a | 162 [128, 206] ab | 6.97 [5.73, 8.47] a | 874 [644, 1190] a | 0.249 [0.209, 0.297] a |
| Legume | 4.24 [4.06, 4.42] c | 271 [176, 417] a | 344 [273, 433] ab | 169 [133, 214] ab | 7.99 [6.57, 9.71] a | 1060 [780, 1440] a | 0.277 [0.232, 0.33] a |
| Ridomil | 4.78 [4.6, 4.96] a | 134 [87.2, 207] bc | 201 [159, 253] c | 182 [143, 231] a | 6.65 [5.48, 8.09] a | 977 [720, 1330] a | 0.1 [0.0838, 0.119] c |
Table 2 analysis
In this table, we have repeated measurements within plot. Thus, we must account for measurements taken on the same experimental unit over time not being independent of one another.
Convert the data into number of months since first measurement, with May 2016 being month 0.
table2 <- table2 %>%
mutate(
Date = fct_recode(Date, `2021.10` = '2021.1') %>% paste0('.01') %>% as.Date(format = '%Y.%m.%d'),
month = lubridate::month(Date) + 12*(lubridate::year(Date) - 2016) - 5,
Plot_Rep = interaction(Plot, Rep)
)Plot the data by month so we can see how the time trend should be modeled. Plot each individual plot’s time trend as a thin line, and the median for each treatment as a thick line.
ptop <- ggplot(table2, aes(x = Date, y = pH, color = Cover, group = interaction(Cover, Rep, Plot))) +
geom_line(alpha = .4, linewidth = .3) +
stat_summary(aes(group = Cover), geom = 'line', fun = 'median', linewidth = 1.5) +
see::scale_color_okabeito(name = 'Soil amendment\ntreatment') +
scale_x_continuous(breaks = unique(table2$Date), labels = substr(as.character(unique(table2$Date)), 1, 7)) +
theme(legend.position = 'inside', legend.position.inside = c(1, 1), legend.justification = c(1, 1), legend.background = element_blank(),
axis.text.x = element_blank(), axis.title.x = element_blank())
pbot <- ggplot(table2, aes(x = Date, y = EC, color = Cover, group = interaction(Cover, Rep, Plot))) +
geom_line(alpha = .4, linewidth = .3) +
stat_summary(aes(group = Cover), geom = 'line', fun = 'median', linewidth = 1.5) +
see::scale_color_okabeito() +
scale_x_continuous(breaks = unique(table2$Date), labels = substr(as.character(unique(table2$Date)), 1, 7)) +
labs(y = 'Electrical conductivity (mS cm<sup>-1</sup>)') +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = 'none', axis.title.y = ggtext::element_markdown())
(fig1 <- ptop / pbot + plot_annotation(tag_levels = 'A'))
# ggsave('fig1.png', fig1, height = 8, width = 4, dpi = 300)
# ggsave('fig1.jpg', fig1, height = 8, width = 4, dpi = 400)Judging from the figures above, the time trends are similar across cover treatments but they are not linear so we will fit a smoothing spline that is allowed to vary by treatment. That will allow us to make predictions to compare treatments on specific dates or averaged over time. The pH variable is assumed to have normally distributed residuals but the EC variable is log-transformed before fitting.
This is a generalized additive mixed model with cover as the fixed effect and a separate smoothing term for time fit to each cover time (allowing for an interaction between cover and time). Replicate and plot within replicate have random intercepts. Thus the repeated measures are accounted for by the time splines rather than a residual covariance structure.
t2_pH_fit <- gamm(pH ~ Cover + s(month, bs = 'tp', k = 5, by = Cover),
random = list(Plot = ~ 1, Plot_Rep = ~ 1), data = table2)
t2_EC_fit <- gamm(log(EC) ~ Cover + s(month, bs = 'tp', k = 5, by = Cover),
random = list(Plot = ~ 1, Plot_Rep = ~ 1), data = table2)Table 2 results: means separation
Because this is a specialized model, it is not straightforward to calculate an ANOVA-like table. Therefore we will proceed directly to estimating the marginal means for cover types within each time point and comparing them with Tukey tests. Separate Tukey adjustments are made within each time point. The procedure is similar to that described for Table 1 above.
t2_pH_EMM <- emmeans(t2_pH_fit, pairwise ~ Cover | month, at = list(month = unique(table2$month)), data = table2)
t2_pH_CLD <- cld(t2_pH_EMM$emmeans, Letters = letters, adjust = 'tukey', decreasing = TRUE)
#Transformation must be specified for EC as the transformation is not detected by emmeans
t2_EC_grid <- update(ref_grid(t2_EC_fit, data = table2, at = list(month = unique(table2$month))), tran = 'log')
t2_EC_EMM <- emmeans(t2_EC_grid, pairwise ~ Cover | month, type = 'response', inv.lbl = 'emmean')
t2_EC_CLD <- cld(t2_EC_EMM$emmeans, Letters = letters, adjust = 'tukey', decreasing = TRUE)
t2_meantables <- bind_rows(
tibble(variable = 'pH', as_tibble(t2_pH_CLD)),
tibble(variable = 'EC', as_tibble(t2_EC_CLD))
) %>%
mutate(month = table2$Date[match(month, table2$month)]) %>%
rename(Date = month)
t2_contrasttables <- bind_rows(
tibble(variable = 'pH', as_tibble(t2_pH_EMM$contrasts)),
tibble(variable = 'EC', as_tibble(t2_EC_EMM$contrasts) %>% mutate(month = as.numeric(as.character(month))))
) %>%
rename(difference = estimate) %>%
relocate(ratio, .after = 'difference') %>%
mutate(null = if_else(is.na(null), 0, 1)) %>%
mutate(month = table2$Date[match(month, table2$month)]) %>%
rename(Date = month)
t2_meantables %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover, Date), names_from = variable, values_from = display) %>%
group_by(Date) %>%
arrange(Cover) %>%
gt() %>%
tabopts()| Cover | pH | EC |
|---|---|---|
| 2016-05-01 | ||
| Control | 5.35 [5, 5.7] a | 0.516 [0.328, 0.811] ab |
| Bark | 5.52 [5.17, 5.87] a | 0.436 [0.277, 0.685] b |
| Brassica | 5.31 [4.96, 5.66] a | 0.828 [0.527, 1.3] ab |
| Grass | 5.42 [5.06, 5.77] a | 0.887 [0.565, 1.39] a |
| Legume | 5.32 [4.97, 5.67] a | 0.805 [0.512, 1.26] ab |
| Ridomil | 5.73 [5.37, 6.08] a | 0.591 [0.376, 0.929] ab |
| 2017-09-01 | ||
| Control | 4.37 [4.05, 4.69] ab | 0.134 [0.0861, 0.207] b |
| Bark | 4.62 [4.32, 4.92] a | 0.0411 [0.0265, 0.0639] c |
| Brassica | 4.17 [3.84, 4.5] ab | 0.192 [0.124, 0.298] ab |
| Grass | 4.03 [3.7, 4.37] b | 0.285 [0.184, 0.441] a |
| Legume | 4.32 [4.01, 4.64] ab | 0.218 [0.14, 0.337] ab |
| Ridomil | 4.37 [4.04, 4.71] ab | 0.142 [0.0913, 0.22] b |
| 2019-02-01 | ||
| Control | 4.44 [4.19, 4.69] a | 0.175 [0.126, 0.241] a |
| Bark | 4.23 [3.99, 4.47] a | 0.163 [0.118, 0.226] a |
| Brassica | 4.24 [3.99, 4.49] a | 0.195 [0.141, 0.27] a |
| Grass | 4.27 [4.02, 4.52] a | 0.19 [0.137, 0.263] a |
| Legume | 4.23 [3.98, 4.48] a | 0.199 [0.144, 0.275] a |
| Ridomil | 4.4 [4.14, 4.65] a | 0.176 [0.127, 0.243] a |
| 2019-09-01 | ||
| Control | 4.59 [4.33, 4.85] a | 0.134 [0.0955, 0.189] a |
| Bark | 4.26 [4.01, 4.51] a | 0.187 [0.133, 0.263] a |
| Brassica | 4.38 [4.12, 4.65] a | 0.141 [0.1, 0.198] a |
| Grass | 4.49 [4.23, 4.75] a | 0.122 [0.087, 0.172] a |
| Legume | 4.33 [4.08, 4.59] a | 0.138 [0.0981, 0.194] a |
| Ridomil | 4.56 [4.3, 4.82] a | 0.132 [0.0941, 0.186] a |
| 2021-05-01 | ||
| Control | 4.72 [4.45, 4.99] a | 0.0662 [0.0456, 0.0961] a |
| Bark | 4.53 [4.28, 4.79] a | 0.0843 [0.0579, 0.123] a |
| Brassica | 4.51 [4.23, 4.78] a | 0.0723 [0.0498, 0.105] a |
| Grass | 4.59 [4.3, 4.87] a | 0.0705 [0.0486, 0.102] a |
| Legume | 4.52 [4.25, 4.78] a | 0.0704 [0.0485, 0.102] a |
| Ridomil | 4.71 [4.43, 4.99] a | 0.0632 [0.0435, 0.0917] a |
| 2021-10-01 | ||
| Control | 4.72 [4.51, 4.93] a | 0.14 [0.107, 0.185] a |
| Bark | 4.57 [4.36, 4.78] a | 0.18 [0.137, 0.237] a |
| Brassica | 4.55 [4.34, 4.76] a | 0.155 [0.118, 0.205] a |
| Grass | 4.59 [4.37, 4.8] a | 0.148 [0.113, 0.195] a |
| Legume | 4.53 [4.32, 4.74] a | 0.15 [0.114, 0.197] a |
| Ridomil | 4.72 [4.51, 4.94] a | 0.136 [0.103, 0.178] a |
| 2022-04-01 | ||
| Control | 4.72 [4.4, 5.04] a | 0.561 [0.366, 0.86] a |
| Bark | 4.6 [4.29, 4.9] a | 0.793 [0.516, 1.22] a |
| Brassica | 4.63 [4.3, 4.95] a | 0.624 [0.407, 0.957] a |
| Grass | 4.61 [4.28, 4.94] a | 0.553 [0.361, 0.847] a |
| Legume | 4.54 [4.22, 4.86] a | 0.587 [0.383, 0.9] a |
| Ridomil | 4.76 [4.43, 5.08] a | 0.554 [0.361, 0.849] a |
Table 3 analysis
For all the statistical models fit here, we set up the model to reflect the split-plot design. Cover is the main plot effect and cultivar is the subplot effect. Cover, cultivar, and their interaction are treated as fixed effects. Replicate and cover nested within replicate (i.e. main plot) are treated as random intercepts.
First, fit binomial mixed models to the binary variables (whether the plant was alive in 9/21 and 6/22). Alive is coded as 1 and dead is coded as 0, so the marginal estimates are the probabilities of being alive conditional on the treatments.
Because there is complete separation in the data (some combinations of cover and cultivar have 100% alive individuals), we include regularizing priors on the fixed effects in the binomial GLMMs to mitigate that problem (prior distributions on the fixed effects are means with normal 0 and standard deviation 5).
t3_live921_fit <- glmmTMB(Live921 ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = binomial(link = 'logit'), data = table3,
priors = data.frame(prior = c('normal(0, 5)'), class = 'fixef'))
t3_live622_fit <- glmmTMB(Live622 ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = binomial(link = 'logit'), data = table3,
priors = data.frame(prior = c('normal(0, 5)'), class = 'fixef'))Model diagnostics look acceptable (not shown).
plot(simulateResiduals(t3_live921_fit))
plot(simulateResiduals(t3_live622_fit))The remaining variables are ordered categorical variables. To fit these, we will fit cumulative logistic mixed models. The response distribution is ordered multinomial with a cumulative logit link function. This models the probability of each of the levels of the score categories, conditional on the treatments and the random effects.
Because of the small dataset where some predictor combinations have little to no variation in the response, we cannot fit this model with “classical” frequentist statistics. We must use a Bayesian approach. We put normal prior distributions with mean 0 and standard deviation 5 on the fixed effects.
t3_leaf_fit <- brm(Leaf ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 114,
control = list(adapt_delta = 0.9),
file = 'table3_leaf_fit')
t3_flower_fit <- brm(Crop ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 115,
control = list(adapt_delta = 0.9),
file = 'table3_flower_fit')
t3_crop_fit <- brm(Crop ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 116,
control = list(adapt_delta = 0.95),
file = 'table3_crop_fit')
t3_ripe_fit <- brm(Ripe ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 117,
control = list(adapt_delta = 0.97),
file = 'table3_ripe_fit')
t3_vig921_fit <- brm(Vig921 ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 118,
control = list(adapt_delta = 0.95),
file = 'table3_vig921_fit')
t3_vig622_fit <- brm(Vig622 ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = cumulative(link = 'logit', threshold = 'flexible'),
data = table3,
prior = c(
prior(normal(0, 5), class = b),
prior(gamma(1, 1), class = sd)
),
chains = 4, iter = 3000, warmup = 2000, seed = 119,
control = list(adapt_delta = 0.9),
file = 'table3_vig622_fit')Model diagnostics indicate convergence and good model fit (not shown).
map(list(t3_leaf_fit, t3_flower_fit, t3_crop_fit, t3_ripe_fit, t3_vig921_fit, t3_vig622_fit), ~ max(rhat(.), na.rm = TRUE))
map(list(t3_leaf_fit, t3_flower_fit, t3_crop_fit, t3_ripe_fit, t3_vig921_fit, t3_vig622_fit), ~ pp_check(., type = 'bars'))Table 3 results: ANOVA tables
We only produce analysis of deviance tables with joint chi-squared tests for the binomial models that modeled probability of the plant being alive. The Bayesian models do not have a meaningful ANOVA table; for those, we will proceed directly to means separation.
t3_freqfits <- tibble(
variable = c('Probability alive 9/21', 'Probability alive 6/22'),
fit = list(t3_live921_fit, t3_live622_fit)
) %>%
mutate(ANOVA = map(fit, joint_tests))
t3_ANOVAtables <- with(t3_freqfits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y)))) %>%
dplyr::select(-c(df2, F.ratio))
t3_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(c(Chisq), decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts %>%
fmt_missing| model term | df1 | Chisq | p.value |
|---|---|---|---|
| Probability alive 9/21 | |||
| Cover | 5 | 1.52 | 0.91 |
| Cultivar | 2 | 2.54 | 0.28 |
| Cover × Cultivar | 10 | 9.70 | 0.47 |
| Probability alive 6/22 | |||
| Cover | 5 | 4.79 | 0.44 |
| Cultivar | 2 | 6.35 | 0.042 |
| Cover × Cultivar | 10 | 12.61 | 0.25 |
Table 3 results: means separation
For the two models where we fit the probability that the plant is alive, we produce means separation tables based on p-values with Tukey post hoc adjustment. For the Bayesian models, we estimate marginal means using the expected value of the posterior predictive, obtaining a mean for each category using the mean of the posterior probabilities for each ordinal category weighted by the score values. Then we generate pairwise contrasts by taking the difference of this weighted mean between each pair of treatments. We produce means separation tables based on the Bayesian maximum a posteriori p-value (\(p_{MAP}\)). This value is analogous to a p-value; it is the ratio of evidence for the null (0 difference) to the highest-probability posterior value. We declare \(p_{MAP} < 0.05\) to indicate a statistically significant difference between two groups and generate the multiple comparison letters accordingly. No multiple comparison adjustment is done because the regularizing priors have a similar effect to that adjustment.
### Frequentist models: two binomial GLMMs
t3_freqfits <- t3_freqfits %>%
mutate(EMM_cover = map(fit, ~ emmeans(., pairwise ~ Cover, type = 'response', inv.lbl = 'emmean')),
CLD_cover = map(EMM_cover, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_cultivar = map(fit, ~ emmeans(., pairwise ~ Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_cultivar = map(EMM_cultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_coverbycultivar = map(fit, ~ emmeans(., pairwise ~ Cover | Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_coverbycultivar = map(EMM_coverbycultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)))
### Bayesian models: 6ordered multinomial GLMMs
# Define function to get mean.class mode predictions
mean_class_pred <- function(fit, pred_grid) {
add_epred_draws(pred_grid, fit, re_formula = ~ 0) %>%
mutate(.category = as.numeric(levels(.category))[.category]) %>%
group_by(Cover, Cultivar, .draw) %>%
summarize(mean_class = weighted.mean(x = .category, w = .epred))
}
# Define function to get multiple comparison letters
bayes_multcompletters <- function(post_cons, post_means, group_col, mean_col = 'emmean', alpha = 0.05, decreasing_letters = FALSE) {
different = with(post_cons, setNames(p_MAP < alpha, contrast))
levels_ordered = as.character(post_means[[group_col]][order(post_means[[mean_col]])])
letter_vec <- multcomp:::insert_absorb(different,
decreasing = decreasing_letters,
separator = ' - ',
lvl_order = levels_ordered)
data.frame(group = names(letter_vec$Letters), letter = letter_vec$Letters) %>% setNames(c(group_col, 'letter'))
}
# Define function to get pairwise differences between mean class predictions
get_pairwise_differences <- function(dt, byvar, xvar, groupvars, op = `-`, reverse = FALSE) {
grps <- unique(dt[[byvar]])
combs <- combn(1:length(grps), 2)
if (reverse) combs <- combs[2:1, ]
dt %>%
group_by(across(all_of(groupvars))) %>%
arrange(!!sym(byvar)) %>%
reframe(grp1 = grps[combs[1, ]], grp2 = grps[combs[2, ]], mean1 = pick(xvar)[[xvar]][combs[1, ]], mean2 = pick(xvar)[[xvar]][combs[2, ]]) %>%
mutate(diff = op(mean1, mean2))
}
t3_predgrid <- table3 %>%
dplyr::select(Cover, Cultivar) %>%
distinct
# Get all draws for the estimated marginal means and contrasts
t3_bayesfits <- tibble(
variable = c('Leaf development', 'Flower development', 'Crop load', 'Fruit ripeness', 'Vigor rating 9/21', 'Vigor rating 6/22'),
fit = list(t3_leaf_fit, t3_flower_fit, t3_crop_fit, t3_ripe_fit, t3_vig921_fit, t3_vig622_fit)
)
t3_bayesfits <- t3_bayesfits %>%
mutate(EMM_draws_coverbycultivar = map(fit, mean_class_pred, pred_grid = t3_predgrid),
EMM_draws_cover = map(EMM_draws_coverbycultivar, ~ group_by(., Cover, .draw) %>% summarize(mean_class = mean(mean_class))),
EMM_draws_cultivar = map(EMM_draws_coverbycultivar, ~ group_by(., Cultivar, .draw) %>% summarize(mean_class = mean(mean_class))),
CON_draws_coverbycultivar = map(EMM_draws_coverbycultivar, ~ get_pairwise_differences(., byvar = 'Cover', xvar = 'mean_class', groupvars = c('.draw', 'Cultivar'))),
CON_draws_cover = map(EMM_draws_cover, ~ get_pairwise_differences(., byvar = 'Cover', xvar = 'mean_class', groupvars = c('.draw'))),
CON_draws_cultivar = map(EMM_draws_cultivar, ~ get_pairwise_differences(., byvar = 'Cultivar', xvar = 'mean_class', groupvars = c('.draw')))
)
# Summarize the draws with mean and 95% CI for the EMMs, mean and 95% CI for the CONs, and pMAP value for the CONs
t3_bayesfits <- t3_bayesfits %>%
mutate(
EMM_coverbycultivar = map(EMM_draws_coverbycultivar, ~ summarize(., describe_posterior(mean_class, ci_method = 'eti', ci = 0.95, test = 'none'))),
EMM_cover = map(EMM_draws_cover, ~ summarize(., describe_posterior(mean_class, ci_method = 'eti', ci = 0.95, test = 'none'))),
EMM_cultivar = map(EMM_draws_cultivar, ~ summarize(., describe_posterior(mean_class, ci_method = 'eti', ci = 0.95, test = 'none'))),
CON_coverbycultivar = map(CON_draws_coverbycultivar, ~ group_by(., Cultivar, grp1, grp2) %>% summarize(describe_posterior(diff, ci_method = 'eti', ci = 0.95, test = 'p_map', null = 0)) %>% mutate(contrast = paste(grp1, grp2, sep = ' - '))),
CON_cover = map(CON_draws_cover, ~ group_by(., grp1, grp2) %>% summarize(describe_posterior(diff, ci_method = 'eti', ci = 0.95, test = 'p_map', null = 0)) %>% mutate(contrast = paste(grp1, grp2, sep = ' - '))),
CON_cultivar = map(CON_draws_cultivar, ~ group_by(., grp1, grp2) %>% summarize(describe_posterior(diff, ci_method = 'eti', ci = 0.95, test = 'p_map', null = 0)) %>% mutate(contrast = paste(grp1, grp2, sep = ' - ')))
)
# Generate multiple comparison letters based on pMAP < 0.05 criterion
t3_bayesfits <- t3_bayesfits %>%
mutate(
CLD_coverbycultivar = map2(EMM_coverbycultivar, CON_coverbycultivar, ~ map_dfr(levels(table3$Cultivar), function(cult) left_join(filter(.x, Cultivar == cult), bayes_multcompletters(filter(.y, Cultivar == cult), filter(.x, Cultivar == cult), group_col = 'Cover', mean_col = 'Median', decreasing_letters = TRUE)))),
CLD_cover = map2(EMM_cover, CON_cover, ~ left_join(.x, bayes_multcompletters(.y, .x, group_col = 'Cover', mean_col = 'Median', decreasing_letters = TRUE))),
CLD_cultivar = map2(EMM_cultivar, CON_cultivar, ~ left_join(.x, bayes_multcompletters(.y, .x, group_col = 'Cultivar', mean_col = 'Median', decreasing_letters = TRUE)))
)
# Combine all CLDs into one data frame and get them ready for printing and exporting
t3_allclds <- bind_rows(
t3_freqfits %>% dplyr::select(variable, CLD_cover, CLD_cultivar, CLD_coverbycultivar),
t3_bayesfits %>% dplyr::select(variable, CLD_cover, CLD_cultivar, CLD_coverbycultivar)
)
t3_meantables_cover <- with(t3_allclds, map2_dfr(variable, CLD_cover, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of(c('asymp.LCL', 'CI_low')), upper.CL = any_of(c('asymp.UCL', 'CI_high')), emmean = any_of('Median'), .group = any_of('letter')))) %>%
dplyr::select(-c(Parameter, CI, df)) %>%
mutate(Cover = factor(Cover, levels = levels(table3$Cover)))
t3_meantables_cultivar <- with(t3_allclds, map2_dfr(variable, CLD_cultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of(c('asymp.LCL', 'CI_low')), upper.CL = any_of(c('asymp.UCL', 'CI_high')), emmean = any_of('Median'), .group = any_of('letter')))) %>%
dplyr::select(-c(Parameter, CI, df))
t3_meantables_coverbycultivar <- with(t3_allclds, map2_dfr(variable, CLD_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of(c('asymp.LCL', 'CI_low')), upper.CL = any_of(c('asymp.UCL', 'CI_high')), emmean = any_of('Median'), .group = any_of('letter')))) %>%
dplyr::select(-c(Parameter, CI, df)) %>%
mutate(Cover = factor(Cover, levels = levels(table3$Cover)))
# Combine all contrasts into one data frame and get them ready for exporting
t3_freqcontrasttables_cover <- with(t3_freqfits, map2_dfr(variable, EMM_cover, ~ tibble(variable = .x, as_tibble(.y$contrasts))))
t3_freqcontrasttables_cultivar <- with(t3_freqfits, map2_dfr(variable, EMM_cultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts))))
t3_freqcontrasttables_coverbycultivar <- with(t3_freqfits, map2_dfr(variable, EMM_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts))))
t3_bayescontrasttables_cover <- with(t3_bayesfits, map2_dfr(variable, CON_cover, ~ tibble(variable = .x, as_tibble(.y)))) %>%
rename(difference = Median) %>% mutate(null = 0)
t3_bayescontrasttables_cultivar <- with(t3_bayesfits, map2_dfr(variable, CON_cultivar, ~ tibble(variable = .x, as_tibble(.y)))) %>%
rename(difference = Median) %>% mutate(null = 0)
t3_bayescontrasttables_coverbycultivar <- with(t3_bayesfits, map2_dfr(variable, CON_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y)))) %>%
rename(difference = Median) %>% mutate(null = 0)
t3_contrasttables_cover <- bind_rows(t3_freqcontrasttables_cover, t3_bayescontrasttables_cover) %>%
dplyr::select(variable, contrast, odds.ratio, difference, SE, CI_low, CI_high, z.ratio, p.value, p_MAP, null)
t3_contrasttables_cultivar <- bind_rows(t3_freqcontrasttables_cultivar, t3_bayescontrasttables_cultivar) %>%
dplyr::select(variable, contrast, odds.ratio, difference, SE, CI_low, CI_high, z.ratio, p.value, p_MAP, null)
t3_contrasttables_coverbycultivar <- bind_rows(t3_freqcontrasttables_coverbycultivar, t3_bayescontrasttables_coverbycultivar) %>%
dplyr::select(variable, Cultivar, contrast, odds.ratio, difference, SE, CI_low, CI_high, z.ratio, p.value, p_MAP, null)
t3_meantables_cover %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover), names_from = variable, values_from = display) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover') %>%
tabopts()| Cover | Probability alive 9/21 | Probability alive 6/22 | Leaf development | Flower development | Crop load | Fruit ripeness | Vigor rating 9/21 | Vigor rating 6/22 |
|---|---|---|---|---|---|---|---|---|
| Control | 0.844 [0.487, 0.969] a | 0.877 [0.531, 0.978] a | 1.16 [0.635, 1.71] b | 1.46 [0.751, 2.46] a | 1.46 [0.731, 2.41] a | 1.23 [0.646, 1.91] a | 2.09 [0.943, 3.17] a | 1.59 [1.21, 2.03] b |
| Bark | 0.802 [0.435, 0.955] a | 0.997 [0.256, 1] a | 3.12 [2.65, 3.54] a | 1.49 [0.712, 2.49] a | 1.5 [0.735, 2.54] a | 1.52 [0.804, 2.28] a | 1.73 [0.488, 2.96] a | 4.33 [4, 4.58] a |
| Brassica | 0.812 [0.436, 0.96] a | 0.932 [0.532, 0.994] a | 1.04 [0.557, 1.63] b | 1.97 [1.08, 3.03] a | 1.97 [1.08, 3.11] a | 1.48 [0.919, 2.14] a | 2.14 [0.841, 3.29] a | 1.7 [1.27, 2.15] b |
| Grass | 0.881 [0.521, 0.981] a | 0.926 [0.504, 0.994] a | 0.884 [0.449, 1.46] b | 1.92 [1.01, 2.98] a | 1.92 [1.02, 3.03] a | 1.37 [0.729, 2.02] a | 1.98 [0.681, 3.13] a | 1.65 [1.23, 2.12] b |
| Legume | 0.929 [0.456, 0.995] a | 0.783 [0.45, 0.941] a | 1.27 [0.708, 1.98] b | 0.801 [0.298, 1.62] a | 0.796 [0.287, 1.63] a | 0.7 [0.255, 1.31] a | 1.89 [0.674, 3.04] a | 1.78 [1.33, 2.22] b |
| Ridomil | 0.887 [0.515, 0.983] a | 0.915 [0.474, 0.992] a | 1.45 [0.916, 2.02] b | 1 [0.393, 1.94] a | 1 [0.381, 1.96] a | 0.923 [0.419, 1.56] a | 3.36 [2.22, 3.93] a | 1.82 [1.43, 2.28] b |
t3_meantables_cultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cultivar), names_from = variable, values_from = display) %>%
arrange(Cultivar) %>%
gt(caption = 'Means separation: cultivar') %>%
tabopts()| Cultivar | Probability alive 9/21 | Probability alive 6/22 | Leaf development | Flower development | Crop load | Fruit ripeness | Vigor rating 9/21 | Vigor rating 6/22 |
|---|---|---|---|---|---|---|---|---|
| Gumbo | 0.797 [0.586, 0.915] a | 0.989 [0.836, 0.999] a | 1.8 [1.28, 2.36] a | 1.64 [1.15, 2.34] a | 1.65 [1.12, 2.39] a | 1.24 [0.863, 1.65] a | 1.91 [1.36, 2.46] b | 2.92 [2.62, 3.24] a |
| Legacy | 0.89 [0.711, 0.964] a | 0.922 [0.609, 0.989] ab | 1.91 [1.48, 2.38] a | 1.29 [0.798, 1.94] a | 1.28 [0.81, 1.95] a | 1.21 [0.804, 1.7] a | 2.35 [1.8, 2.89] a | 1.91 [1.62, 2.22] b |
| Pearl | 0.895 [0.665, 0.973] a | 0.8 [0.48, 0.945] b | 0.76 [0.475, 1.14] b | 1.42 [0.962, 2.05] a | 1.43 [0.982, 2.07] a | 1.17 [0.86, 1.59] a | 2.27 [1.67, 2.82] ab | 1.6 [1.34, 1.91] b |
t3_meantables_coverbycultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover, Cultivar), names_from = variable, values_from = display) %>%
group_by(Cultivar) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover within each cultivar') %>%
tabopts()| Cover | Probability alive 9/21 | Probability alive 6/22 | Leaf development | Flower development | Crop load | Fruit ripeness | Vigor rating 9/21 | Vigor rating 6/22 |
|---|---|---|---|---|---|---|---|---|
| Gumbo | ||||||||
| Control | 0.815 [0.343, 0.974] a | 0.98 [0.563, 0.999] a | 1.35 [0.596, 2.22] b | 1.65 [0.734, 2.85] a | 1.66 [0.758, 2.85] a | 1.29 [0.574, 2.01] a | 1.84 [0.636, 3] ab | 2.32 [1.73, 2.94] b |
| Bark | 0.733 [0.236, 0.961] a | 0.999 [0.0991, 1] a | 2.88 [2.03, 3.56] a | 1.96 [0.73, 3.4] a | 1.99 [0.791, 3.49] a | 1.1 [0.289, 1.98] a | 1.74 [0.394, 2.97] ab | 4.64 [4.26, 4.87] a |
| Brassica | 0.822 [0.303, 0.98] a | 0.993 [0.203, 1] a | 1.8 [0.847, 2.72] ab | 1.75 [0.758, 2.97] a | 1.74 [0.76, 2.97] a | 1.51 [0.767, 2.24] a | 2.01 [0.603, 3.25] ab | 2.76 [2.11, 3.4] b |
| Grass | 0.744 [0.239, 0.964] a | 0.992 [0.205, 1] a | 1.68 [0.785, 2.58] ab | 1.46 [0.586, 2.7] a | 1.46 [0.54, 2.79] a | 1.14 [0.426, 1.9] a | 1.33 [0.196, 2.67] ab | 2.41 [1.77, 3.08] b |
| Legume | 0.648 [0.182, 0.938] a | 0.841 [0.343, 0.982] a | 1.59 [0.761, 2.48] b | 1.42 [0.492, 2.77] a | 1.43 [0.462, 2.82] a | 1.03 [0.34, 1.85] a | 1.18 [0.16, 2.53] b | 2.51 [1.84, 3.2] b |
| Ridomil | 0.923 [0.375, 0.996] a | 0.992 [0.211, 1] a | 1.53 [0.702, 2.44] b | 1.48 [0.487, 2.9] a | 1.49 [0.459, 2.97] a | 1.35 [0.517, 2.19] a | 3.48 [2.19, 3.97] a | 2.89 [2.2, 3.59] b |
| Legacy | ||||||||
| Control | 0.761 [0.259, 0.967] a | 0.906 [0.401, 0.993] a | 1.9 [0.913, 2.87] b | 0.777 [0.155, 2.02] a | 0.767 [0.172, 1.99] a | 1.12 [0.296, 2.07] a | 1.95 [0.624, 3.13] a | 1.78 [1.15, 2.45] b |
| Bark | 0.828 [0.296, 0.982] a | 0.998 [0.00447, 1] a | 3.99 [3.75, 4] a | 0.662 [0.0617, 2.29] a | 0.648 [0.0596, 2.27] a | 1.78 [0.512, 3.2] a | 1.73 [0.389, 2.99] a | 4.06 [3.49, 4.45] a |
| Brassica | 0.914 [0.352, 0.995] a | 0.82 [0.322, 0.978] a | 1.16 [0.375, 2.24] bc | 1.89 [0.783, 3.27] a | 1.9 [0.713, 3.38] a | 0.886 [0.275, 1.7] a | 2.59 [1.11, 3.68] a | 1.21 [0.653, 1.81] b |
| Grass | 0.92 [0.358, 0.996] a | 0.618 [0.201, 0.912] a | 0.33 [0.0517, 1.12] c | 2.65 [0.973, 4.26] a | 2.62 [1.03, 4.26] a | 1.75 [0.708, 2.8] a | 2.31 [0.842, 3.48] a | 0.865 [0.357, 1.53] b |
| Legume | 0.917 [0.358, 0.995] a | 0.902 [0.36, 0.993] a | 1.49 [0.604, 2.53] bc | 0.688 [0.114, 1.92] a | 0.673 [0.117, 1.92] a | 0.796 [0.158, 1.74] a | 2.14 [0.756, 3.33] a | 1.68 [1.08, 2.3] b |
| Ridomil | 0.931 [0.351, 0.997] a | 0.913 [0.36, 0.995] a | 2.54 [1.47, 3.41] b | 0.775 [0.155, 2.09] a | 0.774 [0.151, 2.09] a | 0.885 [0.196, 1.86] a | 3.65 [2.42, 3.98] a | 1.87 [1.28, 2.49] b |
| Pearl | ||||||||
| Control | 0.919 [0.408, 0.995] a | 0.439 [0.117, 0.822] a | 0.163 [0.0161, 0.744] b | 1.9 [0.828, 3.23] a | 1.89 [0.812, 3.2] a | 1.31 [0.536, 2.1] ac | 2.53 [1.25, 3.58] a | 0.648 [0.226, 1.28] b |
| Bark | 0.833 [0.3, 0.983] a | 0.993 [0.0408, 1] a | 2.55 [1.65, 3.28] a | 1.74 [0.723, 2.98] a | 1.76 [0.773, 3.01] a | 1.67 [0.866, 2.45] ab | 1.76 [0.385, 3.05] a | 4.32 [3.83, 4.66] a |
| Brassica | 0.622 [0.168, 0.931] a | 0.81 [0.317, 0.975] a | 0.114 [0.0047, 0.628] b | 2.27 [0.91, 3.78] a | 2.26 [0.943, 3.87] a | 2.03 [1.25, 3.03] a | 1.83 [0.415, 3.22] a | 1.12 [0.576, 1.73] b |
| Grass | 0.924 [0.36, 0.996] a | 0.904 [0.361, 0.994] a | 0.579 [0.154, 1.46] b | 1.63 [0.615, 2.95] ab | 1.62 [0.649, 3.01] ab | 1.23 [0.493, 1.97] ac | 2.28 [0.766, 3.44] a | 1.65 [1.02, 2.34] b |
| Legume | 0.991 [0.124, 1] a | 0.49 [0.135, 0.856] a | 0.694 [0.172, 1.7] b | 0.183 [0.0109, 0.968] b | 0.184 [0.0108, 1] b | 0.191 [0.0112, 0.967] c | 2.3 [0.919, 3.46] a | 1.14 [0.462, 1.89] b |
| Ridomil | 0.747 [0.225, 0.968] a | 0.502 [0.14, 0.862] a | 0.247 [0.0326, 0.92] b | 0.666 [0.146, 1.77] ab | 0.656 [0.126, 1.85] ab | 0.487 [0.0847, 1.33] bc | 3.05 [1.59, 3.9] a | 0.703 [0.257, 1.33] b |
Table 4 analysis
Here, we once again follow the split plot design with the random effects. Stem count is modeled with a Poisson distribution. Percentage live stems is modeled with a beta distribution. Because the beta distribution does not easily accommodate values exactly 1, we adjust the 100% alive values to 99.9%. We include a zero-inflation parameter in the model to model the probability of 0% alive because there are a large number of zeros in the data for that variable. The zero-inflation parameter has fixed effects for the treatments and their interaction. A normal prior distribution with mean 0 and standard deviation 5 is assumed for the fixed effects on the zero inflation parameter, because there is complete separation where no individuals in the bark treatment have a zero. The prior is weakly informative and regularizes the fixed effect coefficients.
table4 <- mutate(table4, Alive = Alive/100)
table4 <- mutate(table4, Alive_adj = if_else(Alive == 1, .999, Alive))
t4_stemcount_fit <- glmmTMB(Shoots ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = poisson(link = 'log'), data = table4)
t4_alive_fit <- glmmTMB(Alive_adj ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
ziformula = ~ Cover * Cultivar,
family = beta_family(link = 'logit'), data = table4,
priors = data.frame(prior = c('normal(0, 5)'), class = 'fixef_zi'))The variables total stem length, cane diameter, and stem weight are positive continuous variables. They only contain a few zeros (5), presumably for dead plants. It is appropriate to exclude that small number of zeros from the analysis. So, the means that we present here are conditional on the plant having survived.
table4 <- table4 %>%
mutate(across(c(`Total length`, `stem dia`, wt), ~ if_else(. == 0, as.numeric(NA), .)))
t4_totallength_fit <- glmmTMB(`Total length` ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = Gamma(link = 'log'), data = table4)
t4_stemdia_fit <- glmmTMB(`stem dia` ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = Gamma(link = 'log'), data = table4)
t4_wt_fit <- glmmTMB(wt ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
family = Gamma(link = 'log'), data = table4)The live length variable is different because it has a large number of zeros (35). This is enough to measurably impact the results. Therefore we fit this variable with a hurdle-gamma distribution, which is a gamma distribution with an extra zero component. The probability of observing zero will be accounted for in the final mean estimate. The hurdle parameter has fixed effects of treatments and their interaction on it. A normal prior distribution with mean 0 and standard deviation 5 is assumed for the fixed effects on the zero inflation parameter, because there is complete separation where no individuals in the bark treatment have a zero. The prior is weakly informative and regularizes the fixed effect coefficients.
t4_livelength_fit <- glmmTMB(`Live length` ~ Cover * Cultivar + (1 | Rep) + (1 | Cover:Rep),
ziformula = ~ Cover * Cultivar,
family = ziGamma(link = 'log'), data = table4,
priors = data.frame(prior = c('normal(0, 5)'), class = 'fixef_zi'))Model diagnostics look acceptable (not shown).
plot(simulateResiduals(t4_stemcount_fit))
plot(simulateResiduals(t4_alive_fit))
plot(simulateResiduals(t4_totallength_fit))
plot(simulateResiduals(t4_stemdia_fit))
plot(simulateResiduals(t4_wt_fit))
plot(simulateResiduals(t4_livelength_fit))Table 4 results: ANOVA tables
Joint Chi-squared tests are done for all variables.
t4_fits <- tibble(variable = c('Shoot count', 'Total stem length (cm)', 'Live stem length (cm)', 'Cane diameter (mm)', 'Stem weight (g)', 'Proportion living stems'),
fit = list(t4_stemcount_fit, t4_totallength_fit, t4_livelength_fit, t4_stemdia_fit, t4_wt_fit, t4_alive_fit)) %>%
mutate(ANOVA = map(fit, joint_tests))
t4_ANOVAtables <- with(t4_fits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y)))) %>%
dplyr::select(-c(df2, F.ratio))
t4_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(c(Chisq), decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts %>%
fmt_missing| model term | df1 | Chisq | p.value |
|---|---|---|---|
| Shoot count | |||
| Cover | 5 | 17.39 | 0.0038 |
| Cultivar | 2 | 12.63 | 0.0018 |
| Cover × Cultivar | 10 | 6.32 | 0.79 |
| Total stem length (cm) | |||
| Cover | 5 | 43.66 | < 0.00001 |
| Cultivar | 2 | 33.48 | < 0.00001 |
| Cover × Cultivar | 10 | 18.69 | 0.044 |
| Live stem length (cm) | |||
| Cover | 5 | 87.69 | < 0.00001 |
| Cultivar | 2 | 39.56 | < 0.00001 |
| Cover × Cultivar | 10 | 33.46 | 0.00023 |
| Cane diameter (mm) | |||
| Cover | 5 | 61.93 | < 0.00001 |
| Cultivar | 2 | 28.98 | < 0.00001 |
| Cover × Cultivar | 10 | 17.06 | 0.073 |
| Stem weight (g) | |||
| Cover | 5 | 82.83 | < 0.00001 |
| Cultivar | 2 | 43.67 | < 0.00001 |
| Cover × Cultivar | 10 | 27.99 | 0.0018 |
| Proportion living stems | |||
| Cover | 5 | 82.59 | < 0.00001 |
| Cultivar | 2 | 26.40 | < 0.00001 |
| Cover × Cultivar | 10 | 30.25 | 0.00078 |
Table 4 results: means separation
Some of the responses have a significant interaction between cover and cultivar. Both main effects and means separation for soil amendment treatment within cultivar are presented.
t4_fits <- t4_fits %>%
mutate(EMM_cover = map(fit, ~ emmeans(., pairwise ~ Cover, type = 'response', inv.lbl = 'emmean')),
CLD_cover = map(EMM_cover, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_cultivar = map(fit, ~ emmeans(., pairwise ~ Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_cultivar = map(EMM_cultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_coverbycultivar = map(fit, ~ emmeans(., pairwise ~ Cover | Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_coverbycultivar = map(EMM_coverbycultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)))
t4_meantables_cover <- with(t4_fits, map2_dfr(variable, CLD_cover, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t4_meantables_cultivar <- with(t4_fits, map2_dfr(variable, CLD_cultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t4_meantables_coverbycultivar <- with(t4_fits, map2_dfr(variable, CLD_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t4_contrasttables_cover <- with(t4_fits, map2_dfr(variable, EMM_cover, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t4_contrasttables_cultivar <- with(t4_fits, map2_dfr(variable, EMM_cultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t4_contrasttables_coverbycultivar <- with(t4_fits, map2_dfr(variable, EMM_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t4_meantables_cover %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover), names_from = variable, values_from = display) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover') %>%
tabopts()| Cover | Shoot count | Total stem length (cm) | Live stem length (cm) | Cane diameter (mm) | Stem weight (g) | Proportion living stems |
|---|---|---|---|---|---|---|
| Control | 2.22 [1.6, 3.08] ab | 156 [120, 202] b | 74.5 [45.7, 122] bc | 18.5 [14.1, 24.2] b | 61.4 [37, 102] b | 0.484 [0.336, 0.635] c |
| Bark | 3.42 [2.63, 4.44] a | 294 [228, 379] a | 283 [195, 411] a | 41.9 [32.1, 54.5] a | 406 [246, 670] a | 0.908 [0.851, 0.945] a |
| Brassica | 2.12 [1.52, 2.96] b | 139 [107, 180] b | 61.5 [41.1, 92.1] c | 17.5 [13.4, 22.8] b | 59.1 [35.7, 97.9] b | 0.548 [0.41, 0.679] bc |
| Grass | 2.3 [1.67, 3.16] ab | 160 [124, 208] b | 81.9 [53.8, 125] bc | 21.5 [16.4, 28.1] b | 81.9 [49.3, 136] b | 0.558 [0.417, 0.69] bc |
| Legume | 2.12 [1.52, 2.95] b | 158 [122, 204] b | 127 [80.1, 201] b | 19.1 [14.6, 24.9] b | 92.3 [55.7, 153] b | 0.714 [0.568, 0.825] b |
| Ridomil | 1.96 [1.38, 2.79] b | 138 [106, 179] b | 62.4 [40, 97.3] c | 17.5 [13.3, 22.9] b | 64.5 [38.7, 107] b | 0.566 [0.412, 0.709] bc |
t4_meantables_cultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cultivar), names_from = variable, values_from = display) %>%
arrange(Cultivar) %>%
gt(caption = 'Means separation: cultivar') %>%
tabopts()| Cultivar | Shoot count | Total stem length (cm) | Live stem length (cm) | Cane diameter (mm) | Stem weight (g) | Proportion living stems |
|---|---|---|---|---|---|---|
| Gumbo | 2.93 [2.44, 3.51] a | 217 [187, 252] a | 158 [122, 204] a | 27.4 [23.3, 32.3] a | 143 [109, 189] a | 0.783 [0.715, 0.838] a |
| Legacy | 2.04 [1.64, 2.54] b | 143 [123, 165] b | 72 [54.7, 94.7] b | 18.4 [15.7, 21.7] b | 75.9 [57.7, 100] b | 0.591 [0.502, 0.673] b |
| Pearl | 2.07 [1.66, 2.58] b | 152 [131, 178] b | 80.7 [59.3, 110] b | 19.4 [16.4, 23] b | 77.8 [58.6, 103] b | 0.567 [0.461, 0.668] b |
t4_meantables_coverbycultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover, Cultivar), names_from = variable, values_from = display) %>%
group_by(Cultivar) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover within each cultivar') %>%
tabopts()| Cover | Shoot count | Total stem length (cm) | Live stem length (cm) | Cane diameter (mm) | Stem weight (g) | Proportion living stems |
|---|---|---|---|---|---|---|
| Gumbo | ||||||
| Control | 2.8 [1.7, 4.6] a | 187 [128, 275] b | 73.9 [40.6, 134] b | 21.5 [14.4, 31.9] b | 71.3 [37.6, 135] b | 0.47 [0.273, 0.677] c |
| Bark | 4 [2.64, 6.06] a | 377 [257, 553] a | 366 [203, 660] a | 51.8 [34.8, 77.1] a | 528 [279, 1000] a | 0.921 [0.825, 0.967] a |
| Brassica | 2.8 [1.7, 4.6] a | 188 [128, 276] b | 147 [81.3, 264] b | 23.3 [15.7, 34.6] b | 115 [60.5, 220] b | 0.821 [0.651, 0.918] ab |
| Grass | 2.8 [1.7, 4.6] a | 190 [129, 279] b | 121 [64.6, 225] b | 24.2 [16.3, 36] b | 110 [57.7, 210] b | 0.586 [0.363, 0.779] bc |
| Legume | 2.4 [1.4, 4.11] a | 193 [131, 286] b | 187 [91.3, 383] ab | 24.8 [16.5, 37.2] b | 143 [74.9, 272] b | 0.86 [0.678, 0.947] ab |
| Ridomil | 3 [1.86, 4.85] a | 213 [145, 314] ab | 174 [96.3, 313] ab | 27.4 [18.4, 40.8] b | 127 [66.3, 242] b | 0.841 [0.682, 0.929] ab |
| Legacy | ||||||
| Control | 2.3 [1.33, 3.98] a | 169 [115, 249] a | 92.5 [46.7, 183] ab | 20.5 [13.8, 30.7] ab | 86.4 [45.1, 165] b | 0.533 [0.315, 0.739] b |
| Bark | 2.7 [1.63, 4.48] a | 200 [136, 294] a | 199 [110, 359] a | 32.6 [21.9, 48.5] a | 289 [153, 548] a | 0.92 [0.823, 0.966] a |
| Brassica | 2 [1.11, 3.6] a | 126 [85.7, 185] a | 51.7 [26.7, 100] b | 15.9 [10.7, 23.7] b | 46.1 [24.3, 87.4] b | 0.452 [0.24, 0.684] b |
| Grass | 1.8 [0.968, 3.35] a | 132 [89.7, 194] a | 60.8 [28.2, 131] b | 15.7 [10.5, 23.4] b | 58.4 [30.3, 112] b | 0.583 [0.314, 0.81] b |
| Legume | 1.8 [0.968, 3.35] a | 122 [82.9, 180] a | 64.3 [32.6, 127] b | 15.7 [10.5, 23.4] b | 55.8 [29.3, 107] b | 0.471 [0.254, 0.7] b |
| Ridomil | 1.8 [0.968, 3.35] a | 123 [83.5, 180] a | 37.3 [20.1, 69.3] b | 15 [10.1, 22.2] b | 51.1 [27, 96.9] b | 0.401 [0.212, 0.625] b |
| Pearl | ||||||
| Control | 1.7 [0.898, 3.22] ab | 120 [79.9, 179] b | 60.6 [23.9, 154] bc | 14.3 [9.41, 21.7] bc | 37.5 [19.4, 72.7] b | 0.449 [0.173, 0.761] bc |
| Bark | 3.7 [2.4, 5.7] a | 336 [229, 493] a | 311 [173, 561] a | 43.4 [29.2, 64.6] a | 436 [230, 827] a | 0.879 [0.746, 0.947] a |
| Brassica | 1.7 [0.898, 3.22] ab | 113 [75.7, 170] b | 30.7 [15.2, 62.1] c | 14.3 [9.47, 21.7] bc | 38.8 [20, 75.2] b | 0.32 [0.142, 0.572] c |
| Grass | 2.4 [1.4, 4.11] ab | 165 [110, 247] b | 74.8 [39.6, 141] bc | 26 [16.9, 39.9] ab | 85.6 [42.6, 172] b | 0.505 [0.291, 0.716] bc |
| Legume | 2.2 [1.26, 3.86] ab | 166 [113, 245] b | 170 [74.3, 387] ab | 17.8 [11.9, 26.5] bc | 98.7 [50.5, 193] b | 0.739 [0.449, 0.908] ab |
| Ridomil | 1.4 [0.693, 2.83] b | 100 [65.5, 153] b | 37.5 [15.2, 92.6] c | 13 [8.39, 20.1] c | 41.3 [20.8, 82.2] b | 0.386 [0.138, 0.712] bc |
Table 5 analysis
In this table, we have three positive right-skewed continuous variables. There are a few zeros (5) in the data, presumably for dead plants. These are excluded from the analysis. So, the means that we present here are conditional on the plant having survived.
For the three continuous variables we log-transform and fit a general linear mixed model with two main fixed effects, cover and cultivar, and their interaction. Replicate gets a random intercept. For the probability of being potbound, we fit a binomial GLMM with the same fixed and random effects. As before, we include a weakly informative prior on the fixed effects in the binomial GLMM to mitigate the issue of complete separation (0 of the individuals in the bark treatment were potbound).
table5subset <- table5 %>% filter(`Root ball diameter` > 0)
t5_rootdiam_fit <- glmmTMB(log(`Root ball diameter`) ~ Cover * Cultivar + (1 | Rep), data = table5subset)
t5_longestroot_fit <- glmmTMB(log(`Longest root`) ~ Cover * Cultivar + (1 | Rep), data = table5subset)
t5_weight_fit <- glmmTMB(log(weight) ~ Cover * Cultivar + (1 | Rep), data = table5subset)
t5_potbound_fit <- glmmTMB(Root_bound ~ Cover * Cultivar + (1 | Rep), family = binomial(link = 'logit'), data = table5subset,
priors = data.frame(prior = 'normal(0, 5)', class = 'fixef'))Model diagnostics look acceptable (not shown).
plot(simulateResiduals(t5_rootdiam_fit))
plot(simulateResiduals(t5_longestroot_fit))
plot(simulateResiduals(t5_weight_fit))
plot(simulateResiduals(t5_potbound_fit))Table 5 analysis: ANOVA tables
As before, joint F-tests for cover, cultivar, and their interaction are done for the three continuous variables. Joint Chi-squared tests are done for the binary variable, probability of being potbound.
t5_fits <- tibble(variable = c('Root ball diameter (cm)', 'Length of longest root (cm)', 'Weight (g)', 'Probability potbound'),
fit = list(t5_rootdiam_fit, t5_longestroot_fit, t5_weight_fit, t5_potbound_fit)) %>%
mutate(ANOVA = map(fit, joint_tests))
t5_ANOVAtables <- with(t5_fits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y)))) %>%
mutate(F.ratio = if_else(is.na(Chisq), F.ratio, as.numeric(NA))) %>%
relocate(Chisq, .after = 'F.ratio')
t5_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(c(F.ratio, Chisq), decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts %>%
fmt_missing| model term | df1 | df2 | F.ratio | Chisq | p.value |
|---|---|---|---|---|---|
| Root ball diameter (cm) | |||||
| Cover | 5 | 155 | 59.12 | — | < 0.00001 |
| Cultivar | 2 | 155 | 34.13 | — | < 0.00001 |
| Cover × Cultivar | 10 | 155 | 1.51 | — | 0.14 |
| Length of longest root (cm) | |||||
| Cover | 5 | 155 | 42.80 | — | < 0.00001 |
| Cultivar | 2 | 155 | 37.02 | — | < 0.00001 |
| Cover × Cultivar | 10 | 155 | 1.81 | — | 0.063 |
| Weight (g) | |||||
| Cover | 5 | 155 | 40.25 | — | < 0.00001 |
| Cultivar | 2 | 155 | 21.33 | — | < 0.00001 |
| Cover × Cultivar | 10 | 155 | 1.61 | — | 0.11 |
| Probability potbound | |||||
| Cover | 5 | Inf | — | 9.36 | 0.096 |
| Cultivar | 2 | Inf | — | 2.79 | 0.25 |
| Cover × Cultivar | 10 | Inf | — | 4.18 | 0.94 |
Table 5 analysis: means separation
We present both the main effects and the interaction effects, comparing cover within each cultivar.
t5_fits <- t5_fits %>%
mutate(EMM_cover = map(fit, ~ emmeans(., pairwise ~ Cover, type = 'response', inv.lbl = 'emmean')),
CLD_cover = map(EMM_cover, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_cultivar = map(fit, ~ emmeans(., pairwise ~ Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_cultivar = map(EMM_cultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)),
EMM_coverbycultivar = map(fit, ~ emmeans(., pairwise ~ Cover | Cultivar, type = 'response', inv.lbl = 'emmean')),
CLD_coverbycultivar = map(EMM_coverbycultivar, ~ cld(.$emmeans, adjust = 'tukey', Letters = letters, decreasing = TRUE)))
t5_meantables_cover <- with(t5_fits, map2_dfr(variable, CLD_cover, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t5_meantables_cultivar <- with(t5_fits, map2_dfr(variable, CLD_cultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t5_meantables_coverbycultivar <- with(t5_fits, map2_dfr(variable, CLD_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t5_contrasttables_cover <- with(t5_fits, map2_dfr(variable, EMM_cover, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t5_contrasttables_cultivar <- with(t5_fits, map2_dfr(variable, EMM_cultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t5_contrasttables_coverbycultivar <- with(t5_fits, map2_dfr(variable, EMM_coverbycultivar, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
dplyr::select(-df) %>%
relocate(odds.ratio, .after = 'ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t5_meantables_cover %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover), names_from = variable, values_from = display) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover') %>%
tabopts()| Cover | Root ball diameter (cm) | Length of longest root (cm) | Weight (g) | Probability potbound |
|---|---|---|---|---|
| Control | 26.4 [21.8, 31.9] b | 16.9 [13.8, 20.7] b | 126 [93.7, 170] b | 0.102 [0.0168, 0.432] a |
| Bark | 85.7 [71.2, 103] a | 47.9 [39.2, 58.5] a | 504 [377, 674] a | 0.00283 [2.8e-06, 0.742] a |
| Brassica | 28 [23.3, 33.7] b | 17.5 [14.3, 21.3] b | 120 [90, 161] b | 0.068 [0.0042, 0.558] a |
| Grass | 27 [22.4, 32.6] b | 16 [13, 19.6] b | 138 [103, 185] b | 0.32 [0.0981, 0.67] a |
| Legume | 32.6 [27.1, 39.4] b | 19 [15.6, 23.3] b | 136 [102, 183] b | 0.108 [0.0188, 0.431] a |
| Ridomil | 26 [21.6, 31.4] b | 15.6 [12.7, 19.1] b | 126 [94.2, 170] b | 0.22 [0.0531, 0.586] a |
t5_meantables_cultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cultivar), names_from = variable, values_from = display) %>%
arrange(Cultivar) %>%
gt(caption = 'Means separation: cultivar') %>%
tabopts()| Cultivar | Root ball diameter (cm) | Length of longest root (cm) | Weight (g) | Probability potbound |
|---|---|---|---|---|
| Gumbo | 44.4 [38.9, 50.6] a | 27.5 [23.9, 31.7] a | 208 [168, 259] a | 0.0628 [0.0105, 0.298] a |
| Legacy | 31.8 [27.8, 36.3] b | 18.9 [16.4, 21.8] b | 174 [140, 216] a | 0.174 [0.0346, 0.553] a |
| Pearl | 27.1 [23.6, 31] c | 15.7 [13.6, 18.2] c | 118 [94.3, 147] b | 0.0417 [0.00384, 0.329] a |
t5_meantables_coverbycultivar %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover, Cultivar), names_from = variable, values_from = display) %>%
group_by(Cultivar) %>%
arrange(Cover) %>%
gt(caption = 'Means separation: cover within each cultivar') %>%
tabopts()| Cover | Root ball diameter (cm) | Length of longest root (cm) | Weight (g) | Probability potbound |
|---|---|---|---|---|
| Gumbo | ||||
| Control | 34.5 [25.7, 46.2] b | 22.3 [16.2, 30.6] b | 139 [89.3, 215] b | 0.0472 [0.00298, 0.451] a |
| Bark | 97.9 [73.1, 131] a | 57.3 [41.7, 78.7] a | 572 [368, 888] a | 0.00251 [7.68e-07, 0.892] a |
| Brassica | 40 [29.8, 53.6] b | 28.2 [20.5, 38.7] b | 175 [113, 272] b | 0.0704 [0.0037, 0.607] a |
| Grass | 38.4 [28.7, 51.5] b | 23.6 [17.1, 32.4] b | 188 [121, 292] b | 0.26 [0.0401, 0.747] a |
| Legume | 36.7 [27.4, 49.2] b | 21.7 [15.8, 29.9] b | 156 [100, 242] b | 0.0746 [0.00421, 0.606] a |
| Ridomil | 40 [29.8, 53.6] b | 23.6 [17.2, 32.4] b | 201 [130, 313] b | 0.254 [0.0385, 0.743] a |
| Legacy | ||||
| Control | 26.8 [19.7, 36.5] b | 15.8 [11.3, 22.1] b | 160 [101, 254] b | 0.268 [0.0403, 0.761] a |
| Bark | 89.1 [66.5, 119] a | 48.1 [35, 66.1] a | 543 [350, 843] a | 0.00436 [4.84e-07, 0.975] a |
| Brassica | 25.8 [19.3, 34.6] b | 15.2 [11.1, 20.9] b | 145 [93.4, 225] b | 0.383 [0.0752, 0.826] a |
| Grass | 22.3 [16.7, 29.9] b | 13.4 [9.74, 18.4] b | 120 [77.5, 187] b | 0.39 [0.0775, 0.83] a |
| Legume | 31 [23.1, 41.5] b | 17.8 [12.9, 24.4] b | 135 [87, 210] b | 0.177 [0.02, 0.693] a |
| Ridomil | 24.1 [18, 32.3] b | 16.7 [12.2, 23] b | 135 [86.8, 209] b | 0.39 [0.0775, 0.83] a |
| Pearl | ||||
| Control | 19.9 [14.6, 27.1] bc | 13.7 [9.79, 19.1] bc | 90.4 [57, 143] b | 0.0755 [0.00383, 0.634] a |
| Bark | 72 [53.8, 96.6] a | 39.9 [29, 54.8] a | 412 [265, 640] a | 0.00208 [2.63e-08, 0.994] a |
| Brassica | 21.3 [15.9, 28.6] bc | 12.4 [9.03, 17] bc | 68.5 [44.1, 106] b | 0.0082 [7.19e-06, 0.905] a |
| Grass | 23 [16.9, 31.3] bc | 12.9 [9.25, 18] bc | 116 [73, 184] b | 0.316 [0.0493, 0.804] a |
| Legume | 30.6 [22.5, 41.6] b | 17.9 [12.8, 25] b | 120 [75.8, 190] b | 0.0919 [0.00508, 0.667] a |
| Ridomil | 18.2 [13.4, 24.8] c | 9.59 [6.87, 13.4] c | 74.5 [47, 118] b | 0.0933 [0.00548, 0.658] a |
Table 6 analysis
Here, we treat assay, cover crop treatment, and their interaction as fixed effects. Replication nested within assay is a random effect. Stem count is analyzed with a Poisson response distribution, and the two weight variables are strictly positive and right-skewed so they are log-transformed before analysis.
t6_count_fit <- glmmTMB(Count ~ Cover * Assay + (1 | Rep:Assay), family = poisson(link = 'log'), data = table6)
t6_wt_fit <- glmmTMB(log(weight_g) ~ Cover * Assay + (1 | Rep:Assay), data = table6)
t6_wtperseedling_fit <- glmmTMB(log(wt_per_seedling) ~ Cover * Assay + (1 | Rep:Assay), data = table6)Model diagnostics look acceptable (not shown) for the most part. The Poisson fit does not have perfect residual behavior but the negative binomial (not shown) is not better, so the Poisson is the best we can do. In any case there are no significant differences between treatments so we are not committing any Type I error with the Poisson model.
plot(simulateResiduals(t6_count_fit))
plot(simulateResiduals(t6_wt_fit))
plot(simulateResiduals(t6_wtperseedling_fit))Table 6 results: ANOVA tables
For the Poisson fit (seedling count), an analysis of deviance table is produced using joint Chi-squared tests for each model term. For the other two variables, joint F-tests for each model term are used to produce an ANOVA table.
t6_fits <- tibble(variable = c('Seedling count', 'Seedling weight (g)', 'Weight/Seedling (g)'),
fit = list(t6_count_fit, t6_wt_fit, t6_wtperseedling_fit)) %>%
mutate(ANOVA = map(fit, joint_tests))
t6_ANOVAtables <- with(t6_fits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y)))) %>%
mutate(F.ratio = if_else(is.na(Chisq), F.ratio, as.numeric(NA)))
t6_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(c(F.ratio, Chisq), decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts %>%
fmt_missing| model term | df1 | df2 | F.ratio | Chisq | p.value |
|---|---|---|---|---|---|
| Seedling count | |||||
| Cover | 5 | Inf | — | 0.76 | 0.98 |
| Assay | 1 | Inf | — | 0.01 | 0.91 |
| Cover × Assay | 5 | Inf | — | 1.80 | 0.88 |
| Seedling weight (g) | |||||
| Cover | 5 | 46 | 0.92 | — | 0.48 |
| Assay | 1 | 46 | 0.18 | — | 0.68 |
| Cover × Assay | 5 | 46 | 0.62 | — | 0.69 |
| Weight/Seedling (g) | |||||
| Cover | 5 | 46 | 0.84 | — | 0.53 |
| Assay | 1 | 46 | 0.43 | — | 0.52 |
| Cover × Assay | 5 | 46 | 0.80 | — | 0.56 |
Table 6 results: means separation
We see no significant effect of assay in any case, so we produce means separation tables averaged over the two levels of assay, and using Tukey HSD to adjust the p-values and Sidak adjustment to adjust the 95% confidence intervals.
There are no significant differences.
t6_fits <- t6_fits %>%
mutate(
EMM = map(fit, ~ emmeans(., pairwise ~ Cover, type = 'response', inv.lbl = 'emmean')),
CLD = map(EMM, ~ cld(.$emmeans, Letters = letters, adjust = 'tukey', decreasing = TRUE))
)
t6_meantables <- with(t6_fits, map2_dfr(variable, CLD, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t6_contrasttables <- with(t6_fits, map2_dfr(variable, EMM, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
relocate(t.ratio, .after = 'z.ratio')
t6_meantables %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover), names_from = variable, values_from = display) %>%
arrange(Cover) %>%
gt() %>%
tabopts()| Cover | Seedling count | Seedling weight (g) | Weight/Seedling (g) |
|---|---|---|---|
| Control | 5.09 [3.52, 7.36] a | 9.54 [5.67, 16] a | 1.89 [1.18, 3.02] a |
| Bark | 4.4 [2.96, 6.54] a | 10.6 [6.3, 17.8] a | 2.48 [1.55, 3.97] a |
| Brassica | 4.89 [3.36, 7.13] a | 11.3 [6.7, 18.9] a | 2.34 [1.46, 3.74] a |
| Grass | 4.85 [3.32, 7.09] a | 10.9 [6.46, 18.3] a | 2.3 [1.44, 3.69] a |
| Legume | 4.58 [3.1, 6.76] a | 7.17 [4.26, 12.1] a | 1.69 [1.06, 2.71] a |
| Ridomil | 4.5 [3.04, 6.66] a | 9.58 [5.7, 16.1] a | 2.17 [1.35, 3.47] a |
Table 7 analysis
For Table 7, number of seedlings 4 and 100 days after seeding are count variables so they are analyzed with a Poisson mixed model. Total weight in grams is a non-negative skewed variable so it is log(x + 0.1) transformed because of the single zero value in the dataset, then modeled with a general linear mixed model. Percentage seedlings germinated is actually a count out of 18 that are germinated, so it is converted back to the count and then modeled with a logistic GLMM, where the modeled response ends up being the percentage seedlings out of 18. In all cases, we have one fixed effect, cover crop treatment, and one random effect, replicate.
table7 <- table7 %>%
rename(total_weight = `Totalwt (gm)`) %>%
mutate(across(c(Count4d, Count100d), as.integer))
t7_totalweight_fit <- glmmTMB(log(total_weight + 0.1) ~ Cover + (1 | Rep), data = table7)
t7_survival4d_fit <- glmmTMB(cbind(Count4d, 18 - Count4d) ~ Cover + (1 | Rep), family = binomial(link = 'logit'), data = table7)
t7_survival100d_fit <- glmmTMB(cbind(Count100d, 18 - Count100d) ~ Cover + (1 | Rep), family = binomial(link = 'logit'), data = table7)Model diagnostics look acceptable (not shown).
plot(simulateResiduals(t7_totalweight_fit))
plot(simulateResiduals(t7_survival4d_fit))
plot(simulateResiduals(t7_survival100d_fit))Table 7 results: ANOVA tables
Analysis of deviance tables are produced using joint Chi-squared tests for each model term (joint F-tests for the ANOVA table in the case of seedling weight). There is only one fixed effect term in each model thus only one test per variable.
t7_fits <- tibble(variable = c('Seedling survival 4 das (%)', 'Seedling survival 100 das (%)', 'Seedling weight 100 das (g)'),
fit = list(t7_survival4d_fit, t7_survival100d_fit, t7_totalweight_fit)) %>%
mutate(ANOVA = map(fit, joint_tests))
t7_ANOVAtables <- with(t7_fits, map2_dfr(variable, ANOVA, ~ tibble(variable = .x, as_tibble(.y)))) %>%
mutate(F.ratio = if_else(is.na(Chisq), F.ratio, as.numeric(NA)))
t7_ANOVAtables %>%
mutate(`model term` = gsub(':', ' × ', `model term`, fixed = TRUE) %>% map(md)) %>%
gt(groupname_col = 'variable') %>%
fmt_number(c(F.ratio, Chisq), decimals = 2) %>%
fmt_number(p.value, n_sigfig = 2) %>%
sub_small_vals(p.value, threshold = 0.00001, small_pattern = '< 0.00001') %>%
tabopts %>%
fmt_missing| model term | df1 | df2 | F.ratio | Chisq | p.value |
|---|---|---|---|---|---|
| Seedling survival 4 das (%) | |||||
| Cover | 8 | Inf | — | 44.22 | < 0.00001 |
| Seedling survival 100 das (%) | |||||
| Cover | 8 | Inf | — | 49.78 | < 0.00001 |
| Seedling weight 100 das (g) | |||||
| Cover | 8 | 25 | 1.06 | — | 0.42 |
Table 7 results: means separation
Means separations are done by levels of cover as before, but with additional treatment levels. (Note that the contrasts for the percent seedlings variable are odds ratios because the model is fit on the log odds scale. The response variable for that model is a proportion between 0 and 1.)
t7_fits <- t7_fits %>%
mutate(
EMM = map(fit, ~ emmeans(., pairwise ~ Cover, type = 'response', inv.lbl = 'emmean')),
CLD = map(EMM, ~ cld(.$emmeans, Letters = letters, adjust = 'tukey', decreasing = TRUE))
)
t7_meantables <- with(t7_fits, map2_dfr(variable, CLD, ~ tibble(variable = .x, as_tibble(.y)) %>% rename(lower.CL = any_of('asymp.LCL'), upper.CL = any_of('asymp.UCL'))))
t7_contrasttables <- with(t7_fits, map2_dfr(variable, EMM, ~ tibble(variable = .x, as_tibble(.y$contrasts)))) %>%
rename(difference = estimate) %>%
relocate(difference, .before = 'odds.ratio') %>%
relocate(z.ratio, .after = 't.ratio') %>%
mutate(null = if_else(is.na(null), 0, 1))
t7_meantables %>%
mutate(across(c(emmean, lower.CL, upper.CL), ~ signif(., 3)),
across(c(emmean, lower.CL, upper.CL), ~ if_else(grepl('survival', variable), . * 100, .)),
display = glue::glue('{emmean} [{lower.CL}, {upper.CL}] {trimws(.group)}')) %>%
pivot_wider(id_cols = c(Cover), names_from = variable, values_from = display) %>%
arrange(Cover) %>%
gt() %>%
tabopts()| Cover | Seedling survival 4 das (%) | Seedling survival 100 das (%) | Seedling weight 100 das (g) |
|---|---|---|---|
| Control | 25.4 [14.3, 40.9] b | 45.4 [27.6, 64.4] abc | 18.4 [5.3, 63.3] a |
| Bark | 25.4 [14.3, 40.9] b | 38.4 [22.1, 57.8] abcd | 14.8 [4.25, 50.9] a |
| Brassica | 12.1 [5.17, 25.7] b | 21.5 [10.3, 39.5] d | 9.41 [2.68, 32.5] a |
| Grass | 16.5 [8, 31] b | 24.8 [12.4, 43.4] cd | 12.7 [3.65, 43.8] a |
| Legume | 29.8 [17.7, 45.7] b | 47.7 [29.5, 66.4] ab | 23.6 [6.82, 81.1] a |
| Ridomil | 15.4 [7.27, 29.7] b | 30.5 [16.3, 49.7] bcd | 6.27 [1.76, 21.7] a |
| Ckpos | 18.1 [6.06, 43.1] b | 22.8 [8.41, 48.8] bcd | 11.6 [1.57, 81.8] a |
| Ckneg | 68 [42.2, 86.1] a | 70.5 [40.3, 89.4] a | 28.3 [3.95, 199] a |
| Ckneut | 37 [17.5, 62] ab | 70.5 [40.3, 89.4] a | 20.4 [2.83, 144] a |
Write output to file
Here we write all output to an Excel spreadsheet with multiple tabs (code not evaluated when notebook is run). The values in the spreadsheet are not formatted for display as the tables in the document above are.
out_list <- list(
'T1 ANOVAs' = t1_ANOVAtables,
'T1 means' = t1_meantables,
'T1 contrasts' = t1_contrasttables,
'T2 means' = t2_meantables,
'T2 contrasts' = t2_contrasttables,
'T3 ANOVAs' = t3_ANOVAtables,
'T3 means cover' = t3_meantables_cover,
'T3 contrasts cover' = t3_contrasttables_cover,
'T3 means cultivar' = t3_meantables_cultivar,
'T3 contrasts cultivar' = t3_contrasttables_cultivar,
'T3 means cover x cultivar' = t3_meantables_coverbycultivar,
'T3 contrasts cover x cultivar'= t3_contrasttables_coverbycultivar,
'T4 ANOVAs' = t4_ANOVAtables,
'T4 means cover' = t4_meantables_cover,
'T4 contrasts cover' = t4_contrasttables_cover,
'T4 means cultivar' = t4_meantables_cultivar,
'T4 contrasts cultivar' = t4_contrasttables_cultivar,
'T4 means cover x cultivar' = t4_meantables_coverbycultivar,
'T4 contrasts cover x cultivar' = t4_contrasttables_coverbycultivar,
'T5 ANOVAs' = t5_ANOVAtables,
'T5 means cover' = t5_meantables_cover,
'T5 contrasts cover' = t5_contrasttables_cover,
'T5 means cultivar' = t5_meantables_cultivar,
'T5 contrasts cultivar' = t5_contrasttables_cultivar,
'T5 means cover x cultivar' = t5_meantables_coverbycultivar,
'T5 contrasts cover x cultivar' = t5_contrasttables_coverbycultivar,
'T6 ANOVAs' = t6_ANOVAtables,
'T6 means' = t6_meantables,
'T6 contrasts' = t6_contrasttables,
'T7 ANOVAs' = t7_ANOVAtables,
'T7 means' = t7_meantables,
'T7 contrasts' = t7_contrasttables
)
# Replace all instances of Inf with NA
replace_infs <- function(dt) mutate(dt, across(where(is.numeric), \(x) if_else(!is.finite(x), as.numeric(NA), x)))
out_list <- map(out_list, replace_infs)
openxlsx::write.xlsx(out_list, file = 'Cover crop blueberry all stat output.xlsx', overwrite = TRUE)