3
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ Thanks for the answer! Can you think of any reason that this phenomenon of a fixed effect being wiped out by a random effect might depend at all on the number or distribution of bouts per monkey? I only have 1-4 bouts per monkey and I'm wondering if that might be playing a role here. The breakdown of numbers of bouts (N=15) per monkey (N=8) is: 1, 1, 1, 1, 2, 2, 3, 4. Does that add anything to your interpretation of whether or not it's best to include the MonkeyID:BoutID term? Thanks again $\endgroup$
    – Jay
    Aug 20 '19 at 13:22
  • $\begingroup$ If you're saying that you have a total of 15 observations among 8 monkeys, I would (alas) say that this is a very small data set for statistical analysis, let alone a multilevel model ... I'd be tempted to do a simple regression on pitch vs hormone, treating it as an exploratory analysis (e.g. not quoting p-values). You could try a power analysis based on this sample size, but again I would expect it to be very small ... $\endgroup$
    – Ben Bolker
    Aug 20 '19 at 21:11
  • $\begingroup$ Sorry, I should have clarified I have a total of 170 calls (observations) spread out among the 15 bouts (among 8 monkeys). I do only have 15 stress hormone measurements. Does that change your advice? The numbers of calls and bouts are on the low end, but not unheard of in the literature especially for field behavioral studies. Thanks again $\endgroup$
    – Jay
    Aug 20 '19 at 21:31
  • 1
    $\begingroup$ I think collapsing the data to the bout level and weighting by number of calls per bout is a great solution. Thank you again! For the record, the effect of hormone remains nonsignificant when I do this, which makes me think that the results of the LMM excluding the BoutID random effect were false positives stemming from inflated sample size. Bad news for me but it feels good to have an idea one way or the other. $\endgroup$
    – Jay
    Aug 21 '19 at 14:19
  • 1
    $\begingroup$ OK. I like your attitude. Don't forget to try the maximal model (random = ~hormone|MonkeyID); if you have a lot of variation among monkeys in their pitch-hormone relationship, this model will actually increase your statistical power ... $\endgroup$
    – Ben Bolker
    Aug 21 '19 at 15:18
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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