Chapter 2 Pollen Consumption

Load libraries

Input data

df <- read_csv("pollen.dates.csv", col_types = cols(treatment = col_factor(levels = c("1","2", "3", "4")), pollen.start = col_date(format = "%m/%d/%Y"),start.time = col_time(format = "%H:%M:%S"),pollen.end = col_date(format = "%m/%d/%Y"),end.time = col_time(format = "%H:%M:%S"),colony.start = col_date(format = "%m/%d/%Y")))
df$colony <- as.factor(df$colony)

Visualize

#Pollen consumption over time by treatment, from https://rpubs.com/alecri/review_longitudinal
ggplot(df, aes(pollen.ball.id, whole_dif)) +
  geom_line(aes(group = factor(colony))) +
  geom_smooth() +
  facet_grid(~ treatment) 

#from https://stats.oarc.ucla.edu/r/seminars/repeated-measures-analysis-with-r/
with(df, interaction.plot(pollen.ball.id,treatment,whole_dif))

#from https://stats.oarc.ucla.edu/r/seminars/repeated-measures-analysis-with-r/

par(cex = .6)
xyplot(whole_dif ~ pollen.ball.id, data = df, groups = colony,
       type = "o", panel = panel.superpose)

xyplot(whole_dif ~ pollen.ball.id | treatment, data = df, groups = colony,
       type = "o", panel = panel.superpose)

Make a time^2 variables for quadratic analysis

# Make a square of time 
df$time2 <- df$pollen.ball.id^2


time.quad <- lme(whole_dif~treatment*pollen.ball.id+time2, random = list(colony=pdDiag(~pollen.ball.id)),data=df)
summary(time.quad)

## Linear mixed-effects model fit by REML
##   Data: df 
##         AIC       BIC   logLik
##   -472.1207 -417.1032 248.0604
## 
## Random effects:
##  Formula: ~pollen.ball.id | colony
##  Structure: Diagonal
##         (Intercept) pollen.ball.id  Residual
## StdDev:   0.0764445     0.01823242 0.1474935
## 
## Fixed effects:  whole_dif ~ treatment * pollen.ball.id + time2 
##                                 Value  Std.Error  DF   t-value p-value
## (Intercept)                0.10158486 0.03884055 693  2.615433  0.0091
## treatment2                -0.04149554 0.04972959  31 -0.834424  0.4104
## treatment3                -0.03582851 0.04951518  31 -0.723586  0.4747
## treatment4                -0.03221396 0.05143537  31 -0.626300  0.5357
## pollen.ball.id             0.05455725 0.00733842 693  7.434473  0.0000
## time2                     -0.00128180 0.00015513 693 -8.262602  0.0000
## treatment2:pollen.ball.id -0.00272373 0.00898256 693 -0.303224  0.7618
## treatment3:pollen.ball.id -0.01265524 0.00896351 693 -1.411861  0.1584
## treatment4:pollen.ball.id -0.00974646 0.00925967 693 -1.052571  0.2929
##  Correlation: 
##                           (Intr) trtmn2 trtmn3 trtmn4 plln.. time2  tr2:..
## treatment2                -0.639                                          
## treatment3                -0.628  0.505                                   
## treatment4                -0.615  0.486  0.489                            
## pollen.ball.id            -0.352  0.106  0.090  0.099                     
## time2                      0.419  0.012  0.045  0.018 -0.499              
## treatment2:pollen.ball.id  0.113 -0.179 -0.092 -0.089 -0.610 -0.009       
## treatment3:pollen.ball.id  0.106 -0.092 -0.173 -0.089 -0.601 -0.028  0.503
## treatment4:pollen.ball.id  0.109 -0.089 -0.090 -0.180 -0.590 -0.011  0.487
##                           tr3:..
## treatment2                      
## treatment3                      
## treatment4                      
## pollen.ball.id                  
## time2                           
## treatment2:pollen.ball.id       
## treatment3:pollen.ball.id       
## treatment4:pollen.ball.id  0.488
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.54359522 -0.46754766 -0.08655237  0.34417022  5.21796278 
## 
## Number of Observations: 733
## Number of Groups: 35

