unbalanced crossed random effects in lme4, lme4, gamlss

crossed_gamlss.R

library(lme4)
library(nlme)
set.seed(101)
## want to sample a *mostly* balanced design?
np <- nrow(Penicillin)
pu <- Penicillin[sample(np, size=round(np/2)),]
pu$dummy <- 1

## WILL NOT WORK without a "grouped data" structure
## set up a dummy grouping variable for the top level (necessary)
## the residual variance and the variance attributed to the dummy will be confounded/unidentifiable,
## but the rest of the model (fixed coefficients, other variance components, etc.) should be OK 
## (but this could complicate getting confidence intervals on variance components ...)

pg <- groupedData(diameter ~ 1| dummy, pu)

## model spec follows example from p. 165 of Pinheiro & Bates:
## https://books.google.ca/books?ei=PBB2T67wN6X20gH7zLzPDQ&id=y54QDUTmvDcC&dq=pinheiro+bates&q=crossed&hl=en#v=snippet&q=crossed&f=false

m1G <- gamlss(diameter ~ 1 + re(random=pdBlocked(list(pdIdent(~1),
                                              pdIdent(~plate-1),
                                              pdIdent(~sample-1)))),
             data=pg)
             
m1L <- lme(diameter ~ 1, random=pdBlocked(list(pdIdent(~1),
                                               pdIdent(~plate-1),
                                               pdIdent(~sample-1))),
             data=pg)

m1L4 <- lmer(diameter ~ 1 + (1|plate) + (1|sample), data = pg)

VarCorr(m1L)
VarCorr(m1L4)

## predictions are close
all.equal(predict(m1G), c(unname(predict(m1L))), tol = 1e-4)
all.equal(unname(predict(m1L4)), c(unname(predict(m1L))), tol = 1e-7)

## original example from P&B

m2 <- lme(logDens ~ sample*dilut, data = Assay,
          random=pdBlocked(list(pdIdent(~1), pdIdent(~sample-1),
                                pdIdent(~dilut-1))))

## groupedData structure required (this fails with "invalid formula
##  for groups")
try(update(m2, data = as.data.frame(Assay)))

## original example balanced, still works when unbalanced
set.seed(101)
a2 <- Assay[sample(nrow(Assay), size=round(0.9*nrow(Assay))),]
m2_un <- update(m2, data = a2)

m4 <- gamlss(logDens ~ sample*dilut +
               re(random=pdBlocked(list(pdIdent(~1), pdIdent(~sample-1),
                                pdIdent(~dilut-1)))),
             data = Assay)
all.equal(predict(m4), c(unname(predict(m2)))) ## TRUE

Hi Ben,

I've tried your code for my unbalanced data with multi-level structure, with one continuous DV, one continuous IV and two categorical random factors.
I grouped my data using dummy value and formatted random effects as shown in your code
When I plotted the predicted values from the result, I got same slope for all groups but with different intercepts.
Is it what is supposed to happen??

Thank you!
Minami