Project Analysis

1.1 Factors: Canopy cover percent (continuous) ~ Ash dieback levels (categories)

#check normality of data 
library(dplyr)
library(ggpubr)
library(car)
library(dunn.test)
library(FSA)
#change the colname of Canopy_cover_percent in Ash data
colnames(Ash)[colnames(Ash)=="Canopy_cover_percent"] <- "Canopy_cover"
colnames(Ash)[colnames(Ash)=="Ash_dieback"] <- "AD"
#the qq plot shows characteristic of binomial data
p <- ggqqplot(Ash$Canopy_cover)
p
#change the factor levels of AD
Ash$AD <- factor(Ash$AD, levels = c("low", "medium", "high"))
#boxplot
p <- ggboxplot(Ash, x="AD", y= "Canopy_cover",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Canopy Cover (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

#GLM of all three factors
library(MASS)
Canopy_AD <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD, data=Ash, family=binomial)

non-integer counts in a binomial glm!

summary(Canopy_AD)

Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD, family = binomial, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8170  -1.0968  -0.1488   1.4751   3.3650  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.41139    0.06808  -6.043 1.52e-09 ***
ADmedium    -0.10631    0.09439  -1.126     0.26    
ADhigh      -1.04340    0.09639 -10.825  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 271.00  on 32  degrees of freedom
Residual deviance: 117.65  on 30  degrees of freedom
AIC: 280.69

Number of Fisher Scoring iterations: 4

step.model1 <- stepAIC(Canopy_AD, direction = "backward", 
                       trace = FALSE)

non-integer #successes in a binomial glm!

summary(step.model1)

Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD, family = binomial, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8170  -1.0968  -0.1488   1.4751   3.3650  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.41139    0.06808  -6.043 1.52e-09 ***
ADmedium    -0.10631    0.09439  -1.126     0.26    
ADhigh      -1.04340    0.09639 -10.825  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 271.00  on 32  degrees of freedom
Residual deviance: 117.65  on 30  degrees of freedom
AIC: 280.69

Number of Fisher Scoring iterations: 4

Anova(step.model1)

non-integer #successes in a binomial glm!

Analysis of Deviance Table (Type II tests)

Response: cbind(Canopy_cover, Canopy_not_cover)
   LR Chisq Df Pr(>Chisq)    
AD   153.35  2  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Canopy_AD_BA <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD * Understorey, data=Ash, family=binomial)

non-integer counts in a binomial glm!

summary(Canopy_AD_BA)

Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * Understorey, 
    family = binomial, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8910  -1.1736  -0.1596   1.4570   3.3494  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)           -0.49178    0.10304  -4.773 1.82e-06 ***
ADmedium              -0.09076    0.17989  -0.504    0.614    
ADhigh                -0.96580    0.13703  -7.048 1.81e-12 ***
Understorey+           0.14375    0.13734   1.047    0.295    
ADmedium:Understorey+ -0.06288    0.21431  -0.293    0.769    
ADhigh:Understorey+   -0.13726    0.19457  -0.705    0.481    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 271.00  on 32  degrees of freedom
Residual deviance: 116.31  on 27  degrees of freedom
AIC: 285.13

Number of Fisher Scoring iterations: 4

Anova(Canopy_AD_BA)

non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!

Analysis of Deviance Table (Type II tests)

Response: cbind(Canopy_cover, Canopy_not_cover)
               LR Chisq Df Pr(>Chisq)    
AD              140.272  2     <2e-16 ***
Understorey       0.843  1     0.3585    
AD:Understorey    0.499  2     0.7791    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

step.model2 <- stepAIC(Canopy_AD_BA, direction = "backward", 
                       trace = FALSE)

non-integer #successes in a binomial glm!non-integer counts in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer counts in a binomial glm!non-integer #successes in a binomial glm!

summary(step.model2)

Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD, family = binomial, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8170  -1.0968  -0.1488   1.4751   3.3650  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.41139    0.06808  -6.043 1.52e-09 ***
ADmedium    -0.10631    0.09439  -1.126     0.26    
ADhigh      -1.04340    0.09639 -10.825  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 271.00  on 32  degrees of freedom
Residual deviance: 117.65  on 30  degrees of freedom
AIC: 280.69

Number of Fisher Scoring iterations: 4

Anova(step.model2)

non-integer #successes in a binomial glm!

Analysis of Deviance Table (Type II tests)

Response: cbind(Canopy_cover, Canopy_not_cover)
   LR Chisq Df Pr(>Chisq)    
AD   153.35  2  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dropterm(step.model2,test = "Chisq")

non-integer #successes in a binomial glm!

Single term deletions

Model:
cbind(Canopy_cover, Canopy_not_cover) ~ AD
       Df Deviance    AIC    LRT   Pr(Chi)    
<none>      117.65 280.69                     
AD      2   271.00 430.05 153.35 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#if dropped AD:BA, AIC increases, so don't drop
Canopy_AD_BA_2 <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD + Understorey, data=Ash, family=binomial)

non-integer counts in a binomial glm!

anova(Canopy_AD_BA, Canopy_AD_BA_2, test="Chisq")

Analysis of Deviance Table

Model 1: cbind(Canopy_cover, Canopy_not_cover) ~ AD * Understorey
Model 2: cbind(Canopy_cover, Canopy_not_cover) ~ AD + Understorey
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        27     116.31                     
2        29     116.81 -2 -0.49918   0.7791

summary(Canopy_AD_BA_2)

Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD + Understorey, 
    family = binomial, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8873  -1.1698  -0.1798   1.6026   3.1812  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.45418    0.08261  -5.498 3.84e-08 ***
ADmedium     -0.12504    0.09658  -1.295    0.195    
ADhigh       -1.03396    0.09691 -10.669  < 2e-16 ***
Understorey+  0.07675    0.08361   0.918    0.359    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 271.00  on 32  degrees of freedom
Residual deviance: 116.81  on 29  degrees of freedom
AIC: 281.82

Number of Fisher Scoring iterations: 4

anova(Canopy_AD_BA,Canopy_AD_BA_2)

Analysis of Deviance Table

Model 1: cbind(Canopy_cover, Canopy_not_cover) ~ AD * Understorey
Model 2: cbind(Canopy_cover, Canopy_not_cover) ~ AD + Understorey
  Resid. Df Resid. Dev Df Deviance
1        27     116.31            
2        29     116.81 -2 -0.49918

#therefore understorey and the interaction between AD and understorey are not important

Design a linear regression plot function

ggplotRegression <- function (fit) {
    
    require(ggplot2)
    
    ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
        geom_point() +
        stat_smooth(method = "lm", col = "red") +
        labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
                           "Intercept =",signif(fit$coef[[1]],5 ),
                           " Slope =",signif(fit$coef[[2]], 5),
                           " P =",signif(summary(fit)$coef[2,4], 5)))
}

Factor 2: temperature fluctuation vs AD & LAI

## Add temperature summaries to Ash data
Ash$Max_temp <- Temp_data$Max_temp[match(Ash$Site, Temp_data$Site)]
Ash$Min_temp <- Temp_data$Min_temp[match(Ash$Site, Temp_data$Site)]
Ash$Temp_range <- Temp_data$Temp_range[match(Ash$Site, Temp_data$Site)] 
#Max temp vs AD
#test normality (not normal, SW test: p = 0.036)
p <- ggqqplot(Ash$Max_temp)
p

shapiro.test(Ash$Max_temp)

    Shapiro-Wilk normality test

data:  Ash$Max_temp
W = 0.93058, p-value = 0.03634

