3
$\begingroup$

Study background

My research question looks at the effect of age group (AgeGr) on gazing. Each infant was observed for 1h and signals with gazing (GazingY) and without (GazingN) were coded. The total number of signals was different for each infant. I amalgamated all events for each infant (e.g., INF_1 = 30 signals, 10 gaze, 20 no gaze) and fitted a binomial GLMM using glmer from the lme4 package. Hypothesis is that the odds (probability?) of GazingY will be higher for the older age groups (G2, G3). I need some advice if I am doing this correctly or not, and if there is anything missing.

Example of data

    ID Sex AgeGr Total GazingY GazingN
 INF_1   M    G1     3       1       2
 INF_2   F    G1     5       2       3
 INF_3   F    G2    31      23       8
 INF_4   M    G2    69      52      17
 INF_5   M    G2     6       2       4
 INF_6   F    G3     8       2       6
 INF_7   M    G3    55      28      27
 INF_8   M    G3     8       2       6
 INF_9   F    G1    20      16       4
 INF_10  M    G2     9       7       2

Model

mod.1 <- glmer(cbind(GazingY,GazingN) ~ AgeGr + Sex + (1|ID), 
               data=mydata,family="binomial")

Output

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
glmerMod]
 Family: binomial  ( logit )
Formula: cbind(GazingY, GazingN) ~ AgeGr + Sex + (1 | ID)
   Data: mydata

     AIC      BIC   logLik deviance df.resid 
   278.2    289.3   -134.1    268.2       62 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3337 -0.4715  0.1102  0.5356  1.2726 

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.5807   0.762   
Number of obs: 67, groups:  ID, 67

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.5050     0.2751   1.836  0.06639 .  
AgeGrG2       0.8301     0.3127   2.655  0.00794 ** 
AgeGrG3       2.3259     0.3672   6.335 2.38e-10 ***
SexM          0.1574     0.2790   0.564  0.57262    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
        (Intr) AgGrG2 AgGrG3
AgeGrG2 -0.606              
AgeGrG3 -0.553  0.414       
SexM    -0.681  0.203  0.234

My main question is whether family = binomial is a correct choice here. I tried to assess the distribution of GazingY using fitdistrplus package. The Collen and Frey graph isn't very informative as binomial disttribution is not included. I also used plotdist as plotdist(mydata$GazingY, histo = TRUE, demp = TRUE), but again was unsure of the output.

Fig.1

Fig.2

I then used DHARMa package to examine the patterns in residuals in mod.1 which didn't pick up any major issues.

enter image description here

$\endgroup$
6
  • $\begingroup$ You might have a question whether the observed variances match the mathematically required ones for the binomial distribution, but that's it. There's no reason to look at a Cullen & Frey graph. All of that is just a distraction. $\endgroup$ Apr 1, 2021 at 14:12
  • 1
    $\begingroup$ If you have the actual ages of the infants, you would probably be better off using those than grouping the infants into age categories. $\endgroup$ Apr 1, 2021 at 14:13
  • $\begingroup$ Is that your entire dataset? $\endgroup$ Apr 1, 2021 at 14:18
  • $\begingroup$ Thank you @gung-ReinstateMonica for your help. My entire data set is based on 67 infants, I just provided a sample for the first 10. I can include the rest if that would help. $\endgroup$ Apr 2, 2021 at 8:38
  • $\begingroup$ Thank you for your suggestions. I'll have to learn a bit more about how to check the observed variances and if they mathematically match the binomial distribution (any sugegstions would be appreciated how to do this). I don't have enough infants per specific age/month (sometimes only one subject) so there were grouped into age categories. $\endgroup$ Apr 2, 2021 at 8:58

1 Answer 1

2
$\begingroup$

You don't have to check the same assumptions for a logistic regression model as you do for a linear regression (OLS) model. A lot of what you're doing seems unnecessary. Even running a GLMM seems like overkill to me. You know the response variable is a binomial, because the infant will either gaze or not with some probability. The main question is whether the probabilities are overdispersed. I would simply run a regular logistic regression, but use the quasibinomial option to allow for overdispersion. The key here is the line "Dispersion parameter ... taken to be 2.371701"; this implies overdispersion (the 'correct' amount of dispersion would be 1).

