Data

library(readr)
library(metafor)
library(tidyverse)
library(here)
rust_yld <- read_csv("data/dat-yld.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(),
##   vi_yld = col_double(),
##   n2 = col_double(),
##   n3 = col_double(),
##   design = col_double()
## )

Absolute yield

library(metafor)

mv_yld <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name,
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld
)

summary(mv_yld)
## 
## Multivariate Meta-Analysis Model (k = 1314; method: ML)
## 
##     logLik    Deviance         AIC         BIC        AICc 
## -9215.4624   4054.2612  18538.9248  18818.6897  18543.6429   
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 175)
## inner factor: brand_name (nlvls = 9)
## 
##                  estim      sqrt  k.lvl  fixed               level 
## tau^2.1    636726.0927  797.9512    175     no             AACHECK 
## tau^2.2    645330.4099  803.3246    145     no         AZOX + BENZ 
## tau^2.3    779999.1995  883.1756    116     no  BIXF + TFLX + PROT 
## tau^2.4    676953.2688  822.7717    118     no         PICO + BENZ 
## tau^2.5    716365.1912  846.3836    175     no         PICO + CYPR 
## tau^2.6    717892.6609  847.2855    154     no         PICO + TEBU 
## tau^2.7    745271.0498  863.2908    116     no  PYRA + EPOX + FLUX 
## tau^2.8    726161.0416  852.1508    146     no         TFLX + CYPR 
## tau^2.9    703865.3637  838.9668    169     no         TFLX + PROT 
## 
##                     rho.AACH  rho.AZ+B  rho.B+T+P  rho.PI+B  rho.PI+C  rho.PI+T 
## AACHECK                    1    0.8246     0.8342    0.8355    0.9335    0.9119 
## AZOX + BENZ           0.8246         1     0.9116    0.9743    0.9132    0.9103 
## BIXF + TFLX + PROT    0.8342    0.9116          1    0.9205    0.9393    0.9322 
## PICO + BENZ           0.8355    0.9743     0.9205         1    0.9281    0.9263 
## PICO + CYPR           0.9335    0.9132     0.9393    0.9281         1    0.9860 
## PICO + TEBU           0.9119    0.9103     0.9322    0.9263    0.9860         1 
## PYRA + EPOX + FLUX    0.8470    0.9540     0.9744    0.9590    0.9357    0.9312 
## TFLX + CYPR           0.9357    0.9216     0.9427    0.9268    0.9958    0.9834 
## TFLX + PROT           0.8834    0.9209     0.9703    0.9306    0.9767    0.9670 
##                     rho.P+E+F  rho.TF+C  rho.TF+P    AACH  AZ+B  B+T+P  PI+B 
## AACHECK                0.8470    0.9357    0.8834       -    no     no    no 
## AZOX + BENZ            0.9540    0.9216    0.9209     145     -     no    no 
## BIXF + TFLX + PROT     0.9744    0.9427    0.9703     116    87      -    no 
## PICO + BENZ            0.9590    0.9268    0.9306     118    88    116     - 
## PICO + CYPR            0.9357    0.9958    0.9767     175   145    116   118 
## PICO + TEBU            0.9312    0.9834    0.9670     154   124     95    97 
## PYRA + EPOX + FLUX          1    0.9455    0.9475     116    87    114   116 
## TFLX + CYPR            0.9455         1    0.9755     146   116    116   118 
## TFLX + PROT            0.9475    0.9755         1     169   145    110   112 
##                     PI+C  PI+T  P+E+F  TF+C  TF+P 
## AACHECK               no    no     no    no    no 
## AZOX + BENZ           no    no     no    no    no 
## BIXF + TFLX + PROT    no    no     no    no    no 
## PICO + BENZ           no    no     no    no    no 
## PICO + CYPR            -    no     no    no    no 
## PICO + TEBU          154     -     no    no    no 
## PYRA + EPOX + FLUX   116    95      -    no    no 
## TFLX + CYPR          146   125    116     -    no 
## TFLX + PROT          169   148    111   140     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 1305) = 173239.6764, p-val < .0001
## 
## Test of Moderators (coefficients 2:9):
## QM(df = 8) = 906.6343, p-val < .0001
## 
## Model Results:
## 
##                                estimate       se     zval    pval      ci.lb 
## intrcpt                       2460.9764  60.8511  40.4426  <.0001  2341.7103 
## brand_nameAZOX + BENZ          910.6261  38.2758  23.7912  <.0001   835.6070 
## brand_nameBIXF + TFLX + PROT  1080.0092  40.9167  26.3953  <.0001   999.8139 
## brand_namePICO + BENZ         1010.9269  38.5883  26.1978  <.0001   935.2953 
## brand_namePICO + CYPR          600.3912  25.2691  23.7599  <.0001   550.8646 
## brand_namePICO + TEBU          682.1111  28.9720  23.5438  <.0001   625.3270 
## brand_namePYRA + EPOX + FLUX   981.5161  38.9349  25.2092  <.0001   905.2050 
## brand_nameTFLX + CYPR          646.4323  25.3750  25.4752  <.0001   596.6983 
## brand_nameTFLX + PROT          891.0381  32.0831  27.7728  <.0001   828.1563 
##                                   ci.ub 
## intrcpt                       2580.2424  *** 
## brand_nameAZOX + BENZ          985.6452  *** 
## brand_nameBIXF + TFLX + PROT  1160.2044  *** 
## brand_namePICO + BENZ         1086.5585  *** 
## brand_namePICO + CYPR          649.9178  *** 
## brand_namePICO + TEBU          738.8952  *** 
## brand_namePYRA + EPOX + FLUX  1057.8271  *** 
## brand_nameTFLX + CYPR          696.1663  *** 
## brand_nameTFLX + PROT          953.9199  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
yield_res<- data.frame(cbind(mv_yld$b, 
                             mv_yld$ci.lb,
                             mv_yld$ci.ub)) %>% 
  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")


names (yield_res) = c("yld", "yld_lower", "yld_upper", "fungicide")
  
yield_res
#write_csv(yield_res, "data/yld_kg.csv")
openxlsx::write.xlsx(yield_res, here("data","yld_kg.xlsx"), colNames = TRUE)

Contrasts

We can set linear contrasts between treatments of interest and get the P-valued using the anova function.

