extending BACI analysis to lme, glm

Question

I am analyzing data from a BACI (before-after, control-impact) design. BACI analyses typically look for changes in the mean response before(B) and after(A) at control(C) and impact(I) sites, but I need to add a continuous predictor variable to explain my response. I am wondering about extensions of BACIs from factor only models (i.e. response ~ period*treatment) to models with other terms (e.g. linear models, glms).

Study design: the impact is removal of large organisms from a population. There are 5 impact sites and 4 control sites. It's unbalanced,there are a varying number of organisms (i.e. responses) per site and period. I am testing the density dependent response of a variety of demographic rates, and some require models require adding terms due to their inherent nature (i.e. maturity is binomial) or sampling considerations (i.e. the organisms available to measure growth rate is affected by the treatment).

Sine growth and maturity are allometric (related to size), and I have removed large individuals, I need to include terms for size of the remaining organisms (i.e. a continuous variable for length). Since small individuals have more potential for growth, and I have removed large individuals from impact sites, I need to include a growth term to avoid type I error that would occur if I only analyzed the distribution of growth rates.

I have been following the wonderful guide Carl Schwarz wrote on BACI analysis and the classics he cites(Underwood, 1991; Roberts, 1993). I was extending the typical method of estimating marginal means for each factor group (i.e. Before-Impact, After-Impact, Before-Control, After-Control) and analyzing the marginal means of the groups with the emmeans() package as suggested by Schwarz. The whole point of this is to test the interaction, or the difference in the differences between periods and groups.

Why use emmeans() and not just summary(BACI_model)? Schwarz says it helps when you have unbalanced data. Further, it seems like it's fine if your BACI model contains a continuous predictor. It essentially becomes an analysis of the difference in intercept.

An example:

library(dplyr)
library(ggplot2)
library(emmeans)
n_level <- c(100, 200) #the number of observations per level combo (levels = control, impact), stays same by period

UC_data <- data.frame(period = factor(c(rep("before", sum(n_level)), 
                                        rep("after", sum(n_level))), levels= c("before", "after"))) %>%
  group_by(period) %>%
  mutate(treatment = factor(c(rep("control", n_level[1]), rep("impact", n_level[2])))) %>%
  ungroup() %>%
  mutate(length = rnorm(sum(n_level*2), mean= 150, sd=80), #simulate lengths
         error = rnorm(sum(n_level*2), mean=0, sd=2)) %>% #sim.error 
  mutate(response = ifelse(treatment=="control", (length*-0.02) + 2 + error, #sim. growth at control 
                           ifelse(period=="before", (length*-0.02) + 2 + error, #sim. gr. at impact before
                                  (length*-0.02)+ 3 + error))) #sim. gr. at impact after - note intercept increase 

#plot 'n see the data
ggplot(UC_data, aes(x=length, y=response, color=treatment, lty=period)) +
  geom_point(alpha=0.5) +
  stat_smooth(method="lm", aes(), se=FALSE) + #lm fit to each factor combo
  labs(title="unbalanced BACI w/continuous predictor", y="growth rate", x= "length + period * treatment")

#check that Schwarz's emmeans() method is same as summary() estimate of coefficient
UC_mod <- lm(response~ length + period*treatment, data=UC_data)
UC_emm <- emmeans(UC_mod,  specs= ~period:treatment)
summary(UC_emm) #look at marginal means for each level

confint(contrast(UC_emm, method= list(baci=c(1,-1,-1,1)))) #method argument flips signs to follow math - i think...
summary(UC_mod)$coefficients[5,1] #estimate of interaction coefficient is same as emmeans... 

I ran some caveman for loops (i.e. Ctrl+A then Ctrl+Enter) and it appears the distribution of contrast/effect estimates is centred on 1, which is what I expect since that what I specified the intercept difference to be.

In my paper, I am using my experimental design (BACI) to test a variety of responses of various rates to an impact, and certain rates need to be modeled differently. I hope to look at the effect size and CI to talk about which rates were significant (difference in effect size will have to be interpreted ecologically). Eventually, I would like to be able to defend this type of analysis for unbalanced data (different number of observations across sites, periods, and treatments) at multiple sites. My (limited) stats intuition tells me it makes sense to test for the significance of the period:treatment interaction in models like lmer(growth ~length +period*treatment + (1|site)) to allow the growth to vary by site or lmer(growth ~ length + period*treatment + (1+length|site)) to allow the allometric relationship (i.e. length's effect on growth) to vary at the site level. The latter results in a much better fit of residuals at the site level and better AIC.

I also want to extend this to binomial maturity models of a form like glm(mat ~ length + period*treatment, family = binomial(link = "logit"))

In summary, why are BACI analyses so heavily based in ANOVA? Is it okay to estimate the BACI from the effect size and SE from summary()? From models that are not factor only (e.g. lme, glm)?

Disclaimer: First post. Happy to edit if asked!