I'm still learning about mixed effects models, so bear with me here. I'm interested in modelling a binary response using a generalized additive mixed effects model with "year" as a covariate and random effect, but I've seen far smarter people than I argue for and against it (a covariate can't be both). For example:

Absolutely. In fact, the vast majority of the time, you absolutely should include a fixed effect. The reason for this is that random effects are restrained to ∑γ=0 , or always centered around 0. Thus, the random effect is the individual's estimated deviation from the group average for that individual. By leaving out the fixed effect, you would imply that the average effect of time must be 0.

I've also heard from colleagues:

"In short, no, a variable can't be both fixed and random. In Frequentist statistics, fixed effects are assumed to represent an actual "true" effect of the variable on the target mean, while random effects are assumed to have an effect on the mean that's randomly drawn from some distribution of possible values."

I'd like to know if there's a way around this problem. Please let me know if I'm misinterpreting these points of view.

My experimental design:

Fish are collected and there stomachs examined to see if they ate something, or not (0/1), at the same 12 sites, every 9 months, every year. Many sites = a zone, and some zones have many more sites than others. My repeated measure is the length of the fish as another covariate.

If I'm interested in differences between years, differences between zones, and their interaction, but I also want to capture similarity between observations taken in the same site, zone, and year, is converting year to a continuous variable in the random effects and a factor in the fixed effects one possible solution? Where every zone has it's own trend through time (s(fZone, CYR, bs='re')) and the fixed effects shows where those differences are (fZone*fCYR)? Or, if it must be one or the other, can fixed effects also capture correlation structures without random effects?

  • fCYR = factor calendar year
  • fZone = factor zone
  • CYR = continuous year

1 Answer 1


I think the way you are thinking about this is conflating different kinds random effects into a single entity.

It wouldn't make sense to fit a model like

y ~ fZone + (1 | fZone)

where we have a parametric fixed effect for zone and a random effect for zone. What one is trying to do with either of these terms is to account for the mean of Y for each zone. The random intercept achieves this with some shrinkage towards the overall model intercept/constant term. But as shown, the pseudo-code model would be including the Zone means twice.

However, in the context of your time covariate for year, you would want to do something like

y ~ cyr + (1 + cyr | fZone)

for example, where the cyr term is referred to as the "population" level effect of cyr and represents that average change in y for a unit change in cyr (having accounted for variation in the effects of cyr on each zone). The (1 + cyr | fZone) term models a separate mean value of y for each zone, plus a separate "random slope" of year for each zone, which captures the variation in the change in y for a unit change in cyr over the different zones (each zone gets it's own slope/trend).

But you shouldn't think that factor terms can't be both fixed and random. If you had more levels of variation in the data, it could very much make sense to have, say, an overall treatment effect, plus different treatment effects for each level of a grouping variable.

y ~ treatment + (1 + treatment | field)

where treatment is a factor and the experiment was conducted over $f$ fields, with replication at the field level. The fixed treatment effect would capture the average effect of the treatment over the whole data set, but the random effect of treatment would allow the effects of the treatment levels to vary between the fields. For example, the effect of a fertilizer treatment over the control would be to increase yields (y) in all samples, but that effect could very much vary between fields because of natural (or human-induced) variation in the organic content of the soils; all else equal, we might expect less of an effect of an application of nitrogen fertilizer to a soil that is already rich in N content, compared with the same application of N fertilizer in a soil that is deplete in N.

The main take away here is to think about how your covariates vary in the data and at what levels in the data hierarchy the covariates may have effects. It is too simplistic to distill "random effects" down to the sorts of blanket statement your colleague provided.

In your specific example of a model (GAM), the term s(fZone, CYR, bs = "re") doesn't include the different group means; one would typically do s(fZone, bs = "re") + s(fZone, CYR, bs = "re") to get the equivalent of (1 | fZone) + (0 + CYR | fZone) (because mgcv isn't fitting correlated random effects).

I personally don't think it makes sense to include both Year as a continuous (linear) trend and as a parametric factor effect as you intend - you are basically including the trend part twice. The factor based time effect is essentially a non-linear trend, but that could easily model the linear trend you are setting up in the random effect smooth; at that point your model is likely getting unidentifiable.

In general you could model the trend(s) as:

  1. simply a factor time fixed effect
  2. a random effect of the time factor
  3. a linear trend via a continuous time variable
  4. a smooth non-linear trend via a spline of the continuous time variable
  5. a MRF representation of an AR(1)

You could have combinations of these; a low-EDF smooth of time plus a random intercept of factor time, for example. All of these could be nested within zones. Which approach you take is highly dependent upon the data you to hand (how many years, how many zones, etc) and the questions you wish to ask.

  • 1
    $\begingroup$ Thank you for the detailed response! I think one thing I forgot to mention that could be important (since you mentioned sample size) is there are only 4 zones and 12 years, so zone may have to just stay a fixed effect, regardless of the nested structure. Station and year could then be random effects. Like this: CYR + fZone + s(fStation, bs = "re") + s(fStation, CYR, bs = "re"). $\endgroup$
    – Nate
    Commented May 4, 2023 at 13:24
  • $\begingroup$ The old argument about the number of levels in a grouping variable being greater than 5 or 6 before you consider using a random effect to represent it it largely motivated from "old style" random effects that are motivated from the view of estimating a variance parameter. "New style" random effects look, walk, and quack like random effects but aren't motivated from the same argument of a variance component. 1/2 $\endgroup$ Commented May 5, 2023 at 15:56
  • 1
    $\begingroup$ So it's less clear to me that it actually matters that you only have four levels; we can subject the estimated $\beta$ to shrinkage through the ridge penalty we place on this "smooth" & choose the smoothing parameter for that "smooth" using the same approaches as we choose the smoothing parameter(s) for other smooth(s). There's nothing here that needs these terms to be related to a Gaussian with a variance term that we estimate. It just so happens that this penalised spline form can be written in random effect form where the smoothing parameter is the inverse of the variance term. 2/2 $\endgroup$ Commented May 5, 2023 at 15:59

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.