1.1 Factors: Canopy cover percent (continuous) ~ Ash dieback levels (categories)
#check normality of data
library(dplyr)
library(ggpubr)
#select out the investigating two factors: Ash dieback level & canopy cover percent
LAI_AD <- select(Ash, Ash_dieback,Canopy_cover_percent)
#check the data
dplyr::sample_n(LAI_AD, 10)
#from central limit theorm, if the data size is larger than 30, no normality test is needed. but here just to make sure i know what's going on, I will still check the normality. plus our sample size is only a bit larger than 30 (n=33)
p <- ggqqplot(LAI_AD$Canopy_cover_percent)
p
#using Shapiro-Wilk’s method to test normality (p-value = 0.05575)
shapiro.test(LAI_AD$Canopy_cover_percent)
Shapiro-Wilk normality test
data: LAI_AD$Canopy_cover_percent
W = 0.93704, p-value = 0.05575
#change the factor levels of Ash_dieback
LAI_AD$Ash_dieback <- factor(LAI_AD$Ash_dieback, levels = c("low", "medium", "high"))
#boxplot
p <- ggboxplot(LAI_AD, x="Ash_dieback", y= "Canopy_cover_percent",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Canopy Cover (%)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
#to be safe, I decided to use non-parametric test (Kruskal-Wallis) p-value = 0.0001066 (very significant difference!)
kruskal.test(Canopy_cover_percent ~ Ash_dieback, data = LAI_AD)
Kruskal-Wallis rank sum test
data: Canopy_cover_percent by Ash_dieback
Kruskal-Wallis chi-squared = 18.293, df = 2, p-value = 0.0001066
#Comparions between groups
pairwise.wilcox.test(LAI_AD$Canopy_cover_percent, LAI_AD$Ash_dieback, p.adjust.method = "BH")
Pairwise comparisons using Wilcoxon rank sum test
data: LAI_AD$Canopy_cover_percent and LAI_AD$Ash_dieback
low medium
medium 0.71970 -
high 0.00014 0.00030
P value adjustment method: BH
1.2 Factors: Canopy cover pecent (continuous) ~ Basal area of all ash trees (categories) We measured the diameter at breast height (DBH) of all the alive ash trees in a 36m^2 square for each site, using the location of insect trap as the central of the square we want to compare the canopy cover present with the basal area of all alive ash trees in 36m^2.
library(ggpmisc)
#Formula for Basal area: BA (m) = 0.0001 * Pi *(Diameter at breast height/2)^2
#total basal area in 36 m^2: sum of all BA
BA_tree_density <- read.csv("tree_density_n_BA.csv")
View(BA_tree_density)
#Add data of Canopy cover precent and understorey
colnames(BA_tree_density)[colnames(BA_tree_density)=="sum_of_BA.of.all.trees"] <- "BA"
BA_tree_density$Canopy_cover <- Ash$Canopy_cover_percent[match(Ash$Site,BA_tree_density$Site)]
BA_tree_density$understorey <- Ash$Understorey[match(Ash$Site,BA_tree_density$Site)]
lm(formula = BA_tree_density$BA ~ BA_tree_density$Canopy_cover)
Call:
lm(formula = BA_tree_density$BA ~ BA_tree_density$Canopy_cover)
Coefficients:
(Intercept) BA_tree_density$Canopy_cover
0.74323 0.02565
#Basal area of alive ash trees vs canopy cover (w/ & w/o understorey)
p <- ggplot(data = BA_tree_density, aes(x = BA, y = Canopy_cover)) +
geom_smooth(method = "lm", se=FALSE, color="black") +
stat_poly_eq(formula = BA_tree_density$Canopy_cover ~ BA_tree_density$BA,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) +
geom_point()
p
#excluded trees that have understorey, compare canopy cover to basal area
no_understorey <- BA_tree_density[BA_tree_density$understorey == "-",]
p <- ggplot(data = no_understorey, aes(x = BA, y = Canopy_cover)) +
geom_smooth(method = "lm", se=FALSE, color="black") +
stat_poly_eq(formula = no_understorey$Canopy_cover ~ no_understorey$BA,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE) + geom_point()
p
#in both case, no correlation was found.
#same as tree density vs canopy cover
GLM of all three factors
library(MASS)
Canopy_AD <- glm(Canopy_cover_digital ~ AD, data=BA_tree_density, family=binomial)
non-integer #successes in a binomial glm!
summary(Canopy_AD)
Call:
glm(formula = Canopy_cover_digital ~ AD, family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.38170 -0.10968 -0.01488 0.14751 0.33650
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4114 0.6808 -0.604 0.546
ADmedium -0.1063 0.9439 -0.113 0.910
ADhigh -1.0434 0.9639 -1.083 0.279
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2.7100 on 32 degrees of freedom
Residual deviance: 1.1765 on 30 degrees of freedom
AIC: 30.376
Number of Fisher Scoring iterations: 4
Canopy_AD_BA <- glm(Canopy_cover_digital ~ AD * BA, data=BA_tree_density, family=binomial)
non-integer #successes in a binomial glm!
summary(Canopy_AD_BA)
Call:
glm(formula = Canopy_cover_digital ~ AD * BA, family = binomial,
data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.41791 -0.09278 0.03795 0.12261 0.28923
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.58136 1.09930 -0.529 0.597
ADmedium -0.19437 1.69695 -0.115 0.909
ADhigh -0.70209 1.66101 -0.423 0.673
BA 0.08047 0.40478 0.199 0.842
ADmedium:BA 0.10574 0.89327 0.118 0.906
ADhigh:BA -0.22029 0.95940 -0.230 0.818
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2.7100 on 32 degrees of freedom
Residual deviance: 1.0561 on 27 degrees of freedom
AIC: 36.443
Number of Fisher Scoring iterations: 4
Canopy_AD_BA_Tree_Dens <- glm(Canopy_cover_digital ~ AD * BA * Tree_density , data=BA_tree_density, family=binomial)
non-integer #successes in a binomial glm!
summary(Canopy_AD_BA_Tree_Dens)
Call:
glm(formula = Canopy_cover_digital ~ AD * BA * Tree_density,
family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.33295 -0.12780 0.01957 0.09279 0.28204
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02576 3.24023 -0.008 0.994
ADmedium -2.23713 5.84936 -0.382 0.702
ADhigh -1.27396 3.88475 -0.328 0.743
BA 0.01454 1.15299 0.013 0.990
Tree_density -1.65815 9.09261 -0.182 0.855
ADmedium:BA 0.59760 2.76878 0.216 0.829
ADhigh:BA -0.23721 2.32475 -0.102 0.919
ADmedium:Tree_density 5.12777 14.24482 0.360 0.719
ADhigh:Tree_density 1.85616 10.58344 0.175 0.861
BA:Tree_density 0.19393 3.18503 0.061 0.951
ADmedium:BA:Tree_density -1.31681 6.47186 -0.203 0.839
ADhigh:BA:Tree_density -0.13146 4.50075 -0.029 0.977
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2.7100 on 32 degrees of freedom
Residual deviance: 0.7605 on 21 degrees of freedom
AIC: 48.623
Number of Fisher Scoring iterations: 4
#NOTE: got very different results when trying to fit the binomial result through dividing canopy cover by 100, rather than using cbind() to get a limit
Canopy_AD <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD, data=BA_tree_density, 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 = BA_tree_density)
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 = BA_tree_density)
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
Canopy_AD_BA <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD * BA, data=BA_tree_density, family=binomial)
non-integer counts in a binomial glm!
summary(Canopy_AD_BA)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA,
family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.1791 -0.9278 0.3795 1.2261 2.8923
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.58136 0.10993 -5.288 1.23e-07 ***
ADmedium -0.19437 0.16969 -1.145 0.2520
ADhigh -0.70209 0.16610 -4.227 2.37e-05 ***
BA 0.08047 0.04048 1.988 0.0468 *
ADmedium:BA 0.10574 0.08933 1.184 0.2365
ADhigh:BA -0.22029 0.09594 -2.296 0.0217 *
---
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: 105.61 on 27 degrees of freedom
AIC: 274.11
Number of Fisher Scoring iterations: 4
step.model2 <- stepAIC(Canopy_AD_BA, direction = "backward",
trace = FALSE)
non-integer #successes in a binomial glm!
summary(step.model2)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA,
family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.1791 -0.9278 0.3795 1.2261 2.8923
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.58136 0.10993 -5.288 1.23e-07 ***
ADmedium -0.19437 0.16969 -1.145 0.2520
ADhigh -0.70209 0.16610 -4.227 2.37e-05 ***
BA 0.08047 0.04048 1.988 0.0468 *
ADmedium:BA 0.10574 0.08933 1.184 0.2365
ADhigh:BA -0.22029 0.09594 -2.296 0.0217 *
---
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: 105.61 on 27 degrees of freedom
AIC: 274.11
Number of Fisher Scoring iterations: 4
dropterm(step.model2,test = "Chisq")
non-integer #successes in a binomial glm!
Single term deletions
Model:
cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA
Df Deviance AIC LRT Pr(Chi)
<none> 105.61 274.11
AD:BA 2 113.83 278.33 8.2252 0.01637 *
---
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_Tree_Dens <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD * BA * Tree_density , data=BA_tree_density, family=binomial)
non-integer counts in a binomial glm!
summary(Canopy_AD_BA_Tree_Dens)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA *
Tree_density, family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3295 -1.2780 0.1957 0.9279 2.8204
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02576 0.32402 -0.080 0.936629
ADmedium -2.23713 0.58494 -3.825 0.000131 ***
ADhigh -1.27396 0.38848 -3.279 0.001040 **
BA 0.01454 0.11530 0.126 0.899646
Tree_density -1.65815 0.90926 -1.824 0.068209 .
ADmedium:BA 0.59760 0.27688 2.158 0.030902 *
ADhigh:BA -0.23721 0.23248 -1.020 0.307547
ADmedium:Tree_density 5.12777 1.42448 3.600 0.000319 ***
ADhigh:Tree_density 1.85616 1.05834 1.754 0.079458 .
BA:Tree_density 0.19393 0.31850 0.609 0.542614
ADmedium:BA:Tree_density -1.31681 0.64719 -2.035 0.041884 *
ADhigh:BA:Tree_density -0.13146 0.45008 -0.292 0.770223
---
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: 76.05 on 21 degrees of freedom
AIC: 257.16
Number of Fisher Scoring iterations: 4
step.model3 <- stepAIC(Canopy_AD_BA_Tree_Dens, direction = "backward", trace = FALSE)
non-integer #successes in a binomial glm!
summary(step.model3)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA *
Tree_density, family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.3295 -1.2780 0.1957 0.9279 2.8204
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02576 0.32402 -0.080 0.936629
ADmedium -2.23713 0.58494 -3.825 0.000131 ***
ADhigh -1.27396 0.38848 -3.279 0.001040 **
BA 0.01454 0.11530 0.126 0.899646
Tree_density -1.65815 0.90926 -1.824 0.068209 .
ADmedium:BA 0.59760 0.27688 2.158 0.030902 *
ADhigh:BA -0.23721 0.23248 -1.020 0.307547
ADmedium:Tree_density 5.12777 1.42448 3.600 0.000319 ***
ADhigh:Tree_density 1.85616 1.05834 1.754 0.079458 .
BA:Tree_density 0.19393 0.31850 0.609 0.542614
ADmedium:BA:Tree_density -1.31681 0.64719 -2.035 0.041884 *
ADhigh:BA:Tree_density -0.13146 0.45008 -0.292 0.770223
---
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: 76.05 on 21 degrees of freedom
AIC: 257.16
Number of Fisher Scoring iterations: 4
dropterm(step.model3,test = "Chisq")
non-integer #successes in a binomial glm!
Single term deletions
Model:
cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA * Tree_density
Df Deviance AIC LRT Pr(Chi)
<none> 76.050 257.16
AD:BA:Tree_density 2 80.391 257.50 4.3411 0.1141
step.model31 <- update(step.model3, . ~ .-AD:BA:Tree_density)
non-integer counts in a binomial glm!
dropterm(step.model31,test = "Chisq")
non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!
Single term deletions
Model:
cbind(Canopy_cover, Canopy_not_cover) ~ AD + BA + Tree_density +
AD:BA + AD:Tree_density + BA:Tree_density
Df Deviance AIC LRT Pr(Chi)
<none> 80.391 257.34
AD:BA 2 85.986 258.93 5.5947 0.06097 .
AD:Tree_density 2 100.504 273.45 20.1130 4.291e-05 ***
BA:Tree_density 1 80.436 255.38 0.0452 0.83159
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
step.model32 <- update(step.model31, . ~ .-BA:Tree_density)
non-integer counts in a binomial glm!
Anova(step.model32)
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!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 138.161 2 < 2.2e-16 ***
BA 3.462 1 0.06279 .
Tree_density 4.405 1 0.03584 *
AD:BA 6.126 2 0.04674 *
AD:Tree_density 20.766 2 3.095e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Canopy_AD_BA_Tree_Dens_Understorey <- glm(cbind(Canopy_cover,Canopy_not_cover) ~ AD * BA * Tree_density * understorey , data=BA_tree_density, family=binomial)
non-integer counts in a binomial glm!
summary(Canopy_AD_BA_Tree_Dens_Understorey)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD * BA *
Tree_density * understorey, family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3765 -0.7984 0.0000 0.5993 2.5233
Coefficients: (3 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.63989 0.71361 -0.897 0.369882
ADmedium -3.76506 5.80376 -0.649 0.516514
ADhigh -0.72826 0.57671 -1.263 0.206666
BA 0.49626 0.44672 1.111 0.266609
Tree_density -0.96781 2.82476 -0.343 0.731887
understorey+ 0.53800 1.03176 0.521 0.602061
ADmedium:BA 3.23276 7.11447 0.454 0.649547
ADhigh:BA -1.49287 0.45194 -3.303 0.000956 ***
ADmedium:Tree_density 5.04788 2.89529 1.743 0.081250 .
ADhigh:Tree_density 1.67027 2.16226 0.772 0.439838
BA:Tree_density -0.51166 1.64771 -0.311 0.756160
ADmedium:understorey+ 1.32257 5.28517 0.250 0.802401
ADhigh:understorey+ 1.75988 0.99406 1.770 0.076662 .
BA:understorey+ -0.48291 0.43721 -1.105 0.269360
Tree_density:understorey+ -0.23314 3.13743 -0.074 0.940765
ADmedium:BA:Tree_density -1.42421 0.84542 -1.685 0.092063 .
ADhigh:BA:Tree_density 1.91763 0.86299 2.222 0.026278 *
ADmedium:BA:understorey+ -2.52360 6.99356 -0.361 0.718214
ADhigh:BA:understorey+ 0.08938 0.49984 0.179 0.858076
ADmedium:Tree_density:understorey+ NA NA NA NA
ADhigh:Tree_density:understorey+ -4.76878 2.41008 -1.979 0.047852 *
BA:Tree_density:understorey+ 0.66038 1.52835 0.432 0.665677
ADmedium:BA:Tree_density:understorey+ NA NA NA NA
ADhigh:BA:Tree_density:understorey+ NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 271.000 on 32 degrees of freedom
Residual deviance: 48.566 on 12 degrees of freedom
AIC: 247.21
Number of Fisher Scoring iterations: 4
step.model4 <- stepAIC(Canopy_AD_BA_Tree_Dens_Understorey, 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 #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!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!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!non-integer #successes in a binomial glm!
summary(step.model4)
Call:
glm(formula = cbind(Canopy_cover, Canopy_not_cover) ~ AD + BA +
Tree_density + understorey + AD:BA + AD:Tree_density + BA:Tree_density +
AD:understorey + Tree_density:understorey + AD:BA:Tree_density +
AD:Tree_density:understorey, family = binomial, data = BA_tree_density)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4860 -0.9217 0.0000 0.8513 2.5211
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1518 0.3837 -0.396 0.692424
ADmedium -1.3479 1.4997 -0.899 0.368778
ADhigh -1.1104 0.4433 -2.505 0.012250 *
BA 0.1149 0.2181 0.527 0.598477
Tree_density -1.4203 1.0904 -1.303 0.192700
understorey+ -0.3619 0.6786 -0.533 0.593811
ADmedium:BA 0.6076 0.3431 1.771 0.076556 .
ADhigh:BA -1.3596 0.3681 -3.693 0.000221 ***
ADmedium:Tree_density 3.3109 3.7300 0.888 0.374744
ADhigh:Tree_density 1.9200 1.2211 1.572 0.115853
BA:Tree_density -0.1013 0.5829 -0.174 0.862060
ADmedium:understorey+ -0.6829 1.6246 -0.420 0.674239
ADhigh:understorey+ 2.3630 0.8035 2.941 0.003273 **
Tree_density:understorey+ 1.2721 1.8020 0.706 0.480251
ADmedium:BA:Tree_density -1.1742 0.8218 -1.429 0.153063
ADhigh:BA:Tree_density 1.9746 0.7737 2.552 0.010707 *
ADmedium:Tree_density:understorey+ 0.6843 4.0400 0.169 0.865490
ADhigh:Tree_density:understorey+ -5.7881 2.0542 -2.818 0.004837 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 271.000 on 32 degrees of freedom
Residual deviance: 50.404 on 15 degrees of freedom
AIC: 242.99
Number of Fisher Scoring iterations: 4
Factor 2: temperature fluctuation vs AD & LAI
#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="Ash_dieback", y= "Max_temp",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Maximum temerature (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
#Wilcoxon rank sum test:
kruskal.test(Max_temp ~ Ash_dieback, data = Ash)
Kruskal-Wallis rank sum test
data: Max_temp by Ash_dieback
Kruskal-Wallis chi-squared = 13.316, df = 2, p-value = 0.001283
pairwise.wilcox.test(Ash$Max_temp, Ash$Ash_dieback, p.adjust.method = "BH")
cannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with ties
Pairwise comparisons using Wilcoxon rank sum test
data: Ash$Max_temp and Ash$Ash_dieback
low medium
medium 0.9020 -
high 0.0045 0.0045
P value adjustment method: BH
#significance is between medium and high levels of AD (p=0.0045)
#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="Ash_dieback", y= "Min_temp",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Minimum temerature (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
#Wilcoxon rank sum test:
kruskal.test(Min_temp ~ Ash_dieback, data = Ash)
Kruskal-Wallis rank sum test
data: Min_temp by Ash_dieback
Kruskal-Wallis chi-squared = 4.9176, df = 2, p-value = 0.08554
pairwise.wilcox.test(Ash$Min_temp, Ash$Ash_dieback, p.adjust.method = "BH")
cannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with ties
Pairwise comparisons using Wilcoxon rank sum test
data: Ash$Min_temp and Ash$Ash_dieback
low medium
medium 0.62 -
high 0.16 0.16
P value adjustment method: BH
#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="Ash_dieback", y= "Temp_range",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Temerature range (°C)") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
kruskal.test(Temp_range ~ Ash_dieback, data = Ash)
Kruskal-Wallis rank sum test
data: Temp_range by Ash_dieback
Kruskal-Wallis chi-squared = 15.28, df = 2, p-value = 0.0004809
pairwise.wilcox.test(Ash$Temp_range, Ash$Ash_dieback, p.adjust.method = "BH")
cannot compute exact p-value with tiescannot compute exact p-value with tiescannot compute exact p-value with ties
Pairwise comparisons using Wilcoxon rank sum test
data: Ash$Temp_range and Ash$Ash_dieback
low medium
medium 0.8371 -
high 0.0021 0.0021
P value adjustment method: BH
#significance between medium and high (0.0021)
#temp ~ BA
#ggplot linear regression
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)))
}
ggplotRegression(lm(Temp_range ~ BA, data = Ash))
ggplotRegression(lm(Temp_range ~ Tree_density, data = Ash))
#temp ~ LAI
ggplotRegression(lm(Canopy_cover_percent ~ Temp_range, data = Ash))
ggplotRegression(lm(Canopy_cover_percent ~ Min_temp, data = Ash))
ggplotRegression(lm(Canopy_cover_percent ~ Max_temp, data = Ash))
#GLM for temperature change ~ canopy cover * AD * BA* tree density
GLM_temp <- glm(Temp_range ~ Canopy_cover_percent * Ash_dieback * BA * Tree_density, data = Ash)
summary(GLM_temp)
Call:
glm(formula = Temp_range ~ Canopy_cover_percent * Ash_dieback *
BA * Tree_density, data = Ash)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.27335 -0.13813 0.00032 0.09401 2.02072
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.9384 40.8115 2.032 0.0727 .
Canopy_cover_percent -1.4899 0.9244 -1.612 0.1415
Ash_diebackmedium -133.6513 55.3879 -2.413 0.0391 *
Ash_diebackhigh -60.2098 41.1296 -1.464 0.1773
BA -37.5738 24.7889 -1.516 0.1639
Tree_density -275.4060 135.3355 -2.035 0.0723 .
Canopy_cover_percent:Ash_diebackmedium 3.0346 1.2666 2.396 0.0402 *
Canopy_cover_percent:Ash_diebackhigh 0.9876 0.9742 1.014 0.3372
Canopy_cover_percent:BA 0.7912 0.5378 1.471 0.1753
Ash_diebackmedium:BA 68.1478 31.8993 2.136 0.0614 .
Ash_diebackhigh:BA 28.0950 25.7757 1.090 0.3040
Canopy_cover_percent:Tree_density 5.7586 3.0270 1.902 0.0895 .
Ash_diebackmedium:Tree_density 430.9651 164.6840 2.617 0.0280 *
Ash_diebackhigh:Tree_density 260.9985 136.3133 1.915 0.0878 .
BA:Tree_density 163.0052 77.8763 2.093 0.0658 .
Canopy_cover_percent:Ash_diebackmedium:BA -1.5497 0.7000 -2.214 0.0541 .
Canopy_cover_percent:Ash_diebackhigh:BA -0.2393 0.6864 -0.349 0.7354
Canopy_cover_percent:Ash_diebackmedium:Tree_density -9.5463 3.6839 -2.591 0.0291 *
Canopy_cover_percent:Ash_diebackhigh:Tree_density -4.8981 3.1757 -1.542 0.1574
Canopy_cover_percent:BA:Tree_density -3.4950 1.7001 -2.056 0.0700 .
Ash_diebackmedium:BA:Tree_density -245.7428 92.4031 -2.659 0.0261 *
Ash_diebackhigh:BA:Tree_density -144.1904 79.6957 -1.809 0.1039
Canopy_cover_percent:Ash_diebackmedium:BA:Tree_density 5.4946 2.0188 2.722 0.0235 *
Canopy_cover_percent:Ash_diebackhigh:BA:Tree_density 2.3554 1.9735 1.194 0.2632
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 1.73008)
Null deviance: 139.320 on 32 degrees of freedom
Residual deviance: 15.571 on 9 degrees of freedom
AIC: 118.86
Number of Fisher Scoring iterations: 2
temp_model <- stepAIC(GLM_temp, direction = "backward", trace = FALSE)
summary(temp_model)
Call:
glm(formula = Temp_range ~ Canopy_cover_percent * Ash_dieback *
BA * Tree_density, data = Ash)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.27335 -0.13813 0.00032 0.09401 2.02072
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.9384 40.8115 2.032 0.0727 .
Canopy_cover_percent -1.4899 0.9244 -1.612 0.1415
Ash_diebackmedium -133.6513 55.3879 -2.413 0.0391 *
Ash_diebackhigh -60.2098 41.1296 -1.464 0.1773
BA -37.5738 24.7889 -1.516 0.1639
Tree_density -275.4060 135.3355 -2.035 0.0723 .
Canopy_cover_percent:Ash_diebackmedium 3.0346 1.2666 2.396 0.0402 *
Canopy_cover_percent:Ash_diebackhigh 0.9876 0.9742 1.014 0.3372
Canopy_cover_percent:BA 0.7912 0.5378 1.471 0.1753
Ash_diebackmedium:BA 68.1478 31.8993 2.136 0.0614 .
Ash_diebackhigh:BA 28.0950 25.7757 1.090 0.3040
Canopy_cover_percent:Tree_density 5.7586 3.0270 1.902 0.0895 .
Ash_diebackmedium:Tree_density 430.9651 164.6840 2.617 0.0280 *
Ash_diebackhigh:Tree_density 260.9985 136.3133 1.915 0.0878 .
BA:Tree_density 163.0052 77.8763 2.093 0.0658 .
Canopy_cover_percent:Ash_diebackmedium:BA -1.5497 0.7000 -2.214 0.0541 .
Canopy_cover_percent:Ash_diebackhigh:BA -0.2393 0.6864 -0.349 0.7354
Canopy_cover_percent:Ash_diebackmedium:Tree_density -9.5463 3.6839 -2.591 0.0291 *
Canopy_cover_percent:Ash_diebackhigh:Tree_density -4.8981 3.1757 -1.542 0.1574
Canopy_cover_percent:BA:Tree_density -3.4950 1.7001 -2.056 0.0700 .
Ash_diebackmedium:BA:Tree_density -245.7428 92.4031 -2.659 0.0261 *
Ash_diebackhigh:BA:Tree_density -144.1904 79.6957 -1.809 0.1039
Canopy_cover_percent:Ash_diebackmedium:BA:Tree_density 5.4946 2.0188 2.722 0.0235 *
Canopy_cover_percent:Ash_diebackhigh:BA:Tree_density 2.3554 1.9735 1.194 0.2632
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 1.73008)
Null deviance: 139.320 on 32 degrees of freedom
Residual deviance: 15.571 on 9 degrees of freedom
AIC: 118.86
Number of Fisher Scoring iterations: 2
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="Ash_dieback", y= "Soil_moisture_ave",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", 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 ~ Ash_dieback,data = Ash)
summary(SM_AD)
Df Sum Sq Mean Sq F value Pr(>F)
Ash_dieback 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
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Soil_moisture_ave ~ Ash_dieback, data = Ash)
$Ash_dieback
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
ggplotRegression(lm(Canopy_cover_percent ~ Soil_moisture_ave, data = Ash))
ggplotRegression(lm(BA ~ Soil_moisture_ave, data = Ash))
ggplotRegression(lm(Tree_density ~ Soil_moisture_ave, data = Ash))
ggplotRegression(lm(Temp_range ~ Soil_moisture_ave, data = Ash))
#GLM for AD ~ Soil moisture * BA * Tree density * Temp
GLM_soil <- glm(Soil_moisture_ave ~ Canopy_cover_percent * Ash_dieback * BA * Tree_density * Temp_range,data = Ash, family = "binomial")
non-integer #successes in a binomial glm!
summary(GLM_soil)
Call:
glm(formula = Soil_moisture_ave ~ Canopy_cover_percent * Ash_dieback *
BA * Tree_density * Temp_range, family = "binomial", data = Ash)
Deviance Residuals:
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Coefficients: (15 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.996e+03 1.633e+04 0.245 0.807
Canopy_cover_percent -8.657e+01 3.557e+02 -0.243 0.808
Ash_diebackmedium -1.683e+04 6.821e+04 -0.247 0.805
Ash_diebackhigh -7.561e+03 3.071e+04 -0.246 0.806
BA 5.718e+03 2.316e+04 0.247 0.805
Tree_density -6.048e+03 2.451e+04 -0.247 0.805
Temp_range -2.123e+02 8.715e+02 -0.244 0.808
Canopy_cover_percent:Ash_diebackmedium 4.883e+02 1.978e+03 0.247 0.805
Canopy_cover_percent:Ash_diebackhigh 3.146e+02 1.276e+03 0.246 0.805
Canopy_cover_percent:BA -2.086e+02 8.462e+02 -0.247 0.805
Ash_diebackmedium:BA 4.641e+03 1.885e+04 0.246 0.805
Ash_diebackhigh:BA 3.240e+03 1.319e+04 0.246 0.806
Canopy_cover_percent:Tree_density 2.025e+02 8.276e+02 0.245 0.807
Ash_diebackmedium:Tree_density 7.650e+03 3.094e+04 0.247 0.805
Ash_diebackhigh:Tree_density -7.576e+03 3.092e+04 -0.245 0.806
BA:Tree_density -5.401e+03 2.192e+04 -0.246 0.805
Canopy_cover_percent:Temp_range 5.279e-01 2.727e+00 0.194 0.846
Ash_diebackmedium:Temp_range 9.726e+02 3.936e+03 0.247 0.805
Ash_diebackhigh:Temp_range 4.683e+02 1.904e+03 0.246 0.806
BA:Temp_range -6.298e+02 2.555e+03 -0.246 0.805
Tree_density:Temp_range 1.079e+03 4.387e+03 0.246 0.806
Canopy_cover_percent:Ash_diebackmedium:BA -6.606e+01 2.689e+02 -0.246 0.806
Canopy_cover_percent:Ash_diebackhigh:BA -1.090e+02 4.439e+02 -0.246 0.806
Canopy_cover_percent:Ash_diebackmedium:Tree_density -2.875e+02 1.165e+03 -0.247 0.805
Canopy_cover_percent:Ash_diebackhigh:Tree_density -1.345e+00 1.622e+01 -0.083 0.934
Canopy_cover_percent:BA:Tree_density 9.864e+01 4.004e+02 0.246 0.805
Ash_diebackmedium:BA:Tree_density 1.006e+03 4.086e+03 0.246 0.806
Ash_diebackhigh:BA:Tree_density 3.405e+03 1.383e+04 0.246 0.805
Canopy_cover_percent:Ash_diebackmedium:Temp_range -2.733e+01 1.106e+02 -0.247 0.805
Canopy_cover_percent:Ash_diebackhigh:Temp_range -1.734e+01 7.035e+01 -0.246 0.805
Canopy_cover_percent:BA:Temp_range 2.311e+01 9.382e+01 0.246 0.805
Ash_diebackmedium:BA:Temp_range -2.304e+02 9.327e+02 -0.247 0.805
Ash_diebackhigh:BA:Temp_range NA NA NA NA
Canopy_cover_percent:Tree_density:Temp_range -2.101e+01 8.609e+01 -0.244 0.807
Ash_diebackmedium:Tree_density:Temp_range NA NA NA NA
Ash_diebackhigh:Tree_density:Temp_range NA NA NA NA
BA:Tree_density:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackmedium:BA:Tree_density NA NA NA NA
Canopy_cover_percent:Ash_diebackhigh:BA:Tree_density NA NA NA NA
Canopy_cover_percent:Ash_diebackmedium:BA:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackhigh:BA:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackmedium:Tree_density:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackhigh:Tree_density:Temp_range NA NA NA NA
Canopy_cover_percent:BA:Tree_density:Temp_range NA NA NA NA
Ash_diebackmedium:BA:Tree_density:Temp_range NA NA NA NA
Ash_diebackhigh:BA:Tree_density:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackmedium:BA:Tree_density:Temp_range NA NA NA NA
Canopy_cover_percent:Ash_diebackhigh:BA:Tree_density:Temp_range NA NA NA NA
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 4.0847e-01 on 32 degrees of freedom
Residual deviance: 9.0661e-16 on 0 degrees of freedom
AIC: 77.845
Number of Fisher Scoring iterations: 5
soil_model <- stepAIC(GLM_soil, direction = "backward", trace = FALSE)
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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!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!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!non-integer #successes in a binomial glm!
summary(soil_model)
Call:
glm(formula = Soil_moisture_ave ~ 1, family = "binomial", data = Ash)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.24904 -0.06700 -0.01971 0.05264 0.28338
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.634 0.471 -3.47 0.000521 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 0.40847 on 32 degrees of freedom
Residual deviance: 0.40847 on 32 degrees of freedom
AIC: 13.763
Number of Fisher Scoring iterations: 4
Factor 4: herb layer
#species richness GLM
##shall we still include BA and tree density?
Richness_mod <- glm(Richness ~ Ash_dieback * Tree_dens * Canopy_cover_percent * BA, 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)
Ash_dieback 0.5610 2 0.75542
Tree_dens 0.0870 1 0.76801
Canopy_cover_percent 4.8879 1 0.02704 *
BA 0.6675 1 0.41391
Ash_dieback:Tree_dens 2.0571 2 0.35752
Ash_dieback:Canopy_cover_percent 2.5442 2 0.28024
Tree_dens:Canopy_cover_percent 5.2365 1 0.02212 *
Ash_dieback:BA 0.7278 2 0.69496
Tree_dens:BA 3.9472 1 0.04695 *
Canopy_cover_percent:BA 2.9479 1 0.08599 .
Ash_dieback:Tree_dens:Canopy_cover_percent 5.3569 2 0.06867 .
Ash_dieback:Tree_dens:BA 5.5832 2 0.06132 .
Ash_dieback:Canopy_cover_percent:BA 3.2940 2 0.19263
Tree_dens:Canopy_cover_percent:BA 3.5941 1 0.05799 .
Ash_dieback:Tree_dens:Canopy_cover_percent:BA 3.7045 2 0.15689
---
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)
ggplotRegression(lm(Temp_range ~ Moss_percent, data = Ash))
ggplotRegression(lm(Tree_dens ~ Moss_percent, data = Ash))
ggplotRegression(lm(Canopy_cover_percent ~ Moss_percent, data = Ash))
ggplotRegression(lm(BA ~ Moss_percent, data = Ash))
ggplotRegression(lm(Tree_dens ~ Moss_percent, data = Ash))
ggplotRegression(lm(Soil_moisture_ave ~ Moss_percent, data = Ash))
#wish AD (no significance)
p <- ggboxplot(Ash, x="Ash_dieback", y= "Moss_percent",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Moss_percent") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
kruskal.test(Moss_percent ~ Ash_dieback, data = Ash)
Kruskal-Wallis rank sum test
data: Moss_percent by Ash_dieback
Kruskal-Wallis chi-squared = 4.0124, df = 2, p-value = 0.1345
Moss_glm <- glm(cbind(Moss_percent,(100-Moss_percent)) ~ Ash_dieback * Tree_dens * Canopy_cover_percent * BA * Temp_range * Soil_moisture_ave, data=Ash, family = "binomial")
Moss_model <- stepAIC(Moss_glm,direction = "backward", trace = FALSE)
glm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurred
Anova(Moss_model)
glm.fit: fitted probabilities numerically 0 or 1 occurred
Analysis of Deviance Table (Type II tests)
Response: cbind(Moss_percent, (100 - Moss_percent))
LR Chisq Df Pr(>Chisq)
Ash_dieback 641.79 2 < 2.2e-16 ***
Tree_dens 101.59 1 < 2.2e-16 ***
Canopy_cover_percent 31.27 1 2.240e-08 ***
BA 17.62 1 2.701e-05 ***
Temp_range 194.60 1 < 2.2e-16 ***
Soil_moisture_ave 69.81 1 < 2.2e-16 ***
Ash_dieback:Tree_dens 50.56 2 1.050e-11 ***
Ash_dieback:Canopy_cover_percent 16.72 2 0.0002343 ***
Tree_dens:Canopy_cover_percent 28.65 1 8.651e-08 ***
Ash_dieback:BA 6.78 2 0.0337894 *
Tree_dens:BA 4.10 1 0.0428482 *
Tree_dens:Temp_range 39.63 1 3.071e-10 ***
BA:Temp_range 84.03 1 < 2.2e-16 ***
Tree_dens:Soil_moisture_ave 110.99 1 < 2.2e-16 ***
Canopy_cover_percent:Soil_moisture_ave 45.87 1 1.263e-11 ***
BA:Soil_moisture_ave 132.46 1 < 2.2e-16 ***
Ash_dieback:Tree_dens:Canopy_cover_percent 98.95 2 < 2.2e-16 ***
Ash_dieback:Tree_dens:BA 98.97 2 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dropterm(Moss_model,test = "Chisq")
glm.fit: fitted probabilities numerically 0 or 1 occurred
Single term deletions
Model:
cbind(Moss_percent, (100 - Moss_percent)) ~ Ash_dieback + Tree_dens +
Canopy_cover_percent + BA + Temp_range + Soil_moisture_ave +
Ash_dieback:Tree_dens + Ash_dieback:Canopy_cover_percent +
Tree_dens:Canopy_cover_percent + Ash_dieback:BA + Tree_dens:BA +
Tree_dens:Temp_range + BA:Temp_range + Tree_dens:Soil_moisture_ave +
Canopy_cover_percent:Soil_moisture_ave + BA:Soil_moisture_ave +
Ash_dieback:Tree_dens:Canopy_cover_percent + Ash_dieback:Tree_dens:BA
Df Deviance AIC LRT Pr(Chi)
<none> 0.000 141.00
Tree_dens:Temp_range 1 39.629 178.63 39.629 3.071e-10 ***
BA:Temp_range 1 84.027 223.03 84.027 < 2.2e-16 ***
Tree_dens:Soil_moisture_ave 1 110.990 250.00 110.990 < 2.2e-16 ***
Canopy_cover_percent:Soil_moisture_ave 1 45.870 184.88 45.870 1.263e-11 ***
BA:Soil_moisture_ave 1 132.462 271.47 132.462 < 2.2e-16 ***
Ash_dieback:Tree_dens:Canopy_cover_percent 2 98.951 235.96 98.951 < 2.2e-16 ***
Ash_dieback:Tree_dens:BA 2 98.969 235.97 98.969 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#ASH SEEDLINGS
#UNIVARIABLE (significant: canopy cover, temp range)
ggplotRegression(lm(Temp_range ~ Ash_percent, data = Ash))
ggplotRegression(lm(Tree_dens ~ Ash_percent, data = Ash))
ggplotRegression(lm(Canopy_cover_percent ~ Ash_percent, data = Ash))
ggplotRegression(lm(BA ~ Ash_percent, data = Ash))
ggplotRegression(lm(Tree_dens ~ Ash_percent, data = Ash))
ggplotRegression(lm(Soil_moisture_ave ~ Ash_percent, data = Ash))
ash_seedlings_mod <- glm(Ash_percent ~ Ash_dieback * Tree_dens * Canopy_cover_percent, data=Ash)
#with ash dieback (no significance)
p <- ggboxplot(Ash, x="Ash_dieback", y= "Ash_percent",color = "Ash_dieback", palette = c("#00AFBB", "#E7B800", "#FC4E07")) + labs(x="Extent of ash tree disback", y = "Ash_percent") + geom_jitter(shape=19, position=position_jitter(0.2)) + geom_point() + theme(legend.position = "none")
p
kruskal.test(Ash_percent ~ Ash_dieback, data = Ash)
Kruskal-Wallis rank sum test
data: Ash_percent by Ash_dieback
Kruskal-Wallis chi-squared = 5.3715, df = 2, p-value = 0.06817
Ash_glm <- glm(cbind(Ash_percent,(100-Ash_percent)) ~ Ash_dieback * Tree_dens * Canopy_cover_percent * BA * Temp_range * Soil_moisture_ave, data=Ash, family = "binomial")
Ash_model <- stepAIC(Ash_glm,direction = "backward", trace = FALSE)
glm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: algorithm did not convergeglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurred
Anova(Ash_model)
glm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurred
Analysis of Deviance Table (Type II tests)
Response: cbind(Ash_percent, (100 - Ash_percent))
LR Chisq Df Pr(>Chisq)
Ash_dieback 17.269 2 0.0001779 ***
Tree_dens 6.149 1 0.0131459 *
Canopy_cover_percent 7.871 1 0.0050225 **
BA 26.262 1 2.981e-07 ***
Temp_range 14.386 1 0.0001489 ***
Soil_moisture_ave 0.003 1 0.9549460
Ash_dieback:Tree_dens 52.521 2 3.938e-12 ***
Ash_dieback:Canopy_cover_percent 25.164 2 3.434e-06 ***
Tree_dens:Canopy_cover_percent 7.696 1 0.0055349 **
Ash_dieback:BA 58.239 2 2.257e-13 ***
Tree_dens:BA 25.255 1 5.023e-07 ***
Canopy_cover_percent:BA 22.840 1 1.761e-06 ***
Ash_dieback:Temp_range 89.444 2 < 2.2e-16 ***
Tree_dens:Temp_range 49.677 1 1.813e-12 ***
BA:Temp_range 76.944 1 < 2.2e-16 ***
Tree_dens:Soil_moisture_ave 31.462 1 2.034e-08 ***
Canopy_cover_percent:Soil_moisture_ave 66.153 1 4.172e-16 ***
BA:Soil_moisture_ave 42.456 1 7.228e-11 ***
Ash_dieback:Tree_dens:Canopy_cover_percent 26.045 2 2.210e-06 ***
Ash_dieback:Canopy_cover_percent:BA 40.621 2 1.511e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
dropterm(Ash_model,test = "Chisq")
glm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurredglm.fit: fitted probabilities numerically 0 or 1 occurred
Single term deletions
Model:
cbind(Ash_percent, (100 - Ash_percent)) ~ Ash_dieback + Tree_dens +
Canopy_cover_percent + BA + Temp_range + Soil_moisture_ave +
Ash_dieback:Tree_dens + Ash_dieback:Canopy_cover_percent +
Tree_dens:Canopy_cover_percent + Ash_dieback:BA + Tree_dens:BA +
Canopy_cover_percent:BA + Ash_dieback:Temp_range + Tree_dens:Temp_range +
BA:Temp_range + Tree_dens:Soil_moisture_ave + Canopy_cover_percent:Soil_moisture_ave +
BA:Soil_moisture_ave + Ash_dieback:Tree_dens:Canopy_cover_percent +
Ash_dieback:Canopy_cover_percent:BA
Df Deviance AIC LRT Pr(Chi)
<none> 0.000 148.29
Tree_dens:BA 1 25.255 171.54 25.255 5.023e-07 ***
Ash_dieback:Temp_range 2 89.444 233.73 89.444 < 2.2e-16 ***
Tree_dens:Temp_range 1 49.677 195.97 49.677 1.813e-12 ***
BA:Temp_range 1 76.944 223.24 76.944 < 2.2e-16 ***
Tree_dens:Soil_moisture_ave 1 31.462 177.75 31.462 2.034e-08 ***
Canopy_cover_percent:Soil_moisture_ave 1 66.153 212.44 66.153 4.172e-16 ***
BA:Soil_moisture_ave 1 42.456 188.75 42.456 7.228e-11 ***
Ash_dieback:Tree_dens:Canopy_cover_percent 2 26.045 170.34 26.045 2.210e-06 ***
Ash_dieback:Canopy_cover_percent:BA 2 40.621 184.91 40.621 1.511e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#the following is from previous analysis during the fieldtrip
par(mfrow=c(2,2))
plot(ash_seedlings_mod)
summary(ash_seedlings_mod)
Call:
glm(formula = Ash_percent ~ Ash_dieback * Tree_dens * Canopy_cover_percent,
data = Ash)
Deviance Residuals:
Min 1Q Median 3Q Max
-15.874 -4.363 -1.117 5.233 16.080
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.32308 54.10089 0.117 0.908
Ash_diebackmedium 37.66678 59.56791 0.632 0.534
Ash_diebackhigh -4.61589 56.55273 -0.082 0.936
Tree_dens -10.00647 115.13476 -0.087 0.932
Canopy_cover_percent -0.08259 1.29809 -0.064 0.950
Ash_diebackmedium:Tree_dens -66.13009 135.15687 -0.489 0.630
Ash_diebackhigh:Tree_dens 47.70873 120.21742 0.397 0.695
Ash_diebackmedium:Canopy_cover_percent -0.75444 1.46566 -0.515 0.612
Ash_diebackhigh:Canopy_cover_percent 0.71861 1.55687 0.462 0.649
Tree_dens:Canopy_cover_percent 0.27787 2.95281 0.094 0.926
Ash_diebackmedium:Tree_dens:Canopy_cover_percent 1.34158 3.34664 0.401 0.693
Ash_diebackhigh:Tree_dens:Canopy_cover_percent -2.36391 3.49806 -0.676 0.507
(Dispersion parameter for gaussian family taken to be 94.41399)
Null deviance: 2972.9 on 32 degrees of freedom
Residual deviance: 1982.7 on 21 degrees of freedom
AIC: 254.81
Number of Fisher Scoring iterations: 2
Anova(ash_seedlings_mod)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.26881 2 0.5303
Tree_dens 0.00119 1 0.9725
Canopy_cover_percent 0.59310 1 0.4412
Ash_dieback:Tree_dens 0.09621 2 0.9530
Ash_dieback:Canopy_cover_percent 0.27100 2 0.8733
Tree_dens:Canopy_cover_percent 0.01045 1 0.9186
Ash_dieback:Tree_dens:Canopy_cover_percent 2.29272 2 0.3178
#remove 3 interaction factor
ash_seedlings_mod2 <- glm(Ash_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Tree_dens + Ash_dieback:Canopy_cover_percent + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod2)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.25287 2 0.5345
Tree_dens 0.00117 1 0.9727
Canopy_cover_percent 0.58565 1 0.4441
Ash_dieback:Tree_dens 0.09500 2 0.9536
Ash_dieback:Canopy_cover_percent 0.26760 2 0.8748
Tree_dens:Canopy_cover_percent 0.01032 1 0.9191
#remove 2 interaction factors
ash_seedlings_mod3 <- glm(Ash_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Canopy_cover_percent + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod3)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.35621 2 0.5076
Tree_dens 0.00127 1 0.9716
Canopy_cover_percent 1.30394 1 0.2535
Ash_dieback:Canopy_cover_percent 0.41227 2 0.8137
Tree_dens:Canopy_cover_percent 0.02942 1 0.8638
ash_seedlings_mod4 <- glm(Ash_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod4)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.23640 2 0.5389
Tree_dens 0.00132 1 0.9710
Canopy_cover_percent 1.35451 1 0.2445
Ash_dieback:Canopy_cover_percent 0.57341 2 0.7507
ash_seedlings_mod4 <- glm(Ash_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod4)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.30277 2 0.5213
Tree_dens 0.13567 1 0.7126
Canopy_cover_percent 1.42722 1 0.2322
ash_seedlings_mod5 <- glm(Ash_percent ~ Ash_dieback + Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod5)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.2090 2 0.5463
Canopy_cover_percent 1.6127 1 0.2041
ash_seedlings_mod6 <- glm(Ash_percent ~ Canopy_cover_percent, data=Ash)
Anova(ash_seedlings_mod6)
Analysis of Deviance Table (Type II tests)
Response: Ash_percent
LR Chisq Df Pr(>Chisq)
Canopy_cover_percent 7.9649 1 0.004769 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#result: ash seedlings percentage correlates to canopy cover percentage (p-value=0.0048)
#BRAMBLE
bramble_mod <- glm(Bramble_percent ~ Ash_dieback * Tree_dens * Canopy_cover_percent, data=Ash)
par(mfrow=c(2,2))
plot(bramble_mod)
summary(bramble_mod)
Call:
glm(formula = Bramble_percent ~ Ash_dieback * Tree_dens * Canopy_cover_percent,
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
Ash_diebackmedium 114.054 159.045 0.717 0.481
Ash_diebackhigh 65.660 150.995 0.435 0.668
Tree_dens 124.803 307.407 0.406 0.689
Canopy_cover_percent 1.896 3.466 0.547 0.590
Ash_diebackmedium:Tree_dens -315.669 360.866 -0.875 0.392
Ash_diebackhigh:Tree_dens -180.187 320.978 -0.561 0.580
Ash_diebackmedium:Canopy_cover_percent -3.257 3.913 -0.832 0.415
Ash_diebackhigh:Canopy_cover_percent -3.756 4.157 -0.904 0.376
Tree_dens:Canopy_cover_percent -6.203 7.884 -0.787 0.440
Ash_diebackmedium:Tree_dens:Canopy_cover_percent 9.641 8.935 1.079 0.293
Ash_diebackhigh:Tree_dens:Canopy_cover_percent 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)
Ash_dieback 1.32272 2 0.5161
Tree_dens 0.95608 1 0.3282
Canopy_cover_percent 0.55259 1 0.4573
Ash_dieback:Tree_dens 2.36412 2 0.3066
Ash_dieback:Canopy_cover_percent 0.16290 2 0.9218
Tree_dens:Canopy_cover_percent 0.42000 1 0.5169
Ash_dieback:Tree_dens:Canopy_cover_percent 1.24544 2 0.5365
#remove the 3 interaction factor
bramble_mod2 <- glm(Bramble_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Tree_dens + Ash_dieback:Canopy_cover_percent + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(bramble_mod2)
Analysis of Deviance Table (Type II tests)
Response: Bramble_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.36759 2 0.5047
Tree_dens 0.98851 1 0.3201
Canopy_cover_percent 0.57133 1 0.4497
Ash_dieback:Tree_dens 2.44431 2 0.2946
Ash_dieback:Canopy_cover_percent 0.16842 2 0.9192
Tree_dens:Canopy_cover_percent 0.43425 1 0.5099
#remove 2 interaction factors
bramble_mod3 <- glm(Bramble_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Tree_dens + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(bramble_mod3)
Analysis of Deviance Table (Type II tests)
Response: Bramble_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.47570 2 0.4781
Tree_dens 1.62449 1 0.2025
Canopy_cover_percent 0.61650 1 0.4324
Ash_dieback:Tree_dens 2.69046 2 0.2605
Tree_dens:Canopy_cover_percent 0.39162 1 0.5314
bramble_mod4 <- glm(Bramble_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Tree_dens , data=Ash)
Anova(bramble_mod4)
Analysis of Deviance Table (Type II tests)
Response: Bramble_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 1.27696 2 0.5281
Tree_dens 1.66342 1 0.1971
Canopy_cover_percent 0.63127 1 0.4269
Ash_dieback:Tree_dens 2.85746 2 0.2396
bramble_mod5 <- glm(Bramble_percent ~ Ash_dieback + Tree_dens + Ash_dieback:Tree_dens , data=Ash)
Anova(bramble_mod5)
Analysis of Deviance Table (Type II tests)
Response: Bramble_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 0.85516 2 0.6521
Tree_dens 1.94704 1 0.1629
Ash_dieback:Tree_dens 2.72217 2 0.2564
bramble_mod6 <- glm(Bramble_percent ~ Ash_dieback + Tree_dens , data=Ash)
Anova(bramble_mod6)
Analysis of Deviance Table (Type II tests)
Response: Bramble_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 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
dogs_mercury_mod <- glm(Dogs_Mercury_percent ~ Ash_dieback * Tree_dens * Canopy_cover_percent, data=Ash)
par(mfrow=c(2,2))
plot(dogs_mercury_mod)
summary(dogs_mercury_mod)
Call:
glm(formula = Dogs_Mercury_percent ~ Ash_dieback * Tree_dens *
Canopy_cover_percent, 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
Ash_diebackmedium 18.223883 102.468323 0.178 0.861
Ash_diebackhigh 76.334845 97.281625 0.785 0.441
Tree_dens -33.235373 198.054048 -0.168 0.868
Canopy_cover_percent 0.563493 2.232967 0.252 0.803
Ash_diebackmedium:Tree_dens 33.014481 232.495951 0.142 0.888
Ash_diebackhigh:Tree_dens -48.478043 206.797212 -0.234 0.817
Ash_diebackmedium:Canopy_cover_percent -0.001679 2.521215 -0.001 0.999
Ash_diebackhigh:Canopy_cover_percent -2.302356 2.678109 -0.860 0.400
Tree_dens:Canopy_cover_percent 2.286326 5.079412 0.450 0.657
Ash_diebackmedium:Tree_dens:Canopy_cover_percent -2.502278 5.756872 -0.435 0.668
Ash_diebackhigh:Tree_dens:Canopy_cover_percent 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)
Ash_dieback 2.34408 2 0.3097
Tree_dens 0.01950 1 0.8889
Canopy_cover_percent 1.55846 1 0.2119
Ash_dieback:Tree_dens 1.53114 2 0.4651
Ash_dieback:Canopy_cover_percent 2.31893 2 0.3137
Tree_dens:Canopy_cover_percent 0.74472 1 0.3882
Ash_dieback:Tree_dens:Canopy_cover_percent 1.04803 2 0.5921
#remove factors
dogs_mercury_mod2 <- glm(Dogs_Mercury_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Tree_dens + Ash_dieback:Canopy_cover_percent + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod2)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 2.44530 2 0.2944
Tree_dens 0.02034 1 0.8866
Canopy_cover_percent 1.62575 1 0.2023
Ash_dieback:Tree_dens 1.59725 2 0.4499
Ash_dieback:Canopy_cover_percent 2.41906 2 0.2983
Tree_dens:Canopy_cover_percent 0.77687 1 0.3781
dogs_mercury_mod3 <- glm(Dogs_Mercury_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Ash_dieback:Canopy_cover_percent + Tree_dens:Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod3)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 2.48533 2 0.2886
Tree_dens 0.02067 1 0.8857
Canopy_cover_percent 1.42746 1 0.2322
Ash_dieback:Canopy_cover_percent 1.63665 2 0.4412
Tree_dens:Canopy_cover_percent 1.11788 1 0.2904
dogs_mercury_mod4 <- glm(Dogs_Mercury_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent
+ Tree_dens:Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod4)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 2.51924 2 0.2838
Tree_dens 0.08684 1 0.7682
Canopy_cover_percent 1.44694 1 0.2290
Tree_dens:Canopy_cover_percent 0.67523 1 0.4112
dogs_mercury_mod5 <- glm(Dogs_Mercury_percent ~ Ash_dieback + Tree_dens + Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod5)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 2.13112 2 0.3445
Tree_dens 0.08786 1 0.7669
Canopy_cover_percent 1.46391 1 0.2263
dogs_mercury_mod6 <- glm(Dogs_Mercury_percent ~ Ash_dieback+ Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod6)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Ash_dieback 2.2475 2 0.3251
Canopy_cover_percent 1.4426 1 0.2297
dogs_mercury_mod7 <- glm(Dogs_Mercury_percent ~ Canopy_cover_percent, data=Ash)
Anova(dogs_mercury_mod7)
Analysis of Deviance Table (Type II tests)
Response: Dogs_Mercury_percent
LR Chisq Df Pr(>Chisq)
Canopy_cover_percent 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
#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$total_invert <- Ash$Invert_sp_3006_0107 + Ash$Invert_sp_0107_0207
ggplotRegression(lm(Invert_sp_3006_0107 ~ Invert_sp_0107_0207, data = Ash))
#Univariable
ggplotRegression(lm(Temp_range ~ total_invert, data = Ash))
ggplotRegression(lm(Tree_dens ~ total_invert, data = Ash))
ggplotRegression(lm(Canopy_cover_percent ~ total_invert, data = Ash))
ggplotRegression(lm(BA ~ total_invert, data = Ash))
ggplotRegression(lm(Tree_dens ~ total_invert, data = Ash))
ggplotRegression(lm(Soil_moisture_ave ~ total_invert, data = Ash))
#no correlation at all!
optional: maybe try PCA to reduce the number of variables? (i am not sure about how to do this)
library("FactoMineR")
library("factoextra")
pca <- Ash[,c(8,9,13,15,23,25,28,29,31)]#excluded Ash dieback levels bc data needs to be numeric
View(pca)
pca.scaled <- scale(pca, center = TRUE, scale = TRUE)
res.cor <- cor(pca.scaled)
round(res.cor, 2)
Soil_moisture_ave Canopy_cover_percent Ash_percent Moss_percent Richness Tree_dens
Soil_moisture_ave 1.00 -0.28 0.24 0.08 0.47 0.34
Canopy_cover_percent -0.28 1.00 -0.45 -0.39 -0.28 0.07
Ash_percent 0.24 -0.45 1.00 0.02 0.25 -0.04
Moss_percent 0.08 -0.39 0.02 1.00 0.10 0.27
Richness 0.47 -0.28 0.25 0.10 1.00 0.00
Tree_dens 0.34 0.07 -0.04 0.27 0.00 1.00
Temp_range 0.25 -0.58 0.37 0.58 0.26 0.04
BA 0.01 0.28 -0.24 -0.31 -0.13 0.18
total_invert 0.16 0.01 -0.03 -0.33 -0.04 -0.06
Temp_range BA total_invert
Soil_moisture_ave 0.25 0.01 0.16
Canopy_cover_percent -0.58 0.28 0.01
Ash_percent 0.37 -0.24 -0.03
Moss_percent 0.58 -0.31 -0.33
Richness 0.26 -0.13 -0.04
Tree_dens 0.04 0.18 -0.06
Temp_range 1.00 -0.25 -0.12
BA -0.25 1.00 0.25
total_invert -0.12 0.25 1.00
#calculate the eigenvalues and eigenvectors of the matrix
res.eig <- eigen(res.cor)
res.eig
eigen() decomposition
$values
[1] 2.7729872 1.5495091 1.3591782 0.9200915 0.7367377 0.6488546 0.4088140 0.3651987 0.2386289
$vectors
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] -0.28714600 0.56465622 -0.15091034 -0.194131993 -0.04566296 0.22966331 0.38310298 0.53030135
[2,] 0.46705328 -0.01568447 -0.19190126 -0.321201120 0.03715849 0.08472012 -0.63586643 0.43719350
[3,] -0.34402411 0.13003576 0.36619532 0.003958477 0.73767807 -0.10540846 -0.27728858 -0.04730679
[4,] -0.37913741 -0.31121479 -0.42876286 0.198742018 -0.24694865 0.03298010 -0.14943999 -0.14236857
[5,] -0.31285613 0.32347461 0.11288795 -0.574933939 -0.40101238 -0.23572001 -0.31115402 -0.36903446
[6,] -0.07278756 0.25470004 -0.70420998 0.008226557 0.38024559 0.21687492 -0.12794225 -0.32127524
[7,] -0.48333121 -0.07906573 -0.04465779 0.335552105 -0.13297882 -0.23245168 -0.35589770 0.48385991
[8,] 0.28688784 0.38897260 -0.21674062 0.279448926 0.01003977 -0.77705019 0.06970127 -0.01113233
[9,] 0.13498136 0.48999428 0.25463145 0.546191287 -0.26115085 0.40919335 -0.32400192 -0.17744582
[,9]
[1,] -0.23543917
[2,] -0.19389806
[3,] -0.31018812
[4,] -0.65680420
[5,] 0.06938574
[6,] 0.35020287
[7,] 0.46157262
[8,] -0.18010409
[9,] -0.08012687
# Transpose eigeinvectors
eigenvectors.t <- t(res.eig$vectors)
# Transpose the adjusted data
pca.scaled.t <- t(pca.scaled)
# The new dataset
pca.new <- eigenvectors.t %*% pca.scaled.t
# Transpose new data ad rename columns
pca.new <- t(pca.new)
colnames(pca.new) <- c("PC1", "PC2", "PC3", "PC4","PC5","PC6","PC7","PC8","PC9")
head(pca.new)
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
[1,] 0.69963110 -2.5329545 0.61381320 1.1321918 -0.6462949 0.34883220 -0.03251557 -0.18859798
[2,] 0.16885521 0.3516004 -1.11053951 -0.4966349 -0.1151903 0.41070352 0.88368197 -0.06282776
[3,] 1.24550288 -2.4366722 -0.13005310 0.3505708 -0.2690557 -0.13428352 -0.27418130 0.47356471
[4,] 0.27517446 -1.5355109 -1.08561168 0.2275626 -0.0507482 0.60604883 0.30323121 0.79352037
[5,] 1.58429526 -0.3929382 -0.92958742 -0.3049923 -0.6361114 -0.95618024 -0.32256721 0.09974079
[6,] 0.02196995 -1.1508229 0.08607259 0.3915251 -1.2635745 0.09715244 -0.58031411 0.39205656
PC9
[1,] 0.253012504
[2,] 0.063804336
[3,] -0.023837012
[4,] 0.003078894
[5,] -0.988437921
[6,] -0.029799317
PCA(pca.new, scale.unit = TRUE, ncp = 5, graph = TRUE)
zero-length arrow is of indeterminate angle and so skipped
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 33 individuals, described by 9 variables
*The results are available in the following objects:
name description
1 "$eig" "eigenvalues"
2 "$var" "results for the variables"
3 "$var$coord" "coord. for the variables"
4 "$var$cor" "correlations variables - dimensions"
5 "$var$cos2" "cos2 for the variables"
6 "$var$contrib" "contributions of the variables"
7 "$ind" "results for the individuals"
8 "$ind$coord" "coord. for the individuals"
9 "$ind$cos2" "cos2 for the individuals"
10 "$ind$contrib" "contributions of the individuals"
11 "$call" "summary statistics"
12 "$call$centre" "mean of the variables"
13 "$call$ecart.type" "standard error of the variables"
14 "$call$row.w" "weights for the individuals"
15 "$call$col.w" "weights for the variables"