anova(mv_yld, 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:  -169.3830 33.3077  -5.0854 <.0001 
## 2:  -100.3007 22.8112  -4.3970 <.0001 
## 3:   310.2349 29.1941  10.6266 <.0001 
## 4:   228.5151 30.0648   7.6008 <.0001 
## 5:   -70.8899 27.2375  -2.6027 0.0093 
## 6:   264.1939 28.4203   9.2959 <.0001 
## 7:    19.5880 28.1047   0.6970 0.4858 
## 8:    69.0823 32.6511   2.1158 0.0344 
## 9:   479.6179 28.4152  16.8789 <.0001 
## 10:  397.8981 30.0046  13.2612 <.0001 
## 11:   98.4931 21.0235   4.6849 <.0001 
## 12:  433.5769 28.1232  15.4171 <.0001 
## 13:  188.9710 23.3893   8.0794 <.0001 
## 14:  410.5357 28.4939  14.4079 <.0001 
## 15:  328.8158 29.3148  11.2167 <.0001 
## 16:   29.4108 25.5630   1.1505 0.2499 
## 17:  364.4946 29.2291  12.4703 <.0001 
## 18:  119.8888 28.2063   4.2504 <.0001 
## 19:  -81.7199 15.4180  -5.3003 <.0001 
## 20: -381.1248 28.5673 -13.3413 <.0001 
## 21:  -46.0411 11.8314  -3.8914 <.0001 
## 22: -290.6469 17.2112 -16.8871 <.0001 
## 23: -299.4050 29.7364 -10.0686 <.0001 
## 24:   35.6788 16.7421   2.1311 0.0331 
## 25: -208.9271 20.1956 -10.3452 <.0001 
## 26:  335.0838 27.1044  12.3627 <.0001 
## 27:   90.4779 27.0003   3.3510 0.0008 
## 28: -244.6059 18.0530 -13.5493 <.0001

Moderator analysis

Mods: year

mv_yld_year <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.numeric(year),
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld %>% mutate(year= year - 2015))

mv_yld_year
## 
## Multivariate Meta-Analysis Model (k = 1314; method: ML)
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 175)
## inner factor: brand_name (nlvls = 9)
## 
##                  estim      sqrt  k.lvl  fixed               level 
## tau^2.1    620332.6378  787.6120    175     no             AACHECK 
## tau^2.2    644840.3574  803.0195    145     no         AZOX + BENZ 
## tau^2.3    766684.4458  875.6052    116     no  BIXF + TFLX + PROT 
## tau^2.4    675455.2431  821.8608    118     no         PICO + BENZ 
## tau^2.5    701904.5140  837.7974    175     no         PICO + CYPR 
## tau^2.6    700466.0081  836.9385    154     no         PICO + TEBU 
## tau^2.7    741471.5668  861.0874    116     no  PYRA + EPOX + FLUX 
## tau^2.8    710636.6290  842.9927    146     no         TFLX + CYPR 
## tau^2.9    693032.8659  832.4860    169     no         TFLX + PROT 
## 
##                     rho.AACH  rho.AZ+B  rho.B+T+P  rho.PI+B  rho.PI+C  rho.PI+T 
## AACHECK                    1    0.8504     0.8305    0.8509    0.9325    0.9139 
## AZOX + BENZ           0.8504         1     0.9240    0.9770    0.9318    0.9383 
## BIXF + TFLX + PROT    0.8305    0.9240          1    0.9318    0.9391    0.9304 
## PICO + BENZ           0.8509    0.9770     0.9318         1    0.9395    0.9477 
## PICO + CYPR           0.9325    0.9318     0.9391    0.9395         1    0.9878 
## PICO + TEBU           0.9139    0.9383     0.9304    0.9477    0.9878         1 
## PYRA + EPOX + FLUX    0.8491    0.9592     0.9759    0.9635    0.9381    0.9365 
## TFLX + CYPR           0.9345    0.9413     0.9423    0.9389    0.9960    0.9843 
## TFLX + PROT           0.8818    0.9362     0.9701    0.9415    0.9763    0.9691 
##                     rho.P+E+F  rho.TF+C  rho.TF+P    AACH  AZ+B  B+T+P  PI+B 
## AACHECK                0.8491    0.9345    0.8818       -    no     no    no 
## AZOX + BENZ            0.9592    0.9413    0.9362     145     -     no    no 
## BIXF + TFLX + PROT     0.9759    0.9423    0.9701     116    87      -    no 
## PICO + BENZ            0.9635    0.9389    0.9415     118    88    116     - 
## PICO + CYPR            0.9381    0.9960    0.9763     175   145    116   118 
## PICO + TEBU            0.9365    0.9843    0.9691     154   124     95    97 
## PYRA + EPOX + FLUX          1    0.9481    0.9490     116    87    114   116 
## TFLX + CYPR            0.9481         1    0.9752     146   116    116   118 
## TFLX + PROT            0.9490    0.9752         1     169   145    110   112 
##                     PI+C  PI+T  P+E+F  TF+C  TF+P 
## AACHECK               no    no     no    no    no 
## AZOX + BENZ           no    no     no    no    no 
## BIXF + TFLX + PROT    no    no     no    no    no 
## PICO + BENZ           no    no     no    no    no 
## PICO + CYPR            -    no     no    no    no 
## PICO + TEBU          154     -     no    no    no 
## PYRA + EPOX + FLUX   116    95      -    no    no 
## TFLX + CYPR          146   125    116     -    no 
## TFLX + PROT          169   148    111   140     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 1296) = 172423.1537, p-val < .0001
## 
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 1016.1004, p-val < .0001
## 
## Model Results:
## 
##                                                 estimate        se     zval 
## intrcpt                                        2270.5064  107.4955  21.1219 
## brand_nameAZOX + BENZ                          1161.6923   62.4000  18.6169 
## brand_nameBIXF + TFLX + PROT                   1025.1549   94.4270  10.8566 
## brand_namePICO + BENZ                          1336.6621   82.2547  16.2503 
## brand_namePICO + CYPR                           619.3209   45.1394  13.7202 
## brand_namePICO + TEBU                           627.3148   50.7278  12.3663 
## brand_namePYRA + EPOX + FLUX                   1055.0337   91.0230  11.5909 
## brand_nameTFLX + CYPR                           655.3253   47.8614  13.6922 
## brand_nameTFLX + PROT                           927.5425   57.2384  16.2049 
## as.numeric(year)                                 80.0760   37.4522   2.1381 
## brand_nameAZOX + BENZ:as.numeric(year)         -109.8643   22.3741  -4.9103 
## brand_nameBIXF + TFLX + PROT:as.numeric(year)     8.8058   29.9421   0.2941 
## brand_namePICO + BENZ:as.numeric(year)         -119.9367   26.4036  -4.5424 
## brand_namePICO + CYPR:as.numeric(year)           -7.9894   15.8947  -0.5026 
## brand_namePICO + TEBU:as.numeric(year)           30.4794   18.5042   1.6472 
## brand_namePYRA + EPOX + FLUX:as.numeric(year)   -37.8911   28.6793  -1.3212 
## brand_nameTFLX + CYPR:as.numeric(year)           -4.8123   16.5700  -0.2904 
## brand_nameTFLX + PROT:as.numeric(year)          -15.6598   20.1060  -0.7789 
##                                                  pval      ci.lb      ci.ub 
## intrcpt                                        <.0001  2059.8191  2481.1937 
## brand_nameAZOX + BENZ                          <.0001  1039.3905  1283.9941 
## brand_nameBIXF + TFLX + PROT                   <.0001   840.0814  1210.2284 
## brand_namePICO + BENZ                          <.0001  1175.4459  1497.8783 
## brand_namePICO + CYPR                          <.0001   530.8492   707.7925 
## brand_namePICO + TEBU                          <.0001   527.8901   726.7395 
## brand_namePYRA + EPOX + FLUX                   <.0001   876.6319  1233.4354 
## brand_nameTFLX + CYPR                          <.0001   561.5187   749.1319 
## brand_nameTFLX + PROT                          <.0001   815.3573  1039.7278 
## as.numeric(year)                               0.0325     6.6711   153.4809 
## brand_nameAZOX + BENZ:as.numeric(year)         <.0001  -153.7168   -66.0119 
## brand_nameBIXF + TFLX + PROT:as.numeric(year)  0.7687   -49.8797    67.4912 
## brand_namePICO + BENZ:as.numeric(year)         <.0001  -171.6868   -68.1867 
## brand_namePICO + CYPR:as.numeric(year)         0.6152   -39.1424    23.1636 
## brand_namePICO + TEBU:as.numeric(year)         0.0995    -5.7881    66.7470 
## brand_namePYRA + EPOX + FLUX:as.numeric(year)  0.1864   -94.1015    18.3193 
## brand_nameTFLX + CYPR:as.numeric(year)         0.7715   -37.2889    27.6643 
## brand_nameTFLX + PROT:as.numeric(year)         0.4361   -55.0668    23.7473 
##  
## 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.numeric(year)                                 * 
## brand_nameAZOX + BENZ:as.numeric(year)         *** 
## brand_nameBIXF + TFLX + PROT:as.numeric(year) 
## brand_namePICO + BENZ:as.numeric(year)         *** 
## brand_namePICO + CYPR:as.numeric(year) 
## brand_namePICO + TEBU:as.numeric(year)           . 
## brand_namePYRA + EPOX + FLUX:as.numeric(year) 
## brand_nameTFLX + CYPR:as.numeric(year) 
## brand_nameTFLX + PROT:as.numeric(year) 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mv_yld_year, btt = 11:18)
## 
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 85.4079, p-val < .0001