mydata = read.table(text="    ID Sex AgeGr Total GazingY GazingN
                           INF_1   M    G1     3       1       2
                           INF_2   F    G1     5       2       3
                           INF_3   F    G2    31      23       8
                           INF_4   M    G2    69      52      17
                           INF_5   M    G2     6       2       4
                           INF_6   F    G3     8       2       6
                           INF_7   M    G3    55      28      27
                           INF_8   M    G3     8       2       6
                           INF_9   F    G1    20      16       4
                           INF_10  M    G2     9       7       2", header=TRUE)

## here I make up actual ages
set.seed(5998)  # this makes the example exactly reproducible
for(i in 1:10){
  mydata$age[i] = with(mydata, 
                       ifelse(AgeGr[i]=="G1", runif(1, min=3, max=6),
                       ifelse(AgeGr[i]=="G2", runif(1, min=6, max=9),
                                           runif(1, min=9, max=12))))
}
mydata$age = floor(mydata$age)
print(mydata, row.names=FALSE)
#     ID Sex AgeGr Total GazingY GazingN age
#  INF_1   M    G1     3       1       2   5
#  INF_2   F    G1     5       2       3   5
#  INF_3   F    G2    31      23       8   8
#  INF_4   M    G2    69      52      17   6
#  INF_5   M    G2     6       2       4   7
#  INF_6   F    G3     8       2       6  10
#  INF_7   M    G3    55      28      27  10
#  INF_8   M    G3     8       2       6  10
#  INF_9   F    G1    20      16       4   3
# INF_10   M    G2     9       7       2   7

## let's fit the model & see what we find
mod.1 = glm(cbind(GazingY, GazingN)~age+Sex, mydata, family=quasibinomial)
summary(mod.1)
#Call:
#glm(formula = cbind(GazingY, GazingN) ~ age + Sex, family = quasibinomial, 
#    data = mydata)
#
#Deviance Residuals: 
#    Min       1Q   Median       3Q      Max  
#-1.6638  -1.4959  -0.8585   0.6467   1.6010  
#
#Coefficients:
#              Estimate Std. Error t value Pr(>|t|)  
#(Intercept)  2.2395126  0.8667995   2.584   0.0363 *
#age         -0.2264747  0.1079546  -2.098   0.0741 .
#SexM        -0.0002036  0.5093888   0.000   0.9997  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#(Dispersion parameter for quasibinomial family taken to be 2.371701)
#
#    Null deviance: 27.428  on 9  degrees of freedom
#Residual deviance: 15.773  on 7  degrees of freedom
#AIC: NA
#
#Number of Fisher Scoring iterations: 4

## to understand the model, let's plot it
library(binom)  # we'll use this package to get CI's (below)
cis   = with(mydata, binom.confint(GazingY, Total, method="exact"))
print(cis, row.names=FALSE)
# method  x  n      mean       lower     upper
#  exact  1  3 0.3333333 0.008403759 0.9057007
#  exact  2  5 0.4000000 0.052744951 0.8533672
#  exact 23 31 0.7419355 0.553866066 0.8814360
#  exact 52 69 0.7536232 0.635095088 0.8494503
#  exact  2  6 0.3333333 0.043271868 0.7772219
#  exact  2  8 0.2500000 0.031854026 0.6508558
#  exact 28 55 0.5090909 0.370705345 0.6464638
#  exact  2  8 0.2500000 0.031854026 0.6508558
#  exact 16 20 0.8000000 0.563385997 0.9426660
#  exact  7  9 0.7777778 0.399906426 0.9718550

xseq  = seq(from=0, to=12, by=.1)  
mpred = predict(mod.1, data.frame(Sex="M", age=xseq), type="response")
fpred = predict(mod.1, data.frame(Sex="F", age=xseq), type="response")
windows()
  with(mydata, plot(age, GazingY/Total, ylim=c(0,1), xlim=c(0,12),
                    pch=ifelse(Sex=="M",15,16), 
                    col=ifelse(Sex=="M","blue","red")))
  with(cis, arrows(x0=mydata$age, y0=lower, y1=upper, 
                   code=3, angle=90, length=.1))
  lines(xseq, mpred, col="blue", lwd=3)
  lines(xseq, fpred, col="red",  lwd=3, lty="dashed")

First, I read in the dataset. Then I made up ages for the infants (I assume you have the real ones, if they've been thrown away, you'll have to stick with the age groupings.) Then I fit a logistic regression model with the odds of gazing (number of gazing trials to non-gazing trials) as the response—in a way, you really have only $1$ data point per infant (i.e., the odds of gazing), which has been estimated from multiple trials. I use the quasibinomial family for the model. That means the model will estimate the amount the dispersion of the observed odds relative to what is mathematically expected. In fact, the observations vary around the fitted line by more than twice ($2.4\times$) as much as they 'should'. As a result, the standard errors are enlarged, and the variables are less significant to take the additional uncertainty into account.

To understand what the model is showing, I plot it over the data. I make a scatterplot of the observed proportion gazing by the infant's age in months. I use the binom package to compute exact 95% confidence intervals for each proportion, that I plot with error bars. Then I plot the model over it. I plot two lines, one (solid blue) for males and one (dashed red) for females. In this case, the lines are on top of each other, because the model finds essentially no difference between males and females.

scatterplot of proportion gazing by age with model overlaid

$\endgroup$
3
  • $\begingroup$ Thank you @gung-ReinstateMonica, you've been a tremendous help. With your support, I was able to get the results that support my hypothesis without unnecessary complicating the analysis. If it's not too much to ask, could you please briefly explain what the plot is depicting? Thank you in advance. $\endgroup$ Apr 3, 2021 at 8:01
  • $\begingroup$ You're welcome, @VioletGibson. $\endgroup$ Apr 3, 2021 at 12:24
  • $\begingroup$ Note that the OP was already using an OLRE (1|ID) that accounts for the overdispersion, which is likely why the DHARMa tests didn't flag overdispersion. Whether this is preferable of the quasi distribution is a matter of debate ... $\endgroup$ Apr 6, 2021 at 21:22

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.