ES3307 Assignment 4
Misha Rashanah
2026-11-04
Question 1: Model selection - Information Theoretic (IT) Approach
In the IT Approach computer tutorial, you have learnt how to use the MuMIn package for model selection.
Using the same roadkill example, use the IT approach to determine the best suite of models that explains the occurrence of roadkills in subsistence agricultural land.
We are interested in these habitat characteristics:
Distance to nearest water body
Distance to nearest non-fruiting tree
Distance to nearest shrub
Distance to nearest human dwelling
Distance to nearest fruit tree
Our analyses are also performed on the assumption that the variables have no interaction with each other.
1A. Which model(s) have you included in the final candidate set?
# extract data from working directory and save it as a data frame.
sal<-read.csv("C:/Users/mishr/Desktop/Y4S1/ES3307/Assignment 4/Assignment 4 data/subsistence.csv")
head(sal)## name kill habitat water tree shrub human fruit cable transectID
## 1 control 0 subsistence 40 20 20 60 10 NA 0
## 2 control 0 subsistence 30 10 15 30 20 NA 0
## 3 control 0 subsistence 30 10 5 140 20 NA 0
## 4 control 0 subsistence 25 15 10 190 10 NA 0
## 5 control 0 subsistence 20 20 15 180 10 NA 0
## 6 control 0 subsistence 30 20 15 120 10 NA 0summary(sal)## name kill habitat water tree
## Length :87 Min. :0.0000 Length :87 Min. : 5.00 Min. : 5.00
## N.unique : 5 1st Qu.:0.0000 N.unique : 1 1st Qu.: 5.00 1st Qu.: 5.00
## N.blank : 0 Median :1.0000 N.blank : 0 Median :10.00 Median :10.00
## Min.nchar: 4 Mean :0.5402 Min.nchar:11 Mean :14.25 Mean :11.44
## Max.nchar: 9 3rd Qu.:1.0000 Max.nchar:11 3rd Qu.:20.00 3rd Qu.:15.00
## Max. :1.0000 Max. :50.00 Max. :30.00
## shrub human fruit cable
## Min. : 5.00 Min. : 10.0 Min. : 5.000 Mode:logical
## 1st Qu.: 5.00 1st Qu.: 50.0 1st Qu.: 5.000 NAs :87
## Median :10.00 Median :100.0 Median :10.000
## Mean :11.67 Mean :130.1 Mean : 9.828
## 3rd Qu.:15.00 3rd Qu.:157.5 3rd Qu.:10.000
## Max. :30.00 Max. :580.0 Max. :50.000
## transectID
## Min. :0.000
## 1st Qu.:0.000
## Median :3.000
## Mean :1.586
## 3rd Qu.:3.000
## Max. :3.000str(sal) # check data structure## 'data.frame': 87 obs. of 10 variables:
## $ name : chr "control" "control" "control" "control" ...
## $ kill : int 0 0 0 0 0 0 0 0 0 0 ...
## $ habitat : chr "subsistence" "subsistence" "subsistence" "subsistence" ...
## $ water : int 40 30 30 25 20 30 10 20 5 10 ...
## $ tree : int 20 10 10 15 20 20 25 10 5 10 ...
## $ shrub : int 20 15 5 10 15 15 20 25 30 15 ...
## $ human : int 60 30 140 190 180 120 120 20 310 50 ...
## $ fruit : int 10 20 20 10 10 10 10 5 20 30 ...
## $ cable : logi NA NA NA NA NA NA ...
## $ transectID: int 0 0 0 0 0 0 0 0 0 0 ...# several issues
# firstly, kill is binary either 0 (no kill) or 1 (kill)
# secondly, cable is NA for all so need to exclude from models
# not sure how big of a problem this is but there is only 2 values in transectID - 0 and 3
# missing? but stored as integer so should be okay?
# variables different from Swiss dataset example
################################## DECISIONS ####################################
# omit cable as it is not necessary and unlikely to give us relevant insights
# exclude transectID too as it is just a record metrics rather than an explanatory variable/predictor
# since our response variable KILL is binary we use glm with family=bionomial instead
# as it is not continuous and unlikely to be linear
# name also has no effect on model, its more of info? so we exclude
# write out possible models (referring from Tharka script 1)
library(car) # for Anova test
library(MuMIn) # for model averaging
# Set missing value handling
# options(na.action = "na.fail")
# MODEL
Q1 <- glm(
kill ~ water + tree + shrub + human + fruit,
family = binomial, na.action = "na.fail",
data = sal
)
mmQ1 <- dredge(Q1, rank = "AICc") # all possible models
mmQ1## Global model call: glm(formula = kill ~ water + tree + shrub + human + fruit, family = binomial,
## data = sal, na.action = "na.fail")
## ---
## Model selection table
## (Intrc) fruit human shrub tree water df logLik AICc delta
## 16 3.9640 -0.2340 0.007149 -0.09202 -0.1346 5 -40.864 92.5 0.00
## 32 4.2160 -0.2255 0.007210 -0.09276 -0.1361 -0.02289 6 -40.679 94.4 1.94
## 12 2.9620 -0.2309 0.006524 -0.1379 4 -43.052 94.6 2.12
## 28 3.1790 -0.2237 0.006506 -0.1379 -0.02035 5 -42.879 96.5 4.03
## 8 2.6430 -0.2444 0.006595 -0.09817 4 -45.053 98.6 6.13
## 14 4.0530 -0.1931 -0.07290 -0.1099 4 -45.656 99.8 7.33
## 24 2.7710 -0.2371 0.006649 -0.09732 -0.01528 5 -44.957 100.7 8.19
## 10 3.1990 -0.1914 -0.1123 3 -47.190 100.7 8.20
## 4 1.4880 -0.2299 0.005798 3 -47.624 101.5 9.07
## 30 4.2360 -0.1838 -0.07157 -0.1114 -0.01930 5 -45.504 101.7 9.28
## 26 3.4190 -0.1821 -0.1132 -0.02113 4 -46.981 102.5 9.98
## 20 1.6710 -0.2229 0.005847 -0.01834 4 -47.459 103.4 10.94
## 6 2.9090 -0.2040 -0.07782 3 -49.071 104.4 11.96
## 2 1.9410 -0.1965 2 -50.851 105.8 13.38
## 22 2.9980 -0.1970 -0.07652 -0.01206 4 -49.002 106.5 14.02
## 31 2.7880 0.005995 -0.09172 -0.1324 -0.05648 5 -48.006 106.8 14.28
## 18 2.0950 -0.1881 -0.01660 3 -50.706 107.7 15.23
## 15 1.9900 0.005832 -0.08785 -0.1351 4 -49.907 108.3 15.83
## 27 1.7490 0.005282 -0.1293 -0.05454 4 -50.502 109.5 17.02
## 11 1.0190 0.005143 -0.1328 3 -52.421 111.1 18.66
## 29 3.0870 -0.07137 -0.1148 -0.05406 4 -51.790 112.1 19.60
## 25 2.1870 -0.1116 -0.05259 3 -53.565 113.4 20.95
## 13 2.3100 -0.06960 -0.1160 3 -53.763 113.8 21.35
## 23 1.2870 0.005258 -0.08570 -0.05509 4 -53.108 114.7 22.24
## 9 1.4660 -0.1145 2 -55.541 115.2 22.76
## 7 0.4986 0.004955 -0.08192 3 -55.332 117.0 24.49
## 19 0.3807 0.004417 -0.05371 3 -55.506 117.3 24.83
## 21 1.6860 -0.06750 -0.05185 3 -56.103 118.5 26.03
## 3 -0.3582 0.004182 2 -57.765 119.7 27.20
## 17 0.8962 -0.05168 2 -57.775 119.7 27.23
## 5 0.9397 -0.06663 2 -58.280 120.7 28.23
## 1 0.1613 1 -60.022 122.1 29.62
## weight
## 16 0.500
## 32 0.189
## 12 0.173
## 28 0.067
## 8 0.023
## 14 0.013
## 24 0.008
## 10 0.008
## 4 0.005
## 30 0.005
## 26 0.003
## 20 0.002
## 6 0.001
## 2 0.001
## 22 0.000
## 31 0.000
## 18 0.000
## 15 0.000
## 27 0.000
## 11 0.000
## 29 0.000
## 25 0.000
## 13 0.000
## 23 0.000
## 9 0.000
## 7 0.000
## 19 0.000
## 21 0.000
## 3 0.000
## 17 0.000
## 5 0.000
## 1 0.000
## Models ranked by AICc(x)mmQ2_avg <- model.avg(mmQ1, subset = delta < 4)
summary(mmQ2_avg) # not very relevant for this question ##
## Call:
## model.avg(object = mmQ1, subset = delta < 4)
##
## Component model call:
## glm(formula = kill ~ <3 unique rhs>, family = binomial, data = sal,
## na.action = na.fail)
##
## Component models:
## df logLik AICc delta weight
## 1234 5 -40.86 92.47 0.00 0.58
## 12345 6 -40.68 94.41 1.94 0.22
## 124 4 -43.05 94.59 2.12 0.20
##
## Term codes:
## fruit human shrub tree water
## 1 2 3 4 5
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 3.818637 1.109965 1.124008 3.397 0.00068 ***
## fruit -0.231520 0.069638 0.070679 3.276 0.00105 **
## human 0.007037 0.002615 0.002654 2.652 0.00801 **
## shrub -0.073732 0.055804 0.056277 1.310 0.19014
## tree -0.135623 0.051748 0.052523 2.582 0.00982 **
## water -0.005031 0.020108 0.020345 0.247 0.80470
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 3.818637 1.109965 1.124008 3.397 0.00068 ***
## fruit -0.231520 0.069638 0.070679 3.276 0.00105 **
## human 0.007037 0.002615 0.002654 2.652 0.00801 **
## shrub -0.092222 0.046795 0.047499 1.942 0.05219 .
## tree -0.135623 0.051748 0.052523 2.582 0.00982 **
## water -0.022893 0.037829 0.038403 0.596 0.55109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# models included in the final candidate set
subset_mmQ1<-subset(mmQ1, delta < 4)
print(subset_mmQ1, abbrev.names = FALSE)## Global model call: glm(formula = kill ~ water + tree + shrub + human + fruit, family = binomial,
## data = sal, na.action = "na.fail")
## ---
## Model selection table
## (Intercept) fruit human shrub tree water df logLik AICc delta
## 16 3.964 -0.2340 0.007149 -0.09202 -0.1346 5 -40.864 92.5 0.00
## 32 4.216 -0.2255 0.007210 -0.09276 -0.1361 -0.02289 6 -40.679 94.4 1.94
## 12 2.962 -0.2309 0.006524 -0.1379 4 -43.052 94.6 2.12
## weight
## 16 0.58
## 32 0.22
## 12 0.20
## Models ranked by AICc(x)Based on the mmQ1 results where models with delta AICc < 4 are selected, the models included in the final candidate set are:
Model 16: glm(kill ~ fruit + human + shrub + tree, family = binomial, data = sal, na.action = “na.fail”)
Model 32: glm(kill ~ fruit + human + shrub + tree + water, family = binomial, data = sal, na.action = “na.fail”)
Model 12: glm(kill ~ fruit + human + tree, family = binomial, data = sal, na.action = “na.fail”)
The lower the AICc, the better the model is supported compared to the other models.
1B. Using the best model from your candidate set, obtain your model predicted lines and confidence intervals using ggpredict() and plot them with your original data.
(Hint: Think of the most suitable way to show how the probability of roadkills changes with each variable.)
The best model from the candidate set is Model 16 which has the lowest AICc value (92.5) compared to the other 2 models, Model 32 (94.4) and Model 12 (94.6). Additionally, Model 16 has an Akaike weight of 0.58 (58%). This means that it is the best supported model among the candidate models. However, Model 16 does not include all explanatory variables it excludes water. This means that the inclusion of distance to water body in Model 32 may not have provided sufficient significance to be a justifiable parameter since AICc value of Model 16 is the lowest. On the other hand, by excluding both water and shrub in Model 12, the AICc value is not lowered, suggesting that the inclusion of shrub may be more justifiable. Hence, among the 3 top options, Model 16 is the best option.
# extract model 16 from the mmQ1 set
best_model <- get.models(mmQ1, subset = 1)[[1]] # top 1 from the subset which is Model 16
summary(best_model)##
## Call:
## glm(formula = kill ~ fruit + human + shrub + tree + 1, family = binomial,
## data = sal, na.action = "na.fail")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.964419 1.030820 3.846 0.000120 ***
## fruit -0.234031 0.069852 -3.350 0.000807 ***
## human 0.007149 0.002613 2.736 0.006217 **
## shrub -0.092020 0.046696 -1.971 0.048771 *
## tree -0.134638 0.051944 -2.592 0.009542 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 120.044 on 86 degrees of freedom
## Residual deviance: 81.727 on 82 degrees of freedom
## AIC: 91.727
##
## Number of Fisher Scoring iterations: 5# + 1 is included in the model as it in the intercept by default there is a + 1 even if we do not explicitly write
# relevant libraries needed for plotting
library(ggeffects)
library(patchwork)
# model predicted line for the 4 explanatory variables
pred_fruit <- ggpredict(best_model, terms = "fruit")
pred_human <- ggpredict(best_model, terms = "human")
pred_shrub <- ggpredict(best_model, terms = "shrub")
pred_tree <- ggpredict(best_model, terms = "tree")
make_plot <- function(pred, variable, label, colour) {
ggplot() +
# Original observations
geom_jitter(
data = sal,
aes(x = .data[[variable]], y = kill),
height = 0.025,
width = 0,
alpha = 0.5,
colour = "black",
size = 2
) +
# 95% confidence interval
geom_ribbon(
data = pred,
aes(
x = x,
ymin = conf.low,
ymax = conf.high,
fill = "95% CI"
),
alpha = 0.25
) +
# Predicted probability line
geom_line(
data = pred,
aes(
x = x,
y = predicted,
colour = "Predicted probability"
),
linewidth = 1.2
) +
scale_colour_manual(
name = NULL,
values = c("Predicted probability" = colour)
) +
scale_fill_manual(
name = NULL,
values = c("95% CI" = colour)
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, 0.2)
) +
labs(
x = label,
y = "Predicted P(Roadkill)"
) +
guides(
colour = guide_legend(order = 1),
fill = guide_legend(order = 2)
) +
theme_classic() +
theme(
axis.title = element_text(size = 11),
axis.text = element_text(size = 10),
legend.position = "bottom",
legend.background = element_rect(
fill = "white",
colour = "black"
),
legend.box.background = element_rect(
colour = "black"
)
)
}
# 4 plots for the 4 explanatory variables
p1 <- make_plot(
pred_fruit, "fruit",
"Distance to nearest fruit tree",
"lightblue"
)
p2 <- make_plot(
pred_human, "human",
"Distance to nearest human dwelling",
"pink"
)
p3 <- make_plot(
pred_shrub, "shrub",
"Distance to nearest shrub",
"purple"
)
p4 <- make_plot(
pred_tree, "tree",
"Distance to nearest non-fruiting tree",
"lightgreen"
)
# combine all 4 plots into 1 figure for easier visualization
(p1 + p2) / (p3 + p4)
1C. Using the best model from your candidate set, report and interpret your results.
summary(best_model)##
## Call:
## glm(formula = kill ~ fruit + human + shrub + tree + 1, family = binomial,
## data = sal, na.action = "na.fail")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.964419 1.030820 3.846 0.000120 ***
## fruit -0.234031 0.069852 -3.350 0.000807 ***
## human 0.007149 0.002613 2.736 0.006217 **
## shrub -0.092020 0.046696 -1.971 0.048771 *
## tree -0.134638 0.051944 -2.592 0.009542 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 120.044 on 86 degrees of freedom
## Residual deviance: 81.727 on 82 degrees of freedom
## AIC: 91.727
##
## Number of Fisher Scoring iterations: 5coef_table <- as.data.frame(
summary(best_model)$coefficients
)
coef_table <- tibble::rownames_to_column(
coef_table,
var = "Variable"
)
# Reorder columns
coef_table <- coef_table[, c(
"Variable", "Estimate", "Std. Error",
"z value", "Pr(>|z|)"
)]
# present in table
knitr::kable(
coef_table,
digits = 4,
caption = "Results for roadkill Model 16",
col.names = c(
"Variable", "Estimate", "Standard Error",
"z-value", "p-value"
)
)| Variable | Estimate | Standard Error | z-value | p-value |
|---|---|---|---|---|
| (Intercept) | 3.9644 | 1.0308 | 3.8459 | 0.0001 |
| fruit | -0.2340 | 0.0699 | -3.3504 | 0.0008 |
| human | 0.0071 | 0.0026 | 2.7361 | 0.0062 |
| shrub | -0.0920 | 0.0467 | -1.9706 | 0.0488 |
| tree | -0.1346 | 0.0519 | -2.5920 | 0.0095 |
# report model fit statistics too
fit_table <- data.frame(
Statistic = c("Null deviance", "Residual deviance", "AIC"),
Value = c(120.044, 81.727, 91.727),
Degrees_of_freedom = c(86, 82, NA)
)
knitr::kable(
fit_table,
digits = 3,
caption = "Model fit statistics",
col.names = c("Statistic", "Value", "Degrees of freedom")
)| Statistic | Value | Degrees of freedom |
|---|---|---|
| Null deviance | 120.044 | 86 |
| Residual deviance | 81.727 | 82 |
| AIC | 91.727 | NA |
From the results, all four explanatory variables had significant associations with occurrence of roadkill. The habitat characteristic or parameter, distance to the nearest fruit tree had a significant negative relationship with roadkill occurrence (β = -0.23 (2 d.p.), p < 0.001), suggesting that roadkill was more likely to occur closer to fruit trees. Similarly, distance to the nearest shrub (β = -0.10 (2 d.p.), p = 0.049 < 0.05) and non-fruiting tree (β = -0.14, p = 0.0010 < 0.01) also has a significant negative association with roadkill occurrence. However, the variable distance to the nearest human dwelling had a significant positive relationship with roadkill occurrence (β = 0.01, p = 0.006 < 0.01), implying that roadkill was more likely to occur further away from human dwellings.
In addition, by including this four parameters, the model becomes a better fit as there is a reduction in deviance from 120.04 (null model, df = 86) to 81.73 (residual model, df = 82), with a deviance difference of 120.04 - 81.73 = 38.32 (2 d.p.) (χ^2 = 38.32, p < 0.001).
Overall, these model results suggest that roadkill occurrence in subsistence agricultural land is associated with the surrounding habitat characteristics. Roadkill occurrence was associated with shorter distances to fruit trees, shrubs and non-fruiting trees in particular, but longer distances from human dwellings.
Question 2: Interactions and Non-Linear data
2A. Plot the raw data
Plot the data to show the effect of DBH and SEX on the number of flowers in dioecious trees.
What does the plot tell you about the relationship between DBH and SEX?
# load assignment 1 data - dioecious tree
trees<-read.csv("C:/Users/mishr/Desktop/Y4S1/ES3307/Assignment 1/dioecious_tree.csv")
# check data structure
str(trees)## 'data.frame': 50 obs. of 3 variables:
## $ SEX : int 2 1 2 1 2 1 2 2 1 1 ...
## $ DBH : int 174 171 212 187 227 274 86 161 285 182 ...
## $ FLOWERS: int 234 110 512 133 1107 624 8 333 637 225 ...# DBH and SEX is integer
# since SEX is just 2 values (binary) FEMALE or MALE, need to change to factor
# 1 is MALE and 2 is FEMALE
trees$SEX <- factor(trees$SEX)
trees## SEX DBH FLOWERS
## 1 2 174 234
## 2 1 171 110
## 3 2 212 512
## 4 1 187 133
## 5 2 227 1107
## 6 1 274 624
## 7 2 86 8
## 8 2 161 333
## 9 1 285 637
## 10 1 182 225
## 11 2 225 512
## 12 1 133 47
## 13 2 154 164
## 14 2 71 1
## 15 1 133 36
## 16 1 267 524
## 17 2 166 198
## 18 2 145 225
## 19 2 207 548
## 20 1 216 272
## 21 2 157 164
## 22 2 87 5
## 23 1 185 133
## 24 1 183 133
## 25 2 193 419
## 26 1 264 442
## 27 1 83 3
## 28 2 300 1482
## 29 1 242 333
## 30 2 316 1138
## 31 2 217 561
## 32 1 235 488
## 33 2 189 282
## 34 1 242 536
## 35 2 84 11
## 36 2 135 148
## 37 1 267 512
## 38 2 122 70
## 39 2 280 1107
## 40 2 295 941
## 41 2 139 89
## 42 2 165 181
## 43 2 81 8
## 44 1 190 173
## 45 2 68 1
## 46 2 314 1091
## 47 1 284 561
## 48 2 77 3
## 49 1 228 312
## 50 2 282 1398summary(trees)## SEX DBH FLOWERS
## 1:20 Min. : 68.0 Min. : 1.0
## 2:30 1st Qu.:140.5 1st Qu.: 115.8
## Median :188.0 Median : 253.0
## Mean :191.6 Mean : 383.5
## 3rd Qu.:242.0 3rd Qu.: 533.0
## Max. :316.0 Max. :1482.0# From Assignment 1
# plot histogram #
ggplot(trees, aes(x = FLOWERS, fill = factor(SEX))) +
geom_histogram(
bins = 5,
color = "black"
) +
facet_wrap(~ factor(SEX)) +
labs(
x = "Sex",
y = "Number of Flowers",
title = "Distribution of Flowers by Sex"
) +
scale_fill_manual(
values = c("skyblue", "lightpink")
) +
theme_classic() +
theme(legend.position = "TRUE")
# plot data using boxplots #
ggplot(trees,
aes(x = factor(SEX), y = FLOWERS, fill = factor(SEX))) +
geom_boxplot() +
labs(
x = "Sex",
y = "Number of Flowers",
title = "Distribution of Flowers by Sex"
) +
scale_fill_manual(
values = c("skyblue", "lightpink")
) +
theme_classic() +
theme(legend.position = "TRUE")
# Plot raw data using scatter plot #
ggplot(trees,
aes(x = DBH, y = FLOWERS, colour = factor(SEX))) +
geom_point(size = 3) +
labs(
x = "Diameter at Breast Height (mm)",
y = "Number of Flowers",
title = "Plotting the effect of DBH and SEX on the number of flowers in dioecious trees",
colour = "Sex"
) +
scale_colour_manual(
values = c("skyblue", "lightpink"),
labels = c("Male", "Female")
) +
theme_classic()
From the scatter plot diagram, we can clearly observe that the range of number of flower for female dioecious trees is greater than that of the male and also more spread out in terms of DBH, ranging from less than 0 all the way to more than 300m. Whereas for male dioecious trees, the DBH range is slightly clustered starting from 130mm and ending before 300m. This suggests that there is an association or relationship present between SEX and DBH of diecious trees, indicating a potential interaction.
2B. Fit linear model
- Based on the plot in Question 2A, fit an appropriate linear model to determine the effect of DBH and SEX on the number of flowers.
- Run diagnostic plots to check your model fit.
- Report and interpret your results.
# with interaction
model_trees<-lm(FLOWERS ~ SEX*DBH, data=trees)
summary(model_trees)##
## Call:
## lm(formula = FLOWERS ~ SEX * DBH, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -212.09 -88.09 -46.77 91.64 403.24
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -440.7002 130.2774 -3.383 0.00147 **
## SEX2 -108.0718 146.7008 -0.737 0.46506
## DBH 3.5399 0.5935 5.964 3.28e-07 ***
## SEX2:DBH 1.9779 0.6885 2.873 0.00614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 145.5 on 46 degrees of freedom
## Multiple R-squared: 0.8646, Adjusted R-squared: 0.8557
## F-statistic: 97.88 on 3 and 46 DF, p-value: < 2.2e-16# without interaction, additive effect
model_trees_1 <-lm(FLOWERS ~ SEX + DBH, data = trees)
summary(model_trees_1)##
## Call:
## lm(formula = FLOWERS ~ SEX + DBH, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -206.31 -87.68 -34.16 63.06 443.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -753.0820 77.0812 -9.770 6.75e-13 ***
## SEX2 294.5838 46.5125 6.333 8.41e-08 ***
## DBH 5.0096 0.3232 15.499 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 156.3 on 47 degrees of freedom
## Multiple R-squared: 0.8403, Adjusted R-squared: 0.8335
## F-statistic: 123.6 on 2 and 47 DF, p-value: < 2.2e-16anova(model_trees, model_trees_1)## Analysis of Variance Table
##
## Model 1: FLOWERS ~ SEX * DBH
## Model 2: FLOWERS ~ SEX + DBH
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 46 973657
## 2 47 1148328 -1 -174671 8.2523 0.006138 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# from anova results
# run diagnostic plots
par(mfrow=c(2, 2)) # to set up a side by side layout
fit<-lm(FLOWERS ~ SEX*DBH, data=trees)# Run a regression analysis
plot(fit)
Based on the linear model FLOWERS ~ SEX*DBH, both tree sex and DBH have a significant interaction effect on the number of flowers (F1,46 = 8.25, p = 0.006). The interaction model is preferred over the additive model (FLOWERS ~ SEX + DBH), as it provides a significantly better fit and also statistically significant (F3,46 = 97.88, p < 0.001). This implies that there is strong evidence to support that the effect of DBH on number of flowers is different between male and female trees.
Interpretation of Diagnostic Plots
Residuals vs Fitted: The residuals show a slight curved pattern and not evenly distributed, suggesting possible non-linearity.
Q-Q Residuals: The points deviate from the qqnorm line at the tails, suggesting that the residuals may not be normally distributed.
Scale-Location: The residuals are not equally spread along the range of predictors, suggesting some heteroscedasticity.
Residuals vs Leverage: There is no significant outliers.
2C. Transform your data
Based on the diagnostic plots in Question 1B, transform your y variable using an appropriate transformation and fit the linear model again.
Explain why you have chosen this transformation.
Plot the model predicted lines and confidence intervals with your original data. (Hint: Plot your predicted lines and confidence interval on the same scale as your original data.)
Run diagnostic plots to check your model fit.
Report and interpret your results.
# log transform number of flowers variable
model_log <- lm(log(FLOWERS) ~ SEX * DBH,
data = trees)
summary(model_log)##
## Call:
## lm(formula = log(FLOWERS) ~ SEX * DBH, data = trees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1843 -0.3926 0.1704 0.6193 1.3565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.704945 0.752965 0.936 0.354
## SEX2 -0.319963 0.847887 -0.377 0.708
## DBH 0.021599 0.003430 6.297 1.04e-07 ***
## SEX2:DBH 0.003743 0.003979 0.941 0.352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8409 on 46 degrees of freedom
## Multiple R-squared: 0.8133, Adjusted R-squared: 0.8011
## F-statistic: 66.78 on 3 and 46 DF, p-value: < 2.2e-16# Compare models
model_log_add <- lm(log(FLOWERS) ~ SEX + DBH,
data = trees)
anova(model_log_add, model_log)## Analysis of Variance Table
##
## Model 1: log(FLOWERS) ~ SEX + DBH
## Model 2: log(FLOWERS) ~ SEX * DBH
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 33.151
## 2 46 32.525 1 0.62552 0.8847 0.3518library(ggplot2)
newdata <- expand.grid(
DBH = seq(min(trees$DBH),
max(trees$DBH),
length.out = 100),
SEX = levels(trees$SEX)
)
pred <- predict(model_log_add,
newdata = newdata,
interval = "confidence")
newdata$fit <- exp(pred[, "fit"])
newdata$lwr <- exp(pred[, "lwr"])
newdata$upr <- exp(pred[, "upr"])
# plot
ggplot(trees, aes(x = DBH, y = FLOWERS, colour = SEX)) +
geom_point(alpha = 0.7) +
geom_ribbon(
data = newdata,
aes(x = DBH, ymin = lwr, ymax = upr,
fill = SEX, group = SEX),
alpha = 0.2,
inherit.aes = FALSE
) +
geom_line(
data = newdata,
aes(x = DBH, y = fit,
colour = SEX, group = SEX),
linewidth = 1
) +
scale_colour_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
scale_fill_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
labs(
title = "Predicted Number of Flowers by DBH and Sex",
x = "DBH (mm)",
y = "Number of Flowers"
) +
theme_minimal()
# run diagnostic plot
par(mfrow=c(2, 2)) # to set up a side by side layout
plot(model_log_add)
From the diagnostic plots in 2B, we are going to log transform the number of flowers (y variable) to make it more normally distributed and overcome the heavily right skewed issue (see summary(trees)).
As such a log transformation, log(FLOWERS), was applied.
Based on the linear model log(FLOWERS) ~ SEX × DBH, the overall model is statistically significant (F3,46 = 66.78, p < 0.001). DBH has a significant positive effect on number of flowers (p < 0.001), while the interaction between SEX and DBH is not significant (F1,46= 0.885, p = 0.35 > 0.05). Therefore, the additive model (SEX + DBH) is preferred, as including the interaction term does not significantly improve the model fit. This suggests that flower production increases with DBH, but the relationship does not differ significantly between male and female trees.
Furthermore, the plot shows that the predicted number of flowers increases with DBH for both sexes, with females generally having higher predicted flower counts than males, although the confidence intervals widen at larger DBH values, indicating greater uncertainty.
Interpretation of Diagnostic Plots
Residuals vs Fitted: The residuals show a curved pattern suggesting likely non-linearity.
Normal Q-Q: Most points follow the diagonal line, although some deviations are observed at the lower tail. This suggests that the residuals are likely normally distributed.
Scale-Location: The residuals are not evenly spread across fitted values. This indicates possible heteroscedasticity.
Residuals vs Leverage: A few potentially influential observations, such as 14, 30 and 45. However, none appear to exceed the Cook’s distance line.
2D. Fit a generalised linear model (GLM) - Poisson
As the number of flowers is count data, we can instead fit a GLM to the data. Now, fit a GLM using a Poisson distribution to determine the effect of DBH and SEX on the number of flowers.
Plot the model predicted lines and confidence intervals with your original data. (Hint: Plot your predicted lines and confidence interval on the same scale as your original data.)
Run diagnostic plots to check your model fit.
Report and interpret your results.
What do you notice about the dispersion parameter in the summary() output?
glm_trees <- glm(FLOWERS ~ DBH*SEX, family = poisson, data = trees)
summary(glm_trees)##
## Call:
## glm(formula = FLOWERS ~ DBH * SEX, family = poisson, data = trees)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2567217 0.0798972 28.245 < 2e-16 ***
## DBH 0.0150589 0.0003198 47.087 < 2e-16 ***
## SEX2 0.9453136 0.0869851 10.868 < 2e-16 ***
## DBH:SEX2 -0.0018137 0.0003460 -5.242 1.59e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17921.6 on 49 degrees of freedom
## Residual deviance: 2458.7 on 46 degrees of freedom
## AIC: 2811.6
##
## Number of Fisher Scoring iterations: 5library(ggplot2)
# Create prediction data
newdata <- expand.grid(
DBH = seq(min(trees$DBH),
max(trees$DBH),
length.out = 100),
SEX = levels(trees$SEX)
)
# Predict on the link (log) scale
pred <- predict(
glm_trees,
newdata = newdata,
type = "link",
se.fit = TRUE
)
# Calculate confidence intervals on log scale
newdata$fit <- exp(pred$fit)
newdata$lwr <- exp(pred$fit - 1.96 * pred$se.fit)
newdata$upr <- exp(pred$fit + 1.96 * pred$se.fit)
# Plot
ggplot(trees, aes(x = DBH, y = FLOWERS, colour = SEX)) +
geom_point(alpha = 0.7) +
geom_ribbon(
data = newdata,
aes(x = DBH, ymin = lwr, ymax = upr,
fill = SEX, group = SEX),
alpha = 0.2,
inherit.aes = FALSE
) +
geom_line(
data = newdata,
aes(x = DBH, y = fit,
colour = SEX, group = SEX),
linewidth = 1
) +
scale_colour_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
scale_fill_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
labs(
title = "Poisson GLM: Predicted Number of Flowers",
x = "DBH (mm)",
y = "Number of Flowers"
) +
theme_minimal()
# run diagnostic plot (use DHARMa instead cos mixed model)
library(DHARMa)
# Simulate residuals
sim_res <- simulateResiduals(
fittedModel = glm_trees,
n = 1000
)
# Combine diagnostic plots into one figure
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(sim_res, quantreg = TRUE)
Based on the Poisson GLM (FLOWERS ~ DBH × SEX), DBH, SEX and their interaction all have significant effects on the number of flowers. DBH has a significant positive effect (p < 0.001), implying that the number of flowers increases with tree size. SEX also has a significant positive effect (p < 0.001), while the negative interaction suggests that the increase in the number of flowers with DBH is weaker for females. The plot shows that predicted flower values increase with DBH for both sexes, with females generally having higher predicted values than males.
However, the DHARMa diagnostic plots indicate significant deviations from the expected residual distribution and significant overdispersion (p < 0.001), suggesting that the Poisson model does not adequately fit the data. The dispersion parameter in the summary is fixed at 1, meaning the Poisson model assumes equal mean and variance, but the diagnostic results suggest this assumption is violated.
2E. Fit a generalised linear model (GLM) - Quasipoisson
Now, fit a GLM with a Quasipoisson distribution. What has changed in your Quasipoisson model results?
Plot the model predicted lines and confidence intervals with your original data. (Hint: Plot your predicted lines and confidence interval on the same scale as your original data.)
Run diagnostic plots to check your model fit.
Report and interpret your results.
glm_quasi_trees <- glm(FLOWERS ~ DBH*SEX, family = quasipoisson, data = trees)
summary(glm_quasi_trees)##
## Call:
## glm(formula = FLOWERS ~ DBH * SEX, family = quasipoisson, data = trees)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.256722 0.569334 3.964 0.000255 ***
## DBH 0.015059 0.002279 6.608 3.54e-08 ***
## SEX2 0.945314 0.619840 1.525 0.134083
## DBH:SEX2 -0.001814 0.002466 -0.736 0.465714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 50.77736)
##
## Null deviance: 17921.6 on 49 degrees of freedom
## Residual deviance: 2458.7 on 46 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5library(ggplot2)
# Create prediction data
newdata <- expand.grid(
DBH = seq(min(trees$DBH),
max(trees$DBH),
length.out = 100),
SEX = levels(trees$SEX)
)
# Predictions on log scale
pred <- predict(
glm_quasi_trees,
newdata = newdata,
type = "link",
se.fit = TRUE
)
# Back-transform to original scale
newdata$fit <- exp(pred$fit)
newdata$lwr <- exp(
pred$fit - 1.96 * pred$se.fit
)
newdata$upr <- exp(
pred$fit + 1.96 * pred$se.fit
)
# Plot original data and predicted lines
ggplot(trees, aes(x = DBH, y = FLOWERS, colour = SEX)) +
geom_point(alpha = 0.7) +
geom_ribbon(
data = newdata,
aes(x = DBH, ymin = lwr, ymax = upr,
fill = SEX, group = SEX),
alpha = 0.2,
inherit.aes = FALSE
) +
geom_line(
data = newdata,
aes(x = DBH, y = fit,
colour = SEX, group = SEX),
linewidth = 1
) +
scale_colour_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
scale_fill_manual(
name = "Sex",
values = c("1" = "skyblue", "2" = "pink"),
labels = c("Male", "Female")
) +
labs(
title = "Quasi-Poisson GLM: Predicted Number of Flowers",
x = "DBH (mm)",
y = "Number of Flowers"
) +
theme_minimal()
# run diagnostic plots
# Simulate residuals
par(mfrow = c(2, 2))
plot(glm_quasi_trees)
Based on the Quasi-Poisson GLM (FLOWERS ~ DBH × SEX), DBH has a significant positive effect on the number of flowers (p < 0.001), while SEX (p = 0.134 > 0.05) and the interaction between SEX and DBH (p = 0.466) are not significant. This suggests that the number of flower increases with DBH, but the relationship does not differ significantly between male and female trees. Compared to the Poisson model, the Quasi-Poisson model produces larger standard errors and p-values due to the adjustment for overdispersion. The dispersion parameter is estimated at 50.78 (2 d.p.), indicating substantial overdispersion, implying that the variance is much greater than the mean. Hence, the Quasi-Poisson model is more appropriate than the Poisson model for accounting for the overdispersion in the data.
Interpretation of Diagnostic Plots
Residuals vs Fitted: The residuals show a curved pattern suggesting likely non-linearity.
Normal Q-Q: Most points follow the diagonal line, although some deviations are observed at the upper tail. This suggests that the residuals are likely normally distributed.
Scale-Location: The residuals are not evenly spread across fitted values. This indicates possible heteroscedasticity.
Residuals vs Leverage: No signficiant outliers
2F.
- From the above models, which model would you select as the final model? Explain your reasoning.
I will select the Quasi-Poisson Model as the final model because the response variable is count data (number of flowers), and the data may exhibit overdispersion. This model is more appropriate than the Poisson model when the variance of the response variable is greater than its mean. As for the other model options, there are also not the most suitable because from the diagnostic plots, they violate the linear model rules such as normality and homoscedasticity.
What are the advantages and disadvantages associated with running a linear model or a generalised linear model?
Comparing GLM and LM GLM LM Advantages Works well with non-linear and non normal data
Variety of family structure including poisson and binomial
Apt for count data
Efficient method, faster to calculate
Each regression coefficient shows exactly how much the response variable changes when a predictor (explanatory variable) changes by one unit
Disadvantages GLMs assume each record is completely independent Must be normally distributed and have homoscedasticity and be linear
Question 3: Multivariate Analysis In this study, the authors investigated the difference in dung beetle community composition between two commonly used trapping methods: (1) dung-baited pitfall traps and (2) flight interception traps.
For more information about this study, you may check out: https://onlinelibrary.wiley.com/doi/abs/10.1111/aec.13124
3A. Based on the analysis you have selected, plot and interpret your ordination plot.
Finding the most suitable k value.
# load dataset
beetles <- read.csv("C:/Users/mishr/Desktop/Y4S1/ES3307/Assignment 4/Assignment 4 data/dungbeetletraptype.csv")
# str(beetles)
# 2 different traps (treatments)
unique(beetles$habitat)## [1] "Old-growth Forest" "Logged Forest" "Riparian Reserve"# 3 unique habitat types: old-growth forest, logged forest and riparian reserve
unique(beetles$year)## [1] 2005 2012# 2 years when experiment was conducted: 2005 and 2012
# 42 dung beetles species
table(beetles$habitat, beetles$year)##
## 2005 2012
## Logged Forest 80 36
## Old-growth Forest 40 0
## Riparian Reserve 0 148# Old-growth Forest was only sampled in 2005.
# Riparian Reserve was only sampled in 2012.
# Logged Forest was sampled in both years.
# USE NMDS to understand ecological community composition of these dung beetles
library(permute)
library(vegan)
# Species abundance data (42 species)
beetle_species <- beetles[, 5:46]
beetle_metadata <- beetles[, 1:4]
# calculate Bray-Curtis dissimilarity
# Bray-Curtis is used for abudance
beetle_brayC <- vegdist(
beetle_species,
method = "bray"
)
# dissimilarity matrix
as.matrix(beetle_brayC)[1:5, 1:5]## 1 2 3 4 5
## 1 0.0000000 0.3424658 0.3333333 0.3442623 0.4062500
## 2 0.3424658 0.0000000 0.2873563 0.3684211 0.2658228
## 3 0.3333333 0.2873563 0.0000000 0.3333333 0.2820513
## 4 0.3442623 0.3684211 0.3333333 0.0000000 0.3134328
## 5 0.4062500 0.2658228 0.2820513 0.3134328 0.0000000# for reproducability
set.seed(42)
nmds_results <- list()
stress_values <- numeric(10)
for (k in 1:10) {
nmds_results[[k]] <- metaMDS(
beetle_species,
distance = "bray",
k = k,
trymax = 50,
autotransform = FALSE,
trace = FALSE
)
stress_values[k] <- nmds_results[[k]]$stress
}
# VISUALIZE STRESS VALUE
stress_table <- data.frame(
Dimensions = 1:10, # k: 1 to 10
Stress = stress_values
)
stress_table## Dimensions Stress
## 1 1 0.42841081
## 2 2 0.20935230
## 3 3 0.15083214
## 4 4 0.12028320
## 5 5 0.10197783
## 6 6 0.08824873
## 7 7 0.07813874
## 8 8 0.06974174
## 9 9 0.06378083
## 10 10 0.05823967ggplot(stress_table,
aes(x = Dimensions, y = Stress)) +
geom_line(linewidth = 0.8) +
geom_point(size = 3) +
scale_x_continuous(breaks = 1:10) +
theme_classic() +
labs(
title = "Visualizing the Stress Value against the Number of NMDS Dimensions",
x = "Number of dimensions (k)",
y = "Stress"
)
# from the plot, we can observe that the stress value starts to decrease logarithmatically from k=3 onward
# Select the final NMDS solution
nmds_final <- nmds_results[[3]]
nmds_final##
## Call:
## metaMDS(comm = beetle_species, distance = "bray", k = k, trymax = 50, autotransform = FALSE, trace = FALSE)
##
## global Multidimensional Scaling using monoMDS
##
## Data: beetle_species
## Distance: bray
##
## Dimensions: 3
## Stress: 0.1508321
## Stress type 1, weak ties
## Best solution was not repeated after 50 tries
## The best solution was from try 46 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'beetle_species'nmds_final$stress## [1] 0.1508321# Shepard plot #
shepard_results <- stressplot(nmds_final)
# Calculate statistics
k_value <- nmds_final$ndim
stress_value <- nmds_final$stress
nonmetric_R2 <- 1 - stress_value^2
linear_R2 <- cor(
shepard_results$y,
shepard_results$yf
)^2
# Create results table
nmds_summary <- data.frame(
k = k_value,
Stress = stress_value,
Nonmetric_R2 = nonmetric_R2,
Linear_fit_R2 = linear_R2
)
library(knitr)
# NMDS Model Fit Statistics
kable(
nmds_summary,
digits = 3,
caption = "Summary of NMDS Model Fit Statistics"
)| k | Stress | Nonmetric_R2 | Linear_fit_R2 |
|---|---|---|---|
| 3 | 0.151 | 0.977 | 0.842 |
Draw out the NMDS Plot.
library(ggplot2)
library(vegan)
library(grid)
# Extract NMDS site scores
nmds_scores <- as.data.frame(
scores(nmds_final, display = "sites")
)
# Add sampling metadata
nmds_scores$habitat <- as.factor(beetle_metadata$habitat)
nmds_scores$year <- as.factor(beetle_metadata$year)
nmds_scores$trap <- as.factor(beetle_metadata$trap)
# Check trapping methods
table(nmds_scores$trap)##
## BPT FIT
## 166 138# Check sample numbers
table(nmds_scores$habitat,
nmds_scores$trap,
nmds_scores$year)## , , = 2005
##
##
## BPT FIT
## Logged Forest 20 60
## Old-growth Forest 10 30
## Riparian Reserve 0 0
##
## , , = 2012
##
##
## BPT FIT
## Logged Forest 24 12
## Old-growth Forest 0 0
## Riparian Reserve 112 36set.seed(42)
# NMDS ordination plot
plot_nmds_clean <- ggplot(
nmds_scores,
aes(
x = NMDS1,
y = NMDS2,
colour = trap,
shape = habitat,
fill = trap
)
) +
# 95% confidence ellipses
stat_ellipse(
# aes(group = interaction(habitat, trap)),
aes(group = trap),
geom = "polygon",
level = 0.95,
alpha = 0.12,
linewidth = 0.7
) +
# Individual sampling points
geom_point(
size = 2.8,
alpha = 0.8,
stroke = 0.8
) +
# Separate sampling years
facet_wrap(~year) +
# Trapping method colours
scale_colour_manual(
values = c(
"BPT" = "blue",
"FIT" = "red"
)
) +
scale_fill_manual(
values = c(
"BPT" = "blue",
"FIT" = "red"
)
) +
# Habitat shapes
scale_shape_manual(
values = c(
"Old-growth Forest" = 21,
"Logged Forest" = 22,
"Riparian Reserve" = 24
)
) +
theme_classic() +
labs(
title = "NMDS Ordination of Dung Beetle Communities",
x = "NMDS1",
y = "NMDS2",
colour = "Trapping Method",
fill = "Trapping Method",
shape = "Habitat"
) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 10),
legend.position = "right",
strip.background = element_rect(fill = "grey90"),
strip.text = element_text(face = "bold"),
plot.caption = element_text(size = 8)
)
plot_nmds_clean
The NMDS ordination was run in 3 dimensions with a stress value of 0.151 (3 d.p.), indicating a reasonable representation of the community dissimilarities although stress is slight above the 0.1 threshold for a good fit. k=3 was preferred over k=2 because the stress value at k=2 is at 0.209 (3 d.p.) which indicates a poor fit but when k increased to 3, the stress value dropped significantly to 0.151, reducing the stress to a more acceptable level, making k=3 a more apt and representative choice. From the NMDS plot (visualized in 2 dimensions rather than 3 for easier understanding) which is separated by the years, we can observe that there are differences in dung beetle community between 2005 and 2012, with respect to the trapping method.
In 2005, BPT and FIT samples clearly exhibit distinct clustering, with no overlap between their respective 95% confidence ellipses, implying that there are significant differences in the community composition captured by the two methods. Additionally, the variability in habitat type overlap covers a larger range or more spread out for FIT than BPT. However, for 2012, greater overlap was observed, indicating less noticeable differences between the two trapping methods. In addition, the distribution of samples in 2012 was wider than the ones in 2005, also suggesting greater variability in community composition. Habitat-related differences were also less distinct, with considerable overlap among the three habitat types.
3B. Using appropriate statistical tests, report and interpret your results.
# PERMANOVA
set.seed(42)
permanova_results <- adonis2(
beetle_brayC ~ habitat + year + trap,
data = beetle_metadata,
permutations = 999,
method = "bray",
by = "margin"
)
permanova_results## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = beetle_brayC ~ habitat + year + trap, data = beetle_metadata, permutations = 999, method = "bray", by = "margin")
## Df SumOfSqs R2 F Pr(>F)
## habitat 2 2.824 0.03096 6.4019 0.001 ***
## year 1 2.666 0.02924 12.0896 0.001 ***
## trap 1 9.174 0.10059 41.5933 0.001 ***
## Residual 299 65.947 0.72309
## Total 303 91.201 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(knitr)
permanova_table <- as.data.frame(permanova_results)
permanova_table$Factor <- rownames(permanova_table)
permanova_table <- permanova_table[, c(
"Df",
"SumOfSqs",
"R2",
"F",
"Pr(>F)"
)]
kable(
permanova_table,
digits = 3,
caption = "PERMANOVA Results for Dung Beetle Community Composition"
)| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| habitat | 2 | 2.824 | 0.031 | 6.402 | 0.001 |
| year | 1 | 2.666 | 0.029 | 12.090 | 0.001 |
| trap | 1 | 9.174 | 0.101 | 41.593 | 0.001 |
| Residual | 299 | 65.947 | 0.723 | NA | NA |
| Total | 303 | 91.201 | 1.000 | NA | NA |
# PERMDISP
# Habitat
disp_habitat <- betadisper(
beetle_brayC,
group = beetle_metadata$habitat
)
test_habitat <- permutest(
disp_habitat,
permutations = 999
)
# Year
disp_year <- betadisper(
beetle_brayC,
group = beetle_metadata$year
)
test_year <- permutest(
disp_year,
permutations = 999
)
# Trapping method
disp_trap <- betadisper(
beetle_brayC,
group = beetle_metadata$trap
)
test_trap <- permutest(
disp_trap,
permutations = 999
)
# Extract results
extract_disp <- function(test, factor_name) {
result <- as.data.frame(test$tab)
result$Factor <- factor_name
result
}
dispersion_table <- rbind(
extract_disp(test_habitat, "Habitat"),
extract_disp(test_year, "Year"),
extract_disp(test_trap, "Trapping method")
)
kable(
dispersion_table,
digits = 3,
caption = "PERMDISP Results for Habitat, Year and Trapping Method"
)| Df | Sum Sq | Mean Sq | F | N.Perm | Pr(>F) | Factor | |
|---|---|---|---|---|---|---|---|
| Groups | 2 | 0.008 | 0.004 | 0.431 | 999 | 0.667 | Habitat |
| Residuals | 301 | 2.932 | 0.010 | NA | NA | NA | Habitat |
| Groups1 | 1 | 0.034 | 0.034 | 3.537 | 999 | 0.064 | Year |
| Residuals1 | 302 | 2.944 | 0.010 | NA | NA | NA | Year |
| Groups2 | 1 | 0.047 | 0.047 | 4.615 | 999 | 0.032 | Trapping method |
| Residuals2 | 302 | 3.048 | 0.010 | NA | NA | NA | Trapping method |
From the PERMANOVA results, there are significant differences in dung beetle community composition associated with habitat (F2,299 = 6.40, R^2= 0.031, p = 0.001), year (F1,299 = 12.10, R^2 = 0.029, p = 0.001) and trapping method (F1,299 = 41.60, R^2 = 0.101, p = 0.001). Trapping method accounted for the largest proportion of variation in community composition (0.10059 x 100 = 10.06%), followed by habitat (3.10%) and year (2.92%) (2 d.p.).
On the other hand, PERMDISP results reveal that there are no significant differences in multivariate dispersion among habitats (F2,301 = 0.43, p = 0.667) or between years (F1,301 = 3.54, p = 0.064). However, dispersion differed significantly between trapping methods (F1,301 = 4.62, p = 0.032). This means that the significant PERMANOVA result for trapping method may partially reflect differences in within-group variability. Therefore, differences associated with trapping method should be interpreted with caution.
Overall, as this a subset of the actual dataset, and with only 2 years being taken into account for and also there is unequal representation across habitats and trapping methods. For example, FIT samples from old-growth forest were only available in 2005, potentially confounding temporal and trapping-method differences with habitat-related variation.
