4
$\begingroup$

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).

enter image description here

$\endgroup$
12
  • 3
    $\begingroup$ I don't know why you say you don't have repeated measures when you clearly do have repeated measures at the site level. You do need to account for between site differences in your model even if you aren't interested in those differences because otherwise you'll be violating the assumptions of the model, messing up any inference you might do. This is because the data conditional upon the model won't be independent. $\endgroup$ Dec 31, 2023 at 12:02
  • 2
    $\begingroup$ How you model the site effects depends on how they are distributed in space; a random intercept of fSite would be less useful if there is spatial variation, which you could model with s(x,y) say. $\endgroup$ Dec 31, 2023 at 12:02
  • 1
    $\begingroup$ Gavin is right by the way. Your data necessarily has repeated measures so you will need to model that directly in mgcv. How that matters with spatial data is less clear to me given that is not my topic expertise. $\endgroup$ Dec 31, 2023 at 16:42
  • 1
    $\begingroup$ First, don't let statistical significance be the barometer of choice here. For the record, the te function constructs a tensor product smooth between latitude and longitude. That is not the same as modeling s(fSite, bs = "re"), which constructs a random intercept based on the grouping of each site. So the first case is simply showing the relation between latitude and longitude while the second case is estimating the differences in conditional means for each site. Again I'm not an expert at spatial data but I would imagine one makes more sense than the other depending on the context. $\endgroup$ Jan 1 at 1:17
  • 1
    $\begingroup$ Again, I can't comment specifically in your case since I don't know the specifics of spatial models, but I would advise looking through this page and this page to figure out which works for you. $\endgroup$ Jan 1 at 1:18

1 Answer 1

1
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.