Replicating the Aquaculture Impact Analysis on Minnow Communities
Authors
Kylie Newcomer
Joaquin Sandoval
Grace Bryan
Published
March 20, 2026
Photo from Canadian Aquaculture Industry Alliance
Introduction
Freshwater aquaculture can influence lake ecosystems through nutrient enrichment, increased primary productivity, and shifts in food web structure. These ecosystem-level changes may cascade to affect fish communities, including small-bodied species such as minnows.
This report replicates and extends part of the analysis presented in “Impacts of freshwater aquaculture on fish communities: A whole-ecosystem experimental approach” by Rennie et al. (2019), which examined ecosystem responses to an experimental rainbow trout aquaculture operation conducted at the Experimental Lakes Area in Ontario, Canada. While the original study evaluates broad ecosystem responses, including nutrient dynamics and food web interactions, the present analysis focuses specifically on whether aquaculture activity influenced the abundance of minnows measured as catch-per-unit-effort (CPUE).
The study uses a Before-After Control-Impact (BACI) design comparing a treated lake, where aquaculture was implemented, with a nearby control lake. Minnow CPUE is observed across three experimental periods: before aquaculture began, during the period of active aquaculture operations, and after operations ended.
While BACI designs are widely used in ecological impact studies, they rely on the assumption that the control site provides a valid counterfactual for the treated site. In this replication, we extend the standard BACI framework by explicitly evaluating this assumption using a Difference-in-Differences (DiD) analysis. This approach allows for a more transparent assessment of whether observed changes in minnow abundance can plausibly be attributed to aquaculture activity. Because the original study focuses on ecosystem-level responses, changes in minnow abundance should be interpreted as one component of broader ecological dynamics rather than a complete measure of aquaculture impacts.
Data Preparation
The dataset contains fish sampling observations collected from lakes within the Experimental Lakes Area. Each observation corresponds to a specific lake, year, and sampling season, and includes CPUE measurements for several fish species.
Two lakes are used in this analysis:
Lake 375, which served as the treatment lake where aquaculture operations were introduced
Lake 373, which served as a control lake
Minnow CPUE values are aggregated by lake, year, and season to produce an average abundance estimate for each sampling period. Because abundance data are often right-skewed, the CPUE values are square-root transformed prior to analysis to stabilize variance and better satisfy model assumptions.
The following code imports the dataset and constructs the variables used in the analysis.
Data Wrangling
Key variables include:
treat - This variable was our predictor variable, assigned as a binary value indicating whether the lake experience aquaculture
period - A categorical value indicating the period of time within the study: before, during, and after aquaculture
season - Another catagorical variable used to differentiate sampling season between spring and fall
The combination of these variables allows us to analyze how treatment during different periods of the study influenced minnow CPUE. Additionally, separating the regressions by seasons allows for seasonal patterns to be present and not mask overall results.
Code
library(tidyverse)library(here)library(janitor)library(performance)library(jtools)library(gt)library(gtsummary)library(interactions)library(car)library(patchwork)minnows <-read_csv(here::here("posts", "eds241-final-project", "data", "Minnows_Aquaculture_Rennie19.csv")) %>% janitor::clean_names()aquaculture <- minnows %>%# Adding fish sp names to ambiguous columns rename(redbell_dace = f182,finescale_dace = f183,fathead_minnow = f209,pearl_dace = f214,slimy_sculpin = f382 ) %>%# Drop na values from minnow column drop_na(minnows) %>%# Adding treatment columnmutate(treat =case_when(lake ==375~1, TRUE~0),# Add period column period =case_when( year <2003~"Before", year <=2007~"During",TRUE~"After"),# Reordering levels; default is alphabetical period =fct_relevel(period, "Before", "During", "After"), # Renaming season row names to be cleaner season =case_when(season =="aSpring"~"Spring", season =="zFall"~"Fall"),# Always order from Fall to Spring season =fct_relevel(season, "Fall", "Spring"),sqrt_cpue =sqrt(minnows))# Mean cpue per lake/year/season, sqrt-transformedmean_cpue <- aquaculture %>%group_by(treat, year, season, period) %>%summarise(mean_cpue =mean(minnows),sqrt_cpue =mean(sqrt_cpue),.groups ="drop" )# Prepare plot dataplot_data <- aquaculture %>%group_by(treat, year, season) %>%summarise(mean_cpue =mean(minnows),se_cpue =sd(minnows) /sqrt(n()),.groups ="drop" )
Exploratory Time Series
Before estimating statistical models, it is useful to visualize how the minnow abundance changes over time in both lakes. This provides a preliminary view of whether the two lakes appear to respond differently during the aquaculture period.
Code
ggplot(plot_data, aes(x = year, y = mean_cpue, color =factor(treat))) +annotate("rect", xmin =2003, xmax =2007, # Highlight aquaculture periodymin =-Inf, ymax =Inf, fill ="grey80", alpha =0.5) +annotate("text", x =1996, label ="Aquaculture period",y =400, size =3) +annotate("rect", ymin =399, ymax =401, xmin =2001, xmax =2004) +geom_line(aes(linetype =factor(treat))) +geom_point(size =2) +geom_errorbar(aes(ymin = mean_cpue - se_cpue,ymax = mean_cpue + se_cpue), width =0.3) +scale_linetype_manual(values =c("dashed", "solid"),labels =c("Control", "Treated")) +scale_color_manual(values =c("black", "red"),labels =c("Control", "Treated")) +facet_wrap(~ season) +labs(x ="Year", y ="Average catch per unit effort",linetype ="", color ="",caption ="Figure 1: Time series of average minnow CPUE for the treatment and control lakes.") +theme_classic()
In spring sampling, the treated and control lakes follow broadly similar trajectories over time, suggesting limited divergence associated with aquaculture. In contrast, fall observations exhibit greater variability with the treated lake showing more pronounced fluctuations during the aquaculture period.
These patterns suggest that potential treatment effects may be season-dependent and motivate separate analyses for spring and fall observations. They also highlight the importance of formally evaluating whether post-treatment differences reflect aquaculture impacts or pre-existing difference between lakes.
BACI Analysis
A standard Before-After Control-Impact (BACI) analysis evaluates whether the treated site changes differently over time relative to a control site. Statistically, this is implemented using an interaction between treatment status and time period. The interaction term captures whether the change in the treated lake differs from the change observed in the control lake following the intervention.
In this framework, a significant interaction term is often interpreted as evidence of a treatment effect. However, this interpretation depends on the assumption that the control lake provides a valid counterfactual trajectory for the treated lake.
Spring BACI Model
Code
spring <- mean_cpue %>%filter(season =="Spring") %>%# Select only spring samplingmutate(treat =as.factor(treat), period =as.factor(period))mod_spring <-aov(sqrt_cpue ~ treat * period, data = spring)summary(mod_spring)
Df Sum Sq Mean Sq F value Pr(>F)
treat 1 0.03 0.029 0.009 0.924
period 2 13.81 6.906 2.238 0.123
treat:period 2 0.72 0.361 0.117 0.890
Residuals 33 101.82 3.085
The spring BACI model shows no statistically significant effects of treatment, period, or their interaction. This indicates that, during spring sampling, minnow abundance in the treated lake does not change in a way that differs meaningfully from the control lake. Post-hoc comparisons from the Tukey test also fail to detect significant differences across periods. From an ecological perspective, these results provide little evidence that aquaculture activity influenced minnow abundance during the spring season.
Fall BACI Model
Code
fall <- mean_cpue %>%filter(season =="Fall") %>%# Select only fall samplingmutate(treat =as.factor(treat), period =as.factor(period))mod_fall <-aov(sqrt_cpue ~ treat * period, data = fall)summary(mod_fall)
Df Sum Sq Mean Sq F value Pr(>F)
treat 1 65.06 65.06 6.930 0.01280 *
period 2 102.46 51.23 5.457 0.00896 **
treat:period 2 79.63 39.82 4.241 0.02295 *
Residuals 33 309.81 9.39
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In contrast, the fall BACI model yields statistically significant effects for treatment, period, and their interaction. The significant interaction indicates that minnow abundance in the treated lake changed differently over time relative to the control lake. Post-hoc comparisons suggest that this divergence is driven primarily by increases in CPUE during the aquaculture period in the treated lake. This pattern is consistent with a potential response of minnow populations to aquaculture-related changes in the lake environment, such as increased nutrient availability or altered food resources. However, this interpretation depends critically on the assumption that the treated and control lakes were following similar trajectories prior to treatment.
The fall ANOVA results can be summarized as follows:
treat: The treated lake had significantly higher fall minnow CPUE overall (treated lake was 2.60 sqrt CPUE units higher than control, p = 0.013).
period: Fall minnow abundance changed significantly over time across both lakes. During was significantly higher than Before (p = 0.006); After did not differ significantly from Before (p = 0.439).
treat:period: Lake 375 responded differently than Lake 373 across periods. The treated lake during aquaculture was significantly higher than: the control lake before aquaculture (1:During-0:Before, p < 0.001), the treated lake before aquaculture (1:During-1:Before, p = 0.002), and the control lake during aquaculture (1:During-0:During, p = 0.009). The comparison 0:After-1:During was also significant (p = 0.012), indicating a decline in the treated lake once operations ended.
Limitations of BACI
One limitation of the traditional BACI framework is that it assumes the control site provides a valid counterfactual for the treated site. In other words, the design assumes that the treated and control lakes would have followed similar trajectories in the absence of treatment. If this assumption is violated, for example, if the lakes were already trending differently before aquaculture began, then the estimated BACI interaction could reflect those pre-existing differences rather than the effect of the intervention. Difference-in-Differences approaches make this assumption explicit through the parallel trends assumption, which requires that treated and control units exhibit similar trends prior to treatment.
DiD: Testing Assumptions
In order to complete a Difference-in-differences analysis, a few assumptions must be met - No spillover: Because these lakes are separated, there is no concern of aquaculture in the treated lake affecting the control lake. - Parallel trends: The control and treatment lakes must have a similar trend of minnow CPUE prior to treatment.
Parallel Trends Test
Restricting to the Before period only, we test whether the two lakes were already trending differently prior to aquaculture. A significant treat:year interaction would indicate diverging pre-treatment trends, violating the parallel trends assumption.
Call:
lm(formula = sqrt_cpue ~ treat * factor(year) + season, data = pre_only)
Residuals:
Min 1Q Median 3Q Max
-4.6090 -0.6913 0.0000 0.6913 4.6090
Coefficients: (4 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.60869 1.73399 3.811 0.00154 **
treat 2.83381 2.38315 1.189 0.25175
factor(year)1989 -2.18335 2.38315 -0.916 0.37318
factor(year)1990 -3.78133 2.38315 -1.587 0.13214
factor(year)1991 -3.87029 2.38315 -1.624 0.12390
factor(year)1992 -4.37006 2.38315 -1.834 0.08536 .
factor(year)1993 -5.01627 2.38315 -2.105 0.05146 .
factor(year)1994 -0.70081 2.38315 -0.294 0.77248
factor(year)1995 -1.15006 2.38315 -0.483 0.63593
factor(year)1996 -5.63358 2.94722 -1.911 0.07402 .
factor(year)1998 0.54971 2.38315 0.231 0.82050
factor(year)1999 -5.35024 2.94722 -1.815 0.08826 .
factor(year)2001 -5.51362 2.94722 -1.871 0.07978 .
factor(year)2002 -5.09343 2.38315 -2.137 0.04836 *
seasonSpring -0.34266 0.81741 -0.419 0.68064
treat:factor(year)1989 -1.90928 3.37028 -0.567 0.57891
treat:factor(year)1990 -3.28652 3.37028 -0.975 0.34400
treat:factor(year)1991 -0.48383 3.37028 -0.144 0.88764
treat:factor(year)1992 -1.31633 3.37028 -0.391 0.70127
treat:factor(year)1993 -2.24063 3.37028 -0.665 0.51563
treat:factor(year)1994 NA NA NA NA
treat:factor(year)1995 NA NA NA NA
treat:factor(year)1996 0.59276 4.20789 0.141 0.88973
treat:factor(year)1998 NA NA NA NA
treat:factor(year)1999 -0.02259 3.79019 -0.006 0.99532
treat:factor(year)2001 -1.59834 4.12773 -0.387 0.70369
treat:factor(year)2002 NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.383 on 16 degrees of freedom
Multiple R-squared: 0.5844, Adjusted R-squared: 0.013
F-statistic: 1.023 on 22 and 16 DF, p-value: 0.4908
Code
pre_only %>%ggplot(aes(x = year, y = mean_cpue,color =factor(treat), linetype =factor(treat))) +geom_line() +geom_point() +scale_color_manual(values =c("black", "red"),labels =c("Control", "Treated")) +scale_linetype_manual(values =c("dashed", "solid"),labels =c("Control", "Treated")) +facet_wrap(~ season, scales ="free_y") +labs(x ="Year", y ="Mean cpue (fish per net days)",color ="Lake", linetype ="Lake",title ="Pre-treatment trends (parallel trends check)",caption ="Figure 2: Time series of pre-treatment CPUE trends between the treated and control lakes.") +theme_classic()
The coefficient of interest is the treatment-by-year interaction, which tests whether the treated lake was already trending differently before aquaculture began. In this context, the parallel trends assumption states that, in the absence of treatment, the treated lake would have followed the same trajectory as the control lake. This assumption underlies that causal interpretation of both BACI and Difference-in-Differences analysis.
The interaction between treatment and year is not statistically significant. This indicates that we fail to reject the parallel trends assumption. The accompanying time series plot (Figure 2) provides a visual check of this assumption, as the treated and control lakes appear to move in broadly similar directions before aquaculture begins. Taken together, these results provide some support for the use of the control lake as a counterfactual.
Difference-in-Differences Analysis
Following the evaluation of parallel trends, treatment effects are estimated using a Difference-in-Differences framework. Models are estimated separately by season to account for seasonal variation in ecological responses. Unlike BACI, which compares average differences before and after treatment, the DiD approach explicitly models changes over time between treated and control units. The interaction between treatment and period represents the estimated treatment effect, capturing whether the treated lake experienced a different change in CPUE relative to the control. Separating the analysis into Before-During and Before-During-After comparisons allows us to distinguish between immediate effects of aquaculture and potential legacy effects after operations cease.
Before and During
Fall
Code
pre_aqua_fall <- aquaculture %>%filter(period !="After", season =="Fall")mod_did_pre_fall <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_fall)summary(mod_did_pre_fall)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_fall)
Residuals:
Min 1Q Median 3Q Max
-13.5857 -2.0914 -0.4574 0.9579 23.7357
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.99967 3.57848 0.559 0.5777
treat 2.47479 1.93544 1.279 0.2042
periodDuring 4.90326 4.01048 1.223 0.2246
factor(year)1989 2.09314 4.59219 0.456 0.6496
factor(year)1990 -0.73068 5.10169 -0.143 0.8864
factor(year)1991 1.05524 5.10169 0.207 0.8366
factor(year)1992 -0.96435 4.48829 -0.215 0.8304
factor(year)1993 -2.02694 4.26968 -0.475 0.6361
factor(year)1994 6.14593 4.73388 1.298 0.1974
factor(year)1995 3.52301 7.15696 0.492 0.6237
factor(year)1996 -0.07277 5.97994 -0.012 0.9903
factor(year)1998 7.59200 7.15696 1.061 0.2916
factor(year)1999 0.60837 4.32884 0.141 0.8885
factor(year)2002 0.10463 4.30628 0.024 0.9807
factor(year)2003 -6.95595 2.80880 -2.476 0.0151 *
factor(year)2004 -5.84793 2.43429 -2.402 0.0183 *
factor(year)2005 -0.75653 2.70196 -0.280 0.7801
factor(year)2006 -3.26386 2.24396 -1.455 0.1492
factor(year)2007 NA NA NA NA
treat:periodDuring 6.11921 2.55321 2.397 0.0186 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.198 on 92 degrees of freedom
Multiple R-squared: 0.3891, Adjusted R-squared: 0.2695
F-statistic: 3.255 on 18 and 92 DF, p-value: 0.000105
Code
par(mfrow =c(1, 2))plot(mod_did_pre_fall, which =1) # Residuals vs fittedplot(mod_did_pre_fall, which =2) # QQ-plot
lag Autocorrelation D-W Statistic p-value
1 -0.1203524 2.230921 0.964
Alternative hypothesis: rho != 0
Diagnostics: The residuals vs. fitted plot shows increasing spread at higher fitted values, driven by the extreme CPUE spike during aquaculture — the sqrt transformation did not fully stabilize variance. The QQ plot confirms a heavy right tail from these outlier years. The Durbin-Watson test indicates no autocorrelation (DW = 2.23, p = 0.984).
The interaction coefficient in this model is positive and marginally significant, suggesting that minnow CPUE increased in the treated lake during aquaculture operations relative to the control lake.
Spring
Code
pre_aqua_spring <- aquaculture %>%filter(period !="After", season =="Spring")mod_did_pre_spring <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_spring)summary(mod_did_pre_spring)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_spring)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.0945 -0.4564 0.9555 6.0714
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.87506 1.34479 8.087 3.47e-13 ***
treat 0.83532 0.54539 1.532 0.128016
periodDuring -8.17150 1.43414 -5.698 7.55e-08 ***
factor(year)1989 -7.04426 1.60607 -4.386 2.34e-05 ***
factor(year)1990 -8.24017 1.47534 -5.585 1.28e-07 ***
factor(year)1991 -8.37095 1.55638 -5.378 3.31e-07 ***
factor(year)1992 -7.88904 1.60607 -4.912 2.62e-06 ***
factor(year)1993 -9.16325 1.49818 -6.116 1.01e-08 ***
factor(year)1994 -7.54756 1.59118 -4.743 5.38e-06 ***
factor(year)1995 -5.82314 1.59118 -3.660 0.000364 ***
factor(year)1996 -10.24260 2.32925 -4.397 2.23e-05 ***
factor(year)1998 -6.49257 1.59118 -4.080 7.74e-05 ***
factor(year)1999 -9.67131 1.61945 -5.972 2.04e-08 ***
factor(year)2001 -9.92257 1.66945 -5.944 2.34e-08 ***
factor(year)2002 -9.11274 1.52806 -5.964 2.13e-08 ***
factor(year)2003 -0.04582 0.70970 -0.065 0.948625
factor(year)2004 -0.27698 0.72811 -0.380 0.704250
factor(year)2005 -1.48769 0.65345 -2.277 0.024416 *
factor(year)2006 -0.87121 0.61827 -1.409 0.161159
factor(year)2007 NA NA NA NA
treat:periodDuring -0.95195 0.69777 -1.364 0.174803
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.902 on 132 degrees of freedom
Multiple R-squared: 0.3654, Adjusted R-squared: 0.2741
F-statistic: 4.001 on 19 and 132 DF, p-value: 9.902e-07
Code
par(mfrow =c(1, 2))plot(mod_did_pre_spring, which =1) # Residuals vs fittedplot(mod_did_pre_spring, which =2) # QQ-plot
lag Autocorrelation D-W Statistic p-value
1 0.09870169 1.794782 0.014
Alternative hypothesis: rho != 0
Diagnostics: The residuals vs. fitted plot shows reasonably even spread. The QQ plot is approximately normal with minor right-tail deviation from a few high-leverage observations. The Durbin-Watson test detects mild positive autocorrelation (DW = 1.79, p = 0.01), which may slightly inflate standard error precision; however, given the non-significant treatment effect, this does not change the substantive conclusion.
The spring DiD model shows no significant treatment effect during the aquaculture period (treat:periodDuring, p = 0.175), consistent with the spring BACI result. This supports the conclusion that aquaculture did not alter minnow abundance in spring, likely due to high overwinter mortality resetting abundance before spring sampling each year.
Robustness Check
Adjusting the window
The robustness check shifts the period boundaries (Before: pre-1996, During: 1996–2004) to assess sensitivity to the period definition.
Code
pre_aqua_fall_mod <- aquaculture %>%mutate(period =case_when( year <1996~"Before", year <=2004~"During",TRUE~"After")) %>%filter(period !="After", season =="Fall") mod_did_pre_fall_mod <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_fall_mod)summary(mod_did_pre_fall_mod)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_fall_mod)
Residuals:
Min 1Q Median 3Q Max
-7.7676 -1.5664 -0.6158 0.4748 18.5138
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.9997 2.7198 0.735 0.4653
treat 0.9902 1.8835 0.526 0.6012
periodDuring 1.3070 3.1557 0.414 0.6803
factor(year)1989 2.6870 3.5219 0.763 0.4487
factor(year)1990 -0.2358 3.8973 -0.061 0.9520
factor(year)1991 1.5501 3.8973 0.398 0.6923
factor(year)1992 -0.2220 3.4617 -0.064 0.9491
factor(year)1993 -1.2022 3.3103 -0.363 0.7179
factor(year)1994 6.1459 3.5980 1.708 0.0931 .
factor(year)1995 3.5230 5.4397 0.648 0.5199
factor(year)1996 -2.5453 3.7277 -0.683 0.4975
factor(year)1998 6.2850 4.9753 1.263 0.2117
factor(year)1999 -1.4756 2.1451 -0.688 0.4944
factor(year)2002 -1.7851 2.1902 -0.815 0.4185
factor(year)2003 -0.5290 2.2856 -0.231 0.8178
factor(year)2004 NA NA NA NA
treat:periodDuring 2.6501 2.4853 1.066 0.2909
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.711 on 56 degrees of freedom
Multiple R-squared: 0.2197, Adjusted R-squared: 0.01071
F-statistic: 1.051 on 15 and 56 DF, p-value: 0.4205
Code
pre_aqua_spring_mod <- aquaculture %>%mutate(period =case_when( year <1996~"Before", year <=2004~"During",TRUE~"After")) %>%filter(period !="After", season =="Spring") mod_did_pre_spring_mod <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_spring_mod)summary(mod_did_pre_spring_mod)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_spring_mod)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.1512 -0.4247 1.1120 5.4117
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.8751 1.4165 7.677 3.23e-11 ***
treat 0.8488 0.6619 1.282 0.203388
periodDuring -8.9034 1.5841 -5.620 2.62e-07 ***
factor(year)1989 -7.0496 1.6968 -4.155 8.02e-05 ***
factor(year)1990 -8.2451 1.5586 -5.290 1.02e-06 ***
factor(year)1991 -8.3786 1.6501 -5.078 2.39e-06 ***
factor(year)1992 -7.8944 1.6968 -4.653 1.26e-05 ***
factor(year)1993 -9.1700 1.5866 -5.780 1.34e-07 ***
factor(year)1994 -7.5476 1.6760 -4.503 2.22e-05 ***
factor(year)1995 -5.8231 1.6760 -3.474 0.000825 ***
factor(year)1996 -1.3392 2.1251 -0.630 0.530350
factor(year)1998 2.4109 1.1426 2.110 0.037946 *
factor(year)1999 -0.6500 1.0169 -0.639 0.524474
factor(year)2001 -0.9602 1.1700 -0.821 0.414252
factor(year)2002 -0.1504 0.9313 -0.161 0.872130
factor(year)2003 0.1996 0.8366 0.239 0.812053
factor(year)2004 NA NA NA NA
treat:periodDuring -0.1313 0.9499 -0.138 0.890401
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.003 on 81 degrees of freedom
Multiple R-squared: 0.3989, Adjusted R-squared: 0.2801
F-statistic: 3.359 on 16 and 81 DF, p-value: 0.0001612
Under this alternative window, the fall treatment effect (treat:periodDuring) is no longer significant (p = 0.291), and the spring effect remains non-significant (p = 0.890). Results are therefore sensitive to how the aquaculture window is defined, which warrants caution in interpretation.
Buffer
Remove observations of transition years between starting and stopping aquaculture.
Code
pre_aqua_fall_buff <- aquaculture %>%filter(period !="After", season =="Fall", year !=2003& year !=2008) # Remove year aquaculture began and endedmod_did_pre_fall_buff <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_fall_buff)summary(mod_did_pre_fall_buff)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_fall_buff)
Residuals:
Min 1Q Median 3Q Max
-14.5942 -1.8385 -0.2952 1.2572 22.7675
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.99967 3.45622 0.579 0.56439
treat 2.47479 1.86931 1.324 0.18904
periodDuring 3.99557 3.88835 1.028 0.30703
factor(year)1989 2.09314 4.43529 0.472 0.63817
factor(year)1990 -0.73068 4.92739 -0.148 0.88246
factor(year)1991 1.05524 4.92739 0.214 0.83093
factor(year)1992 -0.96435 4.33495 -0.222 0.82448
factor(year)1993 -2.02694 4.12381 -0.492 0.62431
factor(year)1994 6.14593 4.57215 1.344 0.18242
factor(year)1995 3.52301 6.91244 0.510 0.61159
factor(year)1996 -0.07277 5.77564 -0.013 0.98998
factor(year)1998 7.59200 6.91244 1.098 0.27513
factor(year)1999 0.60837 4.18095 0.146 0.88465
factor(year)2002 0.10463 4.15916 0.025 0.97999
factor(year)2004 -5.67373 2.35203 -2.412 0.01798 *
factor(year)2005 -0.45396 2.61211 -0.174 0.86244
factor(year)2006 -3.00164 2.16952 -1.384 0.17008
factor(year)2007 NA NA NA NA
treat:periodDuring 7.73287 2.53896 3.046 0.00308 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.986 on 86 degrees of freedom
Multiple R-squared: 0.4463, Adjusted R-squared: 0.3369
F-statistic: 4.078 on 17 and 86 DF, p-value: 6.897e-06
Code
pre_aqua_spring_buff <- aquaculture %>%filter(period !="After", season =="Spring", year !=2003& year !=2008)mod_did_pre_spring_buff <-lm(sqrt_cpue ~ treat * period +factor(year), data = pre_aqua_spring_buff)summary(mod_did_pre_spring_buff)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = pre_aqua_spring_buff)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.0906 -0.4911 1.0705 6.1234
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.8751 1.3528 8.039 6.89e-13 ***
treat 0.8353 0.5486 1.523 0.130477
periodDuring -8.2235 1.4456 -5.688 9.09e-08 ***
factor(year)1989 -7.0443 1.6156 -4.360 2.75e-05 ***
factor(year)1990 -8.2402 1.4841 -5.552 1.70e-07 ***
factor(year)1991 -8.3709 1.5656 -5.347 4.29e-07 ***
factor(year)1992 -7.8890 1.6156 -4.883 3.23e-06 ***
factor(year)1993 -9.1633 1.5071 -6.080 1.44e-08 ***
factor(year)1994 -7.5476 1.6006 -4.715 6.52e-06 ***
factor(year)1995 -5.8231 1.6006 -3.638 0.000405 ***
factor(year)1996 -10.2426 2.3431 -4.371 2.63e-05 ***
factor(year)1998 -6.4926 1.6006 -4.056 8.87e-05 ***
factor(year)1999 -9.6713 1.6291 -5.937 2.85e-08 ***
factor(year)2001 -9.9226 1.6794 -5.909 3.26e-08 ***
factor(year)2002 -9.1127 1.5371 -5.928 2.96e-08 ***
factor(year)2004 -0.2817 0.7325 -0.385 0.701213
factor(year)2005 -1.4877 0.6573 -2.263 0.025404 *
factor(year)2006 -0.8764 0.6220 -1.409 0.161405
factor(year)2007 NA NA NA NA
treat:periodDuring -0.8480 0.7261 -1.168 0.245142
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.913 on 121 degrees of freedom
Multiple R-squared: 0.3828, Adjusted R-squared: 0.291
F-statistic: 4.169 on 18 and 121 DF, p-value: 9.516e-07
Excluding the transition years 2003 and 2008, the buffer specification strengthens the fall treatment effect (treat:periodDuring, p = 0.003), while the spring result remains non-significant (p = 0.245). This suggests the fall effect is robust to the inclusion of transition-year noise.
Call:
lm(formula = sqrt_cpue ~ treat * period, data = all_aqua_fall)
Residuals:
Min 1Q Median 3Q Max
-12.574 -3.073 -1.380 2.557 23.395
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.1287 0.9993 3.131 0.001993 **
treat 1.1544 1.4989 0.770 0.442086
periodDuring 0.8038 1.4132 0.569 0.570103
periodAfter 1.8677 1.2081 1.546 0.123629
treat:periodDuring 7.4872 2.0868 3.588 0.000416 ***
treat:periodAfter -0.7597 1.8697 -0.406 0.684925
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.473 on 208 degrees of freedom
Multiple R-squared: 0.2089, Adjusted R-squared: 0.1899
F-statistic: 10.99 on 5 and 208 DF, p-value: 2.12e-09
Code
par(mfrow =c(1, 2)) # 1 row, 2 columnsplot(mod_did_all_fall, which =1)plot(mod_did_all_fall, which =2)
Code
durbinWatsonTest(mod_did_all_fall)
lag Autocorrelation D-W Statistic p-value
1 0.01530469 1.964652 0.546
Alternative hypothesis: rho != 0
Diagnostics: The residuals vs. fitted plot shows heteroscedasticity at higher fitted values, consistent with the extreme aquaculture-period CPUE values. The QQ plot has a heavy right tail from these same observations. The Durbin-Watson test indicates no autocorrelation (DW = 1.96, p = 0.488).
The fall Before and After model shows a significant treatment effect during aquaculture (treat:periodDuring, p < 0.001), consistent with the Before and During results. The post-aquaculture interaction (treat:periodAfter, p = 0.685) is not significant, suggesting that any treatment effect did not persist after operations ended.
Spring
Code
all_aqua_spring <- aquaculture %>%filter(season =="Spring")mod_did_all_spring <-lm(sqrt_cpue ~ treat * period +factor(year), data = all_aqua_spring)summary(mod_did_all_spring)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = all_aqua_spring)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.0652 -0.3948 0.7467 6.0714
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.87506 1.27869 8.505 4.23e-15 ***
treat 0.83532 0.51858 1.611 0.10881
periodDuring -8.17150 1.36365 -5.992 9.48e-09 ***
periodAfter -9.11371 1.38114 -6.599 3.64e-10 ***
factor(year)1989 -7.04426 1.52712 -4.613 7.08e-06 ***
factor(year)1990 -8.24017 1.40282 -5.874 1.75e-08 ***
factor(year)1991 -8.37095 1.47987 -5.657 5.29e-08 ***
factor(year)1992 -7.88904 1.52712 -5.166 5.77e-07 ***
factor(year)1993 -9.16325 1.42453 -6.432 9.07e-10 ***
factor(year)1994 -7.54756 1.51297 -4.989 1.32e-06 ***
factor(year)1995 -5.82314 1.51297 -3.849 0.00016 ***
factor(year)1996 -10.24260 2.21476 -4.625 6.73e-06 ***
factor(year)1998 -6.49257 1.51297 -4.291 2.77e-05 ***
factor(year)1999 -9.67131 1.53985 -6.281 2.06e-09 ***
factor(year)2001 -9.92257 1.58739 -6.251 2.42e-09 ***
factor(year)2002 -9.11274 1.45294 -6.272 2.16e-09 ***
factor(year)2003 -0.04582 0.67481 -0.068 0.94594
factor(year)2004 -0.27698 0.69232 -0.400 0.68953
factor(year)2005 -1.48769 0.62133 -2.394 0.01757 *
factor(year)2006 -0.87121 0.58788 -1.482 0.13993
factor(year)2007 NA NA NA NA
factor(year)2008 0.25096 0.77315 0.325 0.74583
factor(year)2009 -0.76346 0.70680 -1.080 0.28137
factor(year)2010 0.93905 0.72790 1.290 0.19851
factor(year)2011 -0.66631 0.96256 -0.692 0.48960
factor(year)2012 1.07331 0.90417 1.187 0.23661
factor(year)2013 NA NA NA NA
treat:periodDuring -0.95195 0.66347 -1.435 0.15291
treat:periodAfter -0.43854 0.72330 -0.606 0.54500
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.808 on 200 degrees of freedom
Multiple R-squared: 0.3365, Adjusted R-squared: 0.2503
F-statistic: 3.902 on 26 and 200 DF, p-value: 1.833e-08
Code
par(mfrow =c(1, 2))plot(mod_did_all_spring, which =1)plot(mod_did_all_spring, which =2)
Code
durbinWatsonTest(mod_did_all_spring)
lag Autocorrelation D-W Statistic p-value
1 0.1460512 1.699108 0
Alternative hypothesis: rho != 0
Diagnostics: The residuals vs. fitted plot shows roughly constant spread. The QQ plot shows slight right-tail deviation. The Durbin-Watson test detects positive autocorrelation (DW = 1.70, p < 0.001), which may modestly inflate standard error precision; this is noted further in the Critical Evaluation, and does not alter the non-significant results.
Neither the during-aquaculture (treat:periodDuring, p = 0.153) nor post-aquaculture (treat:periodAfter, p = 0.545) interaction terms are significant in spring, consistent with all prior spring results.
Robustness Check
Adjusting the window
Code
all_aqua_fall_mod <- aquaculture %>%filter(season =="Fall") %>%mutate(period =case_when( year <1996~"Before", year <=2004~"During",TRUE~"After"))mod_did_all_fall_mod <-lm(sqrt_cpue ~ treat * period +factor(year), data = all_aqua_fall_mod)summary(mod_did_all_fall_mod)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = all_aqua_fall_mod)
Residuals:
Min 1Q Median 3Q Max
-11.035 -2.801 -1.005 1.681 26.184
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.6301 1.1954 4.710 4.81e-06 ***
treat 4.5136 1.0562 4.273 3.06e-05 ***
periodBefore -3.6305 3.3812 -1.074 0.28433
periodDuring -2.3235 2.2118 -1.051 0.29483
factor(year)1989 2.6870 4.0955 0.656 0.51258
factor(year)1990 -0.2358 4.5321 -0.052 0.95856
factor(year)1991 1.5501 4.5321 0.342 0.73271
factor(year)1992 -0.2220 4.0255 -0.055 0.95607
factor(year)1993 -1.2022 3.8495 -0.312 0.75517
factor(year)1994 6.1459 4.1841 1.469 0.14353
factor(year)1995 3.5230 6.3257 0.557 0.57823
factor(year)1996 -2.5453 4.3349 -0.587 0.55780
factor(year)1998 6.2850 5.7857 1.086 0.27873
factor(year)1999 -1.4756 2.4946 -0.592 0.55486
factor(year)2002 -1.7851 2.5470 -0.701 0.48424
factor(year)2003 -0.5290 2.6578 -0.199 0.84245
factor(year)2004 NA NA NA NA
factor(year)2005 2.0464 2.3103 0.886 0.37686
factor(year)2006 -0.3589 1.8996 -0.189 0.85034
factor(year)2007 3.5680 1.9125 1.866 0.06365 .
factor(year)2008 -3.4467 1.9340 -1.782 0.07634 .
factor(year)2009 -5.4904 1.7121 -3.207 0.00158 **
factor(year)2010 -2.1425 1.9626 -1.092 0.27636
factor(year)2011 0.6586 1.8931 0.348 0.72830
factor(year)2012 -0.9035 2.0390 -0.443 0.65818
factor(year)2013 NA NA NA NA
treat:periodBefore -3.5234 2.4317 -1.449 0.14901
treat:periodDuring -0.8733 2.1612 -0.404 0.68661
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.478 on 188 degrees of freedom
Multiple R-squared: 0.2837, Adjusted R-squared: 0.1885
F-statistic: 2.979 on 25 and 188 DF, p-value: 1.291e-05
Code
all_aqua_spring_mod <- aquaculture %>%filter(season =="Spring") %>%mutate(period =case_when( year <1996~"Before", year <=2004~"During",TRUE~"After"))mod_did_all_spring_mod <-lm(sqrt_cpue ~ treat * period +factor(year), data = all_aqua_spring_mod)summary(mod_did_all_spring_mod)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = all_aqua_spring_mod)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.0808 -0.4441 0.8936 6.1118
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.76135 0.52190 3.375 0.000887 ***
treat -0.03581 0.35245 -0.102 0.919177
periodBefore 9.11371 1.38081 6.600 3.61e-10 ***
periodDuring 0.21028 0.82584 0.255 0.799275
factor(year)1989 -7.04964 1.53136 -4.604 7.38e-06 ***
factor(year)1990 -8.24506 1.40662 -5.862 1.87e-08 ***
factor(year)1991 -8.37863 1.48920 -5.626 6.16e-08 ***
factor(year)1992 -7.89442 1.53136 -5.155 6.07e-07 ***
factor(year)1993 -9.16998 1.43190 -6.404 1.06e-09 ***
factor(year)1994 -7.54756 1.51261 -4.990 1.31e-06 ***
factor(year)1995 -5.82314 1.51261 -3.850 0.000159 ***
factor(year)1996 -1.33917 1.91786 -0.698 0.485823
factor(year)1998 2.41086 1.03119 2.338 0.020378 *
factor(year)1999 -0.65001 0.91771 -0.708 0.479586
factor(year)2001 -0.96020 1.05596 -0.909 0.364279
factor(year)2002 -0.15038 0.84053 -0.179 0.858191
factor(year)2003 0.19957 0.75502 0.264 0.791800
factor(year)2004 NA NA NA NA
factor(year)2005 -0.58588 0.71254 -0.822 0.411920
factor(year)2006 0.02656 0.68803 0.039 0.969245
factor(year)2007 0.90181 0.69643 1.295 0.196850
factor(year)2008 0.48389 0.74821 0.647 0.518551
factor(year)2009 -0.54717 0.68327 -0.801 0.424192
factor(year)2010 1.18949 0.69717 1.706 0.089528 .
factor(year)2011 -0.66631 0.96233 -0.692 0.489494
factor(year)2012 1.07331 0.90395 1.187 0.236498
factor(year)2013 NA NA NA NA
treat:periodBefore 0.88457 0.69358 1.275 0.203660
treat:periodDuring 0.75326 0.70877 1.063 0.289165
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.808 on 200 degrees of freedom
Multiple R-squared: 0.3369, Adjusted R-squared: 0.2507
F-statistic: 3.908 on 26 and 200 DF, p-value: 1.766e-08
Under the alternative period window, neither the fall nor spring treatment effects are significant across any period comparison, consistent with the Before and During robustness check.
Buffer
Code
all_aqua_fall_buff <- aquaculture %>%filter(period !="After", season =="Fall", year !=2003& year !=2008) # Remove year aquaculture began and endedmod_did_all_fall_buff <-lm(sqrt_cpue ~ treat * period +factor(year), data = all_aqua_fall_buff)summary(mod_did_all_fall_buff)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = all_aqua_fall_buff)
Residuals:
Min 1Q Median 3Q Max
-14.5942 -1.8385 -0.2952 1.2572 22.7675
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.99967 3.45622 0.579 0.56439
treat 2.47479 1.86931 1.324 0.18904
periodDuring 3.99557 3.88835 1.028 0.30703
factor(year)1989 2.09314 4.43529 0.472 0.63817
factor(year)1990 -0.73068 4.92739 -0.148 0.88246
factor(year)1991 1.05524 4.92739 0.214 0.83093
factor(year)1992 -0.96435 4.33495 -0.222 0.82448
factor(year)1993 -2.02694 4.12381 -0.492 0.62431
factor(year)1994 6.14593 4.57215 1.344 0.18242
factor(year)1995 3.52301 6.91244 0.510 0.61159
factor(year)1996 -0.07277 5.77564 -0.013 0.98998
factor(year)1998 7.59200 6.91244 1.098 0.27513
factor(year)1999 0.60837 4.18095 0.146 0.88465
factor(year)2002 0.10463 4.15916 0.025 0.97999
factor(year)2004 -5.67373 2.35203 -2.412 0.01798 *
factor(year)2005 -0.45396 2.61211 -0.174 0.86244
factor(year)2006 -3.00164 2.16952 -1.384 0.17008
factor(year)2007 NA NA NA NA
treat:periodDuring 7.73287 2.53896 3.046 0.00308 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.986 on 86 degrees of freedom
Multiple R-squared: 0.4463, Adjusted R-squared: 0.3369
F-statistic: 4.078 on 17 and 86 DF, p-value: 6.897e-06
Code
all_aqua_spring_buff <- aquaculture %>%filter(period !="After", season =="Spring", year !=2003& year !=2008) # Remove year aquaculture began and endedmod_did_all_spring_buff <-lm(sqrt_cpue ~ treat * period +factor(year), data = all_aqua_spring_buff)summary(mod_did_all_spring_buff)
Call:
lm(formula = sqrt_cpue ~ treat * period + factor(year), data = all_aqua_spring_buff)
Residuals:
Min 1Q Median 3Q Max
-3.7609 -1.0906 -0.4911 1.0705 6.1234
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.8751 1.3528 8.039 6.89e-13 ***
treat 0.8353 0.5486 1.523 0.130477
periodDuring -8.2235 1.4456 -5.688 9.09e-08 ***
factor(year)1989 -7.0443 1.6156 -4.360 2.75e-05 ***
factor(year)1990 -8.2402 1.4841 -5.552 1.70e-07 ***
factor(year)1991 -8.3709 1.5656 -5.347 4.29e-07 ***
factor(year)1992 -7.8890 1.6156 -4.883 3.23e-06 ***
factor(year)1993 -9.1633 1.5071 -6.080 1.44e-08 ***
factor(year)1994 -7.5476 1.6006 -4.715 6.52e-06 ***
factor(year)1995 -5.8231 1.6006 -3.638 0.000405 ***
factor(year)1996 -10.2426 2.3431 -4.371 2.63e-05 ***
factor(year)1998 -6.4926 1.6006 -4.056 8.87e-05 ***
factor(year)1999 -9.6713 1.6291 -5.937 2.85e-08 ***
factor(year)2001 -9.9226 1.6794 -5.909 3.26e-08 ***
factor(year)2002 -9.1127 1.5371 -5.928 2.96e-08 ***
factor(year)2004 -0.2817 0.7325 -0.385 0.701213
factor(year)2005 -1.4877 0.6573 -2.263 0.025404 *
factor(year)2006 -0.8764 0.6220 -1.409 0.161405
factor(year)2007 NA NA NA NA
treat:periodDuring -0.8480 0.7261 -1.168 0.245142
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.913 on 121 degrees of freedom
Multiple R-squared: 0.3828, Adjusted R-squared: 0.291
F-statistic: 4.169 on 18 and 121 DF, p-value: 9.516e-07
Excluding transition years, the fall treatment effect during aquaculture remains significant (treat:periodDuring, p = 0.003), while spring remains non-significant (p = 0.245). The post-aquaculture interaction is non-significant in both seasons, indicating no lasting legacy effect once operations ended.
Rennie et al. 2019 Critical Evaluation
Standard BACI designs test whether a treated site changes more than a control site following an intervention, but they rely on an implicit assumption that the control site provides a valid counterfactual trajectory. If the treated and control lakes were already diverging prior to aquaculture, the BACI interaction term could mistakenly attribute that divergence to the treatment effect.
The results suggest that any treatment effect is both season-specific and temporally limited. Evidence of increased minnow abundance during aquaculture is observed in the fall but not in the spring, and this effect does not persist after aquaculture operations end. This seasonal asymmetry is ecologically plausible and aligns with the broader ecosystem responses documented by Rennie et al. (2019), where increased nutrient inputs during aquaculture altered productivity and food availability, potentially benefiting forage fish during certain time of year.
Addition of Difference-in-Differences
This analysis strengthens the traditional BACI framework by incorporating a Difference-in-Differences perspective and explicitly evaluating the assumptions required for causal interpretation. The parallel trends test makes this assumption explicit and allows it to be evaluated using pre-treatment data. The absence of a statistically significant treatment-by-year interaction suggests that we fail to reject the parallel trends assumption, meaning there is no strong evidence of differential trends prior to aquaculture.
Another important strength of this analysis is the separation of the aquaculture period into distinct phases. Traditional BACI approaches collapse all post-treatment observations into a single “after” category, which can obscure differences between short-term and longer term responses. By estimating separate DiD models for the during-treatment and post-treatment periods, this analysis distinguishes between immediate effects of aquaculture and potential legacy effects after operations cease.
Finally, this analysis added robustness checks which were largely missing from the original analysis.
However, several important limitations remain:
The study includes only a single treated lake. With no replication at the treatment level, it is not possible to fully separate treatment effects from lake-specific dynamics. As a result, the findings should be interpreted as evidence consistent with an aquaculture effect rather than definitive causal proof. In addition, the limited number of observations reduces statistical power, making it more difficult to detect smaller or more gradual treatment effects.
There are model-based limitations. Diagnostic tests indicate heteroscedasticity and occasional deviations from normality, particularly in fall models driven by extreme CPUE values during the aquaculture period. Evidence of mild autocorrelation in the spring models suggest that standard errors may be somewhat underestimated, which could affect inference, although this does not alter the substantive conclusions given the lack of significant effects in spring.
Finally, the robustness checks demonstrate that results are sensitive to modeling choices, including the definition of treatment periods and the inclusion of transition years. While the fall treatment effect remains detectable under some specifications, it weakens or disappears under others, reinforcing the need for cautious interpretation.
Despite these limitations, the combination of visualization, BACI modeling, parallel trends testing, and Difference-in-Differences analysis provides a coherent and transparent evaluation of the aquaculture experiment.
Overall, the results illustrate how combining BACI designs with a DiD framework can strengthen causal interpretation in ecological impact studies. While the analysis is limited by the presence of only a single treatment lake, the parallel trends test and the separation of during- and post-treatment periods provide a more transparent evaluation of how aquaculture activity may influence fish community dynamics.
Conclusion
This analysis provides partial evidence that freshwater aquaculture influenced minnow abundance in the treatment lake, with increases observed during the aquaculture period in fall sampling but not in spring. These effects appear to be short-lived, as there is no strong evidence that difference between the treated and control lakes persist after aquaculture operations end. By extending the BACI framework with a Difference-in-Differences perspective and explicitly evaluating the parallel trends assumption, this study offers a more rigorous and transparent approach to causal inference than a standard BACI analysis alone. The results highlight the importance of testing underlying assumptions and considering temporal structure when interpreting ecological impact studies. At the same time, the findings should be interpreted in the broader context of the ecosystem-level experiment conducted by Rennie et al. (2019). Minnow abundance represents only one component of the lake food web, and observed changes likely reflect underlying shifts in nutrient dynamics and productivity associated with aquaculture activity.
Overall, this analysis illustrates both the potential and the limitations of applying quasi-experimental methods in ecological systems. While the results are consistent with a short-term aquaculture effect on forage fish abundance, the lack of replication and sensitivity to modeling choices limit the strength of causal claims. Future work incorporating additional treated systems or longer time series would improve statistical power and strengthen inference and provide a more complete understanding of how aquaculture influences freshwater ecosystems.