#boxplot of AD ~ Max temp
p <- ggboxplot(Ash, x="AD", y= "Max_temp",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Maximum temperature (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

#Dunn test: 
kruskal.test(Max_temp ~ AD, data = Ash)

    Kruskal-Wallis rank sum test

data:  Max_temp by AD
Kruskal-Wallis chi-squared = 13.316, df = 2, p-value = 0.001283

dunnTest(Ash$Max_temp, Ash$AD, method = "bonferroni")

Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Bonferroni method.

#significance is between medium and high levels of AD (p=0.0049)
#Minimum Temp vs AD
p <- ggqqplot(Ash$Min_temp)
p

shapiro.test(Ash$Min_temp) #not normal again

    Shapiro-Wilk normality test

data:  Ash$Min_temp
W = 0.93681, p-value = 0.05489

p <- ggboxplot(Ash, x="AD", y= "Min_temp",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Minimum temperature (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

#Dunn test: 
kruskal.test(Min_temp ~ AD, data = Ash)

    Kruskal-Wallis rank sum test

data:  Min_temp by AD
Kruskal-Wallis chi-squared = 4.9176, df = 2, p-value = 0.08554

dunnTest(Ash$Min_temp, Ash$AD, method = "bonferroni")

Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Bonferroni method.

#no significance 
#Temp range vs AD
p <- ggqqplot(Ash$Temp_range)
p

shapiro.test(Ash$Temp_range) #not very normal again

    Shapiro-Wilk normality test

data:  Ash$Temp_range
W = 0.93922, p-value = 0.0645

p <- ggboxplot(Ash, x="AD", y= "Temp_range",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Temperature range (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

kruskal.test(Temp_range ~ AD, data = Ash)

    Kruskal-Wallis rank sum test

data:  Temp_range by AD
Kruskal-Wallis chi-squared = 15.28, df = 2, p-value = 0.0004809

dunn.test(Ash$Temp_range, Ash$AD, method = "bonferroni")

  Kruskal-Wallis rank sum test

data: x and group
Kruskal-Wallis chi-squared = 15.2797, df = 2, p-value = 0

                           Comparison of x by group                            
                                 (Bonferroni)                                  
Col Mean-|
Row Mean |       high        low
---------+----------------------
     low |   3.309253
         |    0.0014*
         |
  medium |   3.239603  -0.157886
         |    0.0018*     1.0000

alpha = 0.05
Reject Ho if p <= alpha/2

#significance between medium and high and between low and high (0.001 < p < 0.005)
#temp ~ LAI
m1 <- lm(Temp_range ~ Canopy_cover, data = Ash)
m2 <- lm(Min_temp ~ Canopy_cover  , data = Ash)
m3 <- lm(Max_temp ~ Canopy_cover, data = Ash)
ggplotRegression(lm(Temp_range ~ Canopy_cover, data = Ash)) + xlab("Canopy cover (%)") + ylab("Temperature range (°C)")

tab_model(
    m1,m2,m3,  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
    pred.labels = c("Intercept", "Canopy cover (%)"),
    dv.labels = c("Temperature range (C)", "Minimum temperature (C)", "Maximum temperature (C)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
ggplotRegression(lm(Min_temp ~ Canopy_cover  , data = Ash)) + xlab("Canopy cover (%)") + ylab("Minimum temperature (°C)")

ggplotRegression(lm(Max_temp ~ Canopy_cover, data = Ash)) + xlab("Canopy cover (%)") + ylab("Maximum temperature (°C)")

#GLM for temp range ~ AD + LAI
GLM_range_AD_LAI <- glm(Temp_range ~ AD + Canopy_cover, data = Ash)
summary(GLM_range_AD_LAI)

Call:
glm(formula = Temp_range ~ AD + Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9472  -0.8772  -0.5632   1.0921   3.9905  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.10955    1.44423   8.385 3.06e-09 ***
ADmedium     -0.21535    0.74889  -0.288   0.7757    
ADhigh        2.02649    0.98671   2.054   0.0491 *  
Canopy_cover -0.03109    0.03361  -0.925   0.3626    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.622622)

    Null deviance: 139.320  on 32  degrees of freedom
Residual deviance:  76.056  on 29  degrees of freedom
AIC: 131.2

Number of Fisher Scoring iterations: 2

GLM_range_AD_LAI_model <- stepAIC(GLM_range_AD_LAI, direction = "backward", Trace = FALSE)

Start:  AIC=131.2
Temp_range ~ AD + Canopy_cover

               Df Deviance    AIC
- Canopy_cover  1   78.300 130.16
<none>              76.056 131.20
- AD            2   92.520 133.67

Step:  AIC=130.16
Temp_range ~ AD

       Df Deviance    AIC
<none>       78.30 130.16
- AD    2   139.32 145.18

summary(GLM_range_AD_LAI_model)

Call:
glm(formula = Temp_range ~ AD, data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.881  -1.037  -0.400   1.286   3.600  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.8704     0.5385  20.186  < 2e-16 ***
ADmedium     -0.1370     0.7423  -0.185 0.854775    
ADhigh        2.6772     0.6902   3.879 0.000532 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.610012)

    Null deviance: 139.32  on 32  degrees of freedom
Residual deviance:  78.30  on 30  degrees of freedom
AIC: 130.16

Number of Fisher Scoring iterations: 2

Anova(GLM_range_AD_LAI_model)

Analysis of Deviance Table (Type II tests)

Response: Temp_range
   LR Chisq Df Pr(>Chisq)    
AD   23.379  2  8.381e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#GLM for temp range ~ AD * LAI
GLM_range_AD_LAI <- glm(Temp_range ~ AD + Canopy_cover, data = Ash)
summary(GLM_range_AD_LAI)

Call:
glm(formula = Temp_range ~ AD + Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9472  -0.8772  -0.5632   1.0921   3.9905  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.10955    1.44423   8.385 3.06e-09 ***
ADmedium     -0.21535    0.74889  -0.288   0.7757    
ADhigh        2.02649    0.98671   2.054   0.0491 *  
Canopy_cover -0.03109    0.03361  -0.925   0.3626    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.622622)

    Null deviance: 139.320  on 32  degrees of freedom
Residual deviance:  76.056  on 29  degrees of freedom
AIC: 131.2

Number of Fisher Scoring iterations: 2

GLM_range_AD_LAI_model <- stepAIC(GLM_range_AD_LAI, direction = "backward", Trace = FALSE)

Start:  AIC=131.2
Temp_range ~ AD + Canopy_cover

               Df Deviance    AIC
- Canopy_cover  1   78.300 130.16
<none>              76.056 131.20
- AD            2   92.520 133.67

Step:  AIC=130.16
Temp_range ~ AD

       Df Deviance    AIC
<none>       78.30 130.16
- AD    2   139.32 145.18

summary(GLM_range_AD_LAI_model)

Call:
glm(formula = Temp_range ~ AD, data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.881  -1.037  -0.400   1.286   3.600  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.8704     0.5385  20.186  < 2e-16 ***
ADmedium     -0.1370     0.7423  -0.185 0.854775    
ADhigh        2.6772     0.6902   3.879 0.000532 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.610012)

    Null deviance: 139.32  on 32  degrees of freedom
Residual deviance:  78.30  on 30  degrees of freedom
AIC: 130.16

Number of Fisher Scoring iterations: 2

Anova(GLM_range_AD_LAI_model)

Analysis of Deviance Table (Type II tests)

Response: Temp_range
   LR Chisq Df Pr(>Chisq)    
AD   23.379  2  8.381e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#GLM for temp Max ~ AD + LAI
GLM_Max_AD_LAI <- glm(Max_temp ~ AD + Canopy_cover, data = Ash)
summary(GLM_Max_AD_LAI)

Call:
glm(formula = Max_temp ~ AD + Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6237  -0.9007  -0.4507   1.1257   4.2976  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  23.75651    1.42517  16.669   <2e-16 ***
ADmedium     -0.13895    0.73900  -0.188    0.852    
ADhigh        1.95513    0.97369   2.008    0.054 .  
Canopy_cover -0.02502    0.03316  -0.754    0.457    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.553854)

    Null deviance: 126.682  on 32  degrees of freedom
Residual deviance:  74.062  on 29  degrees of freedom
AIC: 130.33

Number of Fisher Scoring iterations: 2

GLM_Max_AD_LAI_model <- stepAIC(GLM_Max_AD_LAI, direction = "backward", Trace = FALSE)

Start:  AIC=130.33
Max_temp ~ AD + Canopy_cover

               Df Deviance    AIC
- Canopy_cover  1   75.515 128.97
<none>              74.062 130.33
- AD            2   88.657 132.26

Step:  AIC=128.97
Max_temp ~ AD

       Df Deviance    AIC
<none>      75.515 128.97
- AD    2  126.682 142.04

summary(GLM_Max_AD_LAI_model)

Call:
glm(formula = Max_temp ~ AD, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7381  -1.0167  -0.4048   1.0952   3.9833  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22.75926    0.52885  43.035  < 2e-16 ***
ADmedium    -0.07593    0.72897  -0.104 0.917740    
ADhigh       2.47884    0.67785   3.657 0.000971 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.517177)

    Null deviance: 126.682  on 32  degrees of freedom
Residual deviance:  75.515  on 30  degrees of freedom
AIC: 128.97

Number of Fisher Scoring iterations: 2

Anova(GLM_Max_AD_LAI_model)

Analysis of Deviance Table (Type II tests)

Response: Max_temp
   LR Chisq Df Pr(>Chisq)    
AD   20.327  2  3.855e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#GLM for temp Max ~ AD * LAI
GLM_Max_AD_LAI <- glm(Max_temp ~ AD * Canopy_cover, data = Ash)
summary(GLM_Max_AD_LAI)

Call:
glm(formula = Max_temp ~ AD * Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8176  -0.8271  -0.3953   0.7564   2.9415  

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           25.68355    2.33850  10.983 1.83e-11 ***
ADmedium              -6.09723    2.86845  -2.126   0.0428 *  
ADhigh                 2.21083    2.55805   0.864   0.3951    
Canopy_cover          -0.07337    0.05753  -1.275   0.2131    
ADmedium:Canopy_cover  0.15631    0.07179   2.177   0.0384 *  
ADhigh:Canopy_cover   -0.06698    0.07703  -0.870   0.3922    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.894189)

    Null deviance: 126.682  on 32  degrees of freedom
Residual deviance:  51.143  on 27  degrees of freedom
AIC: 122.11

Number of Fisher Scoring iterations: 2

GLM_Max_AD_LAI_model <- stepAIC(GLM_Max_AD_LAI, direction = "backward", Trace = FALSE)

Start:  AIC=122.11
Max_temp ~ AD * Canopy_cover

                  Df Deviance    AIC
<none>                 51.143 122.11
- AD:Canopy_cover  2   74.062 130.33

summary(GLM_Max_AD_LAI_model)

Call:
glm(formula = Max_temp ~ AD * Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8176  -0.8271  -0.3953   0.7564   2.9415  

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           25.68355    2.33850  10.983 1.83e-11 ***
ADmedium              -6.09723    2.86845  -2.126   0.0428 *  
ADhigh                 2.21083    2.55805   0.864   0.3951    
Canopy_cover          -0.07337    0.05753  -1.275   0.2131    
ADmedium:Canopy_cover  0.15631    0.07179   2.177   0.0384 *  
ADhigh:Canopy_cover   -0.06698    0.07703  -0.870   0.3922    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.894189)

    Null deviance: 126.682  on 32  degrees of freedom
Residual deviance:  51.143  on 27  degrees of freedom
AIC: 122.11

Number of Fisher Scoring iterations: 2

Anova(GLM_Max_AD_LAI_model)

Analysis of Deviance Table (Type II tests)

Response: Max_temp
                LR Chisq Df Pr(>Chisq)   
AD                7.7053  2   0.021223 * 
Canopy_cover      0.7674  1   0.381034   
AD:Canopy_cover  12.0995  2   0.002358 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#GLM for temp min ~ AD + LAI
GLM_Min_AD_LAI <- glm(Min_temp ~ AD + Canopy_cover, data = Ash)
summary(GLM_Min_AD_LAI)

Call:
glm(formula = Min_temp ~ AD + Canopy_cover, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5848  -0.1785   0.0308   0.1574   0.5249  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.646961   0.238605  48.813   <2e-16 ***
ADmedium      0.076400   0.123725   0.617    0.542    
ADhigh       -0.071364   0.163017  -0.438    0.665    
Canopy_cover  0.006070   0.005553   1.093    0.283    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.07158495)

    Null deviance: 2.6077  on 32  degrees of freedom
Residual deviance: 2.0760  on 29  degrees of freedom
AIC: 12.369

Number of Fisher Scoring iterations: 2

GLM_Min_AD_LAI_model <- stepAIC(GLM_Min_AD_LAI, direction = "backward", Trace = FALSE)

Start:  AIC=12.37
Min_temp ~ AD + Canopy_cover

               Df Deviance     AIC
- AD            2   2.1527  9.5664
- Canopy_cover  1   2.1615 11.7018
<none>              2.0760 12.3692

Step:  AIC=9.57
Min_temp ~ Canopy_cover

               Df Deviance     AIC
<none>              2.1527  9.5664
- Canopy_cover  1   2.6077 13.8952

summary(GLM_Min_AD_LAI_model)

Call:
glm(formula = Min_temp ~ Canopy_cover, data = Ash)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.50704  -0.15187   0.02257   0.15022   0.54313  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.547142   0.117198   98.53   <2e-16 ***
Canopy_cover  0.009138   0.003569    2.56   0.0156 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.06944061)

    Null deviance: 2.6077  on 32  degrees of freedom
Residual deviance: 2.1527  on 31  degrees of freedom
AIC: 9.5664

Number of Fisher Scoring iterations: 2

Anova(GLM_Min_AD_LAI_model)

Analysis of Deviance Table (Type II tests)

Response: Min_temp
             LR Chisq Df Pr(>Chisq)  
Canopy_cover   6.5536  1    0.01047 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Factor 3: Soil moisture ~ AD

p <- ggqqplot(Ash$Soil_moisture_ave)
p

shapiro.test(Ash$Soil_moisture_ave) #the data is normally distributed

    Shapiro-Wilk normality test

data:  Ash$Soil_moisture_ave
W = 0.94762, p-value = 0.1133

p <- ggboxplot(Ash, x="AD", y= "Soil_moisture_ave",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Soil moisture(m^3/m^3)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

SM_AD <- aov(Soil_moisture_ave ~ AD,data = Ash)
summary(SM_AD)

            Df  Sum Sq  Mean Sq F value Pr(>F)  
AD           2 0.01015 0.005074   3.288 0.0511 .
Residuals   30 0.04629 0.001543                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(SM_AD)#significance between low and high(p=0.042)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Soil_moisture_ave ~ AD, data = Ash)

$AD
                  diff          lwr        upr     p adj
medium-low  0.03041111 -0.014082207 0.07490443 0.2274018
high-low    0.04275397  0.001380884 0.08412705 0.0416947
high-medium 0.01234286 -0.027751278 0.05243699 0.7306348

#soil moisture vs LAI
m4 <- lm(Canopy_cover ~ Soil_moisture_ave, data = Ash)
m5 <- lm(Temp_range ~ Soil_moisture_ave, data = Ash)
tab_model(
  m4,
  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
  pred.labels = c("Intercept", "Canopy cover (%)"),
  dv.labels = c("Soil moisture (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
tab_model(
  m5,
  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
  pred.labels = c("Intercept", "Temperature range (C)"),
  dv.labels = c("Soil moisture (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
ggplotRegression(lm(Canopy_cover ~ Soil_moisture_ave, data = Ash)) + xlab("Canopy cover (%)") + ylab("Soil moisture (m^3/m^3) ")

ggplotRegression(lm(Temp_range ~ Soil_moisture_ave, data = Ash)) + xlab("Temperature range (°C)") + ylab("Soil moisture (m^3/m^3)")

#GLM for soil ~ AD + LAI + temp (with AD: p=0.03)
GLM_soil_AD_LAI_T <- glm(Soil_moisture_ave ~ AD + Canopy_cover + Temp_range, data = Ash)
summary(GLM_soil_AD_LAI_T)

Call:
glm(formula = Soil_moisture_ave ~ AD + Canopy_cover + Temp_range, 
    data = Ash)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.097266  -0.022613  -0.005102   0.023471   0.095765  

Coefficients:
               Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.1246452  0.0669488   1.862   0.0732 .
ADmedium      0.0303401  0.0187869   1.615   0.1175  
ADhigh        0.0367396  0.0264544   1.389   0.1758  
Canopy_cover -0.0001055  0.0008542  -0.124   0.9026  
Temp_range    0.0014213  0.0046518   0.306   0.7622  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.001645804)

    Null deviance: 0.056436  on 32  degrees of freedom
Residual deviance: 0.046083  on 28  degrees of freedom
AIC: -111.29

Number of Fisher Scoring iterations: 2

Anova(GLM_soil_AD_LAI_T)

Analysis of Deviance Table (Type II tests)

Response: Soil_moisture_ave
             LR Chisq Df Pr(>Chisq)
AD             3.2121  2     0.2007
Canopy_cover   0.0153  1     0.9017
Temp_range     0.0934  1     0.7600

GLM_soil_AD_LAI_T_model <- stepAIC(GLM_soil_AD_LAI_T, direction = "backward", Trace = FALSE)

Start:  AIC=-111.29
Soil_moisture_ave ~ AD + Canopy_cover + Temp_range

               Df Deviance     AIC
- Canopy_cover  1 0.046108 -113.27
- Temp_range    1 0.046236 -113.18
- AD            2 0.051369 -111.70
<none>            0.046083 -111.29

Step:  AIC=-113.27
Soil_moisture_ave ~ AD + Temp_range

             Df Deviance     AIC
- Temp_range  1 0.046288 -115.14
<none>          0.046108 -113.27
- AD          2 0.052999 -112.67

Step:  AIC=-115.14
Soil_moisture_ave ~ AD

       Df Deviance     AIC
<none>    0.046288 -115.14
- AD    2 0.056436 -112.60

summary(GLM_soil_AD_LAI_T_model)

Call:
glm(formula = Soil_moisture_ave ~ AD, data = Ash)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.099643  -0.024643  -0.000889   0.020111   0.097357  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.13589    0.01309  10.378 1.91e-11 ***
ADmedium     0.03041    0.01805   1.685   0.1024    
ADhigh       0.04275    0.01678   2.548   0.0162 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.00154294)

    Null deviance: 0.056436  on 32  degrees of freedom
Residual deviance: 0.046288  on 30  degrees of freedom
AIC: -115.14

Number of Fisher Scoring iterations: 2

Anova(GLM_soil_AD_LAI_T_model)

Analysis of Deviance Table (Type II tests)

Response: Soil_moisture_ave
   LR Chisq Df Pr(>Chisq)  
AD    6.577  2    0.03731 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#GLM for soil ~ AD * LAI * temp (p(AD) = 0.037 only appear after elimination of all other terms)
GLM_soil_AD_LAI_T2 <- glm(Soil_moisture_ave ~ AD * Canopy_cover * Temp_range, data = Ash)
summary(GLM_soil_AD_LAI_T2)

Call:
glm(formula = Soil_moisture_ave ~ AD * Canopy_cover * Temp_range, 
    data = Ash)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.081414  -0.015334  -0.000157   0.009398   0.079531  

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)                      -0.921634   0.800870  -1.151   0.2628  
ADmedium                          1.597247   0.866373   1.844   0.0794 .
ADhigh                            1.257761   0.855518   1.470   0.1563  
Canopy_cover                      0.029502   0.021059   1.401   0.1758  
Temp_range                        0.099348   0.074217   1.339   0.1950  
ADmedium:Canopy_cover            -0.042537   0.022391  -1.900   0.0713 .
ADhigh:Canopy_cover              -0.031984   0.024499  -1.306   0.2058  
ADmedium:Temp_range              -0.153826   0.081514  -1.887   0.0730 .
ADhigh:Temp_range                -0.106541   0.077253  -1.379   0.1824  
Canopy_cover:Temp_range          -0.002795   0.001983  -1.409   0.1735  
ADmedium:Canopy_cover:Temp_range  0.004172   0.002123   1.965   0.0628 .
ADhigh:Canopy_cover:Temp_range    0.002742   0.002197   1.248   0.2257  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.001540647)

    Null deviance: 0.056436  on 32  degrees of freedom
Residual deviance: 0.032354  on 21  degrees of freedom
AIC: -108.96

Number of Fisher Scoring iterations: 2

GLM_soil_AD_LAI_T_model2 <- stepAIC(GLM_soil_AD_LAI_T2, direction = "backward", Trace = FALSE)

Start:  AIC=-108.96
Soil_moisture_ave ~ AD * Canopy_cover * Temp_range

                             Df Deviance     AIC
<none>                          0.032354 -108.96
- AD:Canopy_cover:Temp_range  2 0.039195 -106.63

summary(GLM_soil_AD_LAI_T_model2)

Call:
glm(formula = Soil_moisture_ave ~ AD * Canopy_cover * Temp_range, 
    data = Ash)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.081414  -0.015334  -0.000157   0.009398   0.079531  

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)                      -0.921634   0.800870  -1.151   0.2628  
ADmedium                          1.597247   0.866373   1.844   0.0794 .
ADhigh                            1.257761   0.855518   1.470   0.1563  
Canopy_cover                      0.029502   0.021059   1.401   0.1758  
Temp_range                        0.099348   0.074217   1.339   0.1950  
ADmedium:Canopy_cover            -0.042537   0.022391  -1.900   0.0713 .
ADhigh:Canopy_cover              -0.031984   0.024499  -1.306   0.2058  
ADmedium:Temp_range              -0.153826   0.081514  -1.887   0.0730 .
ADhigh:Temp_range                -0.106541   0.077253  -1.379   0.1824  
Canopy_cover:Temp_range          -0.002795   0.001983  -1.409   0.1735  
ADmedium:Canopy_cover:Temp_range  0.004172   0.002123   1.965   0.0628 .
ADhigh:Canopy_cover:Temp_range    0.002742   0.002197   1.248   0.2257  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.001540647)

    Null deviance: 0.056436  on 32  degrees of freedom
Residual deviance: 0.032354  on 21  degrees of freedom
AIC: -108.96

Number of Fisher Scoring iterations: 2

Anova(GLM_soil_AD_LAI_T_model2)

Analysis of Deviance Table (Type II tests)

Response: Soil_moisture_ave
                           LR Chisq Df Pr(>Chisq)
AD                           3.3845  2     0.1841
Canopy_cover                 0.3254  1     0.5684
Temp_range                   0.2573  1     0.6120
AD:Canopy_cover              3.4023  2     0.1825
AD:Temp_range                0.1561  2     0.9249
Canopy_cover:Temp_range      0.8539  1     0.3555
AD:Canopy_cover:Temp_range   4.4405  2     0.1086

dropterm(GLM_soil_AD_LAI_T_model2,test = "Chisq")

Single term deletions

Model:
Soil_moisture_ave ~ AD * Canopy_cover * Temp_range
                           Df Deviance     AIC scaled dev. Pr(Chi)  
<none>                        0.032354 -108.96                      
AD:Canopy_cover:Temp_range  2 0.039195 -106.63        6.33 0.04221 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Factor 4: herb layer

#boxplot: AD ~ species richness (no significance)
p <- ggboxplot(Ash, x="AD", y= "Richness",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Plant species richness") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

#species richness lm (soil moisture, p = 0.006)
m6 <- lm(Richness ~ Temp_range, data = Ash)
m7 <- lm(Richness ~ Canopy_cover, data = Ash)
m8 <- lm(Richness ~ Soil_moisture_ave, data = Ash)
ggplotRegression(lm(Richness ~ Temp_range, data = Ash)) + ylab("Plant species richness") + xlab("Temperature range (°C)")

ggplotRegression(lm(Richness ~ Canopy_cover, data = Ash)) + ylab("Plant species richness") + xlab("Canopy cover (%)")

ggplotRegression(lm(Richness ~ Soil_moisture_ave, data = Ash)) + ylab("Plant species richness") + xlab("Soil moisture (m^3/m^3)")

tab_model(
  m6,m7,m8,
  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
  pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)"),
  dv.labels = c("Plant species richness"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
#species richness GLM 
#soil moisture is significant in the species richness model (p = 0.003369)
Richness_mod <- glm(Richness ~ AD + Temp_range + Canopy_cover + Soil_moisture_ave, data=Ash)
richness_model <- stepAIC(Richness_mod,direction = "backward", trace = FALSE)
Anova(richness_model)

Analysis of Deviance Table (Type II tests)

Response: Richness
                  LR Chisq Df Pr(>Chisq)   
Soil_moisture_ave   8.5962  1   0.003369 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##Now analysis individual species 
#MOSS
#univarable (significance with temperature range and canopy cover)
m9 <- lm(Moss_percent ~ Temp_range, data = Ash)
m10 <- lm(Moss_percent ~ Canopy_cover, data = Ash)
m11 <- lm(Moss_percent ~ Soil_moisture_ave, data = Ash)
tab_model(
  m9,m10,m11,
  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
  pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)"),
  dv.labels = c("Percentage of moss (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
ggplotRegression(lm(Moss_percent ~ Temp_range, data = Ash)) + xlab("Temperature range (°C)") + ylab("Percentage of moss (%)")

ggplotRegression(lm(Moss_percent ~ Canopy_cover, data = Ash))  + xlab("Canopy cover (%)") + ylab("Percentage of moss (%)")

ggplotRegression(lm(Moss_percent ~ Soil_moisture_ave, data = Ash))  + xlab("Soil moisture (m^3/m^3)") + ylab("Percentage of moss (%)")

#with AD (no significance)
p <- ggboxplot(Ash, x="AD", y= "Moss_percent",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Percentage of moss (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

kruskal.test(Moss_percent ~ AD, data = Ash)

    Kruskal-Wallis rank sum test

data:  Moss_percent by AD
Kruskal-Wallis chi-squared = 4.0124, df = 2, p-value = 0.1345

dunn.test(Ash$Moss_percent, Ash$AD, method = "bonferroni")

  Kruskal-Wallis rank sum test

data: x and group
Kruskal-Wallis chi-squared = 4.0124, df = 2, p-value = 0.13

                           Comparison of x by group                            
                                 (Bonferroni)                                  
Col Mean-|
Row Mean |       high        low
---------+----------------------
     low |   1.864456
         |     0.0934
         |
  medium |   0.082926  -1.658978
         |     1.0000     0.1457

alpha = 0.05
Reject Ho if p <= alpha/2

#no significance 
#very over disperse, so use quasibinomial
Moss_glm <- glm(cbind(Moss_percent,(100-Moss_percent)) ~ AD + Canopy_cover + Temp_range + Soil_moisture_ave + Understorey, data=Ash, family = "quasibinomial")
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Canopy_cover + Temp_range + Soil_moisture_ave + Understorey, 
    family = "quasibinomial", data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-9.273  -2.296   1.060   2.838  10.416  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -9.44457    3.86651  -2.443 0.021684 *  
ADmedium           2.87039    0.75093   3.822 0.000742 ***
ADhigh            -0.62644    0.90743  -0.690 0.496098    
Canopy_cover      -0.03538    0.02963  -1.194 0.243275    
Temp_range         1.12687    0.27853   4.046 0.000415 ***
Soil_moisture_ave -9.82194    7.66257  -1.282 0.211223    
Understorey+      -0.07384    0.68440  -0.108 0.914915    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 24.05511)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  605.02  on 26  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

Moss_glm<- update(Moss_glm, . ~ .-Soil_moisture_ave)
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Canopy_cover + Temp_range + Understorey, family = "quasibinomial", 
    data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-8.792  -2.537   1.377   2.741  10.014  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -11.65340    3.65564  -3.188 0.003608 ** 
ADmedium       2.57328    0.70661   3.642 0.001133 ** 
ADhigh        -1.02770    0.85660  -1.200 0.240662    
Canopy_cover  -0.02805    0.02963  -0.947 0.352075    
Temp_range     1.17839    0.28845   4.085 0.000353 ***
Understorey+   0.01984    0.67843   0.029 0.976888    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 25.03009)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  646.26  on 27  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

Moss_glm<- update(Moss_glm, . ~ .-Understorey)
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Canopy_cover + Temp_range, family = "quasibinomial", data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-8.840  -2.523   1.393   2.738   9.947  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -11.59944    3.09030  -3.753 0.000810 ***
ADmedium       2.57625    0.68566   3.757 0.000802 ***
ADhigh        -1.02669    0.83949  -1.223 0.231527    
Canopy_cover  -0.02812    0.02897  -0.970 0.340111    
Temp_range     1.17469    0.25390   4.627 7.69e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 24.05826)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  646.28  on 28  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

Moss_glm<- update(Moss_glm, . ~ .-Canopy_cover)
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Temp_range, family = "quasibinomial", data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-9.595  -3.067   1.280   2.803   9.549  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -13.2021     2.6263  -5.027 2.35e-05 ***
ADmedium      2.7639     0.6755   4.092 0.000312 ***
ADhigh       -0.6188     0.7091  -0.873 0.390015    
Temp_range    1.2191     0.2468   4.940 3.00e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 23.92961)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  668.88  on 29  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Anova(Moss_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Moss_percent, (100 - Moss_percent))
Error estimate based on Pearson residuals 

           Sum Sq Df F value    Pr(>F)    
AD         636.33  2  13.297 7.976e-05 ***
Temp_range 979.79  1  40.949 5.361e-07 ***
Residuals  693.88 29                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Moss_glm <- glm(cbind(Moss_percent,(100-Moss_percent)) ~ AD + Canopy_cover + Temp_range + Soil_moisture_ave , data=Ash, family = "quasibinomial")
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Canopy_cover + Temp_range + Soil_moisture_ave, family = "quasibinomial", 
    data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-9.354  -2.373   1.046   2.746  10.656  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -9.66178    3.27324  -2.952 0.006465 ** 
ADmedium           2.85286    0.72207   3.951 0.000504 ***
ADhigh            -0.63819    0.88759  -0.719 0.478306    
Canopy_cover      -0.03512    0.02917  -1.204 0.238930    
Temp_range         1.14037    0.24709   4.615 8.57e-05 ***
Soil_moisture_ave -9.68068    7.41019  -1.306 0.202436    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 23.46499)

    Null deviance: 2088.0  on 32  degrees of freedom
Residual deviance:  605.3  on 27  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

Moss_glm<- update(Moss_glm, . ~ .-Soil_moisture_ave)
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Canopy_cover + Temp_range, family = "quasibinomial", data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-8.840  -2.523   1.393   2.738   9.947  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -11.59944    3.09030  -3.753 0.000810 ***
ADmedium       2.57625    0.68566   3.757 0.000802 ***
ADhigh        -1.02669    0.83949  -1.223 0.231527    
Canopy_cover  -0.02812    0.02897  -0.970 0.340111    
Temp_range     1.17469    0.25390   4.627 7.69e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 24.05826)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  646.28  on 28  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

Moss_glm<- update(Moss_glm, . ~ .-Canopy_cover)
summary(Moss_glm)

Call:
glm(formula = cbind(Moss_percent, (100 - Moss_percent)) ~ AD + 
    Temp_range, family = "quasibinomial", data = Ash)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-9.595  -3.067   1.280   2.803   9.549  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -13.2021     2.6263  -5.027 2.35e-05 ***
ADmedium      2.7639     0.6755   4.092 0.000312 ***
ADhigh       -0.6188     0.7091  -0.873 0.390015    
Temp_range    1.2191     0.2468   4.940 3.00e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 23.92961)

    Null deviance: 2088.03  on 32  degrees of freedom
Residual deviance:  668.88  on 29  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Anova(Moss_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Moss_percent, (100 - Moss_percent))
Error estimate based on Pearson residuals 

           Sum Sq Df F value    Pr(>F)    
AD         636.33  2  13.297 7.976e-05 ***
Temp_range 979.79  1  40.949 5.361e-07 ***
Residuals  693.88 29                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#ASH SEEDLINGS 
p <- ggboxplot(Ash, x="AD", y= "Ash_percent",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Percentage of ash seedlings (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

#UNIVARIABLE (significant: canopy cover, temp range)
m12 <- lm(Ash_percent ~ Temp_range, data = Ash)
m13 <- lm(Ash_percent ~ Canopy_cover, data = Ash)
m14 <- lm(Ash_percent ~ Soil_moisture_ave, data = Ash)
ggplotRegression(lm(Ash_percent ~ Temp_range, data = Ash)) + xlab("Temperature range (°C)") + ylab("Percentage of ash seedlings (%)")

ggplotRegression(lm(Ash_percent ~ Canopy_cover, data = Ash)) + xlab("Canopy cover (%)") + ylab("Percentage of ash seedlings (%)")

ggplotRegression(lm(Ash_percent ~ Soil_moisture_ave, data = Ash)) + xlab("Soil moisture (m^3/m^3)") + ylab("Percentage of ash seedlings (%)")

tab_model(
  m12,m13,m14,
  show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
  pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)"),
  dv.labels = c("Percentage of ash seedlings (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
#glm very overdispersed (used quasibinomial)
Ash_glm <- glm(cbind(Ash_percent,(100-Ash_percent)) ~ AD + Canopy_cover + Temp_range + Soil_moisture_ave, data=Ash, family = "quasibinomial")
summary(Ash_glm)

Call:
glm(formula = cbind(Ash_percent, (100 - Ash_percent)) ~ AD + 
    Canopy_cover + Temp_range + Soil_moisture_ave, family = "quasibinomial", 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-5.7391  -3.1595  -0.3049   2.0570   4.3474  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)       -3.24131    1.93837  -1.672    0.106
ADmedium           0.82447    0.69207   1.191    0.244
ADhigh             0.70044    0.83805   0.836    0.411
Canopy_cover      -0.02827    0.02359  -1.198    0.241
Temp_range         0.07314    0.12222   0.598    0.555
Soil_moisture_ave  0.98091    4.60348   0.213    0.833

(Dispersion parameter for quasibinomial family taken to be 9.101293)

    Null deviance: 386.56  on 32  degrees of freedom
Residual deviance: 288.27  on 27  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

Anova(Ash_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Ash_percent, (100 - Ash_percent))
Error estimate based on Pearson residuals 

                   Sum Sq Df F value Pr(>F)
AD                 14.437  2  0.7932 0.4627
Canopy_cover       13.054  1  1.4344 0.2415
Temp_range          3.240  1  0.3560 0.5557
Soil_moisture_ave   0.412  1  0.0452 0.8332
Residuals         245.733 27               

Ash_glm<- update(Ash_glm, . ~ .-Soil_moisture_ave)
Anova(Ash_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Ash_percent, (100 - Ash_percent))
Error estimate based on Pearson residuals 

              Sum Sq Df F value Pr(>F)
AD            16.309  2  0.9247 0.4084
Canopy_cover  13.520  1  1.5331 0.2259
Temp_range     3.263  1  0.3700 0.5479
Residuals    246.921 28               

Ash_glm<- update(Ash_glm, . ~ .-Temp_range)
Anova(Ash_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Ash_percent, (100 - Ash_percent))
Error estimate based on Pearson residuals 

              Sum Sq Df F value Pr(>F)
AD            16.863  2  0.9932 0.3826
Canopy_cover  16.829  1  1.9824 0.1698
Residuals    246.190 29               

Ash_glm<- update(Ash_glm, . ~ .-AD)
Anova(Ash_glm, test.statistic = "F")

Analysis of Deviance Table (Type II tests)

Response: cbind(Ash_percent, (100 - Ash_percent))
Error estimate based on Pearson residuals 

              Sum Sq Df F value   Pr(>F)   
Canopy_cover  77.757  1  8.9891 0.005313 **
Residuals    268.157 31                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#BRAMBLE
p <- ggboxplot(Ash, x="AD", y= "Bramble_percent",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Percentage of bramble (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

m15 <- lm(Bramble_percent ~ Temp_range, data = Ash)
m16 <- lm(Bramble_percent ~ Canopy_cover, data = Ash)
m17 <- lm(Bramble_percent ~ Soil_moisture_ave, data = Ash)
ggplotRegression(lm(Bramble_percent ~ Temp_range, data = Ash)) + xlab("Temperature range (°C)") + ylab("Percentage of bramble (%)")

ggplotRegression(lm(Bramble_percent ~ Canopy_cover, data = Ash)) + xlab("Canopy cover (%)") + ylab("Percentage of bramble (%)")

ggplotRegression(lm(Bramble_percent ~ Soil_moisture_ave, data = Ash)) + xlab("Soil moisture (m^3/m^3)") + ylab("Percentage of bramble (%)")
bramble_mod <- glm(Bramble_percent ~ AD * Tree_dens * Canopy_cover, data=Ash)
par(mfrow=c(2,2))

plot(bramble_mod)

summary(bramble_mod)

Call:
glm(formula = Bramble_percent ~ AD * Tree_dens * Canopy_cover, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-33.094  -12.495   -4.108    7.739   79.706  

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)
(Intercept)                      -19.020    144.448  -0.132    0.896
ADmedium                         114.054    159.045   0.717    0.481
ADhigh                            65.660    150.995   0.435    0.668
Tree_dens                        124.803    307.407   0.406    0.689
Canopy_cover                       1.896      3.466   0.547    0.590
ADmedium:Tree_dens              -315.669    360.866  -0.875    0.392
ADhigh:Tree_dens                -180.187    320.978  -0.561    0.580
ADmedium:Canopy_cover             -3.257      3.913  -0.832    0.415
ADhigh:Canopy_cover               -3.756      4.157  -0.904    0.376
Tree_dens:Canopy_cover            -6.203      7.884  -0.787    0.440
ADmedium:Tree_dens:Canopy_cover    9.641      8.935   1.079    0.293
ADhigh:Tree_dens:Canopy_cover      9.282      9.340   0.994    0.332

(Dispersion parameter for gaussian family taken to be 673.0568)

    Null deviance: 18805  on 32  degrees of freedom
Residual deviance: 14134  on 21  degrees of freedom
AIC: 319.62

Number of Fisher Scoring iterations: 2

Anova(bramble_mod)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
                          LR Chisq Df Pr(>Chisq)
AD                         1.32272  2     0.5161
Tree_dens                  0.95608  1     0.3282
Canopy_cover               0.55259  1     0.4573
AD:Tree_dens               2.36412  2     0.3066
AD:Canopy_cover            0.16290  2     0.9218
Tree_dens:Canopy_cover     0.42000  1     0.5169
AD:Tree_dens:Canopy_cover  1.24544  2     0.5365

tab_model(
    m15,m16,m17,
    show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
    pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)"),
    dv.labels = c("Percentage of bramble (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
#remove the 3 interaction factor
bramble_mod2 <- glm(Bramble_percent ~ AD + Tree_dens + Canopy_cover
                          + AD:Tree_dens + AD:Canopy_cover + Tree_dens:Canopy_cover, data=Ash)
Anova(bramble_mod2)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
                       LR Chisq Df Pr(>Chisq)
AD                      1.36759  2     0.5047
Tree_dens               0.98851  1     0.3201
Canopy_cover            0.57133  1     0.4497
AD:Tree_dens            2.44431  2     0.2946
AD:Canopy_cover         0.16842  2     0.9192
Tree_dens:Canopy_cover  0.43425  1     0.5099

#remove 2 interaction factors
bramble_mod3 <- glm(Bramble_percent ~ AD + Tree_dens + Canopy_cover
                    + AD:Tree_dens + Tree_dens:Canopy_cover, data=Ash)
Anova(bramble_mod3)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
                       LR Chisq Df Pr(>Chisq)
AD                      1.47570  2     0.4781
Tree_dens               1.62449  1     0.2025
Canopy_cover            0.61650  1     0.4324
AD:Tree_dens            2.69046  2     0.2605
Tree_dens:Canopy_cover  0.39162  1     0.5314

bramble_mod4 <- glm(Bramble_percent ~ AD + Tree_dens + Canopy_cover
                    + AD:Tree_dens , data=Ash)
Anova(bramble_mod4)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
             LR Chisq Df Pr(>Chisq)
AD            1.27696  2     0.5281
Tree_dens     1.66342  1     0.1971
Canopy_cover  0.63127  1     0.4269
AD:Tree_dens  2.85746  2     0.2396

bramble_mod5 <- glm(Bramble_percent ~ AD + Tree_dens + AD:Tree_dens , data=Ash)
Anova(bramble_mod5)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
             LR Chisq Df Pr(>Chisq)
AD            0.85516  2     0.6521
Tree_dens     1.94704  1     0.1629
AD:Tree_dens  2.72217  2     0.2564

bramble_mod6 <- glm(Bramble_percent ~ AD + Tree_dens , data=Ash)
Anova(bramble_mod6)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
          LR Chisq Df Pr(>Chisq)
AD         0.83439  2     0.6589
Tree_dens  1.89974  1     0.1681

bramble_mod7 <- glm(Bramble_percent ~  Tree_dens , data=Ash)
Anova(bramble_mod7)

Analysis of Deviance Table (Type II tests)

Response: Bramble_percent
          LR Chisq Df Pr(>Chisq)
Tree_dens   1.8064  1     0.1789

#result: bramble cover percentage has no significant correlation with all the factors we investigated 
#DOG MERCURY
p <- ggboxplot(Ash, x="AD", y= "Dogs_Mercury_percent",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Percentage of dogs' mercury (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

m18 <- lm(Dogs_Mercury_percent ~ Temp_range, data = Ash)
m19 <- lm(Dogs_Mercury_percent ~ Canopy_cover, data = Ash)
m20 <- lm(Dogs_Mercury_percent ~ Soil_moisture_ave, data = Ash)
ggplotRegression(lm(Dogs_Mercury_percent ~ Temp_range, data = Ash)) + xlab("Temperature range (°C)") + ylab("Percentage of dogs' mercury (%)")

ggplotRegression(lm(Dogs_Mercury_percent ~ Canopy_cover, data = Ash)) + xlab("Canopy cover (%)") + ylab("Percentage of dogs' mercury (%)")

ggplotRegression(lm(Dogs_Mercury_percent ~ Soil_moisture_ave, data = Ash)) + xlab("Soil moisture (m^3/m^3)") + ylab("Percentage of dogs' mercury (%)")
dogs_mercury_mod <- glm(Dogs_Mercury_percent ~ AD * Tree_dens * Canopy_cover, data=Ash)
par(mfrow=c(2,2))

plot(dogs_mercury_mod)

summary(dogs_mercury_mod)

Call:
glm(formula = Dogs_Mercury_percent ~ AD * Tree_dens * Canopy_cover, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-29.192   -8.382   -4.205    4.178   29.320  

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)
(Intercept)                     -30.658851  93.063995  -0.329    0.745
ADmedium                         18.223883 102.468323   0.178    0.861
ADhigh                           76.334845  97.281625   0.785    0.441
Tree_dens                       -33.235373 198.054048  -0.168    0.868
Canopy_cover                      0.563493   2.232967   0.252    0.803
ADmedium:Tree_dens               33.014481 232.495951   0.142    0.888
ADhigh:Tree_dens                -48.478043 206.797212  -0.234    0.817
ADmedium:Canopy_cover            -0.001679   2.521215  -0.001    0.999
ADhigh:Canopy_cover              -2.302356   2.678109  -0.860    0.400
Tree_dens:Canopy_cover            2.286326   5.079412   0.450    0.657
ADmedium:Tree_dens:Canopy_cover  -2.502278   5.756872  -0.435    0.668
ADhigh:Tree_dens:Canopy_cover     1.773715   6.017339   0.295    0.771

(Dispersion parameter for gaussian family taken to be 279.3772)

    Null deviance: 7813.6  on 32  degrees of freedom
Residual deviance: 5866.9  on 21  degrees of freedom
AIC: 290.61

Number of Fisher Scoring iterations: 2

Anova(dogs_mercury_mod)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
                          LR Chisq Df Pr(>Chisq)
AD                         2.34408  2     0.3097
Tree_dens                  0.01950  1     0.8889
Canopy_cover               1.55846  1     0.2119
AD:Tree_dens               1.53114  2     0.4651
AD:Canopy_cover            2.31893  2     0.3137
Tree_dens:Canopy_cover     0.74472  1     0.3882
AD:Tree_dens:Canopy_cover  1.04803  2     0.5921

tab_model(
    m18,m19,m20,
    show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
    pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)"),
    dv.labels = c("Percentage of dog's mercury (%)"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
#remove factors
dogs_mercury_mod2 <- glm(Dogs_Mercury_percent ~ AD + Tree_dens + Canopy_cover
                    + AD:Tree_dens + AD:Canopy_cover + Tree_dens:Canopy_cover, data=Ash)
Anova(dogs_mercury_mod2)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
                       LR Chisq Df Pr(>Chisq)
AD                      2.44530  2     0.2944
Tree_dens               0.02034  1     0.8866
Canopy_cover            1.62575  1     0.2023
AD:Tree_dens            1.59725  2     0.4499
AD:Canopy_cover         2.41906  2     0.2983
Tree_dens:Canopy_cover  0.77687  1     0.3781

dogs_mercury_mod3 <- glm(Dogs_Mercury_percent ~ AD + Tree_dens + Canopy_cover
                        + AD:Canopy_cover + Tree_dens:Canopy_cover, data=Ash)
Anova(dogs_mercury_mod3)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
                       LR Chisq Df Pr(>Chisq)
AD                      2.48533  2     0.2886
Tree_dens               0.02067  1     0.8857
Canopy_cover            1.42746  1     0.2322
AD:Canopy_cover         1.63665  2     0.4412
Tree_dens:Canopy_cover  1.11788  1     0.2904

dogs_mercury_mod4 <- glm(Dogs_Mercury_percent ~ AD + Tree_dens + Canopy_cover
                         + Tree_dens:Canopy_cover, data=Ash)
Anova(dogs_mercury_mod4)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
                       LR Chisq Df Pr(>Chisq)
AD                      2.51924  2     0.2838
Tree_dens               0.08684  1     0.7682
Canopy_cover            1.44694  1     0.2290
Tree_dens:Canopy_cover  0.67523  1     0.4112

dogs_mercury_mod5 <- glm(Dogs_Mercury_percent ~ AD + Tree_dens + Canopy_cover, data=Ash)
Anova(dogs_mercury_mod5)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
             LR Chisq Df Pr(>Chisq)
AD            2.13112  2     0.3445
Tree_dens     0.08786  1     0.7669
Canopy_cover  1.46391  1     0.2263

dogs_mercury_mod6 <- glm(Dogs_Mercury_percent ~ AD+ Canopy_cover, data=Ash)
Anova(dogs_mercury_mod6)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
             LR Chisq Df Pr(>Chisq)
AD             2.2475  2     0.3251
Canopy_cover   1.4426  1     0.2297

dogs_mercury_mod7 <- glm(Dogs_Mercury_percent ~ Canopy_cover, data=Ash)
Anova(dogs_mercury_mod7)

Analysis of Deviance Table (Type II tests)

Response: Dogs_Mercury_percent
             LR Chisq Df Pr(>Chisq)
Canopy_cover  0.14951  1      0.699

#our result showed that none of the factors we investigated was correlated to the percentage cover of dog's mercury.

Factor 5: invertebrate species

New_ash <- read.csv("Newest_ash.csv")
#add second day data in 
Ash$Invert_sp_0107_0207 <- New_ash$No.of.invertebrate.species..1.7.2.7.[2:34][match(New_ash$Site[2:34],Ash$Site)]
#ash dieback ~ invert
p <- ggboxplot(Ash, x="AD", y= "total_invert",color = "AD", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree dieback", y = "Invertebrate species richness") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p

kruskal.test(total_invert ~ AD, data = Ash)

    Kruskal-Wallis rank sum test

data:  total_invert by AD
Kruskal-Wallis chi-squared = 1.4921, df = 2, p-value = 0.4742

dunnTest(Ash$total_invert, Ash$AD, method = "bonferroni")

Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Bonferroni method.

#Univariable
m21 <- lm(total_invert ~ Temp_range , data = Ash)
m22 <- lm(total_invert ~ Canopy_cover, data = Ash)
m23 <- lm(total_invert ~ Soil_moisture_ave, data = Ash)
m24 <- lm(total_invert ~ Moss_percent, data = Ash)
m25 <- lm(total_invert ~ Ash_percent, data = Ash)
tab_model(
    m21,m22,m23,m24,m25,
    show.se = TRUE, show.ci = FALSE,show.stat = TRUE,
    pred.labels = c("Intercept","Temperature range (C)","Canopy cover (%)","Soil moisture (%)","Percentage of moss (%)","Percentage of ash seedlings(%)"),
    dv.labels = c("Invertebrate species richness"),
    string.pred = "Coefficient",
    string.se = "std.Error",
    string.p = "P-Value", string.stat = "F value"
)
ggplotRegression(lm(total_invert ~ Temp_range , data = Ash)) + xlab("Temperature range (°C)") + ylab("Invertebrate species richness")

ggplotRegression(lm(total_invert ~ Canopy_cover, data = Ash))  + xlab("Canopy cover (%)") + ylab("Invertebrate species richness")

ggplotRegression(lm(total_invert ~ Soil_moisture_ave, data = Ash))  + xlab("Soil moisture (m^3/m^3)") + ylab("Invertebrate species richness")

ggplotRegression(lm(total_invert ~ Moss_percent, data = Ash))  + xlab("Percentage of moss (%)") + ylab("Invertebrate species richness")

ggplotRegression(lm(total_invert ~ Ash_percent, data = Ash))  + xlab("Percentage of ash seedlings (%)") + ylab("Invertebrate species richness")

#no correlation at all!
#GLM for everything (found a weak trend with moss percent)
GLM_invert <- glm(total_invert ~ AD + Canopy_cover + Understorey + Temp_range + Soil_moisture_ave + Moss_percent + Ash_percent, data = Ash)
summary(GLM_invert)

Call:
glm(formula = total_invert ~ AD + Canopy_cover + Understorey + 
    Temp_range + Soil_moisture_ave + Moss_percent + Ash_percent, 
    data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-8.3493  -1.7895  -0.7565   2.9291   9.3105  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)        6.61191    9.01573   0.733   0.4704  
ADmedium          -1.96347    3.25441  -0.603   0.5520  
ADhigh             1.29416    3.25783   0.397   0.6947  
Canopy_cover      -0.02723    0.11599  -0.235   0.8164  
Understorey+       4.02496    2.12309   1.896   0.0701 .
Temp_range         0.17768    0.87890   0.202   0.8415  
Soil_moisture_ave 21.85829   23.32481   0.937   0.3580  
Moss_percent      -0.04367    0.04805  -0.909   0.3725  
Ash_percent       -0.15090    0.13100  -1.152   0.2607  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 23.04334)

    Null deviance: 782.91  on 32  degrees of freedom
Residual deviance: 553.04  on 24  degrees of freedom
AIC: 206.67

Number of Fisher Scoring iterations: 2

GLM_invert_model <- stepAIC(GLM_invert, direction = "backward", Trace = FALSE)

Start:  AIC=206.67
total_invert ~ AD + Canopy_cover + Understorey + Temp_range + 
    Soil_moisture_ave + Moss_percent + Ash_percent

                    Df Deviance    AIC
- AD                 2   570.44 203.70
- Temp_range         1   553.98 204.73
- Canopy_cover       1   554.31 204.75
- Moss_percent       1   572.07 205.79
- Soil_moisture_ave  1   573.28 205.86
- Ash_percent        1   583.62 206.45
<none>                   553.04 206.67
- Understorey        1   635.86 209.28

Step:  AIC=203.7
total_invert ~ Canopy_cover + Understorey + Temp_range + Soil_moisture_ave + 
    Moss_percent + Ash_percent

                    Df Deviance    AIC
- Soil_moisture_ave  1   586.85 202.63
- Canopy_cover       1   593.14 202.98
- Temp_range         1   595.54 203.12
<none>                   570.44 203.70
- Ash_percent        1   625.97 204.76
- Understorey        1   654.18 206.22
- Moss_percent       1   682.40 207.61

Step:  AIC=202.63
total_invert ~ Canopy_cover + Understorey + Temp_range + Moss_percent + 
    Ash_percent

               Df Deviance    AIC
- Canopy_cover  1   616.47 202.26
- Temp_range    1   617.86 202.33
<none>              586.85 202.63
- Ash_percent   1   639.72 203.48
- Understorey   1   675.84 205.29
- Moss_percent  1   705.68 206.72

Step:  AIC=202.26
total_invert ~ Understorey + Temp_range + Moss_percent + Ash_percent

               Df Deviance    AIC
- Ash_percent   1   647.40 201.87
<none>              616.47 202.26
- Temp_range    1   663.15 202.67
- Understorey   1   688.70 203.91
- Moss_percent  1   716.67 205.23

Step:  AIC=201.87
total_invert ~ Understorey + Temp_range + Moss_percent

               Df Deviance    AIC
- Temp_range    1   667.41 200.88
<none>              647.40 201.87
- Understorey   1   692.00 202.07
- Moss_percent  1   725.19 203.62

Step:  AIC=200.88
total_invert ~ Understorey + Moss_percent

               Df Deviance    AIC
- Understorey   1   698.23 200.37
<none>              667.41 200.88
- Moss_percent  1   725.32 201.62

Step:  AIC=200.37
total_invert ~ Moss_percent

               Df Deviance    AIC
<none>              698.23 200.37
- Moss_percent  1   782.91 202.15

summary(GLM_invert_model)

Call:
glm(formula = total_invert ~ Moss_percent, data = Ash)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-7.3813  -3.0946  -0.9513   2.9531  10.5135  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  12.72856    2.00716   6.342 4.66e-07 ***
Moss_percent -0.04777    0.02464  -1.939   0.0617 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 22.5237)

    Null deviance: 782.91  on 32  degrees of freedom
Residual deviance: 698.23  on 31  degrees of freedom
AIC: 200.37

Number of Fisher Scoring iterations: 2

Anova(GLM_invert_model)

Analysis of Deviance Table (Type II tests)

Response: total_invert
             LR Chisq Df Pr(>Chisq)  
Moss_percent   3.7593  1    0.05251 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dropterm(GLM_invert_model,test = "Chisq")

Single term deletions

Model:
total_invert ~ Moss_percent
             Df Deviance    AIC scaled dev. Pr(Chi)  
<none>            698.23 200.37                      
Moss_percent  1   782.91 202.15      3.7772 0.05195 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1