Analysis for Effects of Cormoran Abundance to Avian Breeding Communities

Supplementary Code 3: Analysis for Effects of Cormoran Abundance to Avian Breeding Communities Aare Verliin, Richard Meitern, et al. 2025-11-21 Introduction The below code does statistic analysis using the abundance data of island bird communities. Data Preparation The dataset is the same as in the file “02. . . ” and is loaded from the same files. df_wide <- readRDS(paste0(inputFname, "_wide_fixed.rds")) df_long <- readRDS(paste0(inputFname, "_long_fixed.rds")) Because data before 1960 is scarce, we are using only data since 1960. #using only data since 1960 see the file 02 for plots and more info df_wide <- df_wide[df_wide$year >1959, ] df_long <- df_long[df_long$year >1959, ] Also we are calculating the years from cormorant invasion for each site. #get first year of cormoran first_occurance_cormoran <- aggregate(year ~adjusted_name, df_wide[df_wide$Phalacrocorax_carbo >0, ], FUN=min) cormoran_start <- first_occurance_cormoran$year names(cormoran_start) <- first_occurance_cormoran$adjusted_name rm(first_occurance_cormoran) df_wide$cormoran_start <- cormoran_start[df_wide$adjusted_name] df_wide$years_from_cormoran <- df_wide$year -df_wide$cormoran_start df_wide$cormoran_invaded <- TRUE df_wide[is.na(df_wide$years_from_cormoran),"cormoran_invaded"]<- FALSE df_long <- merge(df_long, df_wide[c("years_from_cormoran","adjusted_name","year")], by=c("adjusted_name","year")) 1 countVars <- c( "Habitat", "Habitat2", "Habitat.Density", "Nest.associated", "Primary.Lifestyle", "Trophic.Niche", "Migration", "Trophic.Level" ) #convert to factors numberVars <- c( "Mass", "Range.Size" ) df_long[countVars] <- lapply(df_long[countVars], factor) df_long$sp_name <- factor(df_long$sp_name) df_long$cormorans <- df_long$Phalacrocorax_carbo df_long$cormoran_present <- factor(df_long$Phalacrocorax_carbo >0,levels = c("No","Yes")) df_long$cormoran_level <- cut(df_long$Phalacrocorax_carbo, breaks = c(0,1,100,Inf), include.lowest = TRUE, labels = c("No","Low (1-10)","High (100+)")) df_long$cormoran_level_ha <- cut(df_long$Phalacrocorax_carbo/df_long$ha, breaks = c(0,1,10,100,Inf), include.lowest = TRUE, labels = c("No","Low (1-10)","Moderate (10-100)","High (100+)")) df_long$place <- factor(df_long$adjusted_name) df_long$year_f <- factor(df_long$year, ordered = T) #get mean counts for each species place before cormoran invasion form1 <- formula(paste0("cbind(mean_year_density=count/ha) ~ ", "sp_name + place + ha + years_from_cormoran")) desity_means <- aggregate(form1, df_long, FUN = mean, na.rm =TRUE) #take 10 years form2 <- formula(paste0("cbind(mean_before_cormoran_density=mean_year_density) ~ ", "sp_name + place")) means_before_cormoran <- aggregate(form2, data=desity_means[desity_means$years_from_cormoran <1& desity_means$years_from_cormoran > -9, ], FUN=mean) df_long <- merge(df_long, means_before_cormoran, by=c("sp_name","place"), all.x = T) 2 #filter data df <- df_long[!df_long$cormoran, ] #keep only open islands df <- df[df$island == 1, ] df <- df[!is.na(df$island), ] In the following section we are going to look at the bird community changes due to cormorant invasion. For this task only islands are included in the analysis. df_is <- df species <- aggregate(cbind(sp_names=as.character(sp_name)) ~adjusted_sp_name, df_is, FUN = function(x){paste0(unique(x), collapse=", ")}) species2 <- aggregate(cbind(count=count) ~adjusted_sp_name, df_is, FUN = sum) species3 <- aggregate(cbind(year_count=year) ~adjusted_sp_name, df_is, FUN = function(x){length(unique(x))}) species4 <- aggregate(cbind(years=year) ~adjusted_sp_name, df_is, FUN = function(x){paste0(unique(x), collapse = ",")}) species5 <- aggregate(cbind(place_count=adjusted_name) ~adjusted_sp_name, df_is, FUN = function(x){length(unique(x))}) species6 <- aggregate(cbind(places=adjusted_name) ~adjusted_sp_name, df_is, FUN = function(x){paste0(unique(x), collapse = ",")}) species <- merge(species, species2, by="adjusted_sp_name") species <- Reduce(function(x, y) merge(x, y, by="adjusted_sp_name",all=TRUE), list(species, species3, species4, species5, species6)) rm(df_is, species2, species3, species4, species5, species6) writexl::write_xlsx(species, paste0(rdir, "species_info_islands_only.xlsx")) Statistical Analysis We are going to fit a series of models to the data to see how the bird community changes due to cormorant invasion. We are going to use the lmerTest package to fit the models. The dependent variable is the normalized density of the species. This calculated by dividing the count of the species by the area of the island an then divide by the 10 year mean density before first occurrence of cormorants. Data Since Cormorant Invasion The predictor is current cormorant density. We are going to look all the years where a place has been a part of breeding pair monitoring including only measurements since the introduction of cormorant. df$cormoran_density <- df$cormorans/df$ha df$density <- df$count /df$ha df$density_norm <- df$density /df$mean_before_cormoran_density dfl <- df[df$years_from_cormoran > -1, ] baseModel <- lmerTest::lmer(log10(density_norm) ~log10(cormoran_density +1)+ (1|place) +(1|sp_name), data = dfl) summary(baseModel) 3 ## Linear mixed model fit by REML. t-tests use Satterthwaite’s method [ ## lmerModLmerTest] ## Formula: log10(density_norm) ~ log10(cormoran_density + 1) + (1 | place) + ## (1 | sp_name) ## Data: dfl ## ## REML criterion at convergence: 6776.7 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -4.9079 -0.5030 -0.0174 0.5096 4.1216 ## ## Random effects: ## Groups Name Variance Std.Dev. ## place (Intercept) 0.01822 0.1350 ## sp_name (Intercept) 0.03696 0.1922 ## Residual 0.24575 0.4957 ## Number of obs: 4591, groups: place, 62; sp_name, 56 ## ## Fixed effects: ## Estimate Std. Error df t value Pr(>|t|) ## (Intercept) -6.008e-02 3.907e-02 1.005e+02 -1.538 0.1272 ## log10(cormoran_density + 1) -1.282e-02 7.786e-03 3.040e+03 -1.647 0.0997 . ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 ## ## Correlation of Fixed Effects: ## (Intr) ## lg10(cr_+1) -0.258 dharma_diagnostics(baseModel) 4 0.0 0.4 0.8 0.0 0.2 0.4 0.6 0.8 1.0 QQ plot residuals Expected Observed KS test: p= 0 Deviation significant Outlier test: p= 0 Deviation significant Dispersion test: p= 0.088 Deviation n.s. Model predictions (rank transformed) DHARMa residual 0.0 0.4 0.8 0.00 0.25 0.50 0.75 1.00 DHARMa residual vs. predicted DHARMa residual # Compute R-squared r2_results <- performance::r2(baseModel) print(r2_results) ## # R2 for Mixed Models ## ## Conditional R2: 0.184 ## Marginal R2: 0.001 The base model above indicates that the cormorant density has a significant negative effect on the bird community. The model explains 18% of the variance in the normalized density of the species. source("./R/lmer_utils.R") Let’s look at the results for different bird guilds. We are going to fit a model for each guild separately. The formula used for these models is: log10(density_norm)∼log10(cormoran_density + 1) + (1|place) + (1|sp_name) This formula models the log-transformed normalized density of birds as a function of the log-transformed density of cormorants, while accounting for random effects associated with the place and bird species name. For each level of the guild variable, the model is fit separately using the data subset corresponding to that guild. guildModels <- list() for(var in countVars){ levels <- levels(dfl[[var]]) 5 varModels <- lapply(levels, function(g){ form <- log10(density_norm) ~log10(cormoran_density +1)+(1|place) +(1|sp_name) data_subset <- dfl[dfl[[var]] == g, ] # Fit the model safely model <- fit_lmer_safely(form, data_subset) return(model) }) names(varModels) <- paste0(var, "-",levels) guildModels <- c(guildModels, varModels) } # Use the generalized function for guildModels guild_results_df <- process_model_list(guildModels, signif = 3,adj_method = "bonferroni") knitr::kable(guild_results_df, digits = 3,row.names = F, caption = paste0("Cormoran effect on bird guilds (p-values", "adjusted with Bonferroni correction) Model", "with species and location as random factor")) Table 1: Cormoran effect on bird guilds (p-valuesadjusted with Bonferroni correction) Modelwith species and location as random factor model estimate std_error t_value p_value p_stars adj_p adj_p_stars Habitat-Coastal -0.026 0.011 -2.380 0.017 * 0.470 Habitat-Human Modified -0.030 0.010 -2.950 0.003 ** 0.092 . Habitat-Marine 0.088 0.027 3.290 0.001 ** 0.030 * Habitat-Shrubland -0.053 0.045 -1.180 0.260 1.000 Habitat-Wetland -0.010 0.014 -0.715 0.470 1.000 Habitat2-building 0.029 0.095 0.308 0.760 1.000 Habitat2-open -0.013 0.008 -1.590 0.110 1.000 Habitat2-wood -0.004 0.014 -0.308 0.760 1.000 Habitat.Density-Open -0.013 0.008 -1.670 0.094 . 1.000 Nest.associated-buildings 0.027 0.038 0.710 0.480 1.000 Nest.associated-ground -0.012 0.008 -1.520 0.130 1.000 Nest.associated-trees -0.036 0.017 -2.080 0.038 * 1.000 Primary.Lifestyle-Aerial -0.087 0.032 -2.740 0.006 ** 0.170 Primary.Lifestyle-Aquatic -0.007 0.013 -0.508 0.610 1.000 Primary.Lifestyle-Generalist 0.021 0.017 1.230 0.220 1.000 Primary.Lifestyle-Terrestrial -0.012 0.009 -1.390 0.170 1.000 Trophic.Niche-Aquatic predator -0.008 0.011 -0.695 0.490 1.000 Trophic.Niche-Herbivore aquatic -0.029 0.016 -1.880 0.061 . 1.000 Trophic.Niche-Herbivore terrestrial 0.005 0.034 0.137 0.890 1.000 Trophic.Niche-Invertivore -0.033 0.011 -2.860 0.004 ** 0.120 Trophic.Niche-Omnivore -0.001 0.016 -0.043 0.970 1.000 Migration-Migratory -0.038 0.012 -3.180 0.002 ** 0.041 * Migration-Partially migratory 0.018 0.011 1.670 0.095 . 1.000 Migration-Resident -0.032 0.010 -3.150 0.002 ** 0.047 * 6 model estimate std_error t_value p_value p_stars adj_p adj_p_stars Trophic.Level-Carnivore -0.018 0.010 -1.880 0.061 . 1.000 Trophic.Level-Herbivore -0.020 0.014 -1.450 0.150 1.000 Trophic.Level-Omnivore 0.032 0.018 1.790 0.074 . 1.000 We also analyze the results for each bird species separately. For this, we fit a model for each species using the following formula: log10(density_norm)∼log10(cormoran_density + 1) + (1|place) This formula models the log-transformed normalized density of each species as a function of the logtransformed density of cormorants, with random effects to account for the variation associated with different places. For each species, the model is fit on a subset of the data corresponding to that species. Models with warnings or errors are excluded to ensure robust results. # Get unique species names species <- unique(dfl$sp_name) # Fit models for each species, excluding those with warnings or errors sp_models <- lapply(species, function(sp) { form <- log10(density_norm) ~log10(cormoran_density +1)+(1|place) data_subset <- dfl[dfl$sp_name == sp, ] # Fit the model safely model <- fit_lmer_safely(form, data_subset) return(model) }) # Assign species names to the list names(sp_models) <- species # Use the generalized function for sp_models sp_results_df <- process_model_list(sp_models, signif = 3,adj_method = "bonferroni") knitr::kable(sp_results_df, digits = 3,row.names = F, caption = paste0("Cormoran effect on bird species (p-values", "adjusted with Bonferroni correction) Model with", "location as random factor")) Table 2: Cormoran effect on bird species (p-valuesadjusted with Bonferroni correction) Model withlocation as random factor model estimate std_error t_value p_value p_stars adj_p adj_p_stars Acrocephalus schoenobaenus -0.058 0.114 -0.508 0.62 1 Anas platyrhynchos -0.109 0.022 -4.900 <0.001 *** <0.001 *** Anser anser 0.006 0.035 0.172 0.86 1 Arenaria interpres -0.015 0.022 -0.686 0.5 1 Aythya fuligula -0.107 0.047 -2.280 0.025 * 0.8 Branta leucopsis -0.029 0.094 -0.307 0.77 1 7 model estimate std_error t_value p_value p_stars adj_p adj_p_stars Charadrius hiaticula -0.018 0.017 -1.030 0.31 1 Chroicocephalus ridibundus -0.016 0.091 -0.175 0.86 1 Corvus corone 0.008 0.013 0.640 0.53 1 Curruca curruca -0.048 0.078 -0.621 0.56 1 Cygnus olor -0.001 0.016 -0.044 0.97 1 Haematopus ostralegus -0.036 0.010 -3.600 <0.001 *** 0.011 * Hirundo rustica -0.072 0.116 -0.621 0.54 1 Hydroprogne caspia 0.056 0.057 0.982 0.33 1 Larus argentatus 0.015 0.015 1.050 0.29 1 Larus canus -0.103 0.024 -4.390 <0.001 *** <0.001 *** Larus fuscus -0.085 0.064 -1.340 0.19 1 Larus marinus 0.073 0.018 3.960 <0.001 *** 0.0029 ** Mareca strepera -0.054 0.033 -1.640 0.11 1 Melanitta fusca 0.016 0.041 0.399 0.69 1 Mergus merganser -0.032 0.022 -1.430 0.15 1 Mergus serrator -0.062 0.033 -1.890 0.063 . 1 Motacilla alba -0.037 0.012 -3.220 0.0015 ** 0.048 * Motacilla flava -0.035 0.054 -0.663 0.52 1 Oenanthe oenanthe -0.196 0.031 -6.280 0.0028 ** 0.09 . Podiceps cristatus 0.098 0.269 0.363 0.73 1 Somateria mollissima 0.049 0.024 2.010 0.045 * 1 Spatula clypeata -0.020 0.047 -0.420 0.68 1 Sterna par s hir -0.197 0.034 -5.700 <0.001 *** <0.001 *** Tadorna tadorna 0.008 0.036 0.228 0.82 1 Thalasseus sandvicensis 0.341 0.126 2.720 0.015 * 0.49 Tringa totanus -0.038 0.015 -2.500 0.014 * 0.45 Data Before and After Cormorant Invasion Do increase power we can add include more measurements from the time before cormorant invasion. However it is dabatable if in such a dataset it would be correct to analyse the effect of cormorants. We will add up to 10 years before cormorant invasion to the dataset. And present similar analysis as above. dfl <- df[df$years_from_cormoran > -11, ] baseModel <- lmerTest::lmer(log10(density_norm) ~log10(cormoran_density +1)+ (1|place) +(1|sp_name), data = dfl) summary(baseModel) ## Linear mixed model fit by REML. t-tests use Satterthwaite’s method [ ## lmerModLmerTest] ## Formula: log10(density_norm) ~ log10(cormoran_density + 1) + (1 | place) + ## (1 | sp_name) ## Data: dfl ## ## REML criterion at convergence: 9554 ## ## Scaled residuals: 8 ## Min 1Q Median 3Q Max ## -5.7713 -0.4452 0.0165 0.4907 5.0934 ## ## Random effects: ## Groups Name Variance Std.Dev. ## sp_name (Intercept) 0.011259 0.10611 ## place (Intercept) 0.009285 0.09636 ## Residual 0.184484 0.42952 ## Number of obs: 8146, groups: sp_name, 64; place, 62 ## ## Fixed effects: ## Estimate Std. Error df t value Pr(>|t|) ## (Intercept) -3.963e-02 2.238e-02 1.080e+02 -1.771 0.07940 ## log10(cormoran_density + 1) -1.838e-02 5.235e-03 6.713e+03 -3.511 0.00045 ## ## (Intercept) . ## log10(cormoran_density + 1) *** ## --- ## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 ## ## Correlation of Fixed Effects: ## (Intr) ## lg10(cr_+1) -0.146 dharma_diagnostics(baseModel) 0.0 0.4 0.8 0.0 0.2 0.4 0.6 0.8 1.0 QQ plot residuals Expected Observed KS test: p= 0 Deviation significant Outlier test: p= 0 Deviation significant Dispersion test: p= 0.2 Deviation n.s. Model predictions (rank transformed) DHARMa residual 0.0 0.4 0.8 0.00 0.25 0.50 0.75 1.00 DHARMa residual vs. predicted DHARMa residual 9