1
$\begingroup$

I am using a GLMM to determine if COVID-19 business closures affected rat activity in the city. The response variable is binomial (no activity/activity), measured in bait stations the council uses to deploy poisoned bait. The model has one continuous predictor DAY, and one categorical predictor Period, and a random intercept for the random factor Station. All predictors are within-Station, but the data is not complete, since stations were not checked every day, but between 5 to every 10 days. The categorical factor has 3 levels (Pre-lockdown, Lockdown, and Post-Lockdown).

Here is what I did:

I fitted a binomial GLMM using glmer from the lme4 package and did model selection using drop1 although no variable where removed. Here is the final model:

M <- glmer(ActivityBi ~  Period + Period:DAY + (1|MBP_ID2),family=binomial(logit), data = (DATA), control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
drop1(M)
summary(M)

This is the output

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: ActivityBi ~ Period + Period:DAY + (1 | MBP_ID2)
   Data: (DATA)
Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))

     AIC      BIC   logLik deviance df.resid 
   15028    15080    -7507    15014    12417 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2887 -0.7699 -0.2530  0.7782  4.7020 

Random effects:
 Groups  Name        Variance Std.Dev.
 MBP_ID2 (Intercept) 1.007    1.004   
Number of obs: 12424, groups:  MBP_ID2, 659

Fixed effects:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               0.1047790  0.0666821   1.571   0.1161    
Period2.Lockdown          8.5385820  1.2521507   6.819 9.16e-12 ***
Period3.Postlockdown     -0.0029200  0.7599282  -0.004   0.9969    
Period1.Prelockdown:DAY   0.0007355  0.0004266   1.724   0.0847 .  
Period2.Lockdown:DAY     -0.0408955  0.0053535  -7.639 2.19e-14 ***
Period3.Postlockdown:DAY -0.0066198  0.0025231  -2.624   0.0087 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Prd2.L Prd3.P P1.P:D P2.L:D
Prd2.Lckdwn -0.033                            
Prd3.Pstlck -0.059  0.005                     
Prd1.Pr:DAY -0.685  0.043  0.069              
Prd2.Lc:DAY -0.002 -0.997 -0.001 -0.002       
Prd3.Ps:DAY -0.002 -0.003 -0.993 -0.002  0.003
convergence code: 0
Model failed to converge with max|grad| = 0.0633013 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?
Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

I then use Anova from car to calculate p-values for the fixe factor and the interactions, and emmeans to calculate the estimates for each level of period

> Anova(M, test="Chisq", type=3)
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: ActivityBi
              Chisq Df Pr(>Chisq)    
(Intercept)  2.4691  1     0.1161    
Period      46.5020  2  7.984e-11 ***
Period:DAY  68.0141  3  1.136e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> marginal = emmeans(M, ~ Period|DAY, type="response") 
NOTE: Results may be misleading due to involvement in interactions
> pairs(marginal2,adjust="tukey")
 contrast                       odds.ratio    SE  df z.ratio p.value
 1.Prelockdown / 2.Lockdown         0.0844  0.04 Inf -5.210  <.0001 
 1.Prelockdown / 3.Postlockdown     2.9291  1.15 Inf  2.728  0.0175 
 2.Lockdown / 3.Postlockdown       34.7132 21.36 Inf  5.766  <.0001 

P value adjustment: tukey method for comparing a family of 3 estimates 
Tests are performed on the log odds ratio scale 
> cld(marginal2,alpha=0.05, Letters=letters, adjust="tukey")
 Period          prob     SE  df asymp.LCL asymp.UCL .group
 3.Postlockdown 0.297 0.0825 Inf     0.141     0.520  a    
 1.Prelockdown  0.553 0.0127 Inf     0.522     0.583   b   
 2.Lockdown     0.936 0.0285 Inf     0.825     0.979    c 

As it can be seen from the model summary and emmeans, rat activity during lockdown was higher than pre-lockdown and post lockdown.

Finally, because there is an effect of DAY I want to calculate estimates for plotting. For that, I used predict() but I keep getting this error:

> EstM<-predict(M, newdata=data, nsim=100, type="response", interval="confidence", re.form=NA, se.fit=TRUE, na.rm = TRUE)
Warning message:
In lme4:::predict.merMod(x, ...) : unused arguments ignored

