LMM: Non-independence of observations sharing a single fixed-effect value

Question

I am interested in the effect of monkey's stress levels on the pitch of their calls. Each stress measurement is associated with a bout of calling (that is, multiple calls), and in some cases I have multiple bouts per monkey.

Each row represents a single call. I plan to use linear mixed modeling with call pitch as the response variable, level of stress hormones as a (continuous) fixed effect, and the ID of the monkey as a random intercept (to account for nonindependence of calls from a single monkey). The different calls from within a single bout, and their pitches, are not independent from one another. Thus, with this model structure, there is a source of nonindependence among observations that is not accounted for by a random effect. However, the fixed effect, stress hormone level, is identical for all calls within a single bout.

My question is, do I need to change my random effects structure to include Bout ID nested within Monkey ID, to account for this nonindependence? Note that since what I am interested in is the relationship between call pitch and stress, the random effect is not really of interest in and of itself; rather its purpose would solely be to account for nonindependence of observations.

The reason I am hesitant to include Bout ID as a random effect is that Bout ID is obviously highly collinear with stress hormone level, due to my study design: all calls from a single bout share a single stress measurement. Thus inclusion of this random effect might "soak up" all the variation due to stress and mask the significance of the fixed effect - see this question. (Indeed that is exactly what happens when I run the models, but I have been trying to keep this question hypothetical because I'm looking for a-priori, theory-based answers.)

I'm not sure this will matter, but FWIW the structure of my data is not at all balanced: some bouts include many calls and some bouts only include a few; and some monkeys are represented in multiple bouts whereas other monkeys only have 1 bout.

How would you recommend I proceed?

Answer 1

Yes, in principle you should use Bout ID nested within Monkey ID as a random effect. If you don't, then you are effectively assuming that the bout-to-bout variation in call pitch within monkeys is entirely explained by variation in stress hormones. This is an optimistic/potentially anti-conservative assumption. The "masking" of the fixed effect of stress would occur because you can't really infer that stress is causing pitch variation; it could be explained by other, unmeasured/uncontrolled variation among bouts. On the other hand, I would expect that if there really is a strong stress effect, it wouldn't be entirely wiped out by the presence of the random effect. (You could do some simulations to test this.)

Just to add some complication, I think the maximal model for this observational design would be

pitch ~ 1 + hormone + (1 + hormone | MonkeyID) + (1 | MonkeyID:BoutID)

— because (as I understand it) it's plausible that the effect of hormone on pitch varies among monkeys, and because you have measured multiple hormone levels for each monkey, you can in principle quantify this variation. This gets you into the usual mess of whether in practice you have sufficient data to quantify it ... see contrasting views between Barr et al. 2013 vs Bates, Kliegl, Matuschek et al ... (or go Bayesian or ...) (see the relevant section of the GLMM FAQ for more discussion)

code dump

With some apologies for the brevity of explanation: here is code that simulates something like your situation, which you could use for informal exploration or power analysis. One thing I realized while setting it up is that you can simplify the analysis slightly by collapsing the data to the bout level (i.e., averaging pitches within bout). With weights proportional to the number of calls per bout (i.e. variance inversely proportional to number of calls), this analysis is equivalent (this is a variation of the argument made by Murtaugh [2007]) to the full, three-level (monkey/bout/call) analysis - it also clarifies that your effective sample size is equivalent to the number of bouts, not the number of calls (although having many calls per bout will reduce the among-bout variation and improve power ...) In one example I tried (with a larger among-monkey variation in hormone-pitch relationship), using the maximal model greatly improved resolution/power ...

nb <- c(1, 1, 1, 1, 2, 2, 3, 4) ## bouts per monkey
b_tot <- sum(nb) ## total bouts
nc <- round(170/b_tot) ## 170 calls total (approx)
dd <- data.frame(MonkeyID=rep(1:length(nb),nb*nc),
                 BoutID=rep(1:b_tot,each=nc),
                 callID=1:(b_tot*nc))
set.seed(101)
## hormone varies by bout
##  (could vary by monkey/bout)
dd$hormone <- rnorm(b_tot,mean=5)[dd$BoutID]
library(lme4)                 
fullform <- pitch ~ 1 + hormone + (1 + hormone | MonkeyID) + (1|MonkeyID:BoutID)
subform <- pitch ~ 1 + hormone + (1 | MonkeyID/BoutID)
dd$pitch <- simulate(fullform[-2],
                     newdata=dd,
                     newparams=list(beta=c(1,2),
                                    ## theta: scaled var-cov
                                    ## pos 1: among bouts
                                    ## pos 2:4: among monkeys
                                    theta=c(1,0.5,0.25,0.5),
                                    sigma=1  ## residual std dev
                                    ),
                     family=gaussian)[[1]]

library(ggplot2)
ggplot(dd,aes(hormone,pitch,colour=factor(MonkeyID)))+
    geom_point()+
    geom_smooth(method="lm",se=FALSE)

fit0 <- lmer(fullform, data=dd) ## singular fit;
## helps figure out order of theta parameters
library(lmerTest)
fit <- lmer(subform, data=dd)
summary(fit)

library(dplyr)
dd0 <- dd %>% group_by(MonkeyID,BoutID) %>%
    summarise(pitch=mean(pitch),hormone=mean(hormone),n=n())

library(nlme)
summary(lme(pitch~hormone, random = ~1|MonkeyID,
            weights=varFixed(~I(1/n)),
            data=dd0))

summary(lme(pitch~hormone, random = ~hormone|MonkeyID,
            weights=varFixed(~I(1/n)),
            data=dd0))

references

Barr, Dale J., Roger Levy, Christoph Scheepers, and Harry J. Tily. “Random Effects Structure for Confirmatory Hypothesis Testing: Keep It Maximal.” Journal of Memory and Language 68, no. 3 (April 2013): 255–78. https://doi.org/10.1016/j.jml.2012.11.001.

Bates, Douglas, Reinhold Kliegl, Shravan Vasishth, and Harald Baayen. “Parsimonious Mixed Models.” ArXiv:1506.04967 [Stat], June 16, 2015. http://arxiv.org/abs/1506.04967.

Matuschek, Hannes, Reinhold Kliegl, Shravan Vasishth, Harald Baayen, and Douglas Bates. “Balancing Type I Error and Power in Linear Mixed Models.” ArXiv:1511.01864 [Stat], November 5, 2015. http://arxiv.org/abs/1511.01864.

Murtaugh, Paul A. “Simplicity and Complexity in Ecological Data Analysis.” Ecology 88, no. 1 (2007): 56–62.

Answer 2

In looking into how to approach data structured in this way, I found this to be helpful:

If one analyzes clustered data with a model that assumes independence of observations, such as the general(ized) linear model, and ignores the clustering of observations, while the regression coefficients will be estimated without bias, the standard errors of time-invariant predictors will be underestimated (see, e.g., Maas and Hox, 2005; Mundfrom and Schultz, 2002).

from: McNeish & Harring, 2015

It seems to me that my version of the model lacking the BoutID random effect probably underestimated the standard error of the effect of stress on pitch (because it assumed a larger number of stress measurements than I actually had), resulting in inflated significance estimates.