Assignment 4_Heng Yan Ying
Heng Yan Ying
2026-11-03
Load packages
library(ggplot2)
library(tidyr)
library(dplyr) # mutate function
library(ggeffects) # ggpredict
library(ggpp) # offset scatter plot on box plot
library(car) # type 2 ANOVA
library(emmeans) # stats
library(lme4) # for lmer
library(lmerTest)
library(readxl) # read excel
library(MuMIn) # dredge()
library(patchwork) # combine graphsQuestion 1: Model selection - Information Theoretic (IT) Approach
Load and review data
oilpalm <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 4/Assignment 4 data/subsistence.csv")
# Check data structure and summary
str(oilpalm)## '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 ...summary(oilpalm)## name kill habitat water
## Length:87 Min. :0.0000 Length:87 Min. : 5.00
## Class :character 1st Qu.:0.0000 Class :character 1st Qu.: 5.00
## Mode :character Median :1.0000 Mode :character Median :10.00
## Mean :0.5402 Mean :14.25
## 3rd Qu.:1.0000 3rd Qu.:20.00
## Max. :1.0000 Max. :50.00
## tree shrub human fruit
## Min. : 5.00 Min. : 5.00 Min. : 10.0 Min. : 5.000
## 1st Qu.: 5.00 1st Qu.: 5.00 1st Qu.: 50.0 1st Qu.: 5.000
## Median :10.00 Median :10.00 Median :100.0 Median :10.000
## Mean :11.44 Mean :11.67 Mean :130.1 Mean : 9.828
## 3rd Qu.:15.00 3rd Qu.:15.00 3rd Qu.:157.5 3rd Qu.:10.000
## Max. :30.00 Max. :30.00 Max. :580.0 Max. :50.000
## cable transectID
## Mode:logical Min. :0.000
## NA's:87 1st Qu.:0.000
## Median :3.000
## Mean :1.586
## 3rd Qu.:3.000
## Max. :3.000Q1A
Which model(s) have you included in the final candidate set?
# Check for multicollinearity
# Fit full model with all terms
full <- glm(kill ~ water + tree + shrub + human + fruit,
family = binomial, data = oilpalm)
# Check VIF values
vif(full)## water tree shrub human fruit
## 1.053694 1.039607 1.044163 1.179552 1.165295As all our variables have a VIF value of less than 3, we can proceed to determine our candidate set of models.
# Fit final model
final <- glm(kill ~ water + tree + shrub + human + fruit,
family = binomial, data = oilpalm,
na.action = "na.fail") # include this to prevent errors when running dredge
# Dredge final model
candidates <- dredge(final)
print(candidates)## Global model call: glm(formula = kill ~ water + tree + shrub + human + fruit, family = binomial,
## data = oilpalm, 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)# Retain retain models where the Akaike’s difference with the top-ranked model was less than 5.0
print(subset(candidates, delta <5))## Global model call: glm(formula = kill ~ water + tree + shrub + human + fruit, family = binomial,
## data = oilpalm, na.action = "na.fail")
## ---
## Model selection table
## (Intrc) 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
## 28 3.179 -0.2237 0.006506 -0.1379 -0.02035 5 -42.879 96.5 4.03
## weight
## 16 0.538
## 32 0.204
## 12 0.186
## 28 0.072
## Models ranked by AICc(x)The models included for the final candidate set are: kill ~ fruit + human + shrub + tree kill ~ fruit + human + shrub + tree + water kill ~ fruit + human + tree kill ~ fruit + human + tree + water
Q1B
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.)
# Fit the best model from candidate set
best <- glm(kill ~ fruit + human + shrub + tree,
family = binomial, data = oilpalm)
summary(best)##
## Call:
## glm(formula = kill ~ fruit + human + shrub + tree, family = binomial,
## data = oilpalm)
##
## 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# Diagnostic plots, AI has been used to generate this bit of code
library(arm)
binnedplot(fitted(best), residuals(best, type = "response"))
Fig 1Ba. Diagnostic plot of best model from candidate set.
# Obtain model predicted lines
best_fruit <- ggpredict(best, terms = c("fruit"))
best_human <- ggpredict(best, terms = c("human"))
best_tree <- ggpredict(best, terms = c("tree"))
best_shrub <- ggpredict(best, terms = c("shrub"))
# Print predicted and CI values of each variable
best_fruit## # Predicted probabilities of kill
##
## fruit | Predicted | 95% CI
## ------------------------------
## 5 | 0.78 | 0.62, 0.88
## 10 | 0.52 | 0.37, 0.67
## 15 | 0.25 | 0.10, 0.50
## 20 | 0.09 | 0.02, 0.36
## 30 | 0.01 | 0.00, 0.17
## 40 | 0.00 | 0.00, 0.07
## 50 | 0.00 | 0.00, 0.03
##
## Adjusted for:
## * human = 100.00
## * shrub = 10.00
## * tree = 10.00best_human## # Predicted probabilities of kill
##
## human | Predicted | 95% CI
## ------------------------------
## 10 | 0.36 | 0.19, 0.57
## 80 | 0.48 | 0.33, 0.64
## 155 | 0.61 | 0.46, 0.75
## 225 | 0.72 | 0.55, 0.85
## 295 | 0.81 | 0.60, 0.93
## 365 | 0.88 | 0.65, 0.97
## 435 | 0.92 | 0.69, 0.98
## 580 | 0.97 | 0.75, 1.00
##
## Adjusted for:
## * fruit = 10.00
## * shrub = 10.00
## * tree = 10.00best_tree## # Predicted probabilities of kill
##
## tree | Predicted | 95% CI
## -----------------------------
## 5 | 0.68 | 0.47, 0.83
## 10 | 0.52 | 0.37, 0.67
## 15 | 0.35 | 0.21, 0.54
## 20 | 0.22 | 0.08, 0.46
## 25 | 0.12 | 0.03, 0.41
## 30 | 0.07 | 0.01, 0.36
##
## Adjusted for:
## * fruit = 10.00
## * human = 100.00
## * shrub = 10.00best_shrub## # Predicted probabilities of kill
##
## shrub | Predicted | 95% CI
## ------------------------------
## 5 | 0.63 | 0.43, 0.80
## 10 | 0.52 | 0.37, 0.67
## 15 | 0.40 | 0.25, 0.58
## 20 | 0.30 | 0.13, 0.54
## 25 | 0.21 | 0.06, 0.52
## 30 | 0.15 | 0.03, 0.51
##
## Adjusted for:
## * fruit = 10.00
## * human = 100.00
## * tree = 10.00# Plot the model predicted lines
best_fruit_graph <-
ggplot() +
geom_point(aes(x = fruit, y = kill),
data = oilpalm,
size = 0.5) +
geom_line(aes(x = x, y = predicted),
data = best_fruit) +
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high),
data = best_fruit,
alpha = 0.3,
fill = "violet") +
labs(x = "Distance to nearest fruit tree (m)",
y = "Estimated probability of roadkills") +
theme_minimal()
best_human_graph <-
ggplot() +
geom_point(aes(x = human, y = kill),
data = oilpalm,
size = 0.5) +
geom_line(aes(x = x, y = predicted),
data = best_human) +
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high),
data = best_human,
alpha = 0.3,
fill = "lightgreen") +
labs(x = "Distance to nearest human dwelling (m)",
y = "Estimated probability of roadkills") +
theme_minimal()
best_tree_graph <-
ggplot() +
geom_point(aes(x = tree, y = kill),
data = oilpalm,
size = 0.5) +
geom_line(aes(x = x, y = predicted),
data = best_tree) +
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high),
data = best_tree,
alpha = 0.3,
fill = "lightpink") +
labs(x = "Distance to nearest non-fruiting tree (m)",
y = "Estimated probability of roadkills") +
theme_minimal()
best_shrub_graph <-
ggplot() +
geom_point(aes(x = shrub, y = kill),
data = oilpalm,
size = 0.5) +
geom_line(aes(x = x, y = predicted),
data = best_shrub) +
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high),
data = best_shrub,
alpha = 0.3,
fill = "lightblue") +
labs(x = "Distance to nearest shrub (m)",
y = "Estimated probability of roadkills") +
theme_minimal()
# Plot all graphs together
(best_fruit_graph + best_human_graph) / (best_tree_graph + best_shrub_graph) +
plot_annotation(tag_levels = "A", title = "Predicted probability of kill") &
theme(axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
Fig 1Bb. Probability of roadkills for each variable with 95% confidence intervals
Q1C
Using the best model from your candidate set, report and interpret your results.
# Use type 2 anova
Anova(best, type = "II")## Analysis of Deviance Table (Type II tests)
##
## Response: kill
## LR Chisq Df Pr(>Chisq)
## fruit 18.0877 1 2.11e-05 ***
## human 9.5854 1 0.001961 **
## shrub 4.3765 1 0.036437 *
## tree 8.3797 1 0.003794 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From type II Anova results, we can see that distance from nearest fruiting tree (χ² = 18.1, df = 1, p-value = <0.01), human dwelling (χ² = 9.58, df = 1, p-value = 0.0019), shrub (χ² = 4.38, df = 1, p-value = 0.0364) and non-fruiting tree (χ² = 8.38, df = 1, p-value = 0.00379) has a statistically significant effect on the probability of roadkill as their p-value is less than 0.05.
Question 2: Interactions and Non-Linear Data
Load and review data
trees <- read_excel("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 1/Assignment 1 data.xls", sheet = "Dioecious trees")
# Check data structure and summary
str(trees)## tibble [50 × 3] (S3: tbl_df/tbl/data.frame)
## $ SEX : num [1:50] 2 1 2 1 2 1 2 2 1 1 ...
## $ DBH : num [1:50] 174 171 212 187 227 274 86 161 285 182 ...
## $ FLOWERS: num [1:50] 234 110 512 133 1107 ...summary(trees)## SEX DBH FLOWERS
## Min. :1.0 Min. : 68.0 Min. : 1.0
## 1st Qu.:1.0 1st Qu.:140.5 1st Qu.: 115.8
## Median :2.0 Median :188.0 Median : 253.0
## Mean :1.6 Mean :191.6 Mean : 383.5
## 3rd Qu.:2.0 3rd Qu.:242.0 3rd Qu.: 533.0
## Max. :2.0 Max. :316.0 Max. :1482.0# Convert SEX to factor variable
trees$SEX <- as.factor(trees$SEX)
str(trees)## tibble [50 × 3] (S3: tbl_df/tbl/data.frame)
## $ SEX : Factor w/ 2 levels "1","2": 2 1 2 1 2 1 2 2 1 1 ...
## $ DBH : num [1:50] 174 171 212 187 227 274 86 161 285 182 ...
## $ FLOWERS: num [1:50] 234 110 512 133 1107 ...Q2A
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? ### Plot the effect of DBH and SEX on the number of flowers
trees %>%
ggplot(aes(x = DBH, y = FLOWERS, colour = SEX)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) + # draw best fit line
labs(x = "DBH", y = "Number of Flowers", colour = "Sex") +
scale_colour_manual(name = "Sex",
values = c("lightblue", "lightpink"),
labels = c("1" = "Male", "2" = "Female")) +
theme_minimal()
Fig 2A. Scatterplot of how SEX and DBH affects the number of flowers
For male and female trees, as DBH increases, number of flowers increases as well. However, the number of flowers produced by female trees increases at a much faster rate than male ones as DBH increases.
Q2B
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.
# Fit a linear model with both explanatory variables
trees_m1 <- lm(FLOWERS ~ SEX * DBH, data = trees)
# Check model fit
anova(trees_m1)## Analysis of Variance Table
##
## Response: FLOWERS
## Df Sum Sq Mean Sq F value Pr(>F)
## SEX 1 171841 171841 8.1186 0.006534 **
## DBH 1 5868927 5868927 277.2748 < 2.2e-16 ***
## SEX:DBH 1 174671 174671 8.2523 0.006138 **
## Residuals 46 973657 21166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Plot diagnostics
par(mfrow = c(2, 2))
plot(trees_m1)
Fig 2B. Diagnostic plots of the linear model: FLOWERS ~ DHB * SEX
From the ANOVA table, interaction between SEX and DBH is statistically significant (F1,46 = 8.25, p = 0.00614) as p-value falls below 0.05. However, from the diagnostic plots, the residuals vs fitted graph shows that the linear model may not be the best fit for the data as the points are not randomly distributed about the y = 0 line and shows funnelling. The residuals mostly follow the QQ line but deviates near the upper tail suggesting non-normal residuals. The scale-location graph shows that the residuals have slight heteroscedasticity as the points increase near the right end.
Q2C
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.
# Transform the data
# Raw y
t1 <- trees %>%
ggplot(aes(x = FLOWERS, fill = SEX)) +
geom_histogram(binwidth = 40) +
labs(x = "Number of Flowers", y = "Count") +
scale_fill_manual(
name = "Sex",
values = c("lightblue", "lightpink"),
labels = c("1" = "Male", "2" = "Female")
)
# sqrt(y)
t2 <- trees %>%
ggplot(aes(x = sqrt(FLOWERS), fill = SEX)) +
geom_histogram(binwidth = 1) +
labs(x = "sqrt(Number of Flowers)", y = "Count") +
scale_fill_manual(
name = "Sex",
values = c("lightblue", "lightpink"),
labels = c("1" = "Male", "2" = "Female")
)
# log(y)
t3 <- trees %>%
ggplot(aes(x = log(FLOWERS), fill = SEX)) +
geom_histogram(binwidth = 0.2) +
labs(x = "log(Number of Flowers)", y = "Count") +
scale_fill_manual(
name = "Sex",
values = c("lightblue", "lightpink"),
labels = c("1" = "Male", "2" = "Female")
)
# Combine the graphs
t1 + t2 + t3 +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = "A")
Fig 2Ca. Histogram of different y transformation data
A square root transformation looks the most normally distributed, while the log-transformed looks left-skewed hence the square root transformation will be selected for subsequent tests.
# Fit a linear model with both explanatory variables and sqrt transformed y
trees_m2 <- lm(sqrt(FLOWERS) ~ DBH * SEX, data = trees)
# Check model fit
anova(trees_m2)## Analysis of Variance Table
##
## Response: sqrt(FLOWERS)
## Df Sum Sq Mean Sq F value Pr(>F)
## DBH 1 4320.2 4320.2 782.7191 < 2.2e-16 ***
## SEX 1 387.2 387.2 70.1451 8.324e-11 ***
## DBH:SEX 1 34.8 34.8 6.3002 0.01565 *
## Residuals 46 253.9 5.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model predicted lines and CI
pred_m2 <- ggpredict(trees_m2, terms = c("DBH", "SEX"))# Plot model predicted lines and CI
ggplot() +
geom_point(aes(x = DBH, y = FLOWERS, colour = SEX),
data = trees) + # original
geom_line(aes(x = x, y = predicted, colour = group),
data = pred_m2) + # predicted lines
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high, fill = group),
alpha = 0.6, # 95% CI
data = pred_m2) +
labs(x = "DBH", y = "Number of Flowers", colour = "Sex") +
guides(fill = FALSE) +
scale_colour_manual(values = c("1" = "blue", "2" = "hotpink"),
labels = c("1" = "Male", "2" = "Female")) +
scale_fill_manual(values = c("1" = "lightblue", "2" = "lightpink"))
Fig 2Cb. Model predicted lines and 95% CI, sqrt transformed y LM
par(mfrow = c(2, 2))
plot(trees_m2)
Fig 2Cc. Diagnostic plots for transformed data with model
Based on the ANOVA table, interaction between DBH and SEX is statistically significant (F1,46 = 6.30, p = 0.0157). The transformed square root y data with its predicted lines look a much better fit (Fig 2Cb) than the linear model (Fig 2A) where the shortest distance between points and the best fit line is reduced. After transformation, the residuals look more normally distributed about the y = 0 line (Fig 2Cc).
Q2D
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 Poisson
trees_m3 <- glm(FLOWERS ~ DBH * SEX, family = poisson, data = trees)
# Model predicted lines and CI
pred_m3 <- ggpredict(trees_m3, terms = c("DBH", "SEX"))
# Plot model predicted lines and CI
ggplot() +
geom_point(aes(x = DBH, y = FLOWERS, colour = SEX),
data = trees) + # original
geom_line(aes(x = x, y = predicted, colour = group),
data = pred_m3) + # predicted lines
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high, fill = group),
alpha = 0.6, # 95% CI
data = pred_m3) +
labs(x = "DBH", y = "Number of Flowers", colour = "Sex") +
guides(fill = FALSE) +
scale_colour_manual(values = c("1" = "blue", "2" = "hotpink"),
labels = c("1" = "Male", "2" = "Female")) +
scale_fill_manual(values = c("1" = "lightblue", "2" = "lightpink"))
Fig 2Da. Model predicted lines and 95% CI, raw y Poisson GLM
# Diagnostic plots
par(mfrow = c(2, 2))
plot(trees_m3)
Fig 2Db. Diagnostic plots
# ANOVA, summary, dispersion parameter
anova(trees_m3) # deviance for GLM## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: FLOWERS
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 49 17921.6
## DBH 1 14343.9 48 3577.6 < 2.2e-16 ***
## SEX 1 1090.9 47 2486.7 < 2.2e-16 ***
## DBH:SEX 1 28.1 46 2458.7 1.168e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(trees_m3)##
## 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: 5The 95% confidence interval for the Poisson distribution GLM is much smaller than that of an LM (Fig 2Da). The residuals vs fitted line is not as horizontal and points are not randomly distributed about the line as well compared to an LM (Fig 2Db). Additionally, there are a few points that line outside of Cook’s distance line hence considered as outliers, likely due to the smaller 95% CI. From ANOVA, the interaction between DBH and SEX is still statistically significant (χ²1,46 = 28.1, p = <0.01).
From summary, the ratio of residual deviance to df is 2458.7/46 = 53.45, suggesting overdispersion as this ratio is much higher than 1.
Q2E
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.
# Quasipoisson model
trees_m4 <- glm(FLOWERS ~ DBH * SEX, family = quasipoisson, data = trees)
# Check for significant interaction, F test for quasi
anova(trees_m4, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: FLOWERS
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 49 17921.6
## DBH 1 14343.9 48 3577.6 282.4870 < 2.2e-16 ***
## SEX 1 1090.9 47 2486.7 21.4837 2.958e-05 ***
## DBH:SEX 1 28.1 46 2458.7 0.5529 0.4609
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(trees_m4)##
## 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: 5Since interaction between DBH and SEX is not statistically significant (F1,46 = 0.553, p = 0.461), the interaction term can be removed before calculating predicted means and confidence intervals.
# Quasipoisson model with no interaction
trees_m5 <- glm(FLOWERS ~ DBH + SEX, family = quasipoisson, data = trees)
anova(trees_m5, test = "F")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: FLOWERS
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 49 17921.6
## DBH 1 14343.9 48 3577.6 284.572 < 2.2e-16 ***
## SEX 1 1090.9 47 2486.7 21.642 2.698e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(trees_m5)##
## Call:
## glm(formula = FLOWERS ~ DBH + SEX, family = quasipoisson, data = trees)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6343031 0.2301817 11.444 3.45e-15 ***
## DBH 0.0135204 0.0008678 15.579 < 2e-16 ***
## SEX2 0.4982782 0.1097203 4.541 3.89e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 50.4053)
##
## Null deviance: 17921.6 on 49 degrees of freedom
## Residual deviance: 2486.7 on 47 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Calculate model predicted lines and CI
pred_m5 <- ggpredict(trees_m5, terms = c("DBH", "SEX"))
# Plot model predicted lines and CI
ggplot() +
geom_point(aes(x = DBH, y = FLOWERS, colour = SEX),
data = trees) + # original
geom_line(aes(x = x, y = predicted, colour = group),
data = pred_m5) + # predicted lines
geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high, fill = group),
alpha = 0.6, # 95% CI
data = pred_m5) +
labs(x = "DBH", y = "Number of Flowers", colour = "Sex") +
guides(fill = FALSE) +
scale_colour_manual(values = c("1" = "blue", "2" = "hotpink"),
labels = c("1" = "Male", "2" = "Female")) +
scale_fill_manual(values = c("1" = "lightblue", "2" = "lightpink"))
Fig 2Ea. Model predicted lines and 95% CI, quasipoisson model
# Diagnostic plots
par(mfrow = c(2, 2))
plot(trees_m5)
Fig 2Eb. Diagnostic plots
For the quasipoisson GLM, the outliers flagged out in Poisson GLM have been addressed and there are no more outliers in this new model. The residuals are generally normally distributed from the QQ plot with slight deviation at the upper tail, but this suggests that the model is a good fit, in addition to the points randomly scattered around the residuals v fitted line (Fig 2Eb). DBH and SEX have a statistically significant effect on number of flowers as their p-values are less than 0.05 (DBH: F1,48 = 285, p = <0.01; SEX: F1,47 = 21.6, p = <0.01).
Q2F
From the above models, which model would you select as the final model? Explain your reasoning. The quasipoisson model as it can handle non-normal data like the number of flowers which is a count data as it is bounded by 0, unlike the LM which involved square root transform of y to get a better fit. It is able to plot the data on the original data without needing to transform the response to fit the model, and it can handle outliers better than Poisson after adjusting for overdispersion.
What are the advantages and disadvantages associated with running a linear model or a generalised linear model? Linear models are easier to interpret and simpler but have stricter assumptions like linearity, normally distributed residuals, constant variance, and independent observations. Count, proportion and binomial data rarely fulfill these assumptions as they are bounded by 0 and or 1 thus linear models are less suited for these data and generalised linear models can handle them better. However, GLMs often are harder to interpret as coefficients are on the link scale thus some manipulation is needed to obtain comparable values.
Question 3: Multivariate Analysis
Load and review data
dung <- read.csv("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 4/Assignment 4 data/dungbeetletraptype.csv")
# Check data structure and summary
str(dung)## 'data.frame': 304 obs. of 46 variables:
## $ trapID : chr "2005_FITs-BRL.2.1" "2005_FITs-BRL.2.2" "2005_FITs-BRL.2.3" "2005_FITs-BRL.2.4" ...
## $ year : int 2005 2005 2005 2005 2005 2005 2005 2005 2005 2005 ...
## $ trap : chr "FIT" "FIT" "FIT" "FIT" ...
## $ habitat : chr "Old-growth Forest" "Old-growth Forest" "Old-growth Forest" "Old-growth Forest" ...
## $ Catharsius.dayacus : int 0 0 1 0 0 0 0 0 0 0 ...
## $ Catharsius.renaudpauliani : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Copris.agnus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Copris.sinicus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Microcopris.doriae : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Microcopris.hidakai : int 0 0 0 0 1 0 1 0 0 0 ...
## $ Ochicanthon.worroae : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Oniticellus.tessellatus : int 2 3 2 5 4 6 2 4 0 2 ...
## $ Onthophagus..Indachorius..aff.aereopictus: int 3 6 2 2 2 0 2 0 3 0 ...
## $ Onthophagus..Indachorius..aff.cheyi : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus..Indachorius..aff.danumensis : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus..Indachorius..aff.woroae : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.aff.pacificus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.aff.phanaeides : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.aff.rutilans : int 0 0 0 0 1 0 0 2 0 0 ...
## $ Onthophagus.aphodiodes : int 0 0 0 1 0 0 0 0 0 0 ...
## $ Onthophagus.borneensis : int 0 0 0 0 0 0 3 0 1 1 ...
## $ Onthophagus.cervicapra : int 0 3 2 2 1 2 1 2 2 1 ...
## $ Onthophagus.deflexicollis : int 1 1 2 0 0 0 1 1 0 0 ...
## $ Onthophagus.hirsutulus : int 1 2 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.incisus : int 0 0 0 1 0 0 2 0 1 1 ...
## $ Onthophagus.johkii : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.mulleri : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.obscurior : int 0 0 0 0 0 0 0 0 0 1 ...
## $ Onthophagus.pacificus : int 0 0 3 2 0 1 2 0 2 2 ...
## $ Onthophagus.parviobscurior : int 0 1 0 0 0 0 0 2 0 0 ...
## $ Onthophagus.pavidus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.rorarius : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.rudis : int 5 7 9 6 8 4 4 2 2 5 ...
## $ Onthophagus.rugicollis : int 5 3 3 2 2 8 5 7 6 7 ...
## $ Onthophagus.sarawacus : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.semiaureus : int 2 4 5 2 3 6 5 6 7 6 ...
## $ Onthophagus.semicupreus : int 1 3 2 1 2 2 2 3 4 0 ...
## $ Onthophagus.sp.08 : int 3 0 4 1 0 2 2 4 3 2 ...
## $ Onthophagus.sp.20 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Onthophagus.vethi : int 0 0 0 1 0 0 0 1 0 0 ...
## $ Onthophagus.vulpes : int 0 2 2 1 3 2 2 3 1 4 ...
## $ Paragymnopleurus.maurus : int 0 0 0 0 0 1 0 0 0 0 ...
## $ Paragymnopleurus.sparsus : int 2 3 1 2 1 0 0 0 0 5 ...
## $ Paragymnopleurus.striatus : int 0 0 0 0 1 0 0 0 0 0 ...
## $ Proagoderus.watanabei : int 0 0 0 0 0 0 0 0 1 0 ...
## $ Sisyphus.thoracicus : int 4 6 5 3 6 7 9 5 8 6 ...summary(dung)## trapID year trap habitat
## Length:304 Min. :2005 Length:304 Length:304
## Class :character 1st Qu.:2005 Class :character Class :character
## Mode :character Median :2012 Mode :character Mode :character
## Mean :2009
## 3rd Qu.:2012
## Max. :2012
## Catharsius.dayacus Catharsius.renaudpauliani Copris.agnus
## Min. : 0.000 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.0000
## Median : 0.000 Median : 0.000 Median : 0.0000
## Mean : 1.217 Mean : 3.191 Mean : 0.2368
## 3rd Qu.: 0.000 3rd Qu.: 2.000 3rd Qu.: 0.0000
## Max. :23.000 Max. :32.000 Max. :13.0000
## Copris.sinicus Microcopris.doriae Microcopris.hidakai Ochicanthon.worroae
## Min. : 0.0000 Min. : 0.0000 Min. : 0.0000 Min. :0.0000
## 1st Qu.: 0.0000 1st Qu.: 0.0000 1st Qu.: 0.0000 1st Qu.:0.0000
## Median : 0.0000 Median : 0.0000 Median : 0.0000 Median :0.0000
## Mean : 0.2928 Mean : 0.4375 Mean : 0.5757 Mean :0.2138
## 3rd Qu.: 0.0000 3rd Qu.: 0.0000 3rd Qu.: 0.0000 3rd Qu.:0.0000
## Max. :20.0000 Max. :36.0000 Max. :33.0000 Max. :4.0000
## Oniticellus.tessellatus Onthophagus..Indachorius..aff.aereopictus
## Min. : 0.000 Min. : 0.0000
## 1st Qu.: 0.000 1st Qu.: 0.0000
## Median : 0.000 Median : 0.0000
## Mean : 1.062 Mean : 0.4868
## 3rd Qu.: 2.000 3rd Qu.: 0.0000
## Max. :14.000 Max. :12.0000
## Onthophagus..Indachorius..aff.cheyi Onthophagus..Indachorius..aff.danumensis
## Min. : 0.0000 Min. : 0.0000
## 1st Qu.: 0.0000 1st Qu.: 0.0000
## Median : 0.0000 Median : 0.0000
## Mean : 0.4737 Mean : 0.3026
## 3rd Qu.: 0.0000 3rd Qu.: 0.0000
## Max. :14.0000 Max. :13.0000
## Onthophagus..Indachorius..aff.woroae Onthophagus.aff.pacificus
## Min. : 0.0000 Min. :0.0000
## 1st Qu.: 0.0000 1st Qu.:0.0000
## Median : 0.0000 Median :0.0000
## Mean : 0.4178 Mean :0.3059
## 3rd Qu.: 0.0000 3rd Qu.:0.0000
## Max. :37.0000 Max. :7.0000
## Onthophagus.aff.phanaeides Onthophagus.aff.rutilans Onthophagus.aphodiodes
## Min. :0.00000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.00000 Median :0.0000 Median :0.0000
## Mean :0.09868 Mean :0.1776 Mean :0.1974
## 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000
## Max. :6.00000 Max. :4.0000 Max. :5.0000
## Onthophagus.borneensis Onthophagus.cervicapra Onthophagus.deflexicollis
## Min. : 0.0000 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 0.0000
## Median : 0.0000 Median : 3.000 Median : 0.0000
## Mean : 0.7303 Mean : 8.599 Mean : 0.7763
## 3rd Qu.: 0.0000 3rd Qu.: 8.000 3rd Qu.: 1.0000
## Max. :14.0000 Max. :125.000 Max. :18.0000
## Onthophagus.hirsutulus Onthophagus.incisus Onthophagus.johkii
## Min. :0.0000 Min. : 0.000 Min. : 0.0000
## 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.: 0.0000
## Median :0.0000 Median : 2.000 Median : 0.0000
## Mean :0.3355 Mean : 3.796 Mean : 0.4342
## 3rd Qu.:0.0000 3rd Qu.: 5.000 3rd Qu.: 0.0000
## Max. :5.0000 Max. :55.000 Max. :17.0000
## Onthophagus.mulleri Onthophagus.obscurior Onthophagus.pacificus
## Min. : 0.0000 Min. : 0.00 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.: 0.00 1st Qu.: 0.000
## Median : 0.0000 Median : 0.00 Median : 0.000
## Mean : 0.8092 Mean : 2.26 Mean : 0.875
## 3rd Qu.: 1.0000 3rd Qu.: 3.00 3rd Qu.: 1.000
## Max. :11.0000 Max. :34.00 Max. :10.000
## Onthophagus.parviobscurior Onthophagus.pavidus Onthophagus.rorarius
## Min. : 0.000 Min. : 0.0000 Min. :0.0000
## 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0.0000
## Median : 0.000 Median : 0.0000 Median :0.0000
## Mean : 2.023 Mean : 0.3651 Mean :0.2697
## 3rd Qu.: 0.000 3rd Qu.: 0.0000 3rd Qu.:0.0000
## Max. :56.000 Max. :32.0000 Max. :5.0000
## Onthophagus.rudis Onthophagus.rugicollis Onthophagus.sarawacus
## Min. :0.000 Min. : 0.000 Min. :0.0000
## 1st Qu.:0.000 1st Qu.: 0.000 1st Qu.:0.0000
## Median :0.000 Median : 1.000 Median :0.0000
## Mean :1.533 Mean : 2.993 Mean :0.3257
## 3rd Qu.:2.000 3rd Qu.: 3.000 3rd Qu.:0.0000
## Max. :9.000 Max. :55.000 Max. :8.0000
## Onthophagus.semiaureus Onthophagus.semicupreus Onthophagus.sp.08
## Min. : 0.000 Min. : 0.0000 Min. :0.000
## 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0.000
## Median : 1.000 Median : 0.0000 Median :0.000
## Mean : 1.618 Mean : 0.8289 Mean :0.273
## 3rd Qu.: 2.250 3rd Qu.: 1.0000 3rd Qu.:0.000
## Max. :10.000 Max. :13.0000 Max. :5.000
## Onthophagus.sp.20 Onthophagus.vethi Onthophagus.vulpes Paragymnopleurus.maurus
## Min. : 0.0000 Min. :0.0000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 0.0000 Median :0.0000 Median : 1.000 Median : 0.000
## Mean : 0.6151 Mean :0.1842 Mean : 3.286 Mean : 2.691
## 3rd Qu.: 0.0000 3rd Qu.:0.0000 3rd Qu.: 3.000 3rd Qu.: 1.000
## Max. :14.0000 Max. :4.0000 Max. :82.000 Max. :93.000
## Paragymnopleurus.sparsus Paragymnopleurus.striatus Proagoderus.watanabei
## Min. : 0.000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 2.000 Median : 0.000 Median : 1.000
## Mean : 4.826 Mean : 1.128 Mean : 4.855
## 3rd Qu.: 5.000 3rd Qu.: 0.000 3rd Qu.: 5.000
## Max. :71.000 Max. :39.000 Max. :62.000
## Sisyphus.thoracicus
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 2.000
## Mean : 5.523
## 3rd Qu.: 6.000
## Max. :92.000# Filter only relevant columns
env <- dung[,1:4]
species <- dung[,5:46]To visualise and determine the community composition of dung beetles NMDS will be used as it is good to simplify complex multivariate data into just 2D in this case. This will help with visualising whether there is overlap of community composition of dung beetles between the 2 traps or not.
Q3A
Based on the analysis you have selected, plot and interpret your ordination plot.
# AI has helped in generating this portion of code, comments have been self-written
library(vegan)
# Large count range of beetles, helps to reduce skewness of count
species_sqrt <- sqrt(species)
# Prevent plot from changing configuration, keeps stress value constant
set.seed(123)
nmds <- metaMDS(species_sqrt, distance = "bray", # bray curtis dissimilarity for sp counts
k = 3, # returns higher linear R^2 fit than 2
trymax = 100,
autotransform = FALSE) ## Run 0 stress 0.1500979
## Run 1 stress 0.1516075
## Run 2 stress 0.1518805
## Run 3 stress 0.1519982
## Run 4 stress 0.1521327
## Run 5 stress 0.1510832
## Run 6 stress 0.1515884
## Run 7 stress 0.1506867
## Run 8 stress 0.1510136
## Run 9 stress 0.1503572
## ... Procrustes: rmse 0.004623916 max resid 0.06803683
## Run 10 stress 0.1528186
## Run 11 stress 0.1501245
## ... Procrustes: rmse 0.001803074 max resid 0.0252112
## Run 12 stress 0.150426
## ... Procrustes: rmse 0.005147621 max resid 0.0678934
## Run 13 stress 0.1506639
## Run 14 stress 0.1513501
## Run 15 stress 0.1503501
## ... Procrustes: rmse 0.004200838 max resid 0.0678212
## Run 16 stress 0.1506543
## Run 17 stress 0.1566428
## Run 18 stress 0.1502198
## ... Procrustes: rmse 0.003505473 max resid 0.03804821
## Run 19 stress 0.1503806
## ... Procrustes: rmse 0.004550322 max resid 0.06795373
## Run 20 stress 0.1518811
## Run 21 stress 0.1523813
## Run 22 stress 0.1524573
## Run 23 stress 0.1503507
## ... Procrustes: rmse 0.004208162 max resid 0.06785644
## Run 24 stress 0.1520601
## Run 25 stress 0.1510834
## Run 26 stress 0.1500982
## ... Procrustes: rmse 0.0003498857 max resid 0.003893322
## ... Similar to previous best
## *** Best solution repeated 1 times# Plot NMDS stress plot
nmds$stress## [1] 0.1500979stressplot(nmds)
Fig 3Aa. NMDS stress plot. Fig 3Ab. NMDS plot of dung beetle composition separated by trap type.
# Separate coordinates by trap type
nmds_df <- as.data.frame(scores(nmds, display = "sites"))
nmds_df$trap <- env$trap
ggplot(nmds_df, aes(NMDS1, NMDS2, colour = trap, fill = trap)) +
geom_point(alpha = 0.6) +
stat_ellipse(geom = "polygon", level = 0.95, alpha = 0.2) + # colour ellipse
labs(colour = "Trap type", fill = "Trap type") +
theme_minimal()
Fig 3Aa. NMDS stress plot. Fig 3Ab. NMDS plot of dung beetle composition separated by trap type.
Stress = 0.150
Q3B
Using appropriate statistical tests, report and interpret your results.
# AI has helped in generating this portion of code, comments have been self-written
dist_bc <- vegdist(species_sqrt, method = "bray")
# Check if community centroids differ between trap types using PERMANOVA
adonis2(dist_bc ~ trap, data = env, permutations = 999)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist_bc ~ trap, data = env, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 13.966 0.19929 75.164 0.001 ***
## Residual 302 56.113 0.80071
## Total 303 70.078 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Check if the spread within each group the same using PERMDISP
bd <- betadisper(dist_bc, env$trap)
anova(bd)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.0006 0.0005802 0.0482 0.8265
## Residuals 302 3.6388 0.0120491permutest(bd, permutations = 999)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0006 0.0005802 0.0482 999 0.832
## Residuals 302 3.6388 0.0120491boxplot(bd)
Fig 3B. Boxplot of distance to centroid for both trap types.
Dung beetle community composition differed significantly between pitfall and flight interception traps (PERMANOVA: pseudo-F1,302 = 75.16, R² = 0.20, p = 0.001), with trap type explaining about 20% of the variation in community composition. Multivariate dispersion did not differ between trap types (PERMDISP: F1,302 = 0.048, p = 0.83), indicating that the difference reflects a shift in community composition and not differences in within-group variability. This is consistent with the NMDS ordination where the two trap types formed distinct clusters with similar spread. From the box plot (Fig 3B), the medians for BPT and FIT are similar and the box sizes overlap almost completely, agreeing with PERMDISP where the two trap types have similar spread.