Decline yield

yld_gain <- rust_yld %>%
  mutate(gain = mean_yld - yld_check) %>% 
  filter(brand_name!= "AACHECK")  
  

yld_gain <- yld_gain %>% 
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" ~ "South",
    state == "SP" ~ "South",
    state == "RS" ~ "South")) 


reg1 = data.frame(mv_yld_year$beta, mv_yld_year$ci.lb, mv_yld_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"))

year = seq(0,6, by = 0.1) 
fungicide = NULL
year_col = NULL
for(i in 1:length(data_model$fungicide)){
data_cache = yld_gain %>% 
    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_gain = intercept_mean + slope_mean*year,
         CIL = intercept_lower + slope_lower*year,
         CIU = intercept_upper + slope_upper*year,
         year = year+2015) %>% 
  mutate(brand_name = fungicide) %>% 
  filter(year <2020.2) %>% 
  dplyr::select(-fungicide)
predicted
openxlsx::write.xlsx(predicted, here("data","predicted_yield.xlsx"), colNames = TRUE)

AZOX + BENZ

library(janitor)
# create the log of the sev variable
rust_yld1 <- rust_yld %>%
  filter(brand_name %in% c("AACHECK", "AZOX + BENZ")) %>% 
  group_by(study) %>% 
  mutate(n2 = n()) %>% 
  filter(n2 != 1)

rust_yld1 = rust_yld1 %>% 
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_yld1 %>%
tabyl(brand_name, region)
rust_yld1 %>%
tabyl(brand_name, state)

NW

library(metafor)
yld_gain <- rust_yld %>%
  mutate(gain = mean_yld - yld_check) %>% 
  filter(brand_name!= "AACHECK")  

rust_yld_nor <- rust_yld1 %>% 
filter(region == "North")

mv_yld_year <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.numeric(year),
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld_nor %>% mutate(year= year - 2015))

mv_yld_year
## 
## Multivariate Meta-Analysis Model (k = 208; method: ML)
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 104)
## inner factor: brand_name (nlvls = 2)
## 
##                  estim      sqrt  k.lvl  fixed        level 
## tau^2.1    404762.2394  636.2093    104     no      AACHECK 
## tau^2.2    407851.7983  638.6328    104     no  AZOX + BENZ 
## 
##              rho.AACH  rho.AZ+B    AACH  AZ+B 
## AACHECK             1    0.7951       -    no 
## AZOX + BENZ    0.7951         1     104     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 204) = 13868.9184, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 463.5621, p-val < .0001
## 
## Model Results:
## 
##                                          estimate        se     zval    pval 
## intrcpt                                 2093.6094  106.3863  19.6793  <.0001 
## brand_nameAZOX + BENZ                   1162.1697   71.7404  16.1997  <.0001 
## as.numeric(year)                         145.3687   39.3894   3.6906  0.0002 
## brand_nameAZOX + BENZ:as.numeric(year)  -126.7211   26.6474  -4.7555  <.0001 
##                                             ci.lb      ci.ub 
## intrcpt                                 1885.0962  2302.1226  *** 
## brand_nameAZOX + BENZ                   1021.5611  1302.7783  *** 
## as.numeric(year)                          68.1669   222.5706  *** 
## brand_nameAZOX + BENZ:as.numeric(year)  -178.9490   -74.4933  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
reg1 = data.frame(mv_yld_year$beta, mv_yld_year$ci.lb, mv_yld_year$ci.ub, mv_yld_year$se) %>%
  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", 1), rep("slope", 1))) %>% 
  mutate(region = rep("North", length(mod))) %>% 
  select(-lado2) 
names(reg1) = c("fungicide", "mean", "ci.lb", "ci.ub", "SE", "mod", "region") 

elatus_norte = reg1

SE

library(metafor)
yld_gain <- rust_yld %>%
  mutate(gain = mean_yld - yld_check) %>% 
  filter(brand_name!= "AACHECK")  

rust_yld_south <- rust_yld1 %>% 
filter(region == "South")

mv_yld_year <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.numeric(year),
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld_south %>% mutate(year= year - 2015))

