0
$\begingroup$

I have a dataset of measurements of "Y" at different locations, and I am trying to determine how variable Y is influenced by variables A, B, and D. I also have another variable, C, that may influence the data and while I am not interested in the variable for the questions I am asking, I can include it in the model if it improves the fit. Importantly, because Y was measured at different locations (sampling areas), I want to include this variable as a random effect in my model. I have tried various models but none seem to have a good fit, and because I'm so new to GLMs/GLMMs, I would appreciate some guidance.

Here is an example of my data:

table <- "   ID location A         B      C     D       Y
1   1       AA 0 0.6181587 -29.67 14.14 168.041
2   2       AA 1 0.5816176 -29.42 14.21 200.991
3   3       AA 2 0.4289670 -28.57 13.55 200.343
4   4       AA 3 0.4158891 -28.59 12.68 215.638
5   5       AA 4 0.3172721 -28.74 12.28 173.299
6   6       AA 5 0.1540603 -27.86 14.01 104.246
7   7       AA 6 0.1219355 -27.18 14.43 128.141
8   8       AA 7 0.1016643 -26.86 13.75 179.330
9   9       BB 0 0.6831649 -28.93 17.03 210.066
10 10       BB 1 0.6796935 -28.54 18.31 280.249
11 11       BB 2 0.5497743 -27.88 17.33 134.023
12 12       BB 3 0.3631052 -27.48 16.79 142.383
13 13       BB 4 0.3875498 -26.98 17.81 136.647
14 14       BB 5 0.3883785 -26.71 17.56 142.179
15 15       BB 6 0.4058061 -26.72 17.71 109.826
16 16       CC 0 0.8647298 -28.53 11.93 220.464
17 17       CC 1 0.8664036 -28.39 11.59 326.868
18 18       CC 2 0.7480748 -27.61 11.75 322.745
19 19       CC 3 0.5959143 -26.81 13.27 170.064
20 20       CC 4 0.4849077 -26.77 14.68 118.092
21 21       CC 5 0.3584687 -26.65 15.65  95.512
22 22       CC 6 0.3018285 -26.33 16.11  71.717
23 23       CC 7 0.2629121 -26.39 16.16  60.052
24 24       DD 0 0.8673077 -27.93 12.09 234.244
25 25       DD 1 0.8226558 -27.96 12.13 244.903
26 26       DD 2 0.7826429 -27.44 12.38 252.485
27 27       DD 3 0.6620447 -27.23 13.84 150.886
28 28       DD 4 0.4453213 -27.03 15.73 102.787
29 29       DD 5 0.3720257 -27.13 16.27 109.201
30 30       DD 6 0.6040217 -27.79 16.41 101.509
31 31       EE 0 0.8770987 -28.62 12.72 239.036
32 32       EE 1 0.8504547 -28.47 12.92 220.600
33 33       EE 2 0.8329484 -28.45 12.94 174.979
34 34       EE 3 0.8181102 -28.37 13.17 138.412
35 35       EE 4 0.7942685 -28.32 13.69 121.330
36 36       EE 5 0.7319724 -28.22 14.62 111.851
37 37       EE 6 0.7014828 -28.24 15.04 110.447
38 38       EE 7 0.7286984 -28.15 15.18 121.831"

#Create a dataframe with the above table
df <- read.table(text=table, header = TRUE)
df

# Make sure location is a factor
df$location<-as.factor(df$location)


# Visualize data
# Single variables at once
plot(Y~A,data=df)
plot(Y~B,data=df)
plot(Y~D,data=df)

# Examine potential differences between locations
ggplot(df, aes(y = Y, x = A)) +
  geom_point() +
  #  scale_x_log10() +
  facet_grid(location ~ .)

Here are the first models that I tried:

# Load libraries
library(ggplot2)
library(pscl)
library(lmtest)
library(lme4)
library(glmmTMB)
library(car)

# Dispersion function
dispersion<-function(model,modeltype='gaussian')
{
  A<-sum(resid(model,type="pearson")^2)
  if(modeltype %in% c("g","p","qp","gaussian","poisson","quasipoisson"))
  {
    B<-length(resid(model))-length(coef(model))
  }
  if(modeltype %in% c("nb","negativebinomial"))
  {
    B<-length(resid(model))-(length(coef(model))+1)
  }
  if(modeltype %in% c("zpoisson","zp"))
  {
    B<-summary(model)$df.resid
  }
  if(modeltype %in% c("znegativebinomial","znb"))
  {
B<-summary(model)$df.resid + 1
  }
  DISP<<-A/B
  return(DISP)
}

# Try a regular GLM first
# Model0
M0 <- glm(Y ~ A + B + D, data = df)
summary(M0)
dispersion(M0)
#1784.75
#High overdispersion

# Plot Residuals vs Fitted Values
E0<-resid(M0)
F0<-fitted(M0)
x0<-data.frame(E0,F0)
ggplot(x0)+
  geom_point(aes(x=F0, y=E0))+
  geom_hline(yintercept=0, linetype='dashed',col='red')+
  theme_bw(22)+ylab("Residuals")+xlab("Fitted Values")


