2
$\begingroup$

I'm using glmms with a negative binomial distribution from the lme4 package. I've run my model (Checked for singularity) and have plotted the fitted values along with the 95% confidence intervals for fixed effects only. However I'm skeptical of my confidence intervals as they don't seem to be overlapping with many of my raw data points. This is data from a repeat-sampling experiment so my raw data does change over time. Still I would like a second opinion on this.

Here is my data (note there's two dataframes, the main one for the model and a secondary one for plotting):

    final<-structure(list(Treatment = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 
2L, 1L), .Label = c("Control", "Treatment"), class = "factor"), 
    Month = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 
    4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
    5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
    6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 
    7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
    9L, 9L, 9L, 9L, 9L, 9L, 9L, 11L, 11L, 11L, 11L, 11L, 11L, 
    11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L), .Label = c("January", 
    "February", "March", "April", "May", "June", "July", "August", 
    "September", "October", "November", "December"), class = c("ordered", 
    "factor")), Occurrence = c(3L, 1L, 6L, 3L, 3L, 5L, 0L, 4L, 
    0L, 3L, 5L, 1L, 2L, 3L, 8L, 3L, 3L, 4L, 2L, 7L, 5L, 15L, 
    4L, 2L, 4L, 1L, 2L, 5L, 4L, 4L, 5L, 10L, 2L, 2L, 11L, 8L, 
    7L, 6L, 7L, 14L, 1L, 4L, 1L, 2L, 4L, 0L, 11L, 2L, 0L, 5L, 
    5L, 4L, 3L, 3L, 12L, 8L, 7L, 7L, 10L, 13L, 5L, 4L, 1L, 4L, 
    2L, 8L, 1L, 2L, 4L, 3L, 6L, 8L, 2L, 2L, 8L, 4L, 8L, 11L, 
    4L, 8L, 4L, 4L, 2L, 1L, 9L, 10L, 2L, 5L, 1L, 5L, 11L, 11L, 
    5L, 11L, 10L, 3L, 9L, 6L, 7L, 8L, 4L, 4L, 2L, 1L, 6L, 11L, 
    4L, 5L, 7L, 2L, 8L, 12L, 10L, 9L, 8L, 8L, 8L, 4L, 9L, 12L, 
    2L, 4L, 1L, 1L, 9L, 5L, 5L, 1L, 9L, 6L, 9L, 10L, 11L, 5L, 
    13L, 3L, 8L, 2L, 4L, 2L, 1L, 2L, 8L, 3L, 11L, 10L, 8L, 9L, 
    3L, 5L, 0L, 4L), Occurrence_scale = c(-2.02998961901208, 
    0.408833764946554, -1.21704849102587, -0.404107363039658, 
    0.00236320095344811, 0.408833764946554, 0.408833764946554, 
    1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 0.00236320095344811, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.00236320095344811, 0.408833764946554, 
    0.408833764946554, 1.22177489293277, -2.02998961901208, 0.00236320095344811, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 2.03471602091898, -0.404107363039658, 
    1.62824545692587, -2.02998961901208, 0.408833764946554, -1.21704849102587, 
    -0.404107363039658, 0.408833764946554, -2.02998961901208, 
    0.00236320095344811, -0.404107363039658, 0.408833764946554, 
    -0.404107363039658, 0.408833764946554, -0.404107363039658, 
    0.408833764946554, 0.00236320095344811, 1.62824545692587)), class = "data.frame", row.names = c(NA, 
-152L))

ref_table<-structure(c(2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
             1L, 1L, 2L, 1L, 2L, 1L, 1L), .Label = c("Control", "Treatment"
             ), class = "factor")

