Raven_GLMM
Grant Rogers
2026-06-22
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the plot. #
==============================================================================
# 1. REQUIRED PACKAGES & DEPENDENCIES #
==============================================================================
Define the complete list of packages used throughout your analysis
data_packages <- c( “tidyverse”, # For data manipulation, pipes, and ggplot2 “readxl”, # For reading the Excel sheets “readr”, # For text and type parsing (type_convert) “janitor”, # For handling mixed Excel date formats “writexl”, “ggplot2”) #for creating plots)
glmm_packages <- c(
“glmmTMB”, # For running the Zero-Inflated GLMM regressions “DHARMa”, # For overdispersion and residual diagnostic tests “emmeans”, “fitdistrplus”, “broom.mixed”, # For tidying GLMM outputs into dataframes “openxlsx” # For exporting results cleanly to Excel tables )
Automatically install missing packages
missing_packages <- data_packages[!(data_packages %in% installed.packages()[, “Package”])] if (length(missing_packages) > 0) { install.packages(missing_packages, dependencies = TRUE) }
missing_glmm_packages <- glmm_packages[!(glmm_packages %in% installed.packages()[, “Package”])] if (length(missing_glmm_packages) > 0) { install.packages(missing_glmm_packages, dependencies = TRUE) }
Load all packages silently
suppressPackageStartupMessages( lapply(data_packages, library, character.only = TRUE) )
suppressPackageStartupMessages( lapply(glmm_packages, library, character.only = TRUE) )
cat(“>>> All required packages successfully verified and loaded! <<<”)
Get all sheet names
raven_sheet_names <- excel_sheets(“Raven_GLMM.xlsx”)
Read and stack them together
raven_all_data <- raven_sheet_names %>% set_names() %>% map_dfr(~ read_excel(“Raven_GLMM.xlsx”, sheet = .x, col_types = “text”), .id = “Feeder”) %>% type.convert()
#Select only columns pertaining to raven treatment raven_all_data <- raven_all_data %>% dplyr::select(Feeder, Date, TarpY_N, RavenPics, DayOfYear, TimeSinceTarp, UntreatFeed)
#format columns into correct class of data, assign na to missing data raven_all_data<-raven_all_data%>%mutate(Date = as.Date(Date, origin = “1899-12-30”)) raven_all_data<-raven_all_data %>% mutate(TarpY_N = na_if(TarpY_N, “N/A”),TarpY_N = factor(TarpY_N)) raven_all_data<-raven_all_data %>% mutate(RavenPics= na_if(RavenPics, “N/A”), RavenPics = as.numeric (RavenPics)) raven_all_data<-raven_all_data %>% mutate(DayOfYear= na_if(DayOfYear, “N/A”), DayOfYear = as.numeric (DayOfYear)) raven_all_data <- raven_all_data %>% mutate( TimeSinceTarp = if_else(TimeSinceTarp %in% c(“N/A”, “365”), NA_character_, TimeSinceTarp), TimeSinceTarp = as.numeric(TimeSinceTarp)) raven_all_data<-raven_all_data %>% mutate(UntreatFeed= na_if(UntreatFeed, “N/A”), UntreatFeed = as.numeric (UntreatFeed))
summary(raven_all_data) glimpse(raven_all_data) write_xlsx(raven_all_data, “raven_all_data.xlsx”)
#remove blank rows raven_clean <- raven_all_data %>% filter(!(is.na(TarpY_N) & is.na(RavenPics))) write_xlsx(raven_clean, “raven_clean.xlsx”)
#linear regression for effect of tarp # 1. Drop missing rows for these specific columns tarp_binary_data <- na.omit(raven_clean[, c(“TarpY_N”, “RavenPics”, “Feeder”)])
2. Make sure TarpY_N is treated as a distinct group (factor)
tarp_binary_data\(TarpY_N <- as.factor(tarp_binary_data\)TarpY_N)
3. Run the GLMM model
TarpGroup_lm <- glmmTMB( RavenPics ~ TarpY_N + (1 | Feeder), data = tarp_binary_data, family = nbinom1 )
Look at the summary table to see if Tarps make a significant difference!
summary(TarpGroup_lm)
1. Get predictions on the link scale (ignoring feeder differences)
pred_binary <- predict(TarpGroup_lm, type = “link”, se.fit = TRUE, re.form = NA)
2. Convert to real count scales
tarp_binary_data\(PredictedPics <- exp(pred_binary\)fit) tarp_binary_data\(LowerCI <- exp(pred_binary\)fit - 1.96 * pred_binary\(se.fit) tarp_binary_data\)UpperCI <- exp(pred_binary\(fit + 1.96 * pred_binary\)se.fit)
ggplot(tarp_binary_data, aes(x = TarpY_N, y = RavenPics)) + # 1. Spread out raw data points sideways so you can see them clearly geom_jitter(alpha = 0.4, color = “darkgray”, width = 0.15) + # 2. Add error bars showing the 95% Confidence Interval for each group geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), color = “red”, width = 0.1, size = 1) + # 3. Add large target dots showing the exact predicted model averages geom_point(aes(y = PredictedPics), color = “red”, size = 3) + # 4. Connect the two group averages with your regression slope line geom_line(aes(y = PredictedPics, group = 1), color = “red”, size = 1, linetype = “dashed”) + labs( title = “Raven Pictures by Tarp Presence”, subtitle = “Model fit using GLMM with random effects”, x = “Tarp Status (TarpY_N)”, y = “Number of Raven Pictures” ) + theme_minimal()+ theme(panel.grid = element_blank())
#linear regression for looking at effect of day of year
1. Create a new dataset without missing values from your model variables
raven_DOY_clean <- raven_clean[raven_clean$TarpY_N == 0, ] DOY_lm_data <- na.omit(raven_DOY_clean[, c(“DayOfYear”, “RavenPics”, “Feeder”)])
DOY_lm_data\(ContinuousDay <- ifelse(DOY_lm_data\)DayOfYear >= 50, DOY_lm_data\(DayOfYear, DOY_lm_data\)DayOfYear + 365)
2. Re-run your model using ContinuousDay as the X variable
DayOfYear_lm <- glmmTMB(RavenPics ~ ContinuousDay + (1|Feeder), data = DOY_lm_data, family = nbinom1) summary(DayOfYear_lm)
Create a perfectly ordered sequence of days from 50 to 415
smooth_timeline <- data.frame(ContinuousDay = seq(50, 415, by = 1))
Predict onto this smooth timeline (ignoring feeder quirks)
pred_smooth <- predict(DayOfYear_lm, newdata = smooth_timeline, type = “link”, se.fit = TRUE, re.form = NA)
3. Calculate your smooth line and ribbon values
smooth_timeline\(PredictedPics <- exp(pred_smooth\)fit) smooth_timeline\(LowerCI <- exp(pred_smooth\)fit - 1.96 * pred_smooth\(se.fit) smooth_timeline\)UpperCI <- exp(pred_smooth\(fit + 1.96 * pred_smooth\)se.fit)
4. Plot using the raw data for points, but the smooth timeline for the line
library(ggplot2) ggplot() + # Use your original data ONLY for the grey scatterplot dots geom_point(data = DOY_lm_data, aes(x = ContinuousDay, y = RavenPics), alpha = 0.5, color = “darkgray”) + # Use the brand-new smooth dataset for the ribbon and the line geom_ribbon(data = smooth_timeline, aes(x = ContinuousDay, ymin = LowerCI, ymax = UpperCI), fill = “blue”, alpha = 0.15) + geom_line(data = smooth_timeline, aes(x = ContinuousDay, y = PredictedPics), color = “blue”, size = 1) + scale_x_continuous( breaks = c(50, 150, 250, 365, 415), labels = c(“50”, “150”, “250”, “365 / 0”, “50 (Next Year)”) ) + labs( title = “Raven Pictures Over the Year (Feb 19 Start Date; Feeders as Random Factor)”, x = “Day of the Year”, y = “Number of Raven Pictures” ) + theme_minimal()
linear regression for testing if tarp loses effectiveness as time passes
tarp_lm_data <- na.omit(raven_clean[, c(“TimeSinceTarp”, “RavenPics”)])
TimeSinceTarp_lm <- glmmTMB( RavenPics ~ TimeSinceTarp, data = tarp_lm_data, family = nbinom1 )
summary(TimeSinceTarp_lm)
1. Get predictions on the link scale (ignoring individual feeder quirks)
pred_tarp <- predict(TimeSinceTarp_lm, type = “link”, se.fit = TRUE, re.form = NA)
2. Convert to real raven count numbers and calculate boundaries
tarp_lm_data\(PredictedPics <- exp(pred_tarp\)fit) tarp_lm_data\(LowerCI <- exp(pred_tarp\)fit - 1.96 * pred_tarp\(se.fit) tarp_lm_data\)UpperCI <- exp(pred_tarp\(fit + 1.96 * pred_tarp\)se.fit)
library(ggplot2) ggplot(tarp_lm_data, aes(x = TimeSinceTarp, y = RavenPics)) + geom_point(alpha = 0.5, color = “darkgray”) + #geom_ribbon(aes(ymin = LowerCI, ymax = UpperCI), fill = “blue”, alpha = 0.15) + geom_line(aes(y = PredictedPics), color = “blue”, size = 1) + labs( title = “Time Since Tarp Placed vs Raven Pics”, x = “Time Since Tarp”, y = “Number of Raven Pictures” ) + theme_minimal()
install.packages(c(“sjPlot”, “performance”))
library(sjPlot)
Export the multiple regression results straight to a Word document
tab_model( raven_simple_nb, transform = “exp”, # Exponentiates log coefficients into clear Rate/Odds Ratios show.zeroinf = TRUE, # Explicitly includes the Zero-Inflation table structure p.style = “numeric_stars”, # Formats p-values neatly alongside significance stars title = “Table 1: Effect of Time of Tarp Treament on Raven Activity”, file = “TimeSinceTarp.doc” # Instantly saves this file to your computer )
#Plot linear model of effect of additional untreated feeders on raven pics
feed_lm_data <- na.omit(raven_clean[, c(“UntreatFeed”, “RavenPics”, “Feeder”)])
UntreatFeed_mixed_lm <- glmmTMB( RavenPics ~ UntreatFeed + (1|Feeder), data = feed_lm_data, family = nbinom1 )
summary(UntreatFeed_mixed_lm)
1. Get predictions on the link scale (ignoring individual feeder quirks)
pred_feed <- predict(UntreatFeed_mixed_lm, type = “link”, se.fit = TRUE, re.form = NA)
2. Convert to real raven counts and calculate boundaries
feed_lm_data\(PredictedPics <- exp(pred_feed\)fit) feed_lm_data\(LowerCI <- exp(pred_feed\)fit - 1.96 * pred_feed\(se.fit) feed_lm_data\)UpperCI <- exp(pred_feed\(fit + 1.96 * pred_feed\)se.fit)
ggplot(feed_lm_data, aes(x = UntreatFeed, y = RavenPics)) + # 1. Add the raw data points geom_point(alpha = 0.4, color = “darkgray”) + # 2. Add the 95% confidence ribbon geom_ribbon(aes(ymin = LowerCI, ymax = UpperCI), fill = “darkorange”, alpha = 0.15) + # 3. Add the model trend line geom_line(aes(y = PredictedPics), color = “darkorange”, size = 1) + # 4. Add clean titles and labels labs( title = “Raven Pictures Relative to Untreated Feeder Activity”, subtitle = “Model fit with a GLMM (Feeder controlled)”, x = “Untreated Feeder Value”, y = “Number of Raven Pictures” ) + theme_minimal()
#zero inflated negative binomial 2 (assumes variance increases quadratically as count increases) raven_glmm <- glmmTMB(RavenPics ~ TarpY_N + DayOfYear + UntreatFeed + (1|Feeder), data=raven_all_data, family=nbinom2(),ziformula = ~TarpY_N) summary(raven_glmm) library(DHARMa)
1. Simulate standardized residuals from your model
sim_res <- simulateResiduals(fittedModel = raven_glmm, n = 250)
2. Run the explicit dispersion test
testDispersion(sim_res)
3. Plot all diagnostics (QQ-plot and Residuals vs. Predicted)
plot(sim_res)
raven_glmm_final <- glmmTMB( RavenPics ~ TarpY_N + scale(DayOfYear) + UntreatFeed + (1|Feeder), data = raven_clean, # Use the version with raw integers family = nbinom2(), ziformula = ~TarpY_N)
summary(raven_glmm_final)
1. Simulate standardized residuals from your model
sim_res <- simulateResiduals(fittedModel = raven_glmm_final, n = 250)
2. Run the explicit dispersion test
testDispersion(sim_res)
3. Plot all diagnostics (QQ-plot and Residuals vs. Predicted)
plot(sim_res)
==============================================================================
FINAL SELECTED MODEL
==============================================================================
Chosen based on Likelihood Ratio Test and verified via DHARMa diagnostics.
Passed dispersion test with p-value = 0.456.
Switch family to nbinom1
raven_glmm_nb1 <- glmmTMB( RavenPics ~ TarpY_N + scale(DayOfYear) + UntreatFeed + (1 | Feeder), data = raven_clean, family = nbinom1, # Linear variance structure ziformula = ~ TarpY_N ) summary(raven_glmm_nb1)
Test the new model’s dispersion
library(DHARMa) sim_res_nb1 <- simulateResiduals(raven_glmm_nb1) testDispersion(sim_res_nb1)
Extract publication-ready estimates and 95% confidence intervals
cc <- confint(raven_glmm_nb1) print(cc)
library(sjPlot)
1. Force R globally to use scientific (E) notation for all output
options(scipen = -999)
2. Run your tab_model code with scientific p-value styling added
tab_model( raven_glmm_nb1, transform = “exp”, # Exponentiates log coefficients into clear Rate/Odds Ratios show.zeroinf = TRUE, # Explicitly includes the Zero-Inflation table structure digits = 2, # Controls how many decimal points appear in your E-notation (e.g., 4.96e+34) p.style = “scientific_stars”,# Swaps numeric p-values to scientific format alongside stars title = “Table 1: Effect of Time of Tarp Treatment on Raven Activity”, file = “RavenGLMM.doc” # Instantly saves this file to your computer )
tab_model( raven_glmm_nb1, transform = “exp”, # Exponentiates log coefficients into clear Rate/Odds Ratios show.zeroinf = TRUE, # Explicitly includes the Zero-Inflation table structure p.style = “numeric_stars”, # Formats p-values neatly alongside significance stars title = “Table 1: Effect of Time of Tarp Treament on Raven Activity”, file = “RavenGLMM.doc” # Instantly saves this file to your computer )
Export the multiple regression results straight to a Word document
install.packages(c(“broom.mixed”, “openxlsx”))
library(tidyverse) library(broom.mixed) library(openxlsx)
1. Extract the raw model results into a clean dataframe
raw_table <- tidy(raven_glmm_nb1, conf.int = TRUE)
2. Format, round, and make it publication-ready
clean_table <- raw_table %>% # Filter to just the fixed effects
we care about filter(component %in% c(“cond”, “zi”)) %>% mutate( #
Exponentiate the estimates and confidence intervals to get Ratios Ratio
= round(exp(estimate), 2), Lower_CI = round(exp(conf.low), 2), Upper_CI
= round(exp(conf.high), 2), # Round the p-value to 3 decimal places
p_value = round(p.value, 3), # Combine the CI columns into a clean
string format 95% CI = paste0(“[”, Lower_CI, ”, ”,
Upper_CI, ”]”) ) %>% # Clean up the row names to read nicely mutate(
Predictor = case_when( term == “(Intercept)” ~ “Baseline Intercept”,
term == “TarpY_N1” ~ “Tarp Present”, term == “scale(DayOfYear)” ~ “Day
of Year (Scaled)”, term == “UntreatFeed” ~ “Untreated Feeders”, TRUE ~
term ), Model_Component = if_else(component == “cond”, “Count Model
(Rate Ratio)”, “Zero-Inflation Model (Odds Ratio)”) ) %>% # Keep only
the columns you want in your final paper select(Model_Component,
Predictor, Ratio, 95% CI, p_value)
3. Save it directly to an Excel file on your computer
write.xlsx(clean_table, “Raven_Final_Table.xlsx”)
#family=poisson raven_glmm2 <- glmmTMB( RavenPics ~ TarpY_N +
scale(DayOfYear) + UntreatFeed + (1|Feeder),
data=raven_all_data, family=poisson(),ziformula = ~TarpY_N)
summary(raven_glmm2)
Compare the Poisson and Negative Binomial models side-by-side
anova(raven_glmm2, raven_glmm_nb1)
library(sjPlot)
Plot the conditional (count) part
plot_model(raven_glmm_nb1, type = “est”, component = “cond”)
Plot the zero-inflation part (probability of extra zeros)
plot_model(raven_glmm_nb1, type = “est”, component = “zi”)
library(ggeffects) library(ggplot2)
1. Plot effect of Tarp being present on the number of Raven pics
pred_tarp <- ggpredict(raven_glmm_nb1, terms = “TarpY_N”) plot(pred_tarp) + labs(title = “Predicted Raven Pics by Tarp Presence”)
1. Generate bias-corrected predictions
pred_tarp <- ggpredict(raven_glmm_nb1, terms = “TarpY_N”, bias_correction = TRUE)
2. Plot a publication-ready bar graph
ggplot(pred_tarp, aes(x = x, y = predicted, fill = x)) + geom_bar(stat = “identity”, width = 0.5, color = “black”, show.legend = FALSE) + geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) + scale_fill_manual(values = c(“0” = “gray70”, “1” = “steelblue”)) + scale_x_discrete(labels = c(“0” = “No Tarp”, “1” = “Tarp Present”)) + labs( x = “Deterrent Treatment”, y = “Predicted Raven Pictures (Count)”, title = “Effect of Tarp Deterrent on Raven Activity” ) + theme_classic(base_size = 14)
2. Plot effect of Temperature (continuous variable)
pred_temp <- ggpredict(raven_glmm_nb1, terms = “MaxTemp”) plot(pred_temp) + labs(title = “Effect of Max Temp on Raven Activity”)
pred_day <- ggpredict(raven_glmm_nb1, terms=“DayOfYear”) plot(pred_day) + labs(title = “Effect of DayOfYear on Raven Activity”)
pred_Untreatfeed <- ggpredict(raven_glmm_nb1, terms = “UntreatFeed”) plot(pred_Untreatfeed) + labs(title = “Effect of Additional Untreated Feeders on Raven Activity”)