[R-sig-ME] Random intercept model with nested effects

[Tibor Kiss]

Dear list members,

I am struggling with random intercept models with nested random effects. To put it simply: I do not understand why I find influential random effects in a random intercept model that cannot occur with the intercept specification. As an illustration, I am using the _dative_ data set, which is found in _languageR_ and also discussed in Baayen (2009: 278-284).

The two relevant factors are SemanticClass (fixed) and Verb (random). Verbs are sampled from the set of ditransitive verbs, but each verb in this set is a member of one of five semantic classes (SemanticClass = c(a, c, f, p, t). Verbs may also appear in more than one semantic class. In Baayen's (2009) analysis, no nesting is used (perhaps for reasons of exposition), but I assume that explicit nesting must be coded.

In order to make the model somewhat easier to interpret (in comparison to the model in Baayen (2009)), I have changed the levels of two predictors (AccessOfRec and AccessOfTheme) from three to two, and also changed the reference level of SemanticClass to _t_ (which is more plausible given the other reference levels). In addition, I am using a function to extract the random structure and calculate levels of the random effect that stay on one side of 0. The code is below.

The crucial result is shown here: 

##     name_class class names     interc   confint
## 5       feed_t     t  feed -2.0099213 1.6863733
## 8       give_t     t  give -1.1068355 0.3574284
## 12      lend_t     t  lend -2.1997060 1.3972652
## 17     offer_t     t offer -2.2969699 1.7727265
## 18       pay_t     t   pay -2.1214353 0.5253131
## 24      sell_t     t  sell  0.8782523 0.5436892
## 32      take_t     t  take  1.9223862 1.8289943
## 35     write_t     t write  1.8757040 1.3270534
## 40     allow_a     a allow -2.5067738 2.2254017
## 44     bring_a     a bring  1.4675516 1.3846133
## 50        do_a     a    do -2.4082479 1.9457716
## 54      give_a     a  give -1.7403999 0.2911642
## 59     issue_a     a issue  2.7544666 2.3742973
## 65       owe_a     a   owe -2.1404955 1.8719220
## 66       pay_a     a   pay  5.2520666 1.0379571
## 92     teach_c     c teach -2.7920566 1.4288249
## 93      tell_c     c  tell -4.1851336 1.1225658
## 101    offer_f     f offer  2.7622290 1.4753894

Recall that the reference level for SemanticClass is set to _t_. But  in addition to verbs belonging to class _t_, verbs belonging to classes _a_, _c_, and _f_ are listed as influential random effects. This strikes me as mysterious because these verbs cannot occur with the intercept.  

Perhaps there is some logical fallacy in my thinking, but I am not able to identify it and would thus appreciate your comments. 

Thanks in advance


Tibor


# libraries

library(languageR)
library(tidyverse)
library(lme4)

# function to extract relevant random effects

relevant.ranefs <- function(model) {
  randoms <- ranef(model, condVar = T)
  ranefs <- data.frame( names = rownames(randoms[[1]]))
  ranefs$interc = randoms[[1]][,1]
  variances <- attr(randoms[[1]], "postVar")
  ranefs$confint <- 1.96 * sqrt(variances[,,1:length(variances)])
  relevant.ranefs <- subset(ranefs, abs(interc)-confint > 0)

# data manipulation

dative.data <- 
  dative %>%
  mutate(AccessRec = ifelse(AccessOfRec == "given", "given", "non_given"),
         AccessTheme = ifelse(AccessOfTheme == "given", "given", "non_given"),
         SemanticClass = relevel(factor(SemanticClass), ref = "t"),
         AccessRec = relevel(factor(AccessRec), ref = "non_given"),
         DefinOfRec = relevel(factor(DefinOfRec), ref = "indefinite")
         )

# model
dative.glmm3 <- glmer(RealizationOfRecipient ~ AccessRec + AccessTheme + SemanticClass + 
                        PronomOfRec + AnimacyOfRec + AnimacyOfTheme + PronomOfTheme + DefinOfTheme + 
                        DefinOfRec + LengthOfRecipient + LengthOfTheme + Modality +
                        (1| SemanticClass:Verb), 
      data = dative.data, family = "binomial", 
      control = glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=2e4)))

## Dealing with a convergence warning, as suggested by Ben Bolker

params <- getME(dative.glmm3,c("theta","fixef"))

dative.glmm3 <- 
  update(dative.glmm3,start = params,
         control = glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=2e4)))

## Convergence warning gone

print(summary(dative.glmm3), corr = FALSE)

# extraction of relevant random effects

dative3.relevant.ranefs <- relevant.ranefs(dative.glmm3)

dative3.relevant.ranefs<-
  dative3.relevant.ranefs %>%
  separate(names, c("class", "names")) %>%
  unite(name_class, c(names, class), remove = FALSE) %>%
  mutate(name_class = factor(name_class, levels = name_class[order(interc, decreasing = TRUE)]))

dative3.relevant.ranefs


—————————
Prof. Dr. Tibor Kiss
Linguistic Data Science Lab
Ruhr-University Bochum