mv_yld_year
## 
## Multivariate Meta-Analysis Model (k = 82; method: ML)
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 41)
## inner factor: brand_name (nlvls = 2)
## 
##                  estim      sqrt  k.lvl  fixed        level 
## tau^2.1    521876.7119  722.4103     41     no      AACHECK 
## tau^2.2    636887.6651  798.0524     41     no  AZOX + BENZ 
## 
##              rho.AACH  rho.AZ+B    AACH  AZ+B 
## AACHECK             1    0.8587       -    no 
## AZOX + BENZ    0.8587         1      41     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 78) = 9033.2665, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 211.6400, p-val < .0001
## 
## Model Results:
## 
##                                          estimate        se     zval    pval 
## intrcpt                                 2443.3541  180.4830  13.5379  <.0001 
## brand_nameAZOX + BENZ                   1217.7822  108.3954  11.2346  <.0001 
## as.numeric(year)                         174.8334   81.1147   2.1554  0.0311 
## brand_nameAZOX + BENZ:as.numeric(year)  -149.9365   49.6502  -3.0199  0.0025 
##                                             ci.lb      ci.ub 
## intrcpt                                 2089.6139  2797.0942  *** 
## brand_nameAZOX + BENZ                   1005.3311  1430.2332  *** 
## as.numeric(year)                          15.8516   333.8152    * 
## brand_nameAZOX + BENZ:as.numeric(year)  -247.2492   -52.6238   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
reg1 = data.frame(mv_yld_year$beta, mv_yld_year$ci.lb, mv_yld_year$ci.ub, mv_yld_year$se) %>%
  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", 1), rep("slope", 1))) %>% 
  mutate(region = rep("South", length(mod))) %>% 
  select(-lado2) 
names(reg1) = c("fungicide", "mean", "ci.lb", "ci.ub", "SE", "mod", "region") 

elatus_sul = reg1

PICO + BENZ

library(janitor)
# create the log of the sev variable
rust_yld1 <- rust_yld %>%
  filter(brand_name %in% c("AACHECK", "PICO + BENZ")) %>% 
  group_by(study) %>% 
  mutate(n2 = n()) %>% 
  filter(n2 != 1)

rust_yld1 = rust_yld1 %>% 
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_yld1 %>%
tabyl(brand_name, region)
rust_yld1 %>%
tabyl(brand_name, state)

NW

library(metafor)
yld_gain <- rust_yld %>%
  mutate(gain = mean_yld - yld_check) %>% 
  filter(brand_name!= "AACHECK")  

rust_yld_nor <- rust_yld1 %>% 
filter(region == "North")

mv_yld_year <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.numeric(year),
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld_nor %>% mutate(year= year - 2015))

mv_yld_year
## 
## Multivariate Meta-Analysis Model (k = 176; method: ML)
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 88)
## inner factor: brand_name (nlvls = 2)
## 
##                  estim      sqrt  k.lvl  fixed        level 
## tau^2.1    694041.3466  833.0914     88     no      AACHECK 
## tau^2.2    798894.3137  893.8089     88     no  PICO + BENZ 
## 
##              rho.AACH  rho.PI+B    AACH  PI+B 
## AACHECK             1    0.8468       -    no 
## PICO + BENZ    0.8468         1      88     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 172) = 22658.9911, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 336.0714, p-val < .0001
## 
## Model Results:
## 
##                                          estimate        se     zval    pval 
## intrcpt                                 1779.3432  290.0553   6.1345  <.0001 
## brand_namePICO + BENZ                   1486.4364  173.8693   8.5492  <.0001 
## as.numeric(year)                         179.3583   82.3555   2.1779  0.0294 
## brand_namePICO + BENZ:as.numeric(year)  -154.9267   49.5106  -3.1292  0.0018 
##                                             ci.lb      ci.ub 
## intrcpt                                 1210.8453  2347.8411  *** 
## brand_namePICO + BENZ                   1145.6589  1827.2139  *** 
## as.numeric(year)                          17.9444   340.7721    * 
## brand_namePICO + BENZ:as.numeric(year)  -251.9657   -57.8877   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
reg1 = data.frame(mv_yld_year$beta, mv_yld_year$ci.lb, mv_yld_year$ci.ub, mv_yld_year$se) %>%
  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", 1), rep("slope", 1))) %>% 
  mutate(region = rep("North", length(mod))) %>% 
  select(-lado2) 
names(reg1) = c("fungicide", "mean", "ci.lb", "ci.ub", "SE", "mod", "region") 

vess_norte = reg1 

SE

library(metafor)
yld_gain <- rust_yld %>%
  mutate(gain = mean_yld - yld_check) %>% 
  filter(brand_name!= "AACHECK")  

rust_yld_south <- rust_yld1 %>% 
filter(region == "South")

mv_yld_year <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.numeric(year),
  random = list(~brand_name | study),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld_south %>% mutate(year= year - 2015))

mv_yld_year
## 
## Multivariate Meta-Analysis Model (k = 60; method: ML)
## 
## Variance Components:
## 
## outer factor: study      (nlvls = 30)
## inner factor: brand_name (nlvls = 2)
## 
##                  estim      sqrt  k.lvl  fixed        level 
## tau^2.1    518346.7061  719.9630     30     no      AACHECK 
## tau^2.2    535819.5912  731.9970     30     no  PICO + BENZ 
## 
##              rho.AACH  rho.PI+B    AACH  PI+B 
## AACHECK             1    0.8879       -    no 
## PICO + BENZ    0.8879         1      30     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 56) = 4877.1263, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 160.2478, p-val < .0001
## 
## Model Results:
## 
##                                          estimate        se     zval    pval 
## intrcpt                                 3528.7061  441.8376   7.9864  <.0001 
## brand_namePICO + BENZ                    699.2209  230.0585   3.0393  0.0024 
## as.numeric(year)                        -167.8161  134.7084  -1.2458  0.2128 
## brand_namePICO + BENZ:as.numeric(year)    52.9786   70.8458   0.7478  0.4546 
##                                             ci.lb      ci.ub 
## intrcpt                                 2662.7202  4394.6919  *** 
## brand_namePICO + BENZ                    248.3145  1150.1274   ** 
## as.numeric(year)                        -431.8396    96.2074      
## brand_namePICO + BENZ:as.numeric(year)   -85.8766   191.8338      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
reg1 = data.frame(mv_yld_year$beta, mv_yld_year$ci.lb, mv_yld_year$ci.ub, mv_yld_year$se) %>%
  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", 1), rep("slope", 1))) %>% 
  mutate(region = rep("South", length(mod))) %>% 
  select(-lado2) 
names(reg1) = c("fungicide", "mean", "ci.lb", "ci.ub", "SE", "mod", "region") 

vess_sul = reg1 

data_two_fung = rbind(elatus_norte, elatus_sul, vess_norte, vess_sul)

openxlsx::write.xlsx(data_two_fung, here("data","Elatus_Vessarya_yld_simulation.xlsx"), colNames = TRUE)

Mods: region

rust_yld1 = rust_yld %>% 
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_yld1 %>%
tabyl(brand_name, region)
mv_yld_region <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.factor(region),
  random = list(~factor(brand_name) | factor(study)),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld1)