I then try to use ```predictInterval`` instead

EstM2 <- predictInterval(M, newdata = data, which="fixed", n.sims=1000, type = "linear.prediction")

But the predicted values don't match the model summary, with rat activity pre-lockdown being higher than during lockdown.

enter image description here

Even if I use "probability"

EstM2 <- predictInterval(M, newdata = data, which="fixed", n.sims=1000, type = "probability")

enter image description here

In both cases, pre-lockdown values appear higher than lockdown. Any ideas why this is happening?

$\endgroup$
0

1 Answer 1

2
$\begingroup$

I think the issue here is that there is a dependency between Period and DAY. Since DAY is a covariate, emmeans() reduces it to its average when it determines the reference grid upon which the EMMs are based.

Here is a simpler (and reproducible) example that illustrates this:

set.seed(1.2021)
fake = data.frame(Day = 1:50, 
                  Period = factor(c(rep("A",20), rep("B",10), rep("C",20))),
                  y = rnorm(50) + 100 - 1:50)
plot(y ~ Day, col = Period, data = fake)

Here we can see that there is a linear trend with Day, and the periods comprise the first 20 days, the next 10 days, and the last 20 days, respectively.

Now let's fit a model and obtain the EMMs:

fake.lm = lm(y ~ Period + Day:Period, data = fake)

library(emmeans)

emmeans(fake.lm, "Period")
#> NOTE: Results may be misleading due to involvement in interactions
#>  Period emmean    SE df lower.CL upper.CL
#>  A        75.0 0.531 44     73.9     76.1
#>  B        74.4 0.270 44     73.8     74.9
#>  C        74.5 0.531 44     73.4     75.6
#> 
#> Confidence level used: 0.95

All the EMMs are about 75. This is because we are predicting everything on the same average value of Day:

ref_grid(fake.lm)@grid
#>   Period  Day .wgt.
#> 1      A 25.5    20
#> 2      B 25.5    10
#> 3      C 25.5    20

But now let us account for the dependency between Period and Day by requiring that we compute the average Day for each Period:

ref_grid(fake.lm, cov.reduce = Day ~ Period)@grid
#>   Period  Day .wgt.
#> 1      A 10.5    20
#> 2      B 25.5    10
#> 3      C 40.5    20

With this change, the EMMs reflect the averages we see in the graph:

emmeans(fake.lm, "Period", cov.reduce = Day ~ Period)
#> NOTE: Results may be misleading due to involvement in interactions
#>  Period emmean     SE df lower.CL upper.CL
#>  A       89.69 0.1906 44    89.31    90.07
#>  B       74.37 0.2696 44    73.82    74.91
#>  C       59.63 0.1906 44    59.24    60.01
#> 
#> Confidence level used: 0.95

Created on 2021-01-25 by the reprex package (v0.3.0)

In this example, the Period and interaction effects are very close to zero because the generated y values are just a linear function of Day. That makes all the EMMs in the first set about equal. In your example, there really are Period and interaction effects, making the results more irregular.

$\endgroup$
5
  • $\begingroup$ Thank you for that explanation Russ, that makes sense. What I still don't understand is why the model summary gives estimates that don't seem to reflect the pattern you explain. The estimate for the lockdown period is 8.538 and highly significant compared to pre-lockdown. At the same time, when I calculate the estimates with emmeans as you explained, I get that activity was higher pre-lockdown. $\endgroup$
    – Miguel
    Commented Jan 26, 2021 at 22:15
  • $\begingroup$ Because those estimates have a complicated interpretation. You cannot just look at the main effects and have a simple interpretation. Look at ref_grid(...)@linfct to see what goes into each estimate. $\endgroup$
    – Russ Lenth
    Commented Jan 26, 2021 at 22:19
  • $\begingroup$ @Miguel it seems like you are comparing the intercepts, but what is the meaning/interpretation of intercept according to you? It is not the average of the specific period. $\endgroup$ Commented Jan 27, 2021 at 17:18
  • $\begingroup$ Hi Sextus, I am probably absolutely wrong, and please keep in mind I am not a statistician. I understand the "estimates" from the model summary are the intercepts. The way I interpret them is that they are calculated in reference to the level that is not shown in the summary (In my case, the pre-lockdown period). As I understand it, positive or negative estimates illustrate the direction in the relationship between the level kept fixed and the one shown. In my case Period: 2.Lockdown is 8.54. Since is positive, to me Activity during lockdown is higher than Pre-lockdown. $\endgroup$
    – Miguel
    Commented Jan 27, 2021 at 23:17
  • $\begingroup$ If you're using the default coding in R, only one of the coefficients is an intercept, only one is a slope, and the rest are differences from that intercept and slope. Trying to interpret regression coefficients is time-consuming and not all that useful. You are better off looking at what the model predicts with carious key predictor combinations. In this model it is also confusing because the predictors don't vary independently. Each day can appear in combination with only one period. $\endgroup$
    – Russ Lenth
    Commented Jan 28, 2021 at 1:43

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.