library(geex) library(lme4)
A user asked a reasonable question about score contributions summing to zero for a logistic-Normal model.
Here we analyze 500 pairs of 2 units generated from:
[ \Pr(Y_{ij} = 1 | X, b) = \mbox{logit}^{-1}(\beta_0 + \beta_1 x_{ij} + b_i) = h(\beta, x_{ij}, b_i) ]
where $X_{ij} \sim N(0, 1)$ and $b_i \sim f(b_i) = N(0, \sigma^2)$. Let $\theta = {\beta_0, \beta_1, \sigma^2 } = {2, 2, 1}$.
The loglihood for this model is:
[ l(\theta; Y, X) = \sum_{i = 1}^m { \int \prod_{j = 1}^{n_i} h(\beta, x_{ij}, b_i)^y{ij}[1 - h(\beta, x_{ij}, b_i)]^{1 - y{ij}} f(b_i) } ]
Take the derivative of $l(\theta; Y, X)$, set equal to zero, and solve for $\theta$ to find $\hat{\theta}$.
set.seed(1) #generate data; 500 pairs n <- 500 m <- 2 beta <- c(2, 2) id <- rep(1:n, each=m) x <- rnorm(m*n) b <- rep(rnorm(n), each=m) y <- rbinom(m*n, 1, plogis(cbind(1, x) %*% beta + b)) d <- data.frame(y,x,id)
$\theta$ can be estimated from the lme4::glmer
function.
#fit logistic-normal random effects model fit <- lme4::glmer(y~x+(1|id), family=binomial, data=d) print(summary(fit))
Note that geex
does most of this work internally.
eefun <- function(data, lme4model){ f <- grab_psiFUN(lme4model, data = data) function(theta) f(theta) } splitdt <- split(d, f = d$id) theta_hat <- unlist(lme4::getME(fit, c('beta', 'theta'))) psi <- lapply(splitdt, function(dt){ eefun(dt, fit) }) # Evaluate psi functions psi_eval <- lapply(psi, function(f){ f(theta_hat) })
Do these evaluated score equations sum to zero?
colSums(matrix(unlist(psi_eval), ncol = 3, byrow = TRUE))
No.
As another check, how does the observed loglihood compare to lme4::glmer
?
# geex internal functions objFun_glmerMod_binomial and binomial integrand copied # here since they are not exported objFun_glmerMod_binomial <- function(parms, response, xmatrix, linkinv) { log(stats::integrate(binomial_integrand, lower = -Inf, upper = Inf, parms = parms, response = response, xmatrix = xmatrix, linkinv = linkinv)$value) } binomial_integrand <- function(b, response, xmatrix, parms, linkinv){ if(class(xmatrix) != 'matrix'){ xmatrix <- as.matrix(xmatrix) } pr <- linkinv( drop(outer(xmatrix %*% parms[-length(parms)], b, '+') ) ) hh <- stats::dbinom(response, 1, prob = pr) hha <- apply(hh, 2, prod) hha * stats::dnorm(b, mean = 0, sd = parms[length(parms)]) } # Compute the per-unit contribution the loglihood objFun <- lapply(splitdt, function(dt){ X <- cbind(1, dt$x) objFun_glmerMod_binomial(parms = theta_hat, response = dt$y, xmatrix = X, linkinv = plogis ) }) sum(unlist(objFun)) logLik(fit)
Not bad.
What happens if instead of using $\hat{\theta}_{glmer}$ to compute the scores we use estimated parameters obtained from the finding the roots of the score functions as they are defined within the grab_eeFUN
function?
NOTE: the evaluation of m_estimate
is slow!
# NOT run to speed creation of documentation attempt <- m_estimate( estFUN = eefun, data = d, units = 'id', root_control = setup_rootFUN(start = theta_hat), # start root finding at lme4's beta hat. compute_vcov = FALSE, outer_args = list(lme4model = fit) ) ## parameter estimates obtained from geex attempt psi_eval2 <- lapply(psi, function(f){ f(roots(attempt)) }) colSums(matrix(unlist(psi_eval2), ncol = 3, byrow = TRUE))
Ah! That looks much better.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.