mv_yld_region
## 
## Multivariate Meta-Analysis Model (k = 1314; method: ML)
## 
## Variance Components:
## 
## outer factor: factor(study)      (nlvls = 175)
## inner factor: factor(brand_name) (nlvls = 9)
## 
##                  estim      sqrt  k.lvl  fixed               level 
## tau^2.1    606007.6711  778.4649    175     no             AACHECK 
## tau^2.2    601937.4814  775.8463    145     no         AZOX + BENZ 
## tau^2.3    721679.6495  849.5173    116     no  BIXF + TFLX + PROT 
## tau^2.4    655626.7851  809.7078    118     no         PICO + BENZ 
## tau^2.5    653264.0342  808.2475    175     no         PICO + CYPR 
## tau^2.6    661692.6217  813.4449    154     no         PICO + TEBU 
## tau^2.7    720031.9783  848.5470    116     no  PYRA + EPOX + FLUX 
## tau^2.8    656849.4796  810.4625    146     no         TFLX + CYPR 
## tau^2.9    628275.7729  792.6385    169     no         TFLX + PROT 
## 
##                     rho.AACH  rho.AZ+B  rho.B+T+P  rho.PI+B  rho.PI+C  rho.PI+T 
## AACHECK                    1    0.8142     0.8273    0.8341    0.9320    0.9082 
## AZOX + BENZ           0.8142         1     0.9057    0.9715    0.9069    0.9038 
## BIXF + TFLX + PROT    0.8273    0.9057          1    0.9198    0.9353    0.9290 
## PICO + BENZ           0.8341    0.9715     0.9198         1    0.9324    0.9317 
## PICO + CYPR           0.9320    0.9069     0.9353    0.9324         1    0.9849 
## PICO + TEBU           0.9082    0.9038     0.9290    0.9317    0.9849         1 
## PYRA + EPOX + FLUX    0.8453    0.9502     0.9749    0.9580    0.9384    0.9355 
## TFLX + CYPR           0.9350    0.9164     0.9391    0.9321    0.9952    0.9821 
## TFLX + PROT           0.8801    0.9169     0.9685    0.9349    0.9747    0.9651 
##                     rho.P+E+F  rho.TF+C  rho.TF+P    AACH  AZ+B  B+T+P  PI+B 
## AACHECK                0.8453    0.9350    0.8801       -    no     no    no 
## AZOX + BENZ            0.9502    0.9164    0.9169     145     -     no    no 
## BIXF + TFLX + PROT     0.9749    0.9391    0.9685     116    87      -    no 
## PICO + BENZ            0.9580    0.9321    0.9349     118    88    116     - 
## PICO + CYPR            0.9384    0.9952    0.9747     175   145    116   118 
## PICO + TEBU            0.9355    0.9821    0.9651     154   124     95    97 
## PYRA + EPOX + FLUX          1    0.9499    0.9504     116    87    114   116 
## TFLX + CYPR            0.9499         1    0.9728     146   116    116   118 
## TFLX + PROT            0.9504    0.9728         1     169   145    110   112 
##                     PI+C  PI+T  P+E+F  TF+C  TF+P 
## AACHECK               no    no     no    no    no 
## AZOX + BENZ           no    no     no    no    no 
## BIXF + TFLX + PROT    no    no     no    no    no 
## PICO + BENZ           no    no     no    no    no 
## PICO + CYPR            -    no     no    no    no 
## PICO + TEBU          154     -     no    no    no 
## PYRA + EPOX + FLUX   116    95      -    no    no 
## TFLX + CYPR          146   125    116     -    no 
## TFLX + PROT          169   148    111   140     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 1296) = 155126.5694, p-val < .0001
## 
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 977.5113, p-val < .0001
## 
## Model Results:
## 
##                                                       estimate        se 
## intrcpt                                              2351.0840   69.9976 
## brand_nameAZOX + BENZ                                 894.3456   45.1105 
## brand_nameBIXF + TFLX + PROT                         1032.7464   47.2686 
## brand_namePICO + BENZ                                1009.8853   44.6194 
## brand_namePICO + CYPR                                 554.8111   29.0685 
## brand_namePICO + TEBU                                 645.6712   33.8301 
## brand_namePYRA + EPOX + FLUX                          967.0207   44.9301 
## brand_nameTFLX + CYPR                                 593.7882   28.9020 
## brand_nameTFLX + PROT                                 832.3037   36.8105 
## as.factor(region)South                                394.8497  132.2577 
## brand_nameAZOX + BENZ:as.factor(region)South           58.9984   85.1539 
## brand_nameBIXF + TFLX + PROT:as.factor(region)South   150.1983   90.9738 
## brand_namePICO + BENZ:as.factor(region)South          -36.6158   85.3171 
## brand_namePICO + CYPR:as.factor(region)South          159.7342   54.6664 
## brand_namePICO + TEBU:as.factor(region)South          129.9113   63.1823 
## brand_namePYRA + EPOX + FLUX:as.factor(region)South    10.4437   86.4463 
## brand_nameTFLX + CYPR:as.factor(region)South          182.2832   54.5530 
## brand_nameTFLX + PROT:as.factor(region)South          202.9433   69.3671 
##                                                         zval    pval      ci.lb 
## intrcpt                                              33.5881  <.0001  2213.8912 
## brand_nameAZOX + BENZ                                19.8257  <.0001   805.9307 
## brand_nameBIXF + TFLX + PROT                         21.8485  <.0001   940.1016 
## brand_namePICO + BENZ                                22.6333  <.0001   922.4328 
## brand_namePICO + CYPR                                19.0863  <.0001   497.8378 
## brand_namePICO + TEBU                                19.0857  <.0001   579.3654 
## brand_namePYRA + EPOX + FLUX                         21.5228  <.0001   878.9594 
## brand_nameTFLX + CYPR                                20.5449  <.0001   537.1414 
## brand_nameTFLX + PROT                                22.6105  <.0001   760.1565 
## as.factor(region)South                                2.9855  0.0028   135.6293 
## brand_nameAZOX + BENZ:as.factor(region)South          0.6928  0.4884  -107.9002 
## brand_nameBIXF + TFLX + PROT:as.factor(region)South   1.6510  0.0987   -28.1072 
## brand_namePICO + BENZ:as.factor(region)South         -0.4292  0.6678  -203.8342 
## brand_namePICO + CYPR:as.factor(region)South          2.9220  0.0035    52.5900 
## brand_namePICO + TEBU:as.factor(region)South          2.0561  0.0398     6.0762 
## brand_namePYRA + EPOX + FLUX:as.factor(region)South   0.1208  0.9038  -158.9880 
## brand_nameTFLX + CYPR:as.factor(region)South          3.3414  0.0008    75.3613 
## brand_nameTFLX + PROT:as.factor(region)South          2.9256  0.0034    66.9863 
##                                                          ci.ub 
## intrcpt                                              2488.2767  *** 
## brand_nameAZOX + BENZ                                 982.7605  *** 
## brand_nameBIXF + TFLX + PROT                         1125.3911  *** 
## brand_namePICO + BENZ                                1097.3378  *** 
## brand_namePICO + CYPR                                 611.7843  *** 
## brand_namePICO + TEBU                                 711.9769  *** 
## brand_namePYRA + EPOX + FLUX                         1055.0819  *** 
## brand_nameTFLX + CYPR                                 650.4350  *** 
## brand_nameTFLX + PROT                                 904.4509  *** 
## as.factor(region)South                                654.0701   ** 
## brand_nameAZOX + BENZ:as.factor(region)South          225.8970      
## brand_nameBIXF + TFLX + PROT:as.factor(region)South   328.5038    . 
## brand_namePICO + BENZ:as.factor(region)South          130.6027      
## brand_namePICO + CYPR:as.factor(region)South          266.8784   ** 
## brand_namePICO + TEBU:as.factor(region)South          253.7463    * 
## brand_namePYRA + EPOX + FLUX:as.factor(region)South   179.8754      
## brand_nameTFLX + CYPR:as.factor(region)South          289.2051  *** 
## brand_nameTFLX + PROT:as.factor(region)South          338.9004   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mv_yld_region, btt = 11:18)
## 
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 34.0610, p-val < .0001
reg1 = data.frame(mv_yld_region$beta, mv_yld_region$ci.lb, mv_yld_region$ci.ub, mv_yld_region$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_yld", "ci.lb_yld", "ci.ub_yld", "SE_yld","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","yld_SE_region.xlsx"), colNames = TRUE)
eff_log = read_excel("data/eff_log_region.xlsx")
yld_reg <- read_excel("data/yld_SE_region.xlsx")

data_simulation = full_join(eff_log, yld_reg, by = c("fungicide", "region"))
openxlsx::write.xlsx(data_simulation, here("data","data_simulation.xlsx"), colNames = TRUE)

Mod effect

reg1 = data.frame(mv_yld_region$beta, mv_yld_region$ci.lb, mv_yld_region$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", "class", "mean", "ci.lb", "ci.ub", "n") 

reg2 = reg1 %>% 
  filter(n < 9) %>% 
  mutate(class = rep("North", length(fungicide)))

reg3 = reg1 %>% 
  filter(n > 8) %>% 
  mutate(class = rep("South", length(fungicide)))

reg4 = rbind(reg2,reg3)  

mean = reg4%>% 
  group_by(fungicide) %>% 
  select(1:3) %>% 
  spread(class, mean) %>% 
  mutate(mean = North + South) %>% 
  select(1,4)

lower = reg4%>% 
  group_by(fungicide) %>% 
  select(1,2,4) %>% 
  spread(class, ci.lb) %>% 
  mutate(lower = North + South) %>%  
  select(1,4)

upper = reg4%>% 
  group_by(fungicide) %>% 
  select(1,2,5) %>% 
  spread(class, ci.ub) %>% 
  mutate(upper = North + South) %>% 
  select(1,4)

reg5 = left_join(mean, lower, by= c("fungicide")) %>% 
  left_join(upper, by = c("fungicide")) %>% 
  mutate(class = rep("South", length("fungicide"))) %>% 
  select("fungicide", "class", "mean", "lower", "upper")


North = reg4 %>% 
  filter(class == "North") %>% 
  select(1:5)
names(North) = c("fungicide", "class", "mean", "lower", "upper") 

reg6 = full_join(North,reg5)
## Joining, by = c("fungicide", "class", "mean", "lower", "upper")
reg6
openxlsx::write.xlsx(reg6, here("data","yld_region.xlsx"), colNames = TRUE)

Mods: severity baseline

library(tidyverse)

rust_yld <- rust_yld %>%
    mutate(
    sev_check_class = case_when(
      sev_check < 70 ~ "Low",
      sev_check >= 70 ~ "High"))

rust_yld %>%
tabyl(brand_name, sev_check_class)
mv_yld_sev_class <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * as.factor(sev_check_class) ,
  random = list(~factor(brand_name) | factor(study)),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld)

mv_yld_sev_class
## 
## Multivariate Meta-Analysis Model (k = 1314; method: ML)
## 
## Variance Components:
## 
## outer factor: factor(study)      (nlvls = 175)
## inner factor: factor(brand_name) (nlvls = 9)
## 
##                  estim      sqrt  k.lvl  fixed               level 
## tau^2.1    588997.6993  767.4619    175     no             AACHECK 
## tau^2.2    634330.7301  796.4488    145     no         AZOX + BENZ 
## tau^2.3    762949.9431  873.4701    116     no  BIXF + TFLX + PROT 
## tau^2.4    668811.0794  817.8087    118     no         PICO + BENZ 
## tau^2.5    691060.7120  831.3006    175     no         PICO + CYPR 
## tau^2.6    691070.0371  831.3062    154     no         PICO + TEBU 
## tau^2.7    735711.1697  857.7361    116     no  PYRA + EPOX + FLUX 
## tau^2.8    704707.8950  839.4688    146     no         TFLX + CYPR 
## tau^2.9    683393.4981  826.6762    169     no         TFLX + PROT 
## 
##                     rho.AACH  rho.AZ+B  rho.B+T+P  rho.PI+B  rho.PI+C  rho.PI+T 
## AACHECK                    1    0.8255     0.8343    0.8413    0.9336    0.9100 
## AZOX + BENZ           0.8255         1     0.9120    0.9728    0.9121    0.9097 
## BIXF + TFLX + PROT    0.8343    0.9120          1    0.9193    0.9381    0.9314 
## PICO + BENZ           0.8413    0.9728     0.9193         1    0.9289    0.9278 
## PICO + CYPR           0.9336    0.9121     0.9381    0.9289         1    0.9855 
## PICO + TEBU           0.9100    0.9097     0.9314    0.9278    0.9855         1 
## PYRA + EPOX + FLUX    0.8525    0.9568     0.9746    0.9585    0.9365    0.9324 
## TFLX + CYPR           0.9371    0.9210     0.9413    0.9268    0.9958    0.9831 
## TFLX + PROT           0.8829    0.9206     0.9696    0.9306    0.9760    0.9663 
##                     rho.P+E+F  rho.TF+C  rho.TF+P    AACH  AZ+B  B+T+P  PI+B 
## AACHECK                0.8525    0.9371    0.8829       -    no     no    no 
## AZOX + BENZ            0.9568    0.9210    0.9206     145     -     no    no 
## BIXF + TFLX + PROT     0.9746    0.9413    0.9696     116    87      -    no 
## PICO + BENZ            0.9585    0.9268    0.9306     118    88    116     - 
## PICO + CYPR            0.9365    0.9958    0.9760     175   145    116   118 
## PICO + TEBU            0.9324    0.9831    0.9663     154   124     95    97 
## PYRA + EPOX + FLUX          1    0.9455    0.9479     116    87    114   116 
## TFLX + CYPR            0.9455         1    0.9745     146   116    116   118 
## TFLX + PROT            0.9479    0.9745         1     169   145    110   112 
##                     PI+C  PI+T  P+E+F  TF+C  TF+P 
## AACHECK               no    no     no    no    no 
## AZOX + BENZ           no    no     no    no    no 
## BIXF + TFLX + PROT    no    no     no    no    no 
## PICO + BENZ           no    no     no    no    no 
## PICO + CYPR            -    no     no    no    no 
## PICO + TEBU          154     -     no    no    no 
## PYRA + EPOX + FLUX   116    95      -    no    no 
## TFLX + CYPR          146   125    116     -    no 
## TFLX + PROT          169   148    111   140     - 
## 
## Test for Residual Heterogeneity:
## QE(df = 1296) = 167547.7250, p-val < .0001
## 
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 983.1316, p-val < .0001
## 
## Model Results:
## 
##                                                              estimate        se 
## intrcpt                                                     2272.7953   77.4679 
## brand_nameAZOX + BENZ                                       1005.3144   49.6013 
## brand_nameBIXF + TFLX + PROT                                1154.1764   53.9432 
## brand_namePICO + BENZ                                       1119.6025   50.1055 
## brand_namePICO + CYPR                                        651.0735   32.9249 
## brand_namePICO + TEBU                                        728.4738   37.7547 
## brand_namePYRA + EPOX + FLUX                                1088.0312   50.7783 
## brand_nameTFLX + CYPR                                        707.2705   33.1093 
## brand_nameTFLX + PROT                                        956.3937   41.7684 
## as.factor(sev_check_class)Low                                438.0983  118.3458 
## brand_nameAZOX + BENZ:as.factor(sev_check_class)Low         -220.7779   75.6310 
## brand_nameBIXF + TFLX + PROT:as.factor(sev_check_class)Low  -171.4304   81.4070 
## brand_namePICO + BENZ:as.factor(sev_check_class)Low         -250.2076   75.5580 
## brand_namePICO + CYPR:as.factor(sev_check_class)Low         -117.3651   50.2694 
## brand_namePICO + TEBU:as.factor(sev_check_class)Low         -107.1294   58.1481 
## brand_namePYRA + EPOX + FLUX:as.factor(sev_check_class)Low  -237.2981   76.3145 
## brand_nameTFLX + CYPR:as.factor(sev_check_class)Low         -137.6340   50.2045 
## brand_nameTFLX + PROT:as.factor(sev_check_class)Low         -151.6628   63.7592 
##                                                                zval    pval 
## intrcpt                                                     29.3385  <.0001 
## brand_nameAZOX + BENZ                                       20.2679  <.0001 
## brand_nameBIXF + TFLX + PROT                                21.3962  <.0001 
## brand_namePICO + BENZ                                       22.3449  <.0001 
## brand_namePICO + CYPR                                       19.7745  <.0001 
## brand_namePICO + TEBU                                       19.2949  <.0001 
## brand_namePYRA + EPOX + FLUX                                21.4271  <.0001 
## brand_nameTFLX + CYPR                                       21.3617  <.0001 
## brand_nameTFLX + PROT                                       22.8976  <.0001 
## as.factor(sev_check_class)Low                                3.7018  0.0002 
## brand_nameAZOX + BENZ:as.factor(sev_check_class)Low         -2.9191  0.0035 
## brand_nameBIXF + TFLX + PROT:as.factor(sev_check_class)Low  -2.1058  0.0352 
## brand_namePICO + BENZ:as.factor(sev_check_class)Low         -3.3115  0.0009 
## brand_namePICO + CYPR:as.factor(sev_check_class)Low         -2.3347  0.0196 
## brand_namePICO + TEBU:as.factor(sev_check_class)Low         -1.8424  0.0654 
## brand_namePYRA + EPOX + FLUX:as.factor(sev_check_class)Low  -3.1095  0.0019 
## brand_nameTFLX + CYPR:as.factor(sev_check_class)Low         -2.7415  0.0061 
## brand_nameTFLX + PROT:as.factor(sev_check_class)Low         -2.3787  0.0174 
##                                                                 ci.lb 
## intrcpt                                                     2120.9609 
## brand_nameAZOX + BENZ                                        908.0976 
## brand_nameBIXF + TFLX + PROT                                1048.4498 
## brand_namePICO + BENZ                                       1021.3974 
## brand_namePICO + CYPR                                        586.5418 
## brand_namePICO + TEBU                                        654.4759 
## brand_namePYRA + EPOX + FLUX                                 988.5077 
## brand_nameTFLX + CYPR                                        642.3774 
## brand_nameTFLX + PROT                                        874.5292 
## as.factor(sev_check_class)Low                                206.1448 
## brand_nameAZOX + BENZ:as.factor(sev_check_class)Low         -369.0118 
## brand_nameBIXF + TFLX + PROT:as.factor(sev_check_class)Low  -330.9852 
## brand_namePICO + BENZ:as.factor(sev_check_class)Low         -398.2986 
## brand_namePICO + CYPR:as.factor(sev_check_class)Low         -215.8914 
## brand_namePICO + TEBU:as.factor(sev_check_class)Low         -221.0976 
## brand_namePYRA + EPOX + FLUX:as.factor(sev_check_class)Low  -386.8717 
## brand_nameTFLX + CYPR:as.factor(sev_check_class)Low         -236.0329 
## brand_nameTFLX + PROT:as.factor(sev_check_class)Low         -276.6287 
##                                                                 ci.ub 
## intrcpt                                                     2424.6296  *** 
## brand_nameAZOX + BENZ                                       1102.5311  *** 
## brand_nameBIXF + TFLX + PROT                                1259.9031  *** 
## brand_namePICO + BENZ                                       1217.8075  *** 
## brand_namePICO + CYPR                                        715.6051  *** 
## brand_namePICO + TEBU                                        802.4717  *** 
## brand_namePYRA + EPOX + FLUX                                1187.5548  *** 
## brand_nameTFLX + CYPR                                        772.1636  *** 
## brand_nameTFLX + PROT                                       1038.2582  *** 
## as.factor(sev_check_class)Low                                670.0518  *** 
## brand_nameAZOX + BENZ:as.factor(sev_check_class)Low          -72.5439   ** 
## brand_nameBIXF + TFLX + PROT:as.factor(sev_check_class)Low   -11.8756    * 
## brand_namePICO + BENZ:as.factor(sev_check_class)Low         -102.1166  *** 
## brand_namePICO + CYPR:as.factor(sev_check_class)Low          -18.8388    * 
## brand_namePICO + TEBU:as.factor(sev_check_class)Low            6.8388    . 
## brand_namePYRA + EPOX + FLUX:as.factor(sev_check_class)Low   -87.7245   ** 
## brand_nameTFLX + CYPR:as.factor(sev_check_class)Low          -39.2350   ** 
## brand_nameTFLX + PROT:as.factor(sev_check_class)Low          -26.6970    * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mv_yld_sev_class, btt = 11:18)
## 
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 15.1418, p-val = 0.0564

Mod effect

reg1 = data.frame(mv_yld_sev_class$beta, mv_yld_sev_class$ci.lb, mv_yld_sev_class$ci.ub) %>%
  rownames_to_column("trat") %>%
  separate(trat, into = c("lado1", "lado2"), sep = ":") %>%
  separate(lado2, into = c("lixo","lado3"),sep = "sev_check_class") %>% 
  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", "class", "mean", "ci.lb", "ci.ub", "n") 

reg2 = reg1 %>%
  filter(n < 9) %>% 
  mutate(class = rep("High", length(fungicide)))

reg3 = reg1 %>% 
  filter(n > 8) %>% 
  mutate(class = rep("Low", length(fungicide)))

reg4 = rbind(reg2,reg3)  

mean = reg4%>% 
  group_by(fungicide) %>% 
  select(1:3) %>% 
  spread(class, mean) %>% 
  mutate(mean = High + Low) %>% 
  select(1,4)

lower = reg4%>% 
  group_by(fungicide) %>% 
  select(1,2,4) %>% 
  spread(class, ci.lb) %>% 
  mutate(lower = High + Low) %>%  
  select(1,4)

upper = reg4%>% 
  group_by(fungicide) %>% 
  select(1,2,5) %>% 
  spread(class, ci.ub) %>% 
  mutate(upper = High + Low) %>% 
  select(1,4)

reg5 = left_join(mean, lower, by= c("fungicide")) %>% 
  left_join(upper, by = c("fungicide")) %>% 
  mutate(class = rep("Low", length("fungicide"))) %>% 
  select("fungicide", "class", "mean", "lower", "upper")


High = reg4 %>% 
  filter(class == "High") %>% 
  select(1:5)
names(High) = c("fungicide", "class", "mean", "lower", "upper") 

reg6 = full_join(High,reg5)
## Joining, by = c("fungicide", "class", "mean", "lower", "upper")
reg6
openxlsx::write.xlsx(reg6, here("data","yld_baseline.xlsx"), colNames = TRUE)

Inconsistency

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_yld_design <- rma.mv(mean_yld, vi_yld,
  mods = ~brand_name * design,
  random = list(~ 1 | study / design / brand_name),
  struct = "UN",
  method = "ML",
  control = list(optimizer = "nlm"),
  data = rust_yld)

mv_yld_design
## 
## Multivariate Meta-Analysis Model (k = 1314; method: ML)
## 
## Variance Components:
## 
##                  estim      sqrt  nlvls  fixed                   factor 
## sigma^2.1  326626.5655  571.5125    175     no                    study 
## sigma^2.2  326626.5651  571.5125    175     no             study/design 
## sigma^2.3   44199.7132  210.2373   1314     no  study/design/brand_name 
## 
## Test for Residual Heterogeneity:
## QE(df = 1296) = 169160.9325, p-val < .0001
## 
## Test of Moderators (coefficients 2:18):
## QM(df = 17) = 2377.4326, p-val < .0001
## 
## Model Results:
## 
##                                       estimate        se     zval    pval 
## intrcpt                              2593.4693  106.3961  24.3756  <.0001 
## brand_nameAZOX + BENZ                 757.4275   42.5418  17.8043  <.0001 
## brand_nameBIXF + TFLX + PROT          987.2565   49.6496  19.8845  <.0001 
## brand_namePICO + BENZ                1028.2639   47.7568  21.5313  <.0001 
## brand_namePICO + CYPR                 587.2696   41.6389  14.1039  <.0001 
## brand_namePICO + TEBU                 678.9830   42.9294  15.8163  <.0001 
## brand_namePYRA + EPOX + FLUX          911.4023   48.1021  18.9473  <.0001 
## brand_nameTFLX + CYPR                 631.5014   42.8403  14.7408  <.0001 
## brand_nameTFLX + PROT                 855.6228   41.6484  20.5439  <.0001 
## design                                -29.2953   19.6717  -1.4892  0.1364 
## brand_nameAZOX + BENZ:design           40.2080    7.7787   5.1690  <.0001 
## brand_nameBIXF + TFLX + PROT:design    26.6692   16.7424   1.5929  0.1112 
## brand_namePICO + BENZ:design          -42.0625   15.0385  -2.7970  0.0052 
## brand_namePICO + CYPR:design            0.1022    7.7117   0.0133  0.9894 
## brand_namePICO + TEBU:design           -0.7833    7.7578  -0.1010  0.9196 
## brand_namePYRA + EPOX + FLUX:design     3.3399   15.3704   0.2173  0.8280 
## brand_nameTFLX + CYPR:design            2.8844    9.1329   0.3158  0.7521 
## brand_nameTFLX + PROT:design            7.3169    7.7320   0.9463  0.3440 
##                                          ci.lb      ci.ub 
## intrcpt                              2384.9368  2802.0019  *** 
## brand_nameAZOX + BENZ                 674.0471   840.8079  *** 
## brand_nameBIXF + TFLX + PROT          889.9450  1084.5679  *** 
## brand_namePICO + BENZ                 934.6624  1121.8655  *** 
## brand_namePICO + CYPR                 505.6588   668.8804  *** 
## brand_namePICO + TEBU                 594.8428   763.1231  *** 
## brand_namePYRA + EPOX + FLUX          817.1240  1005.6806  *** 
## brand_nameTFLX + CYPR                 547.5360   715.4668  *** 
## brand_nameTFLX + PROT                 773.9934   937.2523  *** 
## design                                -67.8511     9.2606      
## brand_nameAZOX + BENZ:design           24.9619    55.4540  *** 
## brand_nameBIXF + TFLX + PROT:design    -6.1454    59.4838      
## brand_namePICO + BENZ:design          -71.5374   -12.5875   ** 
## brand_namePICO + CYPR:design          -15.0124    15.2168      
## brand_namePICO + TEBU:design          -15.9884    14.4217      
## brand_namePYRA + EPOX + FLUX:design   -26.7854    33.4652      
## brand_nameTFLX + CYPR:design          -15.0157    20.7845      
## brand_nameTFLX + PROT:design           -7.8376    22.4714      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mv_yld_design, btt = 11:18)
## 
## Test of Moderators (coefficients 11:18):
## QM(df = 8) = 55.7813, p-val < .0001
Copyright 2017 Emerson Del Ponte