# Try random intercept model without interactions
# Mixed Model1: Random Intercept Model
M1 <- glmer(Y ~ A + B + D + (1|location), data = df) #(1|location) is random intercept
summary(M1)
dispersion(M1)
#1111.311
#High overdispersion

# Plot Residuals vs Fitted Values #
E1<-resid(M1)
F1<-fitted(M1)
x1<-data.frame(E1,F1)
ggplot(x1)+
  geom_point(aes(x=F1, y=E1))+
  geom_hline(yintercept=0, linetype='dashed',col='red')+
  theme_bw(22)+ylab("Residuals")+xlab("Fitted Values")

# Try random intercept model with interactions
# Mixed Model2: Random Intercept Model with interactions
M2 <- glmer(Y ~ A * B * D * (1|location), data = df) #(1|location) is random intercept
summary(M2)
dispersion(M2)
#942.4339
#Less overdispersion than M0 and M1, but still high

# Plot Residuals vs Fitted Values
E2<-resid(M2)
F2<-fitted(M2)
x2<-data.frame(E1,F1)
ggplot(x2)+
  geom_point(aes(x=F2, y=E2))+
  geom_hline(yintercept=0, linetype='dashed',col='red')+
  theme_bw(22)+ylab("Residuals")+xlab("Fitted Values")


# Try negative binomial with interactions
# Mixed Model3: Negative Binomial Random Intercept Model with interactions
M3 <- glmer.nb(Y ~ A * B * D * (1|location), data = df) #(1|location) is random intercept
summary(M3)
dispersion(M3, modeltype = 'nb')
#0.9195665
#Dispersion statistic is closer to 1 than M0, M1 and M2, slight underdispersion

# Plot Residuals vs Fitted Values
E3<-resid(M3)
F3<-fitted(M3)
x3<-data.frame(E3,F3)
ggplot(x3)+
  geom_point(aes(x=F3, y=E3))+
  geom_hline(yintercept=0, linetype='dashed',col='red')+
  theme_bw(22)+ylab("Residuals")+xlab("Fitted Values")


# Try negative binomial without interactions
# Mixed Model4: Negative Binomial Random Intercept Model
M4 <- glmer.nb(Y ~ A + B + D + (1|location), data = df) #(1|location) is random intercept
summary(M4)
dispersion(M4, modeltype = 'nb')
#1.00957
#Good dispersion statistic

# Plot Residuals vs Fitted Values
E4<-resid(M4)
F4<-fitted(M4)
x4<-data.frame(E4,F4)
ggplot(x4)+
  geom_point(aes(x=F4, y=E4))+
  geom_hline(yintercept=0, linetype='dashed',col='red')+
  theme_bw(22)+ylab("Residuals")+xlab("Fitted Values")

