ES3307 Project 2026
Misha Rashanah
2026-11-13
RMARKDOWN FILE FOR ES3307 PROJECT
Far Eastern Curlews & Mudflats
Introduction
Research Gap:
Mudflats in Bako Buntal Bay have been recognised as a significant stopover for Far Eastern Curlews (FEC), who migrate along the East Asian-Australasian Flyway (EAAF)
Yet the drivers behind the local abundance patterns across mudflats remain poorly understood
Identify how prey density and tidal exposure determine FEC numbers so as to allow for targeted site-specific conservation
Research Question:
How do differences in tidal exposure and prey density influence Far Eastern Curlew abundance across mudflat sites during the non-breeding season ?
Hypothesis:
A higher prey density and tidal exposure will lead to more FEC abundance as they indicate more food availability and increase in fitness, where it is crucial for FECs to replenish and store up energy in preparation for their long migration journeys.
Predictions:
Increase in prey density → increase in FEC abundance
Increase in tidal exposure → increase in FEC abundance
Experimental Design Diagram
5 mudflats x 3 Tidal Exposure Levels x 3 Prey Density Levels x 5 samples = 225 samples
Fixed Effects: Tidal Exposure & Prey Density (2 treatments)
Random Effects: Mudflats & Plots
CODE
Load Relevant Packages
Creation of Dummy Dataset
# install.packages("truncnorm")
library(truncnorm)
set.seed(102) # for reproducibility
# 3 levels --> categorical variable
tide_sequence <- rep(c(rep("Low", length.out = 15), rep("Med", length.out = 15), rep("High", length.out = 15)), length.out = 225)
# 3 levels --> categorical variable
prey_sequence <- rep(c(rep("Low", length.out = 5), rep("Med", length.out = 5), rep("High", length.out = 5)), length.out = 225)
df <- data.frame(Abundance = round(rtruncnorm(n = 225, a = 0, b = Inf, mean = 3, sd = 2)), Mudflat = rep(1:5, each = 45), TidalExposure = tide_sequence, PreyDensity = prey_sequence )
# except for
df$Mudflat <- as.factor(df$Mudflat)
df$TidalExposure <- as.factor(df$TidalExposure)
df$PreyDensity <- as.factor(df$PreyDensity)
df$Plot <- interaction(df$Mudflat, df$TidalExposure, df$PreyDensity, drop = TRUE)
# df$Plot1 <- interaction(df$TidalExposure, df$PreyDensity, drop = TRUE)
# plot represents the 9 combinations of tidal exposure and prey density
df ## Abundance Mudflat TidalExposure PreyDensity Plot
## 1 3 1 Low Low 1.Low.Low
## 2 5 1 Low Low 1.Low.Low
## 3 0 1 Low Low 1.Low.Low
## 4 7 1 Low Low 1.Low.Low
## 5 5 1 Low Low 1.Low.Low
## 6 5 1 Low Med 1.Low.Med
## 7 5 1 Low Med 1.Low.Med
## 8 4 1 Low Med 1.Low.Med
## 9 4 1 Low Med 1.Low.Med
## 10 7 1 Low Med 1.Low.Med
## 11 6 1 Low High 1.Low.High
## 12 3 1 Low High 1.Low.High
## 13 2 1 Low High 1.Low.High
## 14 0 1 Low High 1.Low.High
## 15 3 1 Low High 1.Low.High
## 16 2 1 Med Low 1.Med.Low
## 17 7 1 Med Low 1.Med.Low
## 18 4 1 Med Low 1.Med.Low
## 19 5 1 Med Low 1.Med.Low
## 20 3 1 Med Low 1.Med.Low
## 21 0 1 Med Med 1.Med.Med
## 22 1 1 Med Med 1.Med.Med
## 23 2 1 Med Med 1.Med.Med
## 24 3 1 Med Med 1.Med.Med
## 25 3 1 Med Med 1.Med.Med
## 26 3 1 Med High 1.Med.High
## 27 5 1 Med High 1.Med.High
## 28 7 1 Med High 1.Med.High
## 29 6 1 Med High 1.Med.High
## 30 3 1 Med High 1.Med.High
## 31 1 1 High Low 1.High.Low
## 32 3 1 High Low 1.High.Low
## 33 5 1 High Low 1.High.Low
## 34 7 1 High Low 1.High.Low
## 35 6 1 High Low 1.High.Low
## 36 2 1 High Med 1.High.Med
## 37 3 1 High Med 1.High.Med
## 38 9 1 High Med 1.High.Med
## 39 4 1 High Med 1.High.Med
## 40 3 1 High Med 1.High.Med
## 41 1 1 High High 1.High.High
## 42 2 1 High High 1.High.High
## 43 4 1 High High 1.High.High
## 44 1 1 High High 1.High.High
## 45 8 1 High High 1.High.High
## 46 1 2 Low Low 2.Low.Low
## 47 4 2 Low Low 2.Low.Low
## 48 2 2 Low Low 2.Low.Low
## 49 6 2 Low Low 2.Low.Low
## 50 5 2 Low Low 2.Low.Low
## 51 1 2 Low Med 2.Low.Med
## 52 3 2 Low Med 2.Low.Med
## 53 1 2 Low Med 2.Low.Med
## 54 4 2 Low Med 2.Low.Med
## 55 2 2 Low Med 2.Low.Med
## 56 6 2 Low High 2.Low.High
## 57 2 2 Low High 2.Low.High
## 58 3 2 Low High 2.Low.High
## 59 2 2 Low High 2.Low.High
## 60 2 2 Low High 2.Low.High
## 61 3 2 Med Low 2.Med.Low
## 62 7 2 Med Low 2.Med.Low
## 63 7 2 Med Low 2.Med.Low
## 64 5 2 Med Low 2.Med.Low
## 65 3 2 Med Low 2.Med.Low
## 66 1 2 Med Med 2.Med.Med
## 67 3 2 Med Med 2.Med.Med
## 68 4 2 Med Med 2.Med.Med
## 69 2 2 Med Med 2.Med.Med
## 70 4 2 Med Med 2.Med.Med
## 71 2 2 Med High 2.Med.High
## 72 3 2 Med High 2.Med.High
## 73 5 2 Med High 2.Med.High
## 74 7 2 Med High 2.Med.High
## 75 3 2 Med High 2.Med.High
## 76 3 2 High Low 2.High.Low
## 77 3 2 High Low 2.High.Low
## 78 3 2 High Low 2.High.Low
## 79 2 2 High Low 2.High.Low
## 80 5 2 High Low 2.High.Low
## 81 2 2 High Med 2.High.Med
## 82 3 2 High Med 2.High.Med
## 83 2 2 High Med 2.High.Med
## 84 3 2 High Med 2.High.Med
## 85 2 2 High Med 2.High.Med
## 86 0 2 High High 2.High.High
## 87 2 2 High High 2.High.High
## 88 1 2 High High 2.High.High
## 89 6 2 High High 2.High.High
## 90 3 2 High High 2.High.High
## 91 1 3 Low Low 3.Low.Low
## 92 1 3 Low Low 3.Low.Low
## 93 3 3 Low Low 3.Low.Low
## 94 4 3 Low Low 3.Low.Low
## 95 6 3 Low Low 3.Low.Low
## 96 0 3 Low Med 3.Low.Med
## 97 4 3 Low Med 3.Low.Med
## 98 3 3 Low Med 3.Low.Med
## 99 3 3 Low Med 3.Low.Med
## 100 2 3 Low Med 3.Low.Med
## 101 5 3 Low High 3.Low.High
## 102 5 3 Low High 3.Low.High
## 103 3 3 Low High 3.Low.High
## 104 1 3 Low High 3.Low.High
## 105 3 3 Low High 3.Low.High
## 106 4 3 Med Low 3.Med.Low
## 107 3 3 Med Low 3.Med.Low
## 108 1 3 Med Low 3.Med.Low
## 109 5 3 Med Low 3.Med.Low
## 110 2 3 Med Low 3.Med.Low
## 111 4 3 Med Med 3.Med.Med
## 112 6 3 Med Med 3.Med.Med
## 113 4 3 Med Med 3.Med.Med
## 114 4 3 Med Med 3.Med.Med
## 115 2 3 Med Med 3.Med.Med
## 116 2 3 Med High 3.Med.High
## 117 2 3 Med High 3.Med.High
## 118 4 3 Med High 3.Med.High
## 119 4 3 Med High 3.Med.High
## 120 2 3 Med High 3.Med.High
## 121 2 3 High Low 3.High.Low
## 122 3 3 High Low 3.High.Low
## 123 1 3 High Low 3.High.Low
## 124 7 3 High Low 3.High.Low
## 125 2 3 High Low 3.High.Low
## 126 5 3 High Med 3.High.Med
## 127 2 3 High Med 3.High.Med
## 128 2 3 High Med 3.High.Med
## 129 4 3 High Med 3.High.Med
## 130 2 3 High Med 3.High.Med
## 131 2 3 High High 3.High.High
## 132 3 3 High High 3.High.High
## 133 3 3 High High 3.High.High
## 134 6 3 High High 3.High.High
## 135 4 3 High High 3.High.High
## 136 1 4 Low Low 4.Low.Low
## 137 5 4 Low Low 4.Low.Low
## 138 5 4 Low Low 4.Low.Low
## 139 6 4 Low Low 4.Low.Low
## 140 5 4 Low Low 4.Low.Low
## 141 1 4 Low Med 4.Low.Med
## 142 4 4 Low Med 4.Low.Med
## 143 0 4 Low Med 4.Low.Med
## 144 6 4 Low Med 4.Low.Med
## 145 4 4 Low Med 4.Low.Med
## 146 7 4 Low High 4.Low.High
## 147 6 4 Low High 4.Low.High
## 148 3 4 Low High 4.Low.High
## 149 2 4 Low High 4.Low.High
## 150 6 4 Low High 4.Low.High
## 151 3 4 Med Low 4.Med.Low
## 152 4 4 Med Low 4.Med.Low
## 153 3 4 Med Low 4.Med.Low
## 154 1 4 Med Low 4.Med.Low
## 155 4 4 Med Low 4.Med.Low
## 156 0 4 Med Med 4.Med.Med
## 157 3 4 Med Med 4.Med.Med
## 158 5 4 Med Med 4.Med.Med
## 159 1 4 Med Med 4.Med.Med
## 160 1 4 Med Med 4.Med.Med
## 161 7 4 Med High 4.Med.High
## 162 2 4 Med High 4.Med.High
## 163 2 4 Med High 4.Med.High
## 164 7 4 Med High 4.Med.High
## 165 4 4 Med High 4.Med.High
## 166 4 4 High Low 4.High.Low
## 167 3 4 High Low 4.High.Low
## 168 0 4 High Low 4.High.Low
## 169 7 4 High Low 4.High.Low
## 170 6 4 High Low 4.High.Low
## 171 3 4 High Med 4.High.Med
## 172 4 4 High Med 4.High.Med
## 173 6 4 High Med 4.High.Med
## 174 5 4 High Med 4.High.Med
## 175 4 4 High Med 4.High.Med
## 176 0 4 High High 4.High.High
## 177 4 4 High High 4.High.High
## 178 2 4 High High 4.High.High
## 179 2 4 High High 4.High.High
## 180 3 4 High High 4.High.High
## 181 4 5 Low Low 5.Low.Low
## 182 5 5 Low Low 5.Low.Low
## 183 5 5 Low Low 5.Low.Low
## 184 2 5 Low Low 5.Low.Low
## 185 3 5 Low Low 5.Low.Low
## 186 3 5 Low Med 5.Low.Med
## 187 6 5 Low Med 5.Low.Med
## 188 3 5 Low Med 5.Low.Med
## 189 3 5 Low Med 5.Low.Med
## 190 2 5 Low Med 5.Low.Med
## 191 2 5 Low High 5.Low.High
## 192 3 5 Low High 5.Low.High
## 193 6 5 Low High 5.Low.High
## 194 4 5 Low High 5.Low.High
## 195 4 5 Low High 5.Low.High
## 196 3 5 Med Low 5.Med.Low
## 197 3 5 Med Low 5.Med.Low
## 198 1 5 Med Low 5.Med.Low
## 199 4 5 Med Low 5.Med.Low
## 200 2 5 Med Low 5.Med.Low
## 201 3 5 Med Med 5.Med.Med
## 202 6 5 Med Med 5.Med.Med
## 203 2 5 Med Med 5.Med.Med
## 204 2 5 Med Med 5.Med.Med
## 205 5 5 Med Med 5.Med.Med
## 206 0 5 Med High 5.Med.High
## 207 2 5 Med High 5.Med.High
## 208 1 5 Med High 5.Med.High
## 209 3 5 Med High 5.Med.High
## 210 1 5 Med High 5.Med.High
## 211 1 5 High Low 5.High.Low
## 212 1 5 High Low 5.High.Low
## 213 3 5 High Low 5.High.Low
## 214 5 5 High Low 5.High.Low
## 215 0 5 High Low 5.High.Low
## 216 4 5 High Med 5.High.Med
## 217 2 5 High Med 5.High.Med
## 218 5 5 High Med 5.High.Med
## 219 1 5 High Med 5.High.Med
## 220 1 5 High Med 5.High.Med
## 221 4 5 High High 5.High.High
## 222 3 5 High High 5.High.High
## 223 2 5 High High 5.High.High
## 224 2 5 High High 5.High.High
## 225 3 5 High High 5.High.Highstr(df) # check data structure ## 'data.frame': 225 obs. of 5 variables:
## $ Abundance : num 3 5 0 7 5 5 5 4 4 7 ...
## $ Mudflat : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ TidalExposure: Factor w/ 3 levels "High","Low","Med": 2 2 2 2 2 2 2 2 2 2 ...
## $ PreyDensity : Factor w/ 3 levels "High","Low","Med": 2 2 2 2 2 3 3 3 3 3 ...
## $ Plot : Factor w/ 45 levels "1.High.High",..: 21 21 21 21 21 36 36 36 36 36 ...# except for abundance which represents count data and should be numeric, the rest should be stored as factorsOur dummy dataset was created using various functions on R.
The columns headers include
Abundance (referring to the Far Eastern Curlews seen in each sample)
Mudflat (1 to 5 - Random Effect)
Tidal Exposure (Categorical Variable - Low, Medium, High) aka Treatment 1/explanatory variable/fixed effect
Prey Density (Categorical Variable - Low, Medium, High) aka Treatment 2/explanatory variable/fixed effect
Plot (Mudflat| TidalExposureLevel | PreyDensityLevel - Random Effect)
In particular, our abundance data was created using rtruncnorm (truncnorm package) to generate 225 random numbers from the truncated normal distribution, each representing the FEC abundance count, as outlined in our experimental design. This function is appropriate in ensuring there is no negative counts as that is not biologically possible. Only FEC abundance is stored numerically and represents count data, while the rest are stored as factors.
Plots is created using the unique combination of Mudflat X corssed with the respective Tidal Exposure Level further crossed with the respective Prey Density Level. Before adding the Mudflat, there should be 9 unique Tidal Exposure Level x Prey Density combinations (can be seen from the diagram above). Essentially, making plots nested in mudflats.
Visualizing Raw Data
We mostly use histogram/boxplots as our explanatory variables are categorical
library(ggplot2)
# Visualize Abundance Data
hist(df$Abundance, col = "purple", border = "white", main = "Abundance Distribution Histogram ")
# Generated Abundance distribution is right-skewed
# Visualize Abundance (response variable) against Prey Density (explanatory variable)
df %>%
ggplot(aes(x = Mudflat, y = Abundance, fill = PreyDensity)) +
geom_boxplot(
fill = "yellow",
outlier.shape = NA,
width = 0.7,
position = position_dodge(width = 0.8)
) +
geom_jitter(
aes(color = PreyDensity),
position = position_jitterdodge(
jitter.width = 0.15,
dodge.width = 0.8
),
alpha = 0.6,
size = 2
) +
labs(
title = "Far Eastern Curlew Abundance by Prey Density Across Mudflats",
x = "Mudflat",
y = "FEC Abundance",
fill = "Prey Density",
color = "Prey Density"
) +
theme_light()
# Visualize Abundance (response variable) against Tidal Exposure Level (explanatory variable)
df %>%
ggplot(aes(x = Mudflat, y = Abundance, fill = TidalExposure)) +
geom_boxplot(
fill = "orange",
outlier.shape = NA,
width = 0.7,
position = position_dodge(width = 0.8)
) +
geom_jitter(
aes(color = TidalExposure),
position = position_jitterdodge(
jitter.width = 0.15,
dodge.width = 0.8
),
alpha = 0.6,
size = 2
) +
labs(
title = "Far Eastern Curlew Abundance by Tidal Exposure Level Across Mudflats",
x = "Mudflat",
y = "FEC Abundance",
fill = "Tidal Exposure Level",
color = "Tidal Exposure Level"
) +
theme_light()
From the histogram, we can see that the generated abundance is heavily right-skewed, implying that our data is not normally distributed. Next, we plot a boxplot by overlaying the FEC abundance across the 5 mudflats for tidal exposure and prey density respectively. We can observe that for both boxplots, there is no clear linear patterns and there is unequal variance in abundance numbers across the 5 mudflats, implying heteroscedasticity. Thus, we will not use a linear model to analyze our data given that abundance is a discrete and not continuous data.
Interaction Plot
Check for any potential interactions between our two explanatory variables.
ggplot(df, aes(x = PreyDensity, y = Abundance,
colour = TidalExposure, group = TidalExposure)) +
stat_summary(fun = mean, geom = "point", size = 3) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun.data = mean_se, geom = "errorbar",
width = 0.1) +
labs(
title = "Interaction between Tidal Exposure and Prey Density",
x = "Prey Density",
y = "FEC Abundance",
colour = "Tidal Exposure"
) +
theme_minimal()
There could be some biological interactions between the tidal exposure variable and prey density variable which is why we decided to test it out using an interaction plot. From our visualization, there appears to be an interaction as seen from how the lines intersect between the Low and Medium prey density levels.
However, this pattern is not very clear so we still consider interaction as one of our model options.
Model Selection
With respect to our experimental design structure.
library(lme4)
# start with the most complex model including interaction
Model1 <- glmer(Abundance ~ TidalExposure*PreyDensity +(1|Mudflat/Plot), family = 'poisson', data = df)
# without interaction
Model2 <- glmer(Abundance ~ TidalExposure + PreyDensity +(1|Mudflat/Plot), family = 'poisson', data = df)
# try simplifying it further by removing mudflat
Model3<- glmer(Abundance ~ TidalExposure + PreyDensity +(1|Plot), family = 'poisson', data = df)
# Not correct because Mudflat is a random effect in our experiment and needs to be included
# AIC(Model1,Model2) use AIC when models have complicated or many parameters
# generally the interaction is not significant
# anova here is used as a likelihood ratio test to compare model performance
anova(Model1, Model2, test = "Chisq")## Data: df
## Models:
## Model2: Abundance ~ TidalExposure + PreyDensity + (1 | Mudflat/Plot)
## Model1: Abundance ~ TidalExposure * PreyDensity + (1 | Mudflat/Plot)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## Model2 7 919.3 943.22 -452.65 905.3
## Model1 11 924.7 962.27 -451.35 902.7 2.6058 4 0.6258# knit to table
aov_table <- as.data.frame(anova(Model1, Model2, test = "Chisq"))
knitr::kable(
aov_table,
digits = 3,
caption = "Comparing Model 1 and Model 2"
)| npar | AIC | BIC | logLik | -2*log(L) | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| Model2 | 7 | 919.303 | 943.215 | -452.651 | 905.303 | NA | NA | NA |
| Model1 | 11 | 924.697 | 962.274 | -451.348 | 902.697 | 2.606 | 4 | 0.626 |
Given that our abundance data is a count response variable and from our raw data plots we know that our data is not normally distributed, we decided to use Poisson Generalized Linear Mixed Effects Model (glmer) as our model structure, as it accounts for the hierarchical structure of our experiment, with plots nested within mudflats included as random effects.
A likelihood ratio test was conducted to compare the full model, which included the interaction between tidal exposure and prey density, with the reduced additive model. The interaction between tidal exposure and prey density was not statistically significant (\(\chi^2_4\) = 2.61 p = 0.63) (2 d.p.). This implies that there is insufficient evidence to support that the effect of prey density on FEC abundance varied across tidal exposure levels. Therefore, we decided to drop the interaction term and use the simpler additive model.
Therefore, our chosen model is model 2, which we will use in our subsequent statistical analyses.
Model Assumptions & Diagnostic Plots
Evaluate chosen model performance
library(DHARMa)
library(lme4)
# Simulate residuals
set.seed(102) # for reporducability
sim_res <- simulateResiduals(
fittedModel = Model2,
n = 1000
)
# Visual diagnostic assessment
plot(sim_res)
# Residual uniformity
testUniformity(sim_res)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033991, p-value = 0.9573
## alternative hypothesis: two-sided# Overdispersion
testDispersion(sim_res)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0408, p-value = 0.672
## alternative hypothesis: two.sided# Zero inflation
testZeroInflation(sim_res)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.2909, p-value = 0.508
## alternative hypothesis: two.sided# Check model convergence and singularity
summary(Model2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Abundance ~ TidalExposure + PreyDensity + (1 | Mudflat/Plot)
## Data: df
##
## AIC BIC logLik -2*log(L) df.resid
## 919.3 943.2 -452.7 905.3 218
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9762 -0.6591 -0.0979 0.5739 3.2697
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plot:Mudflat (Intercept) 0.005337 0.07305
## Mudflat (Intercept) 0.002832 0.05322
## Number of obs: 225, groups: Plot:Mudflat, 45; Mudflat, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.13835 0.09033 12.602 <2e-16 ***
## TidalExposureLow 0.10407 0.09343 1.114 0.265
## TidalExposureMed 0.04108 0.09475 0.434 0.665
## PreyDensityLow 0.07426 0.09228 0.805 0.421
## PreyDensityMed -0.05371 0.09503 -0.565 0.572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TdlExL TdlExM PryDnL
## TidlExpsrLw -0.543
## TidlExpsrMd -0.536 0.517
## PreyDnstyLw -0.530 0.001 0.002
## PreyDnstyMd -0.516 0.003 0.004 0.502isSingular(Model2, tol = 1e-4)## [1] FALSE# Examine random effects
VarCorr(Model2)## Groups Name Std.Dev.
## Plot:Mudflat (Intercept) 0.073053
## Mudflat (Intercept) 0.053218par(mfrow = c(1, 2))
# Check residuals against tidal exposure
plotResiduals(sim_res, df$TidalExposure)
# Check residuals against prey density
plotResiduals(sim_res, df$PreyDensity)
(AI was used to run the relevant functions from the DHARMa package and to understand the tests results)
For our diagnostic plots, we used the DHARMa package to assess the suitability of our selected poisson glmer model, instead of plot(), which plots standard residuals which are appropriate for linear models and gaussian distribution. But for our case, where it is count data, DHARMa is more appropriate as it uses simulation-based quantile residuals to guarantee an interpretable, uniform (0,1) baseline for any distributional family.
The residual uniformity test indicates that there is no significant deviation from the expected residual distribution (Kolmogorov-Smirnov test, D = 0.034, p = 0.957). The dispersion test shows no evidence of overdispersion (Dispersion = 1.04, p = 0.672), wheras the zero-inflation test suggest that there is no significant excess of zero observations relative to model expectations (ratio of observed to simulated zeros = 1.291, p = 0.508). Furthermore, the model did not exhibit a singular fit (isSingular = FALSE).
Thus, there is no diagnostic evidence suggesting that our selected Poisson GLMM was inappropriate for the dataset.
Anova Test
library(car)
Anova(Model2, type = "2")## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Abundance
## Chisq Df Pr(>Chisq)
## TidalExposure 1.2680 2 0.5305
## PreyDensity 1.9036 2 0.3860Our results indicate that Tidal Exposure has no significant effect on the abundance of the Far Eastern Curlews (\(\chi^2_2\) = 1.27 p = 0.53) (2 d.p.). Similarly, Prey Density does not significantly influence FEC abundance (\(\chi^2_2\) = 1.90, p = 0.39) (2 d.p.). These results support our findings from the boxplots, which also showed no clear patterns or directional trends in FEC abundance across different levels of tidal exposure or prey densities. Hence, there is no strong evidence to conclude that our two treatments have a statistically significant effect on the FEC abundance numbers.
Model Explanatory Power
# calculate R square
library(MuMIn)
r.squaredGLMM(Model2)## R2m R2c
## delta 0.01466709 0.04068479
## lognormal 0.01668142 0.04627231
## trigamma 0.01265789 0.03511150The chosen additive GLMM model explains a relatively small proportion of the variation in FEC abundance, with a marginal R2 of 0.013 and a conditional R2 of 0.035 (3 d.p.). This suggests that the fixed effects of Tidal Exposure and Prey Density explains abous 1.3% of the variation, while the fixed and random effects combined explains approximately 3.5%.
Estimate Differences of Each Treatment (Tidal Exposure & Prey Density)
Additional Post-Hoc Test since we have more than 2 levels for each categorical explanatory variable
library(emmeans)
# Calculate estimated marginal means and pairwise comparisons for Tidal Exposure back-transformed to the original count scale
emmeans(Model2, pairwise ~ TidalExposure, type = "response")## $emmeans
## TidalExposure rate SE df asymp.LCL asymp.UCL
## High 3.14 0.226 Inf 2.73 3.62
## Low 3.49 0.241 Inf 3.05 3.99
## Med 3.27 0.232 Inf 2.85 3.76
##
## Results are averaged over the levels of: PreyDensity
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## High / Low 0.901 0.0842 Inf 1 -1.114 0.5055
## High / Med 0.960 0.0909 Inf 1 -0.434 0.9016
## Low / Med 1.065 0.0985 Inf 1 0.681 0.7747
##
## Results are averaged over the levels of: PreyDensity
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale# Do the same for Prey Density
emmeans(Model2, pairwise ~ PreyDensity, type = "response")## $emmeans
## PreyDensity rate SE df asymp.LCL asymp.UCL
## High 3.28 0.232 Inf 2.85 3.76
## Low 3.53 0.243 Inf 3.08 4.04
## Med 3.11 0.225 Inf 2.69 3.58
##
## Results are averaged over the levels of: TidalExposure
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## High / Low 0.928 0.0857 Inf 1 -0.805 0.7001
## High / Med 1.055 0.1000 Inf 1 0.565 0.8386
## Low / Med 1.137 0.1060 Inf 1 1.369 0.3574
##
## Results are averaged over the levels of: TidalExposure
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scaleThe function emmeans is appropriate for our analysis because our model includes two categorical explanatory variables (Tidal Exposure and Prey Density). We use pairwise, which gives us the predicted means and also allows us to determine whether their different levels (low, medium, high) can influence the FEC abundance rate which is highlighted by the $contrasts feature.
From the emmeans results, we can see that the estimated FEC abundance rate was slightly greater at low tidal exposure (3.49), compared to medium (3.27) and high tidal exposure (3.14) (2 d.p.). However, Tukey-adjusted pairwise comparisons ($contrasts) reveal that there are no significant differences among tidal exposure levels (all p > 0.05). Similarly, estimated FEC abundance rate was slightly graeter at low prey density (3.53), followed by high (3.28) and medium prey density (3.11) (2 d.p.) . Furthermore, pairwise comparisons also show that there are no significant differences among prey density levels (all p > 0.05).
Therefore, there is no strong evidence to support that there are significant differences in the estimated FEC abundance rates among the theree tidal exposure levels and three prey density levels.
Conclusion
In conclusion, to answer our research question, our findings suggest that there is insufficient statistical evidence to support that tidal exposure and prey density significantly influence Far Eastern Curlew (FEC) abundance across mudflat sites in Bako Buntal Bay (Sarawak) during the non-breeding season.