anova(time.quad)

##                          numDF denDF   F-value p-value
## (Intercept)                  1   693 177.61545  <.0001
## treatment                    3    31   0.28572  0.8353
## pollen.ball.id               1   693  28.67272  <.0001
## time2                        1   693  69.02697  <.0001
## treatment:pollen.ball.id     3   693   0.85768  0.4628

fitted2 <- fitted(time.quad, level = 0)
a <- with(df,
          data.frame(pollen.ball.id, fitted2, treatment)[order(treatment, pollen.ball.id), ])

with(a, {
  plot(pollen.ball.id[treatment==1], fitted2[treatment==1],
       xlab = "time", ylab = "predicted", col = "green", type = "b")
  points(pollen.ball.id[treatment==2], fitted2[treatment==2],
         pch = 4, col = "red", type = "b")
  points(pollen.ball.id[treatment==3], fitted2[treatment==3],
         pch = 16, col = "blue", type = "b")
  points(pollen.ball.id[treatment==4], fitted2[treatment==4],
         pch = 16, col = "black", type = "b")
  title("Time Quadratic Effect")})

time.quad2 <- lme(whole_dif ~ treatment * pollen.ball.id + treatment * time2,
                  random = list(colony = pdDiag(~ pollen.ball.id)), data = df)
summary(time.quad2)

## Linear mixed-effects model fit by REML
##   Data: df 
##         AIC       BIC   logLik
##   -434.6329 -365.9233 232.3165
## 
## Random effects:
##  Formula: ~pollen.ball.id | colony
##  Structure: Diagonal
##         (Intercept) pollen.ball.id  Residual
## StdDev:  0.07853685     0.01811885 0.1466662
## 
## Fixed effects:  whole_dif ~ treatment * pollen.ball.id + treatment * time2 
##                                 Value  Std.Error  DF   t-value p-value
## (Intercept)                0.00217397 0.05009561 690  0.043396  0.9654
## treatment2                 0.08403008 0.07020347  31  1.196951  0.2404
## treatment3                 0.09236816 0.06937861  31  1.331364  0.1928
## treatment4                 0.09859801 0.07228950  31  1.363933  0.1824
## pollen.ball.id             0.07689942 0.01011987 690  7.598856  0.0000
## time2                     -0.00222936 0.00033505 690 -6.653781  0.0000
## treatment2:pollen.ball.id -0.03084097 0.01410296 690 -2.186844  0.0291
## treatment3:pollen.ball.id -0.04107491 0.01376972 690 -2.982988  0.0030
## treatment4:pollen.ball.id -0.03895283 0.01444011 690 -2.697545  0.0072
## treatment2:time2           0.00118722 0.00045838 690  2.590026  0.0098
## treatment3:time2           0.00118812 0.00043312 690  2.743128  0.0062
## treatment4:time2           0.00123021 0.00046533 690  2.643719  0.0084
##  Correlation: 
##                           (Intr) trtmn2 trtmn3 trtmn4 plln.. time2  tr2:..
## treatment2                -0.714                                          
## treatment3                -0.722  0.515                                   
## treatment4                -0.693  0.494  0.500                            
## pollen.ball.id            -0.628  0.448  0.453  0.435                     
## time2                      0.702 -0.501 -0.507 -0.486 -0.781              
## treatment2:pollen.ball.id  0.451 -0.620 -0.325 -0.312 -0.718  0.560       
## treatment3:pollen.ball.id  0.461 -0.329 -0.607 -0.320 -0.735  0.574  0.527
## treatment4:pollen.ball.id  0.440 -0.314 -0.318 -0.616 -0.701  0.547  0.503
## treatment2:time2          -0.513  0.697  0.371  0.356  0.571 -0.731 -0.774
## treatment3:time2          -0.543  0.388  0.691  0.376  0.604 -0.774 -0.433
## treatment4:time2          -0.505  0.361  0.365  0.694  0.562 -0.720 -0.403
##                           tr3:.. tr4:.. trt2:2 trt3:2
## treatment2                                           
## treatment3                                           
## treatment4                                           
## pollen.ball.id                                       
## time2                                                
## treatment2:pollen.ball.id                            
## treatment3:pollen.ball.id                            
## treatment4:pollen.ball.id  0.515                     
## treatment2:time2          -0.419 -0.400              
## treatment3:time2          -0.762 -0.423  0.565       
## treatment4:time2          -0.413 -0.771  0.526  0.557
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.25288400 -0.49575835 -0.07801207  0.31343459  5.11349891 
## 
## Number of Observations: 733
## Number of Groups: 35

