3
$\begingroup$

I could use some help interpreting a glmer output? I am unsure how to get the odds of a non-reference level since there is an interaction. I have observations of whether two female sharks (one young and one old) are resting or not (1s and 0s), and want to know whether the presence of a male changes their resting behavior. I included several variables subsequently, checking with AIC whether the fit improved. The variables:

  • rusten: the response variable, 1s and 0s whether the sharks were resting or not; scored every five minutes during several days across time
  • individu: factor with two levels---young female (jv) and old female (ov)
  • locatie_jm: the location of the male (jm)---Oceaan (where the females are, so 1) or Kleine Oceaan (where the females are not, so 0, the reference)
  • voederdag: whether the sharks were fed (voederdag) or not (vastdag, reference)---feeding occurred on Monday, Wednesday and Friday
  • tijd_13: morning (before 13:00) or afternoon (after 13:00)---I checked different hour splits and this got the best fit (lowest AIC)

Output jv as reference (default)

Summary main model

model_rusten_6_13 <- lme4::glmer(rusten ~ locatie_jm * individu + voederdag + (1|datum) + tijd_13,
                                 data = scans_rusten_zonder_verplaatst, family = binomial)

summary(model_rusten_6_13)

# Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
# Family: binomial  ( logit )
# Formula: rusten ~ locatie_jm * individu + voederdag + (1 | datum) + tijd_13
# Data: scans_rusten_zonder_verplaatst
# 
# AIC      BIC   logLik deviance df.resid 
# 923.7    959.0   -454.9    909.7     1132 
# 
# Scaled residuals: 
#   Min       1Q   Median       3Q      Max 
# -10.9737  -0.3573  -0.1892   0.3858  20.5250 
# 
# Random effects:
#   Groups Name        Variance Std.Dev.
# datum  (Intercept) 2.609    1.615   
# Number of obs: 1139, groups:  datum, 34
# 
# Fixed effects:
#                               Estimate Std. Error z value Pr(>|z|)    
# (Intercept)                  -4.8011     0.5811  -8.262  < 2e-16 ***
# locatie_jmOceaan              1.3194     0.6665   1.980   0.0477 *  
# individuov                    3.0289     0.3114   9.725  < 2e-16 ***
# voederdagvoederdag            1.7316     0.6570   2.636   0.0084 ** 
# tijd_13middag                 3.3027     0.2360  13.992  < 2e-16 ***
# locatie_jmOceaan:individuov  -1.7689     0.3835  -4.613 3.97e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#   (Intr) lct_jO indvdv vdrdgv tjd_13
# locat_jmOcn -0.661                            
# individuov  -0.466  0.301                     
# vodrdgvdrdg -0.363 -0.068  0.041              
# tijd_13mddg -0.329  0.025  0.389  0.120       
# lct_jmOcn:n  0.314 -0.371 -0.754  0.007 -0.182

Confidence Intervals

cc <- exp(confint(model_rusten_6_13,parm="beta_"))
cbind(est=exp(lme4::fixef(model_rusten_6_13)),cc)

#                               est        2.5 %      97.5 %
# (Intercept)                  0.008220788  0.002262453  0.02478207
# locatie_jmOceaan             3.741049456  0.997920767 15.50578399
# individuov                  20.674713236 11.423008423 38.98574776
# voederdagvoederdag           5.649618416  1.496438539 22.88038355
# tijd_13middag               27.186527366 17.342727724 44.05745016
# locatie_jmOceaan:individuov  0.170524558  0.079310269  0.35936984

Interpretation

  • The reference is young female while the male is absent, on not-feeding days in the morning. So the intercept of -4.8 tells me the odds of finding the young female resting are exp(-4.8) = 0.008
    • When the male is present (locatie_jmOceaan = 1), the odds of finding her resting are exp(1.32) = 3.74 times higher
    • In the afternoon (tijd_13 = middag), the odds of finding her resting are exp(3.30) = 27.19 times higher
    • On feeding days (voederdag = voederdag), the odds of finding her resting are exp(1.73) = 5.65 times higher

Output ov as reference

Summary main model

scans_rusten_zonder_verplaatst$individu <- relevel(scans_rusten_zonder_verplaatst$individu, ref = "ov")

model_rusten_6_13_ov <- lme4::glmer(rusten ~ locatie_jm * individu + voederdag + (1|datum) + tijd_13,
                                 data = scans_rusten_zonder_verplaatst, family = binomial)

summary(model_rusten_6_13_ov)

# Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
# Family: binomial  ( logit )
# Formula: rusten ~ locatie_jm * individu + voederdag + (1 | datum) + tijd_13
# Data: scans_rusten_zonder_verplaatst
# 
# AIC      BIC   logLik deviance df.resid 
# 923.7    959.0   -454.9    909.7     1132 
# 
# Scaled residuals: 
#   Min       1Q   Median       3Q      Max 
# -10.9734  -0.3573  -0.1892   0.3858  20.5244 
# 
# Random effects:
#   Groups Name        Variance Std.Dev.
# datum  (Intercept) 2.609    1.615   
# Number of obs: 1139, groups:  datum, 34
# 
# Fixed effects:
#   Estimate Std. Error z value Pr(>|z|)    
# (Intercept)                  -1.7722     0.5159  -3.435 0.000592 ***
#   locatie_jmOceaan             -0.4495     0.6338  -0.709 0.478212    
# individujv                   -3.0289     0.3114  -9.725  < 2e-16 ***
#   voederdagvoederdag            1.7315     0.6569   2.636 0.008398 ** 
#   tijd_13middag                 3.3027     0.2360  13.992  < 2e-16 ***
#   locatie_jmOceaan:individujv   1.7688     0.3835   4.613 3.97e-06 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Correlation of Fixed Effects:
#   (Intr) lct_jO indvdj vdrdgv tjd_13
# locat_jmOcn -0.653                            
# individujv  -0.079  0.139                     
# vodrdgvdrdg -0.384 -0.067 -0.041              
# tijd_13mddg -0.135 -0.084 -0.389  0.120       
# lct_jmOcn:n  0.101 -0.215 -0.754 -0.007  0.182

Confidence Intervals

exp(confint(model_rusten_6_13_ov,parm="beta_"))

# Computing profile confidence intervals ...
# Error in zeta(shiftpar, start = opt[seqpar1][-w]) : 
#     profiling detected new, lower deviance

According to this post, this is apparently caused by having fits that are better, and can be overruled with lme4::confint.merMod and setting the devtol parameter manually. I tried that, but no matter how large I make it, the model gives the same problems (and too large gives other problems). Of course, if I could extract what I need from the first model's summary, I don't need this.

lme4::confint.merMod(model_rusten_6_13_ov, parm="beta_", devtol = 1e-3)
# Computing profile confidence intervals ...
# Error in zeta(shiftpar, start = opt[seqpar1][-w]) : 
#     profiling detected new, lower deviance

lme4::confint.merMod(model_rusten_6_13_ov, parm="beta_", devtol = 1e-3)
# Computing profile confidence intervals ...
# Error in approxfun(obj1[, 2], obj1[, 1]) : 
#   need at least two non-NA values to interpolate
# In addition: There were 50 or more warnings (use warnings() to see the first 50)

Interpretation

  • The old female is the reference (without the male, on non-feeding days in the morning), so the intercept of -1.77 tells me the odds of finding her resting are exp(-1.77) = 0.17
  • for tijd_13 and voederdag, the estimates are similar to the young female, because there are no interactions
  • When the male is present (locatie_jmOceaan = 1), the odds of finding the old female resting are exp(-0.45) = 0.64 times higher

Questions

  1. Are my interpretations of the outputs correct?
  2. how do I get the odds of finding the old female from the first summary? It's tedious to run the model again with another reference (ok with only 2 levels, but not with many)

enter image description here

enter image description here

$\endgroup$

1 Answer 1

1
$\begingroup$

Interpretation of the Intercept includes that it's for the situation when all predictors are at reference levels. My sense is that you recognize this, but other visitors to the site might not. Note that exponentiating the locatie_jmOceaan:individuov interaction coefficient in the first model is the extra odds of the old female resting when the male is present, beyond what you would have predicted from the Intercept, the locatie_jmOceaan coefficient, and the individuov coefficient. To get that, there's no need to re-fit the model with individuov as the reference.

For Wald confidence intervals (which seem pretty close to the profile intervals for the first model, given how wide some intervals are) there's no need to re-fit the model to handle different reference levels of a categorical predictor. All of the information is in the initial model, you just have to extract it properly. You can use the formula for the variance of a sum of correlated variables to estimate the standard error for any linear combination of predictor values, based on the coefficient values and the covariance matrix of the coefficient estimates. To avoid errors in calculation, try the linearHypothesis() function of the car package or the tools in the emmeans package.

Profile-likelihood intervals actually require refitting the model over a range of potential coefficient values, although that might be hidden from you in practice. I don't know why the profile likelihood gave an error with the second model but not with the first.

$\endgroup$
1
  • $\begingroup$ Thanks, good idea to look at packages that help with extracting what I need. I did some manual calculations with a glm and got to the same estimates, which helps with understanding. Indeed, I recognize the intercept estimate is for all variables at reference level. I'll change it. $\endgroup$ Aug 29, 2022 at 8:14

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.