lme4 / lme4pureR

A pure R implementation of the numerical steps in lmer and glmer
35 stars 4 forks source link

rBind and family #10

Open alberto-arletti opened 5 months ago

alberto-arletti commented 5 months ago

Dear lme4 authors, thank you for this interesting implementation in pure R. I am trying to "sandbox around" with the estimating function for the multi-level linear model. The first thing I noticed is that the function rBind in mkRanefStructures and mkRanefRepresentation seems to be deprecated for rbind.

Also I was wondering if the pure R implementation allows for estimation of a binomial variable. The useful example code

set.seed(1); n <- 1000; x <- rnorm(n); z <- rnorm(n); X <- cbind(1, x); ZZ <- cbind(1, z); grp <- gl(n/5,5)
RE <- mkRanefStructures(list(grp), list(ZZ)); Z <- t(RE$Zt)
y <- as.numeric(X%*%rnorm(ncol(X)) + Z%*%rnorm(ncol(Z)) + rnorm(n))
m <- lmer.fit(y,X,ZZ,grp)

is for a continuous case. I see the plsform function takes a family argument. Is there a way to run lmer.fit on a binomial target variable y? If not, do you have any knowledge of an implementation in R that can be tinkered with with other link functions too? Thank you

ajinkya-k commented 4 months ago

@alberto-arletti

Fixing rBind

Solving the rBind error is easy, just replace rBind with rbind everywhere in the code.

lmer.fit equivalent for binomial models:

while there doesn't appear to be a single function that does all the steps lmer.fit does for linear models, you could do the following:

  1. get the deviance function for a model with family = binomial using the pirls function
  2. optimize the function using bobyqa similar to the last step of lmer.fit

In fact such an example is included in tests/cbpp.R lines 22-29 which essentially recreate the steps in lmer.fit

# create deviance function with `lme4pureR`  
#devf <- pirls(glmod, ratios, binomial(), weights=weights, eta=eta, nAGQ=1)
devf <- pirls(glmod$X,ratios,glmod$reTrms$Zt,glmod$reTrms$Lambdat,
              glmod$reTrms$thfun,glmod$reTrms$theta,
              weights=weights,eta=eta,family=binomial)

# run `bobyqa`
bobyqa(c(1,beta0), devf, lower=c(0,rep.int(-Inf,length(beta0))))$par

CAUTION: The code in pirls.R suggests that you would have to use the flexlambda branch of lme4 for it to work

Hope this helps!