After this, I compared these models (output commented out so you don't have to run the code):

# Based on dispersion statistic (above), model 4 is likely the best fit

# AIC
AIC(M0, M1, M2, M3, M4)
# df      AIC
# M0  5 398.1200
# M1  6 374.1251
# M2 10 351.8459
# M3 10 391.7005
# M4  6 389.9635
#Based on AIC, model 2 is best fit


# Likelihood ratio tests
#Model 0 and model 1
lrtest(M0, M1)
# Model 1: Y ~ A + B + D
# Model 2: Y ~ A + B + D + (1 | location)
# #Df  LogLik Df  Chisq Pr(>Chisq)    
# 1   5 -194.06                         
# 2   6 -181.06  1 25.995  3.423e-07 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 1 fits better than model 0

#Model 1 and model 2
lrtest(M1, M2)
# Model 1: Y ~ A + B + D + (1 | location)
# Model 2: Y ~ A * B * D * (1 | location)
# #Df  LogLik Df  Chisq Pr(>Chisq)    
# 1   6 -181.06                         
# 2  10 -165.92  4 30.279  4.294e-06 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 2 fits better than model 1

#Model 2 and model 3
lrtest(M2, M3)
# Model 1: Y ~ A * B * D * (1 | location)
# Model 2: Y ~ A * B * D * (1 | location)
# #Df  LogLik Df  Chisq Pr(>Chisq)    
# 1  10 -165.92                         
# 2  10 -185.85  0 39.855  < 2.2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 2 fits better than model 3

#Model 2 and model 4
lrtest(M2, M4)
# Model 1: Y ~ A * B * D * (1 | location)
# Model 2: Y ~ A + B + D + (1 | location)
# #Df  LogLik Df  Chisq Pr(>Chisq)    
# 1  10 -165.92                         
# 2   6 -188.98 -4 46.118  2.328e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 2 fits better than model 4

#Model 3 and model 4
lrtest(M3, M4)
# Model 1: Y ~ A * B * D * (1 | location)
# Model 2: Y ~ A + B + D + (1 | location)
# #Df  LogLik Df Chisq Pr(>Chisq)
# 1  10 -185.85                    
# 2   6 -188.98 -4 6.263     0.1803
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 3 fits better than model 4

# Based on likelihood ratio test, model 2 is best fit

As you can see, model 2 seems to be the best fit based on AIC and Likelihood ratio test, but there is lots of overdispersion, which led me to believe I should keep looking for a better fit. In the negative binomial models, the dispersion was good but the fit was not, so next, I tried:

# Model 2 is still quite overdispersed, so lets try adding another variable, C
# Mixed Model5: Random Intercept Model with interactions ----
M5 <- glmer(Y ~ A * B * D * C * (1|location), data = df) #(1|location) is random intercept
summary(M5)
dispersion(M5)
#716.2611
#Less overdispersion than M2, but still high


# We can also try adding the other variable to the negative binomial, since that had the best dispersion statistic
# Mixed Model6: Negative Binomial Random Intercept Model
M6 <- glmer.nb(Y ~ A + B + D + C + (1|location), data = df) #(1|location) is random intercept
summary(M6)
dispersion(M6, modeltype = 'nb')
#1.004846
#Good dispersion statistic

Then I compared these models with the two strongest ones from before:

AIC(M2, M4, M5, M6)
# M2 10 351.8459
# M4  6 389.9635
# M5 18 295.1666
# M6  7 391.6009
#Based on AIC, model 5 is best fit

lrtest(M5, M6)
# Model 1: Y ~ A * B * D * C * (1 | location)
# Model 2: Y ~ A + B + D + C + (1 | location)
# #Df  LogLik  Df  Chisq Pr(>Chisq)    
# 1  18 -129.58                          
# 2   7 -188.80 -11 118.43  < 2.2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 5 fits better than model 6

lrtest(M2, M5)
# Model 1: Y ~ A * B * D * (1 | location)
# Model 2: Y ~ A * B * D * C * (1 | location)
# #Df  LogLik Df  Chisq Pr(>Chisq)    
# 1  10 -165.92                         
# 2  18 -129.58  8 72.679  1.436e-12 ***
#   ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Model 5 fits better than model 2

So this leads me to believe that Model 5 is the best fit, despite the high overdispersion.

I've tried running Anova() from the package car on some of these, and the results vary between them. For example, Model 5 says that Y is not influenced by A, B, C or D, whereas Model 2 shows that B and D influence Y.

AnovaM5 <- Anova(M5)
AnovaM5

AnovaM2 <- Anova(M2)
AnovaM2

Mixed models and random effects are extremely new to me, and while I've been reading a lot of forums and introductory books/websites, I am finding it difficult to grasp. Any help with this scenario would be greatly appreciated.

$\endgroup$

1 Answer 1

3
$\begingroup$

I think you're getting way ahead of yourself by looking at overdispersion.

The way to choose the distribution of your model is by looking at your data, specifically the distribution of your response variable. Right away we can see that negative binomial isn't appropriate for this data because that's a distribution for count data, but your df$Y can't be counts because it is is not whole numbers.

You've used the glm and glmer functions to estimate them, but your $M0$, $M1$, $M2$, and $M5$ could be estimated with lm and lmer. That's because the default response type is Gaussian, and you haven't specified an alternative with the family argument. There is no concept of "overdispersion" with the Gaussian family, since the dispersion is estimated by the residual variance (and so need not take any particular value).

Since these are linear models and not generalized linear models, the best place to start with diagnostics are the built-in diagnostic plots. I've plotted them here:

layout( matrix(1:4, 2, 2) )
plot(M0)

Model diagnostic plots showing a curved shape to the fitted versus residual line

There is a curved shape to the "Fitted vs Residual" line, and the Q-Q plot shows that the right tail is heavier than expected. So this model isn't ready for prime time. But the curve shape suggests that there is a non-linear relationship with one or more of your input variables. Let's look at the data itself (note that this should always be the first step, but since I am starting from your models rather than from the beginning, I have left it until now.)

plot(df[, 3:7])

scatterplot of all variable combinations of the data

Here you can see that df$Y has a nonlinear relationship with df$D. The pattern of dots looks vaguely like a parabola, so I'll add a second-order polynomial on D to the model for Y and then create the diagnostic plots:

linmod1 = lm(Y ~ A + B + poly(D, 2), data = df)
plot(linmod1)

diagnostic plots for the new model, which has repaired the flaws in the previous version

That looks pretty good. I'd use this model going forward. It can be your starting point for deciding whether to add a random effect for location.

> summary(linmod1)

Call:
lm(formula = Y ~ A + B + poly(D, 2), data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-58.915 -24.022  -1.213  20.535  65.998 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  207.446     24.920   8.324 1.29e-09 ***
A            -11.276      3.309  -3.407 0.001743 ** 
B             -5.412     30.070  -0.180 0.858265    
poly(D, 2)1 -158.561     36.283  -4.370 0.000116 ***
poly(D, 2)2  194.337     35.902   5.413 5.45e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 31.21 on 33 degrees of freedom
Multiple R-squared:  0.8008,    Adjusted R-squared:  0.7767 
F-statistic: 33.17 on 4 and 33 DF,  p-value: 3.892e-11
$\endgroup$
1
  • $\begingroup$ Thank you! This explanation helped make it more clear to me, as you can tell this is still quite new for me so I am a bit lost. This helps, I appreciate it! $\endgroup$
    – cgxytf
    Commented Sep 28, 2021 at 16: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.