Here is the code that I used for the model and plot (the plotting code comes from Ben Bolker's code):

library(lme4)
model_1<-glmer.nb(data=final, Occurrence ~ Occurrence_scale + Treatment + (1|Month))
isSingular(model_1)
summary(model_1)
newdat<-data.frame(Treatment=rep(ref_table, 20), 
                   Occurrence_scale = seq(min(final$Occurrence_scale), max(final$Occurrence_scale),
                                                length.out = 200)) #make new data

ta2mm<-model.matrix(~Occurrence_scale + Treatment,newdat) #cereate a model matrix of you paramters from the dummy data
y<-ta2mm%*%fixef(model_1) #multiply the m2m2 matrix with the model coefficients


pvar1 <- diag(ta2mm %*% tcrossprod(vcov(model_1),ta2mm))  #Take the cross-product of the transpose of a matrix
tvar1 <- pvar1+VarCorr(model_1)$Month[1]

newdat <- data.frame(
  Treatment=newdat$Treatment,Ref = newdat$Occurrence_scale,
  y=exp(y),
  plo = exp(y-1.96*sqrt(pvar1))
  , phi = exp(y+1.96*sqrt(pvar1))
  , tlo = exp(y-1.96*sqrt(tvar1))
  , thi = exp(y+1.96*sqrt(tvar1))
)

###graph regression line
newdat
graph_1<-ggplot(NULL) + 
  geom_ribbon(data = newdat,aes(x = Ref, ymin=plo,ymax = phi, fill = Treatment), alpha = 0.4) +
  geom_point(data=final, aes(x = final$Occurrence_scale, y = final$Occurrence, color = Treatment), position = position_jitter(),
             alpha = 0.5,size=5) +
  geom_line(data = newdat,aes(x = Ref, y = y, lty = Treatment)) +
  theme_classic() + 
  labs(x = expression(bold("Initial sampling occurrence (scaled)")), 
       y = expression(bold("Occurrence"))) +
  theme(axis.text.x = element_text(size=12), 
        axis.text.y = element_text(size=12), axis.title = element_text(size = 14),
        legend.position = 'right')+  scale_color_manual(values=c("#000000", "#2088f7"), 
                                                        name="Treatment Type",
                                                        labels=c("Control", "Treated")) +
  scale_fill_manual(values=c("#000000", "#2088f7"), 
                    name="Treatment Type",
                    labels=c("Control", "Treated")) 
graph_1<-graph_1 + theme(axis.line = element_blank(), panel.grid.major.y = element_line(color='grey'))
graph_1

To add some context. I ran an experiment over the course of nearly a year where I reduced predator occurrences in ten plots while kept the other 10 as controls. I wanted to statistically show that the treatment (reducing predators) plots had less predators than control plots, given that we were using a novel reduction field method. So my response variable was the occurrence of predators as a function of a covariate of initial occurrences of these predators prior to any reduction attempts (this was to account for natural variation in predator occurrence between plots), the fixed effect of treatment (two levels, control or treated), and the random intercept of months to account for time. It's important to note that I scaled the covariate for the model hence the name "occurrence_scale".

I plotted the fixed effects of the model with my x-axis being the covariate and also plotted the raw data. However my 95% Confidence intervals are rather narrow and don't really overlap with the raw data which concerns me and makes me doubt whether or not I'm doing this right. I realize 95% CIs don't necessarily have to overlap with raw data as they are estimated around the mean of the model rather than the mean of the data. However, I'm still rather skeptical of my plot. Any help would be greatly appreciated!

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

In the solution you suggested it was not clear to me why you calculate the sum of the standard errors of the fitted values with the variance of the random intercepts, i.e., this line: tvar1 <- pvar1+VarCorr(model_1)$Month[1].

You could try fitting the model with GLMMadaptive that calculates confidence intervals for the fixed effects relationship using the following code:

library("GLMMadaptive")
model_1 <- mixed_model(Occurrence ~ Occurrence_scale + Treatment, random = ~ 1 | Month,
                       data = final, family = negative.binomial())
summary(model_1)

newdat <- with(final, 
               expand.grid(Treatment = levels(Treatment),
                           Occurrence_scale = seq(min(Occurrence_scale),
                                                  max(Occurrence_scale),
                                                  length.out = 200)))
plot_data <- effectPlotData(model_1, newdat)

library("lattice")

xyplot(exp(pred) + exp(low) + exp(upp) ~ Occurrence_scale | Treatment, data = plot_data,
       lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2, type = "l")
Improve this answer
$\endgroup$
5
  • $\begingroup$ Hi Dimitirs. Thanks for the suggestion. I calculated summed the stand. errors with the variance of the random intercept in case I wanted my confidence intervals to reflect the uncertainty of the fixed effects and the variance of the random effects (this was an option that Ben Bolker offered as a way to visualize both and random and fixed effects). However in my final plot I never end up using it. I tried your alternative solution and came out with the same confidence intervals as my solution. So I guess my question still remains on whether my CIs are right? $\endgroup$
    – antR
    Commented Jun 4, 2019 at 13:05
  • $\begingroup$ These will be indeed the asymptotic 95% confidence intervals for the average occurrence in the grid of values you have specified in newdat. $\endgroup$
    – Dimitris Rizopoulos
    Commented Jun 4, 2019 at 13:08
  • $\begingroup$ Interesting. In the past I have never seen confidence intervals that overlapped very little with the actual raw data. Is this because the confidence intervals are based on off of the mean of the model's predictions and not the means of the actual data? $\endgroup$
    – antR
    Commented Jun 4, 2019 at 13:11
  • $\begingroup$ The confidence interval is for the mean occurrence. In a prediction interval you would expect 95% of the raw data to fall within. $\endgroup$
    – Dimitris Rizopoulos
    Commented Jun 4, 2019 at 13:15
  • $\begingroup$ Understood. Is it odd how less than 50% of the raw data fits within these intervals though? The r-squared (both marginal and conditional) for this model is rather low (~0.2) so that may be another reason for the lack of overlap between raw data and confidence intervals? $\endgroup$
    – antR
    Commented Jun 4, 2019 at 13: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.