anova(time.quad2)

##                          numDF denDF   F-value p-value
## (Intercept)                  1   690 173.34249  <.0001
## treatment                    3    31   0.27883  0.8402
## pollen.ball.id               1   690  29.01589  <.0001
## time2                        1   690  69.98457  <.0001
## treatment:pollen.ball.id     3   690   0.86897  0.4568
## treatment:time2              3   690   3.37306  0.0181

Anova(time.quad2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: whole_dif
##                            Chisq Df Pr(>Chisq)    
## treatment                 1.2331  3    0.74507    
## pollen.ball.id           97.1917  1    < 2e-16 ***
## time2                    69.2173  1    < 2e-16 ***
## treatment:pollen.ball.id 10.8282  3    0.01269 *  
## treatment:time2          10.1192  3    0.01758 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(time.quad2)

Run further quadratic analysis

#https://stackoverflow.com/questions/75580470/emtrends-in-presence-of-a-quadratic-term-difference-in-quadratic-trends
LM.fit <- lmer(whole_dif~crithidia*fungicide*poly(pollen.ball.id,2)+(1|colony), data = df)
LM.fit2 <- lmer(whole_dif~(crithidia+fungicide+poly(pollen.ball.id,2))^2+(1|colony), data = df)
summary(LM.fit2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: whole_dif ~ (crithidia + fungicide + poly(pollen.ball.id, 2))^2 +  
##     (1 | colony)
##    Data: df
## 
## REML criterion at convergence: -338.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6226 -0.5382 -0.1349  0.4351  4.7043 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  colony   (Intercept) 0.04837  0.2199  
##  Residual             0.03089  0.1758  
## Number of obs: 733, groups:  colony, 35
## 
## Fixed effects:
##                                        Estimate Std. Error t value
## (Intercept)                             0.50815    0.07447   6.824
## crithidiaTRUE                          -0.14416    0.10854  -1.328
## fungicideTRUE                          -0.06875    0.10529  -0.653
## poly(pollen.ball.id, 2)1                2.95216    0.31903   9.254
## poly(pollen.ball.id, 2)2               -2.80322    0.32784  -8.550
## crithidiaTRUE:fungicideTRUE             0.03852    0.15118   0.255
## crithidiaTRUE:poly(pollen.ball.id, 2)1 -1.31763    0.35872  -3.673
## crithidiaTRUE:poly(pollen.ball.id, 2)2  0.93547    0.36117   2.590
## fungicideTRUE:poly(pollen.ball.id, 2)1 -0.15999    0.35966  -0.445
## fungicideTRUE:poly(pollen.ball.id, 2)2  0.82807    0.36222   2.286
## 
## Correlation of Fixed Effects:
##               (Intr) crTRUE fnTRUE p(..,2)1 p(..,2)2 cTRUE: cTRUE:(..,2)1
## crithidTRUE   -0.686                                                     
## fungicdTRUE   -0.707  0.485                                              
## ply(p..,2)1    0.014 -0.009 -0.008                                       
## ply(p..,2)2    0.014 -0.008 -0.008  0.123                                
## crTRUE:TRUE    0.493 -0.718 -0.696  0.007    0.007                       
## cTRUE:(..,2)1 -0.008  0.006  0.003 -0.554   -0.055   -0.005              
## cTRUE:(..,2)2 -0.008  0.007  0.003 -0.057   -0.562   -0.006  0.043       
## fTRUE:(..,2)1 -0.008  0.005  0.007 -0.590   -0.053   -0.007 -0.042       
## fTRUE:(..,2)2 -0.008  0.003  0.008 -0.053   -0.585   -0.006 -0.021       
##               cTRUE:(..,2)2 fTRUE:(..,2)1
## crithidTRUE                              
## fungicdTRUE                              
## ply(p..,2)1                              
## ply(p..,2)2                              
## crTRUE:TRUE                              
## cTRUE:(..,2)1                            
## cTRUE:(..,2)2                            
## fTRUE:(..,2)1 -0.020                     
## fTRUE:(..,2)2 -0.059         0.048

anova(LM.fit2, LM.fit, test = "Chi")

## Data: df
## Models:
## LM.fit2: whole_dif ~ (crithidia + fungicide + poly(pollen.ball.id, 2))^2 + (1 | colony)
## LM.fit: whole_dif ~ crithidia * fungicide * poly(pollen.ball.id, 2) + (1 | colony)
##         npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## LM.fit2   12 -331.86 -276.69 177.93  -355.86                     
## LM.fit    14 -329.40 -265.04 178.70  -357.40 1.5418  2     0.4626

#LM.fit2 fits better


emtrends(LM.fit2,pairwise~crithidia, var = "pollen.ball.id", max.degree = 2)

## $emtrends
## degree = linear:
##  crithidia pollen.ball.id.trend       SE  df lower.CL  upper.CL
##  FALSE                  0.01823 0.001517 693  0.01525  0.021208
##   TRUE                  0.00995 0.001496 693  0.00701  0.012885
## 
## degree = quadratic:
##  crithidia pollen.ball.id.trend       SE  df lower.CL  upper.CL
##  FALSE                 -0.00234 0.000261 693 -0.00285 -0.001830
##   TRUE                 -0.00142 0.000241 693 -0.00190 -0.000951
## 
## Results are averaged over the levels of: fungicide 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## degree = linear:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE  0.008282 0.002126 693   3.895  0.0001
## 
## degree = quadratic:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE -0.000917 0.000354 693  -2.590  0.0098
## 
## Results are averaged over the levels of: fungicide 
## Degrees-of-freedom method: kenward-roger

emtrends(LM.fit2,pairwise~fungicide, var = "pollen.ball.id", max.degree = 2)

## $emtrends
## degree = linear:
##  fungicide pollen.ball.id.trend       SE  df lower.CL upper.CL
##  FALSE                  0.01477 0.001566 693  0.01170  0.01785
##   TRUE                  0.01341 0.001448 693  0.01056  0.01625
## 
## degree = quadratic:
##  fungicide pollen.ball.id.trend       SE  df lower.CL upper.CL
##  FALSE                 -0.00229 0.000266 693 -0.00281 -0.00177
##   TRUE                 -0.00148 0.000236 693 -0.00194 -0.00101
## 
## Results are averaged over the levels of: crithidia 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## degree = linear:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE  0.001367 0.002131 693   0.641  0.5215
## 
## degree = quadratic:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE -0.000812 0.000355 693  -2.286  0.0226
## 
## Results are averaged over the levels of: crithidia 
## Degrees-of-freedom method: kenward-roger

##NEW TIME Variable 

df$pollen.time <- as.numeric(df$pollen.time)
df$block <- as.factor(df$block)


LM.fit <- lmer(whole_dif~crithidia*fungicide*poly(pollen.time,2)+(1|colony) + block, data = df)
LM.fit2 <- lmer(whole_dif~(crithidia+fungicide+poly(pollen.time,2))^2+(1|colony) + block, data = df)
LM.fit3 <- lmer(whole_dif~(crithidia+fungicide+poly(pollen.time,2))^2+(1|colony), data = df)
summary(LM.fit2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: whole_dif ~ (crithidia + fungicide + poly(pollen.time, 2))^2 +  
##     (1 | colony) + block
##    Data: df
## 
## REML criterion at convergence: -344.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6514 -0.5244 -0.1313  0.4560  4.6582 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  colony   (Intercept) 0.02156  0.1468  
##  Residual             0.03095  0.1759  
## Number of obs: 733, groups:  colony, 35
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                          0.42954    0.08771   4.897
## crithidiaTRUE                       -0.12779    0.07454  -1.714
## fungicideTRUE                       -0.06971    0.07161  -0.974
## poly(pollen.time, 2)1                2.94132    0.31792   9.252
## poly(pollen.time, 2)2               -2.75860    0.32161  -8.577
## block4                               0.45781    0.10748   4.259
## block6                              -0.15034    0.10713  -1.403
## block7                               0.02902    0.10743   0.270
## block8                               0.13906    0.10737   1.295
## block9                               0.04396    0.10738   0.409
## block10                              0.21771    0.11704   1.860
## block11                             -0.09109    0.10723  -0.849
## block12                              0.07204    0.10725   0.672
## crithidiaTRUE:fungicideTRUE          0.02148    0.10331   0.208
## crithidiaTRUE:poly(pollen.time, 2)1 -1.35371    0.35888  -3.772
## crithidiaTRUE:poly(pollen.time, 2)2  0.81627    0.36106   2.261
## fungicideTRUE:poly(pollen.time, 2)1 -0.10258    0.35990  -0.285
## fungicideTRUE:poly(pollen.time, 2)2  0.87969    0.36215   2.429

anova(LM.fit2, LM.fit, test = "Chi")

## Data: df
## Models:
## LM.fit2: whole_dif ~ (crithidia + fungicide + poly(pollen.time, 2))^2 + (1 | colony) + block
## LM.fit: whole_dif ~ crithidia * fungicide * poly(pollen.time, 2) + (1 | colony) + block
##         npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## LM.fit2   20 -352.13 -260.19 196.06  -392.13                     
## LM.fit    22 -349.80 -248.67 196.90  -393.80 1.6763  2     0.4325

anova(LM.fit2, LM.fit3, test = "Chi")

## Data: df
## Models:
## LM.fit3: whole_dif ~ (crithidia + fungicide + poly(pollen.time, 2))^2 + (1 | colony)
## LM.fit2: whole_dif ~ (crithidia + fungicide + poly(pollen.time, 2))^2 + (1 | colony) + block
##         npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## LM.fit3   12 -330.81 -275.64 177.41  -354.81                         
## LM.fit2   20 -352.13 -260.19 196.06  -392.13 37.318  8  1.006e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(LM.fit2, LM.fit3)

##         df       AIC
## LM.fit2 20 -304.8886
## LM.fit3 12 -313.1797

#LM.fit3 fits better

Anova(LM.fit3)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: whole_dif
##                                   Chisq Df Pr(>Chisq)    
## crithidia                        2.7397  1    0.09788 .  
## fungicide                        0.4755  1    0.49048    
## poly(pollen.time, 2)           279.3323  2  < 2.2e-16 ***
## crithidia:fungicide              0.0649  1    0.79894    
## crithidia:poly(pollen.time, 2)  20.0616  2  4.402e-05 ***
## fungicide:poly(pollen.time, 2)   5.7454  2    0.05654 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

qqnorm(resid(LM.fit3));qqline(resid(LM.fit3))

emtrends(LM.fit3,pairwise~crithidia, var = "pollen.time", max.degree = 2)

## $emtrends
## degree = linear:
##  crithidia pollen.time.trend       SE  df  lower.CL  upper.CL
##  FALSE              0.009076 7.55e-04 693  0.007594  0.010558
##   TRUE              0.004890 7.46e-04 693  0.003425  0.006355
## 
## degree = quadratic:
##  crithidia pollen.time.trend       SE  df  lower.CL  upper.CL
##  FALSE             -0.000565 6.40e-05 694 -0.000691 -0.000440
##   TRUE             -0.000366 6.06e-05 694 -0.000485 -0.000247
## 
## Results are averaged over the levels of: fungicide 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## degree = linear:
##  contrast     estimate       SE  df t.ratio p.value
##  FALSE - TRUE  0.00419 1.06e-03 693   3.951  0.0001
## 
## degree = quadratic:
##  contrast     estimate       SE  df t.ratio p.value
##  FALSE - TRUE -0.00020 8.79e-05 694  -2.269  0.0236
## 
## Results are averaged over the levels of: fungicide 
## Degrees-of-freedom method: kenward-roger

emtrends(LM.fit3,pairwise~fungicide, var = "pollen.time", max.degree = 2)

## $emtrends
## degree = linear:
##  fungicide pollen.time.trend       SE  df  lower.CL  upper.CL
##  FALSE              0.007269 7.78e-04 693  0.005742  0.008797
##   TRUE              0.006697 7.24e-04 693  0.005276  0.008117
## 
## degree = quadratic:
##  fungicide pollen.time.trend       SE  df  lower.CL  upper.CL
##  FALSE             -0.000569 6.53e-05 694 -0.000698 -0.000441
##   TRUE             -0.000362 5.93e-05 693 -0.000478 -0.000245
## 
## Results are averaged over the levels of: crithidia 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## degree = linear:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE  0.000573 1.06e-03 693   0.539  0.5899
## 
## degree = quadratic:
##  contrast      estimate       SE  df t.ratio p.value
##  FALSE - TRUE -0.000208 8.82e-05 694  -2.353  0.0189
## 
## Results are averaged over the levels of: crithidia 
## Degrees-of-freedom method: kenward-roger

Plot

par(mar = c(5, 5, 5, 5))


df$time2 <- df$pollen.time^2


# Calculate fitted values
time.quad2 <- lme(whole_dif ~ treatment * pollen.time + treatment * time2,
                  random = list(colony = pdDiag(~ pollen.time)), data = df)


fitted3 <- fitted(time.quad2, level = 0)

# Create a data frame with ordered values
a <- with(df,
          data.frame(pollen.time, fitted3, treatment)[order(treatment, pollen.time), ])

# Set the font size multiplier and symbol size multiplier
font_multiplier <- 2
symbol_multiplier <- 2  # Adjust this for larger symbols
legend_cex <- 1.5

# Plot
with(a, {
  # Plot data for treatment 1
  plot(pollen.time[treatment==1], fitted3[treatment==1],
       ylab = "Average Consumption (g)", xlab = "Time (days)",
       col = "darkorchid4", type = "b", pch = 1,
       cex.lab = font_multiplier, # Size of axis labels
       cex.axis = font_multiplier, # Size of axis numbers
       cex = symbol_multiplier)    # Size of symbols
  
  # Add points for other treatments with increased symbol size
  points(pollen.time[treatment==2], fitted3[treatment==2],
         pch = 4, col = "dodgerblue", type = "b",
         cex = symbol_multiplier)    # Size of symbols
  points(pollen.time[treatment==3], fitted3[treatment==3],
         pch = 16, col = "orange4", type = "b",
         cex = symbol_multiplier)    # Size of symbols
  points(pollen.time[treatment==4], fitted3[treatment==4],
         pch = 17, col = "orange", type = "b",
         cex = symbol_multiplier)    # Size of symbols
  
  # Add a legend with larger font
  legend("topleft", legend = c("Control", "Fungicide", "Fungicide & Parasite", "Parasite"),
         col = c("darkorchid4", "dodgerblue", "orange4", "orange"), pch = c(1, 4, 16, 17), 
         lty = 1, cex = legend_cex, bty= "n") # Size of legend text
})