knitr::opts_chunk$set(echo = TRUE, warning = FALSE)library(tidyverse)
library(viridis)
library(cowplot)
library(survival)
library(gsheet)
library(ggthemes)
library(ggfortify)
library(survminer)
library(lme4)
library(metafor)
library(janitor)
library(here)
rust_sev <- read_csv("data/dat-sev.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## study = col_double(),
## year = col_double(),
## location = col_character(),
## state = col_character(),
## n_spray = col_double(),
## brand_name = col_character(),
## mean_sev = col_double(),
## mean_yld = col_double(),
## v_sev = col_double(),
## v_yld = col_double(),
## yld_check = col_double(),
## v_yld_check = col_double(),
## sev_check = col_double(),
## v_sev_check = col_double(),
## log_sev = col_double(),
## vi_sev = col_double(),
## n2 = col_double(),
## n3 = col_double(),
## design = col_double()
## )head(rust_sev)summary(rust_sev$vi_sev)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00093 0.00307 0.24779 0.01193 149.30794rust_sev = rust_sev %>%
filter(vi_sev>0.0001)mv_sev_HCS <- rma.mv(log_sev, vi_sev,
mods = ~brand_name,
random = list(~brand_name | factor(study)),
struct = "HCS",
method = "ML",
control = list(optimizer = "nlm"),
data = rust_sev
)
mv_sev_HCS##
## Multivariate Meta-Analysis Model (k = 1266; method: ML)
##
## Variance Components:
##
## outer factor: factor(study) (nlvls = 174)
## inner factor: brand_name (nlvls = 9)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.2171 0.4659 149 no AACHECK
## tau^2.2 0.9952 0.9976 144 no AZOX + BENZ
## tau^2.3 0.8154 0.9030 115 no BIXF + TFLX + PROT
## tau^2.4 0.8687 0.9320 116 no PICO + BENZ
## tau^2.5 0.5489 0.7409 169 no PICO + CYPR
## tau^2.6 0.7062 0.8403 149 no PICO + TEBU
## tau^2.7 0.7445 0.8629 115 no PYRA + EPOX + FLUX
## tau^2.8 0.5990 0.7740 143 no TFLX + CYPR
## tau^2.9 0.8344 0.9135 166 no TFLX + PROT
## rho 0.8495 no
##
## Test for Residual Heterogeneity:
## QE(df = 1257) = 292126.3351, p-val < .0001
##
## Test of Moderators (coefficients 2:9):
## QM(df = 8) = 1233.1513, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 4.1347 0.0360 114.9271 <.0001 4.0642
## brand_nameAZOX + BENZ -1.3017 0.0541 -24.0435 <.0001 -1.4078
## brand_nameBIXF + TFLX + PROT -1.4612 0.0504 -28.9723 <.0001 -1.5601
## brand_namePICO + BENZ -1.3482 0.0520 -25.9251 <.0001 -1.4501
## brand_namePICO + CYPR -0.8268 0.0338 -24.4286 <.0001 -0.8932
## brand_namePICO + TEBU -1.0793 0.0418 -25.8019 <.0001 -1.1613
## brand_namePYRA + EPOX + FLUX -1.2812 0.0468 -27.3721 <.0001 -1.3730
## brand_nameTFLX + CYPR -0.8651 0.0375 -23.0799 <.0001 -0.9385
## brand_nameTFLX + PROT -1.2719 0.0462 -27.5374 <.0001 -1.3624
## ci.ub
## intrcpt 4.2052 ***
## brand_nameAZOX + BENZ -1.1956 ***
## brand_nameBIXF + TFLX + PROT -1.3624 ***
## brand_namePICO + BENZ -1.2463 ***
## brand_namePICO + CYPR -0.7605 ***
## brand_namePICO + TEBU -0.9973 ***
## brand_namePYRA + EPOX + FLUX -1.1895 ***
## brand_nameTFLX + CYPR -0.7916 ***
## brand_nameTFLX + PROT -1.1814 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## HCS
efficacy_sev <- data.frame(cbind(
(1 - exp(mv_sev_HCS$b)) * 100,
(1 - exp(mv_sev_HCS$ci.lb)) * 100,
(1 - exp(mv_sev_HCS$ci.ub)) * 100
))
efficacy_sev#Organize the data.frame
efficacy_sev = efficacy_sev
names (efficacy_sev) = c("efficacy", "efficacy_up", "efficacy_lw")
efficacy_sev = efficacy_sev %>%
mutate(fungicide = c("check", "AZOX + BENZ", "BIXF + TFLX + PROT", "PICO + BENZ", "PICO + CYPR","PICO + TEBU", "PYRA + EPOX + FLUX", "TFLX + CYPR", "TFLX + PROT")) %>%
filter(fungicide != "check")
efficacy_sev#write_csv(efficacy_index, "data/efficacy_index.csv")
openxlsx::write.xlsx(efficacy_sev, here("data","efficacy_sev.xlsx"), colNames = TRUE)We can set linear contrasts between treatments of interest and get the P-value using the anova function.
anova(mv_sev_HCS, L = rbind(
c(0, 1, -1, 0, 0, 0, 0, 0, 0),
c(0, 1, 0, -1, 0, 0, 0, 0, 0),
c(0, 1, 0, 0, -1, 0, 0, 0, 0),
c(0, 1, 0, 0, 0, -1, 0, 0, 0),
c(0, 1, 0, 0, 0, 0, -1, 0, 0),
c(0, 1, 0, 0, 0, 0, 0, -1, 0),
c(0, 1, 0, 0, 0, 0, 0, 0, -1),
c(0, 0, 1, -1, 0, 0, 0, 0, 0),
c(0, 0, 1, 0, -1, 0, 0, 0, 0),
c(0, 0, 1, 0, 0, -1, 0, 0, 0),
c(0, 0, 1, 0, 0, 0, -1, 0, 0),
c(0, 0, 1, 0, 0, 0, 0, -1, 0),
c(0, 0, 1, 0, 0, 0, 0, 0, -1),
c(0, 0, 0, 1, -1, 0, 0, 0, 0),
c(0, 0, 0, 1, 0, -1, 0, 0, 0),
c(0, 0, 0, 1, 0, 0, -1, 0, 0),
c(0, 0, 0, 1, 0, 0, 0, -1, 0),
c(0, 0, 0, 1, 0, 0, 0, 0, -1),
c(0, 0, 0, 0, 1, -1, 0, 0, 0),
c(0, 0, 0, 0, 1, 0, -1, 0, 0),
c(0, 0, 0, 0, 1, 0, 0, -1, 0),
c(0, 0, 0, 0, 1, 0, 0, 0, -1),
c(0, 0, 0, 0, 0, 1, -1, 0, 0),
c(0, 0, 0, 0, 0, 1, 0, -1, 0),
c(0, 0, 0, 0, 0, 1, 0, 0, -1),
c(0, 0, 0, 0, 0, 0, 1, -1, 0),
c(0, 0, 0, 0, 0, 0, 1, 0, -1),
c(0, 0, 0, 0, 0, 0, 0, 1, -1)
))##
## Hypotheses:
## 1: brand_nameAZOX + BENZ - brand_nameBIXF + TFLX + PROT = 0
## 2: brand_nameAZOX + BENZ - brand_namePICO + BENZ = 0
## 3: brand_nameAZOX + BENZ - brand_namePICO + CYPR = 0
## 4: brand_nameAZOX + BENZ - brand_namePICO + TEBU = 0
## 5: brand_nameAZOX + BENZ - brand_namePYRA + EPOX + FLUX = 0
## 6: brand_nameAZOX + BENZ - brand_nameTFLX + CYPR = 0
## 7: brand_nameAZOX + BENZ - brand_nameTFLX + PROT = 0
## 8: brand_nameBIXF + TFLX + PROT - brand_namePICO + BENZ = 0
## 9: brand_nameBIXF + TFLX + PROT - brand_namePICO + CYPR = 0
## 10: brand_nameBIXF + TFLX + PROT - brand_namePICO + TEBU = 0
## 11: brand_nameBIXF + TFLX + PROT - brand_namePYRA + EPOX + FLUX = 0
## 12: brand_nameBIXF + TFLX + PROT - brand_nameTFLX + CYPR = 0
## 13: brand_nameBIXF + TFLX + PROT - brand_nameTFLX + PROT = 0
## 14: brand_namePICO + BENZ - brand_namePICO + CYPR = 0
## 15: brand_namePICO + BENZ - brand_namePICO + TEBU = 0
## 16: brand_namePICO + BENZ - brand_namePYRA + EPOX + FLUX = 0
## 17: brand_namePICO + BENZ - brand_nameTFLX + CYPR = 0
## 18: brand_namePICO + BENZ - brand_nameTFLX + PROT = 0
## 19: brand_namePICO + CYPR - brand_namePICO + TEBU = 0
## 20: brand_namePICO + CYPR - brand_namePYRA + EPOX + FLUX = 0
## 21: brand_namePICO + CYPR - brand_nameTFLX + CYPR = 0
## 22: brand_namePICO + CYPR - brand_nameTFLX + PROT = 0
## 23: brand_namePICO + TEBU - brand_namePYRA + EPOX + FLUX = 0
## 24: brand_namePICO + TEBU - brand_nameTFLX + CYPR = 0
## 25: brand_namePICO + TEBU - brand_nameTFLX + PROT = 0
## 26: brand_namePYRA + EPOX + FLUX - brand_nameTFLX + CYPR = 0
## 27: brand_namePYRA + EPOX + FLUX - brand_nameTFLX + PROT = 0
## 28: brand_nameTFLX + CYPR - brand_nameTFLX + PROT = 0
##
## Results:
## estimate se zval pval
## 1: 0.1595 0.0523 3.0509 0.0023
## 2: 0.0465 0.0523 0.8897 0.3736
## 3: -0.4749 0.0466 -10.1977 <.0001
## 4: -0.2224 0.0475 -4.6805 <.0001
## 5: -0.0205 0.0508 -0.4036 0.6865
## 6: -0.4366 0.0474 -9.2085 <.0001
## 7: -0.0298 0.0472 -0.6320 0.5274
## 8: -0.1130 0.0511 -2.2100 0.0271
## 9: -0.6344 0.0451 -14.0788 <.0001
## 10: -0.3819 0.0471 -8.1070 <.0001
## 11: -0.1800 0.0491 -3.6676 0.0002
## 12: -0.5961 0.0457 -13.0346 <.0001
## 13: -0.1893 0.0476 -3.9808 <.0001
## 14: -0.5214 0.0460 -11.3453 <.0001
## 15: -0.2689 0.0476 -5.6430 <.0001
## 16: -0.0670 0.0495 -1.3534 0.1759
## 17: -0.4831 0.0465 -10.3887 <.0001
## 18: -0.0763 0.0480 -1.5908 0.1117
## 19: 0.2525 0.0375 6.7311 <.0001
## 20: 0.4544 0.0422 10.7548 <.0001
## 21: 0.0382 0.0352 1.0855 0.2777
## 22: 0.4450 0.0399 11.1566 <.0001
## 23: 0.2019 0.0449 4.5001 <.0001
## 24: -0.2143 0.0393 -5.4510 <.0001
## 25: 0.1925 0.0419 4.5916 <.0001
## 26: -0.4162 0.0432 -9.6415 <.0001
## 27: -0.0094 0.0457 -0.2046 0.8379
## 28: 0.4068 0.0412 9.8724 <.0001library(metafor)
mv_sev_year <- rma.mv(log_sev, vi_sev,
mods = ~brand_name*year,
random = list(~brand_name | factor(study)),
struct = "HCS",
method = "ML",
control = list(optimizer = "nlm"),
data = rust_sev %>% mutate(year= year - 2015))
mv_sev_year##
## Multivariate Meta-Analysis Model (k = 1266; method: ML)
##
## Variance Components:
##
## outer factor: factor(study) (nlvls = 174)
## inner factor: brand_name (nlvls = 9)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.2167 0.4655 149 no AACHECK
## tau^2.2 0.8939 0.9455 144 no AZOX + BENZ
## tau^2.3 0.8120 0.9011 115 no BIXF + TFLX + PROT
## tau^2.4 0.8253 0.9085 116 no PICO + BENZ
## tau^2.5 0.5510 0.7423 169 no PICO + CYPR
## tau^2.6 0.7126 0.8441 149 no PICO + TEBU
## tau^2.7 0.7464 0.8639 115 no PYRA + EPOX + FLUX
## tau^2.8 0.6011 0.7753 143 no TFLX + CYPR
## tau^2.9 0.8342 0.9134 166 no TFLX + PROT
## rho 0.8612 no
##
## Test for Residual Heterogeneity:
## QE(df = 1248) = 284314.6611, p-val < .0001
##
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 1369.5916, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 4.1151 0.0643 64.0318 <.0001 3.9892
## brand_nameAZOX + BENZ -1.8080 0.0887 -20.3895 <.0001 -1.9818
## brand_nameBIXF + TFLX + PROT -1.4326 0.1338 -10.7044 <.0001 -1.6949
## brand_namePICO + BENZ -1.9013 0.1355 -14.0350 <.0001 -2.1668
## brand_namePICO + CYPR -0.8418 0.0607 -13.8781 <.0001 -0.9607
## brand_namePICO + TEBU -1.1064 0.0758 -14.6045 <.0001 -1.2549
## brand_namePYRA + EPOX + FLUX -1.4365 0.1251 -11.4825 <.0001 -1.6816
## brand_nameTFLX + CYPR -0.9171 0.0794 -11.5487 <.0001 -1.0727
## brand_nameTFLX + PROT -1.3617 0.0826 -16.4858 <.0001 -1.5235
## year 0.0094 0.0224 0.4167 0.6769 -0.0346
## brand_nameAZOX + BENZ:year 0.2295 0.0319 7.1976 <.0001 0.1670
## brand_nameBIXF + TFLX + PROT:year 0.0039 0.0401 0.0982 0.9217 -0.0747
## brand_namePICO + BENZ:year 0.1812 0.0406 4.4657 <.0001 0.1017
## brand_namePICO + CYPR:year 0.0069 0.0209 0.3292 0.7420 -0.0342
## brand_namePICO + TEBU:year 0.0044 0.0284 0.1531 0.8783 -0.0513
## brand_namePYRA + EPOX + FLUX:year 0.0584 0.0375 1.5572 0.1194 -0.0151
## brand_nameTFLX + CYPR:year 0.0240 0.0259 0.9268 0.3541 -0.0268
## brand_nameTFLX + PROT:year 0.0375 0.0287 1.3092 0.1905 -0.0187
## ci.ub
## intrcpt 4.2411 ***
## brand_nameAZOX + BENZ -1.6342 ***
## brand_nameBIXF + TFLX + PROT -1.1703 ***
## brand_namePICO + BENZ -1.6358 ***
## brand_namePICO + CYPR -0.7229 ***
## brand_namePICO + TEBU -0.9580 ***
## brand_namePYRA + EPOX + FLUX -1.1913 ***
## brand_nameTFLX + CYPR -0.7614 ***
## brand_nameTFLX + PROT -1.1998 ***
## year 0.0533
## brand_nameAZOX + BENZ:year 0.2920 ***
## brand_nameBIXF + TFLX + PROT:year 0.0826
## brand_namePICO + BENZ:year 0.2608 ***
## brand_namePICO + CYPR:year 0.0479
## brand_namePICO + TEBU:year 0.0600
## brand_namePYRA + EPOX + FLUX:year 0.1320
## brand_nameTFLX + CYPR:year 0.0748
## brand_nameTFLX + PROT:year 0.0938
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mv_sev_year, btt = 11:18)##
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 93.6212, p-val < .0001reg1 = data.frame(mv_sev_year$beta, mv_sev_year$ci.lb, mv_sev_year$ci.ub) %>%
rownames_to_column("trat") %>%
separate(trat, into = c("lado1", "lado2"), sep = ":") %>%
separate(lado1, into = c("lixo", "lado3"),sep = "brand_name") %>%
select(-lixo) %>%
filter(lado3 != "NA") %>%
mutate(mod = c(rep("intercept", 8), rep("slope", 8))) %>%
select(-lado2)
names(reg1) = c("fungicide", "mean", "ci.lb", "ci.ub", "mod")
mean = reg1 %>%
group_by(fungicide) %>%
select(1:2,5) %>%
spread(mod, mean)
names(mean) = c("fungicide", "intercept_mean", "slope_mean")
upper = reg1 %>%
group_by(fungicide) %>%
select(1,3,5) %>%
spread(mod, ci.lb)
names(upper) = c("fungicide", "intercept_upper", "slope_upper")
lower = reg1 %>%
group_by(fungicide) %>%
select(1,4:5) %>%
spread(mod, ci.ub)
names(lower) = c("fungicide", "intercept_lower", "slope_lower")
data_model = left_join(mean, lower, by= c("fungicide")) %>%
left_join(upper, by = c("fungicide"))
library(tidyverse)
rust_sev <- read_csv("data/dat-sev.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## study = col_double(),
## year = col_double(),
## location = col_character(),
## state = col_character(),
## n_spray = col_double(),
## brand_name = col_character(),
## mean_sev = col_double(),
## mean_yld = col_double(),
## v_sev = col_double(),
## v_yld = col_double(),
## yld_check = col_double(),
## v_yld_check = col_double(),
## sev_check = col_double(),
## v_sev_check = col_double(),
## log_sev = col_double(),
## vi_sev = col_double(),
## n2 = col_double(),
## n3 = col_double(),
## design = col_double()
## )sbr_effic <- rust_sev %>%
mutate(efficacy = (1-(mean_sev/sev_check))) %>%
mutate(efficacy1 = efficacy*100) %>%
filter(brand_name!= "AACHECK")
year = seq(0,6, by = 0.1)
fungicide = NULL
year_col = NULL
for(i in 1:length(data_model$fungicide)){
data_cache = sbr_effic %>%
filter(brand_name == data_model$fungicide[i])
years = unique(data_cache$year)-2015
year = sort(years)
year = seq(first(year), last(year), by = 0.1)
year_col = c(year_col,year)
fungicide = c(fungicide, rep(data_model$fungicide[i], length(year)))
}
predicted = data.frame(year_col, fungicide) %>%
mutate(year = year_col) %>%
right_join(data_model, by = "fungicide") %>%
mutate(mean_efficacy = (1-exp(intercept_mean + slope_mean*year))*100,
CIL = (1-exp(intercept_lower + slope_lower*year))*100,
CIU = (1-exp(intercept_upper + slope_upper*year))*100,
year = year+2015) %>%
mutate(brand_name = fungicide) %>%
filter(year <2020.2) %>%
dplyr::select(-fungicide)
predictedopenxlsx::write.xlsx(predicted, here("data","predicted_efficacy.xlsx"), colNames = TRUE)library(metafor)
rust_sev1 <- rust_sev %>%
mutate(region = case_when(
state == "MT" ~ "North",
state == "BA" ~ "North",
state == "DF" ~ "North",
state == "MS" ~ "North",
state == "GO" ~ "North",
state == "TO" ~ "North",
state == "PR" ~ "South",
state == "MG" ~ "North",
state == "SP" ~ "South",
state == "RS" ~ "South"))
rust_sev1 %>%
group_by(study,region) %>%
summarise() %>%
tabyl(region) ## `summarise()` regrouping output by 'study' (override with `.groups` argument)rust_sev1 %>%
tabyl(brand_name, region)mv_sev_reg <- rma.mv(log_sev, vi_sev,
mods = ~brand_name * as.factor(region),
random = list(~brand_name | factor(study)),
struct = "HCS",
method = "ML",
#control = list(optimizer = "nlm"),
data = rust_sev1)
mv_sev_reg##
## Multivariate Meta-Analysis Model (k = 1331; method: ML)
##
## Variance Components:
##
## outer factor: factor(study) (nlvls = 177)
## inner factor: brand_name (nlvls = 9)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.2187 0.4676 177 no AACHECK
## tau^2.2 1.0512 1.0253 147 no AZOX + BENZ
## tau^2.3 0.7709 0.8780 118 no BIXF + TFLX + PROT
## tau^2.4 0.9004 0.9489 120 no PICO + BENZ
## tau^2.5 0.5318 0.7293 177 no PICO + CYPR
## tau^2.6 0.7015 0.8376 156 no PICO + TEBU
## tau^2.7 0.7461 0.8638 118 no PYRA + EPOX + FLUX
## tau^2.8 0.5622 0.7498 148 no TFLX + CYPR
## tau^2.9 0.7690 0.8769 170 no TFLX + PROT
## rho 0.8344 no
##
## Test for Residual Heterogeneity:
## QE(df = 1313) = 776986.9790, p-val < .0001
##
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 1443.4894, p-val < .0001
##
## Model Results:
##
## estimate se zval
## intrcpt 4.1704 0.0414 100.6735
## brand_nameAZOX + BENZ -1.3092 0.0658 -19.9061
## brand_nameBIXF + TFLX + PROT -1.3821 0.0557 -24.8045
## brand_namePICO + BENZ -1.3343 0.0618 -21.5902
## brand_namePICO + CYPR -0.7415 0.0385 -19.2592
## brand_namePICO + TEBU -1.0025 0.0487 -20.5784
## brand_namePYRA + EPOX + FLUX -1.2480 0.0542 -23.0389
## brand_nameTFLX + CYPR -0.7779 0.0418 -18.6131
## brand_nameTFLX + PROT -1.1836 0.0509 -23.2557
## as.factor(region)South -0.0055 0.0788 -0.0694
## brand_nameAZOX + BENZ:as.factor(region)South -0.0958 0.1254 -0.7645
## brand_nameBIXF + TFLX + PROT:as.factor(region)South -0.3855 0.1167 -3.3036
## brand_namePICO + BENZ:as.factor(region)South -0.0852 0.1213 -0.7024
## brand_namePICO + CYPR:as.factor(region)South -0.3295 0.0738 -4.4631
## brand_namePICO + TEBU:as.factor(region)South -0.3486 0.0922 -3.7798
## brand_namePYRA + EPOX + FLUX:as.factor(region)South -0.1397 0.1092 -1.2792
## brand_nameTFLX + CYPR:as.factor(region)South -0.3846 0.0815 -4.7187
## brand_nameTFLX + PROT:as.factor(region)South -0.4251 0.0987 -4.3077
## pval ci.lb ci.ub
## intrcpt <.0001 4.0892 4.2516
## brand_nameAZOX + BENZ <.0001 -1.4381 -1.1803
## brand_nameBIXF + TFLX + PROT <.0001 -1.4913 -1.2728
## brand_namePICO + BENZ <.0001 -1.4555 -1.2132
## brand_namePICO + CYPR <.0001 -0.8170 -0.6660
## brand_namePICO + TEBU <.0001 -1.0980 -0.9071
## brand_namePYRA + EPOX + FLUX <.0001 -1.3541 -1.1418
## brand_nameTFLX + CYPR <.0001 -0.8598 -0.6960
## brand_nameTFLX + PROT <.0001 -1.2833 -1.0838
## as.factor(region)South 0.9447 -0.1598 0.1489
## brand_nameAZOX + BENZ:as.factor(region)South 0.4446 -0.3416 0.1499
## brand_nameBIXF + TFLX + PROT:as.factor(region)South 0.0010 -0.6142 -0.1568
## brand_namePICO + BENZ:as.factor(region)South 0.4825 -0.3229 0.1525
## brand_namePICO + CYPR:as.factor(region)South <.0001 -0.4742 -0.1848
## brand_namePICO + TEBU:as.factor(region)South 0.0002 -0.5294 -0.1678
## brand_namePYRA + EPOX + FLUX:as.factor(region)South 0.2008 -0.3538 0.0744
## brand_nameTFLX + CYPR:as.factor(region)South <.0001 -0.5444 -0.2249
## brand_nameTFLX + PROT:as.factor(region)South <.0001 -0.6186 -0.2317
##
## intrcpt ***
## brand_nameAZOX + BENZ ***
## brand_nameBIXF + TFLX + PROT ***
## brand_namePICO + BENZ ***
## brand_namePICO + CYPR ***
## brand_namePICO + TEBU ***
## brand_namePYRA + EPOX + FLUX ***
## brand_nameTFLX + CYPR ***
## brand_nameTFLX + PROT ***
## as.factor(region)South
## brand_nameAZOX + BENZ:as.factor(region)South
## brand_nameBIXF + TFLX + PROT:as.factor(region)South ***
## brand_namePICO + BENZ:as.factor(region)South
## brand_namePICO + CYPR:as.factor(region)South ***
## brand_namePICO + TEBU:as.factor(region)South ***
## brand_namePYRA + EPOX + FLUX:as.factor(region)South
## brand_nameTFLX + CYPR:as.factor(region)South ***
## brand_nameTFLX + PROT:as.factor(region)South ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mv_sev_reg, btt = 11:18)##
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 54.2877, p-val < .0001reg1 = data.frame(mv_sev_reg$beta, mv_sev_reg$ci.lb, mv_sev_reg$ci.ub, mv_sev_reg$se) %>%
rownames_to_column("trat") %>%
separate(trat, into = c("lado1", "lado2"), sep = ":") %>%
separate(lado2, into = c("lixo","lado3"),sep = "region" ) %>%
select(-lixo) %>%
separate(lado1, into = c("lixo","lado1"),sep = "brand_name") %>%
select(-lixo) %>%
filter(lado1 != "NA") %>%
mutate(n = seq(1:16))
names(reg1) = c("fungicide", "region", "mean_log", "ci.lb_log", "ci.ub_log", "SE_eff","n")
reg2 = reg1 %>%
filter(n < 9) %>%
mutate(region = rep("North", length(fungicide)))
reg3 = reg1 %>%
filter(n > 8) %>%
mutate(region = rep("South", length(fungicide)))
reg4 = rbind(reg2,reg3) %>%
select(-"n")
openxlsx::write.xlsx(reg4, here("data","eff_log_region.xlsx"), colNames = TRUE)reg1 = data.frame(mv_sev_reg$beta, mv_sev_reg$ci.lb, mv_sev_reg$ci.ub) %>%
rownames_to_column("trat") %>%
separate(trat, into = c("lado1", "lado2"), sep = ":") %>%
separate(lado2, into = c("lixo","lado3"),sep = "region" ) %>%
select(-lixo) %>%
separate(lado1, into = c("lixo","lado1"),sep = "brand_name") %>%
select(-lixo) %>%
filter(lado1 != "NA") %>%
mutate(n = seq(1:16))
names(reg1) = c("fungicide", "region", "mean", "ci.lb", "ci.ub", "n")
reg2 = reg1 %>%
filter(n < 9) %>%
mutate(region = rep("North", length(fungicide)))
reg3 = reg1 %>%
filter(n > 8) %>%
mutate(region = rep("South", length(fungicide)))
reg4 = rbind(reg2,reg3)
mean = reg4%>%
group_by(fungicide) %>%
select(1:3) %>%
spread(region, mean) %>%
mutate(mean = (1-exp((North + South)))*100) %>%
select(1,4)
upper = reg4%>%
group_by(fungicide) %>%
select(1,2,4) %>%
spread(region, ci.lb) %>%
mutate(upper = (1-exp((North + South)))*100) %>%
select(1,4)
lower = reg4%>%
group_by(fungicide) %>%
select(1,2,5) %>%
spread(region, ci.ub) %>%
mutate(lower = (1-exp((North + South)))*100) %>%
select(1,4)
reg5 = left_join(mean, lower, by= c("fungicide")) %>%
left_join(upper, by = c("fungicide")) %>%
mutate(region = rep("South", length("fungicide"))) %>%
select("fungicide", "region", "mean", "lower", "upper")
North = reg4 %>%
filter(region == "North") %>%
select(1:5) %>%
group_by(fungicide, region) %>%
summarise(mean = (1-exp(mean))*100,
lower = (1-exp(ci.ub))*100,
upper = (1-exp(ci.lb))*100)## `summarise()` regrouping output by 'fungicide' (override with `.groups` argument)names(North) = c("fungicide", "region", "mean", "lower", "upper")
reg6 = full_join(North,reg5)## Joining, by = c("fungicide", "region", "mean", "lower", "upper")reg6openxlsx::write.xlsx(reg6, here("data","eff_region.xlsx"), colNames = TRUE)We used a factorial-type ANOVA model to determine the significance of the treatment x design interaction, evaluated based on the Wald test statistic.
mv_sev_design <- rma.mv(log_sev, vi_sev,
mods = ~brand_name * design,
random = list(~ 1 | study / design / brand_name),
struct = "HCS",
method = "ML",
control = list(optimizer = "nlm"),
data = rust_sev)
mv_sev_design##
## Multivariate Meta-Analysis Model (k = 1331; method: ML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2365 0.4863 177 no study
## sigma^2.2 0.2365 0.4863 177 no study/design
## sigma^2.3 0.1171 0.3423 1331 no study/design/brand_name
##
## Test for Residual Heterogeneity:
## QE(df = 1313) = 967012.5011, p-val < .0001
##
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 2080.1453, p-val < .0001
##
## Model Results:
##
## estimate se zval pval
## intrcpt 4.1825 0.0966 43.2977 <.0001
## brand_nameAZOX + BENZ -1.0430 0.0649 -16.0801 <.0001
## brand_nameBIXF + TFLX + PROT -1.5484 0.0769 -20.1269 <.0001
## brand_namePICO + BENZ -1.6047 0.0764 -21.0079 <.0001
## brand_namePICO + CYPR -0.8693 0.0620 -14.0103 <.0001
## brand_namePICO + TEBU -1.1680 0.0648 -18.0172 <.0001
## brand_namePYRA + EPOX + FLUX -1.4054 0.0760 -18.4929 <.0001
## brand_nameTFLX + CYPR -0.9127 0.0643 -14.2038 <.0001
## brand_nameTFLX + PROT -1.2438 0.0640 -19.4381 <.0001
## design -0.0039 0.0180 -0.2137 0.8308
## brand_nameAZOX + BENZ:design -0.0683 0.0125 -5.4553 <.0001
## brand_nameBIXF + TFLX + PROT:design 0.0475 0.0248 1.9105 0.0561
## brand_namePICO + BENZ:design 0.1242 0.0241 5.1473 <.0001
## brand_namePICO + CYPR:design 0.0092 0.0117 0.7824 0.4340
## brand_namePICO + TEBU:design 0.0205 0.0120 1.7099 0.0873
## brand_namePYRA + EPOX + FLUX:design 0.0622 0.0247 2.5205 0.0117
## brand_nameTFLX + CYPR:design 0.0118 0.0140 0.8441 0.3986
## brand_nameTFLX + PROT:design -0.0080 0.0123 -0.6509 0.5151
## ci.lb ci.ub
## intrcpt 3.9931 4.3718 ***
## brand_nameAZOX + BENZ -1.1701 -0.9159 ***
## brand_nameBIXF + TFLX + PROT -1.6992 -1.3977 ***
## brand_namePICO + BENZ -1.7544 -1.4550 ***
## brand_namePICO + CYPR -0.9909 -0.7477 ***
## brand_namePICO + TEBU -1.2950 -1.0409 ***
## brand_namePYRA + EPOX + FLUX -1.5544 -1.2565 ***
## brand_nameTFLX + CYPR -1.0387 -0.7868 ***
## brand_nameTFLX + PROT -1.3692 -1.1184 ***
## design -0.0392 0.0315
## brand_nameAZOX + BENZ:design -0.0929 -0.0438 ***
## brand_nameBIXF + TFLX + PROT:design -0.0012 0.0961 .
## brand_namePICO + BENZ:design 0.0769 0.1715 ***
## brand_namePICO + CYPR:design -0.0138 0.0321
## brand_namePICO + TEBU:design -0.0030 0.0440 .
## brand_namePYRA + EPOX + FLUX:design 0.0138 0.1106 *
## brand_nameTFLX + CYPR:design -0.0157 0.0393
## brand_nameTFLX + PROT:design -0.0320 0.0161
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mv_sev_design, btt = 11:18)##
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 94.2196, p-val < .0001