Random effects without repeated measures?

Question

As I understand it, random effects are good for repeated measure designs where a random intercept, for example, represents deviations from an over-all trend. Given data with many grouping factors, however, what if the repeated measure aspect isn't at the lowest / observation / subject ID level, but "higher up" in the grouping structure? Since I'm not particularly interested in differences between "subject ID's" (locations, in my case), should I still include them? How does one include them in mgcv? For example, I have a dataset where each observation (each row) is a single count at a particular site. All the same 47 sites (same unique non-randomly selected locations) are surveyed 2x per year (2 seasons), every year. So, in a way, there ARE multiple observations, but it's at the season and year level, not the site level (1 obs. per site, 47 sites * 2 seasons =94 obs. per year, over 15 years).

My fSeason covariate only has 2 levels, so I'm modeling it as a fixed effect. My year effect is also fixed (and can probably vary by season). fSite only has 1 observation per level, but many over the entire study period. However, isn't the repeated aspect taken care of in the season/year fixed effects? So in this instance, is a random effect even necessary?

My guess:

From the list of models below, I think model m2 is the one I need. What I think I'm doing with each term (correct me if I'm saying this wrong) is modeling the mean count over the entire range of the covariate(s). When I exclude the site term, R is just modeling the mean count over the time (and depth) effect, and I'm explicitly saying there are no deviations from these main effects (all sites follow this trend). I need a random effect, despite the lack of a traditional repeated measure per site, because it makes sense that there are deviations from these main effects (?) With a random intercept I'm modeling for a typical site and with a fixed effect I model for every one of these sites specifically.

# fCYR = factor year
# CYR.std = continuous year, rescaled (df$CYR.std <- df$CYR - min(df$CYR))
# --> So, 2008 = year "0", 2009 = year "1", etc...
# fSite = factor site
# fSeason = factor season

List of possible models:

# No fSite term (b/c only 1 obs. per fSite)? 
m1 <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
              offset(log(area_sampled)), # Always 3m^2 at each site
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)

# or IS (fSite, bs='re') important and itself equals many obs per year?
m2 <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            (fSite, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)

# 2 random effects? But then I think I'm modeling the year effect twice...
m3a <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            (fSite, bs='re') + 
            (fCYR, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)

# What about an interaction term?
m3b <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            (fSite, bs='re') + 
            (fSite, fCYR, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)

Data example (rounded number of sites to 50 to make this easier):

set.seed(12345)

# Hypothetical fish counts from negative binomial distribution
df <- as.data.frame(rnbinom(1000, mu = 4, size = 1))

df$year <- rep(2011:2020, each=100)
df$CYR.std <- df$year - min(df$year)
df$fCYR <- as.factor(df$year)

df$site <- seq(1, 50, 1)
df$fSite <- as.factor(df$site)

df$season <- rep(c("DRY", "WET"), each=50)
df$fSeason <- as.factor(df$season)

# Depth (continuous covariate)
df$depth <- sample(30:200, 1000, replace = TRUE)

names(df)[1] <- "count"

UPDATE:

Below is a map of how the sites are distributed. Some are only 300 feet away from each another, while others can be a little under a mile apart. I have the lat/lon of each site, of which the coordinates have a 60ft radius (around the boat we sample from) of space we have the option of surveying (traps are thrown anywhere in this radius).

Answer 1

It seems model m2 has the correct syntax for my situation (plus the te() spatial variation term). This is because although the surveyed counts are only recorded once per season, they ARE "repeatedly measured" every year, and each year after that (a sort of "delayed" repeated measure, but repeated nonetheless).

Also, a variable is only entered into a model once in the same format (as a continuous or factor covariate). So year as a smooth fixed effect, with season, fSeason + s(CYR.std, by=fSeason) is sufficient to capture the long-term trend and there's no need for the additional year random effect (fCYR, bs='re'), for example. So the final model could be...

m2 <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            te(Latitude, Longitude) + 
            (fSite, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)

If, on the other hand, I wanted to test if an annual trend that varied per site was present, I could either use the parametric version fSite * fCYR, or the random effect equivalent: (fSite, fCYR, bs='re'), but not both. This is what the factor interaction term is for.

Lastly, as part of my original post, model comparisons with AIC/BIC should only be done when each model has the same random effect syntax and the same fixed effect variables (+ select = FALSE).

# These models are comparable:
model_a <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            te(Latitude, Longitude) + 
            (fSite, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            family = nb(),
            data = df)

model_b <- bam(count ~ s(depth) + fSeason * CYR.std + 
            te(Latitude, Longitude) + 
            (fSite, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            family = nb(),
            data = df)

# But this one isn't comparable to the two above (te(Latitude, Longitude) is missing)
m2 <- bam(count ~ s(depth) + fSeason + s(CYR.std, by=fSeason) + 
            (fSite, bs='re') + 
              offset(log(area_sampled)),
            method = "fREML",
            discrete = TRUE,
            select = TRUE,
            family = nb(),
            data = df)