How to estimate residuals of subgroup with lme4 in R

Question

I'd like to reproduce the results reported in Hoffman & Rovine's work (Multilevel models for the experimental psychologist: foundations and illustrative examples) with lme4 package in R.

In their first example they compared reaction time between elders and young people. Each of their participants have many task trials. So, in the individual level, participants' reaction time were affected by various variables related to their manipulation of trials. In the second level, participants' age and age group would affect participants' reaction time.

In Hoffman's model 2B, they estimate first level residuals for elders and young people respectively, with two dummy variable for young people and old people.

Hoffman's equation is Level1 equation

I'd like to know how to estimate two residuals in lme4 package.

Hoffman's article and example data could be found in Hoffman's website.

I've successfully replicated their result of model 2A, where the variance of young people and old people were assumed to be the same, with the following code.

lmer(lg_rt ~ c_mean + c_sal + (1|Item) + oldage + yrs65 + (1|id), Ex1, REML = F)

this is clear, but still seems a bit "give me the codes". How far have you gotten in trying to implement this for yourself? — Ben Bolker, Dec 23 '17 at 00:11
I've tried to google it and read several books about mixed model in R from library. Maybe I didn't use the right keywords that I can't find solutions in google results. In books, all of them only mention about the situation that variances across groups were equal. — Ju YuJeng, Dec 23 '17 at 01:31
I've tried to add dummy variables of young group and old group but the result turned out not the same as Hoffman's. summary(lmer(lg_rt ~ c_mean + c_sal + oldage + yrs65 + (1|Item) +(1|Young/id)+(1|Oold/id), Ex1, REML = F)) — Ju YuJeng, Dec 23 '17 at 01:38
Unfortunately, I found that it seems only `nlme` could handle the heteroscedasticity issue for now. But, in this model there are two crossed random effect that there is no easy way to handle crossed random effect in `nlme`, compared to `lme4`. — Ju YuJeng, Dec 25 '17 at 02:31
there are ways to handle heteroscedasticity by group in lme4 : https://stackoverflow.com/questions/21409340/how-to-allow-for-factor-specific-variance-of-random-effect-in-lme or, perhaps more easily, in the newer `glmmTMB` package ... — Ben Bolker, Dec 25 '17 at 02:44

score 0 · Answer 1 · answered Feb 10 '21 at 09:38

You can handle heteroscedasticity in lme4 using the modular fitting functions. Here is an example with two groups, which should be extendable to other types of heteroscedasticity. Note that although the weights are estimated, the uncertainty about the weights is not taken into account in the standard errors of the parameters in the final fit. This issue should be possible to solve using the delta method, see e.g. the first equation in Section 2.3.3 of https://10.3102/1076998611417628.

set.seed(1234)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
library(lme4)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack

n <- 100 # number of level-2 units
m <- 3 # number of repeated observations per unit
sd_b <- .3 # random intercept standard deviation
sd_eps1 <- .1 # residual standard deviation in group 1
sd_eps2 <- .3 # residual standard deviation in group 2

# Simulate data
dat <- tibble(
  # unique ID
  id = seq_len(n),
  # explanatory variable, constant over repetitions
  x = runif(n),
  # random intercept
  b = rnorm(n, sd = sd_b),
  # group membership
  grp = sample(1:2, n, replace = TRUE)
) %>% 
  uncount(3) %>% 
  mutate(
    # residual
    eps = rnorm(nrow(.), sd = c(sd_eps1, sd_eps2)[grp]),
    # response, fixed effect is beta=1
    y = x + b + eps
  )

# now optimize over residual weights, fixing the group 1 weight to 1.
# optimize() would be sufficient, but I show it with optim() because it
# then can be directly extended to a larger number of groups
opt <- optim(
  # initial value for group 2 residual relative to group 1
  par = 1,
  fn = function(weight){
    # Compute weights from group variable
    df <- dat %>% 
      mutate(weight = c(1, weight)[grp])
    ## 1.  Parse the data and formula:
    lmod <- lFormula(y ~ x + (1|id), data = df, weights = df$weight)
    ## 2.  Create the deviance function to be optimized:
    devfun <- do.call(mkLmerDevfun, lmod)
    ## 3.  Optimize the deviance function:
    opt <- optimizeLmer(devfun)
    # return the deviance 
    opt$fval
  },
  # Use a method that allows box constraints
  method = "L-BFGS-B",
  # Weight cannot be negative
  lower = 0.01
)

# The weight estimates the following ratio, and it is pretty close
sd_eps1^2/sd_eps2^2
#> [1] 0.1111111
opt$par
#> [1] 0.1035914
# We can now fit the final model at the chosen weights
df <- dat %>% 
  mutate(weight = c(1, opt$par)[grp])
mod <- lmer(y ~ x + (1|id), data = df, weights = df$weight)

# Our estimate of sd_eps1
sigma(mod)
#> [1] 0.09899687
# True value
sd_eps1
#> [1] 0.1
# Our estimate of sd_eps2
sigma(mod) * sqrt(1/opt$par)
#> [1] 0.307581
# True value
sd_eps2
#> [1] 0.3

^{Created on 2021-02-10 by the reprex package (v1.0.0)}

How to estimate residuals of subgroup with lme4 in R

1 Answers1