Closed mgoplerud closed 3 years ago
Apologies for the slow response. I fixed the specific issue in 2f2d55c, however the best way to fit a model with a fixed random effect covariance is likely to change the optimization problem. By calling bglmer
(or glmer
, even) with devFunOnly = TRUE
, a function of just the parameters in the model is returned. The first part of those parameters are "theta", or the Choleski decompositions of the covariance matrices, while the second part are "beta", or the fixed effects. For example, something like:
library(lme4)
y <- rbinom(100, 1, prob = 0.5)
x <- rnorm(100)
g <- sample(1:10, 100, replace = T)
fixed_theta <- 1.5
dev_fun <- glmer(y ~ x + (1 | g), family = binomial(), devFunOnly = TRUE)
wrapper <- function(beta) {
pars <- c(fixed_theta, beta)
dev_fun(pars)
}
# Get beta for fixed random effect variance using optim:
optim_fit <- optim(c(0, 0), wrapper, hessian = TRUE)
# beta:
optim_fit$par
# s.e.
sqrt(diag(solve(optim_fit$hessian)))
# It is possible to install the values in a glmerMod lme4 object
dummy_fit <- suppressWarnings(
glmer(y ~ x + (1 | g), family = binomial(),
control = glmerControl("Nelder_Mead", optCtrl = list(maxfun = 1L)))
)
rho <- environment(dev_fun)
rho$baseOffset
rho$pp <- dummy_fit@pp
rho$resp <- dummy_fit@resp
# Calling the devfun does most of the work, installing the desired theta
# and updating random effects.
invisible(wrapper(optim_fit$par))
# Beta is handled separately, and must be installed elsewhere. The
# large value 1e16 is simply so that the standard error for the
# random effects covariance inverts to a value close to 0. It is
# not used in any meaningful calculation.
dummy_fit@beta <- optim_fit$par
dummy_fit@optinfo$Hessian <-
as.matrix(Matrix::bdiag(diag(1e16, 1), optim_fit$hessian))
summary(dummy_fit)
Thanks! That makes a lot of sense.
Hello! This response is super helpful for a project I have been working on. However, I had a question regarding the bit where the values are installed back into the lme4 object. I noticed when running your code that the standard errors for the fixed effects as shown in the lme4 object (from: summary(dummy_fit)), don't quite align with the standard errors from above (sqrt(diag(solve(optim_fit$hessian)))). Do you know what is happening there/what I am missing? Thanks!
I forgot a factor of 2, since the objective function is on the deviance scale. Also, lme4 might have changed a tiny bit since I posted that (not sure), and the Hessian goes in a slightly different place. The now-correct version would be:
library(lme4)
y <- rbinom(100, 1, prob = 0.5)
x <- rnorm(100)
g <- sample(1:10, 100, replace = T)
fixed_theta <- 1.5
dev_fun <- glmer(y ~ x + (1 | g), family = binomial(), devFunOnly = TRUE)
wrapper <- function(beta) {
pars <- c(fixed_theta, beta)
dev_fun(pars)
}
# Get beta for fixed random effect variance using optim:
optim_fit <- optim(c(0, 0), wrapper, hessian = TRUE)
# beta:
optim_fit$par
# s.e.
sqrt(diag(solve(0.5 * optim_fit$hessian)))
# It is possible to install the values in a glmerMod lme4 object
dummy_fit <- suppressWarnings(
glmer(y ~ x + (1 | g), family = binomial(),
control = glmerControl("Nelder_Mead", optCtrl = list(maxfun = 1L)))
)
rho <- environment(dev_fun)
rho$baseOffset
rho$pp <- dummy_fit@pp
rho$resp <- dummy_fit@resp
# Calling the devfun does most of the work, installing the desired theta
# and updating random effects.
invisible(wrapper(optim_fit$par))
# Beta is handled separately, and must be installed elsewhere. The
# large value 1e16 is simply so that the standard error for the
# random effects covariance inverts to a value close to 0. It is
# not used in any meaningful calculation.
dummy_fit@beta <- optim_fit$par
dummy_fit@optinfo$derivs$Hessian <-
as.matrix(Matrix::bdiag(diag(1e16, length(dummy_fit@theta)), optim_fit$hessian))
summary(dummy_fit)
Excellent - thank you so much; this is incredibly helpful!
Hello! I am trying to fit a
blme
where the variance of the random effect is fixed and I found a nice solution for the residual variance for the linear model in blme (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q1/021740.html). It seems that there isn't a way to specifypoint
for the variance so I thought about writing a custom function (e.g. strongly penalizing any divergence from the fixed value). I encountered some issues, as documented below.Created on 2021-02-24 by the reprex package (v0.3.0)