First let's simulate some data for a logistic regression with fixed and random parts:
set.seed(1)
n <- 100
x <- runif(n)
z <- sample(c(0,1), n, replace=TRUE)
b <- rnorm(2)
beta <- c(0.4, 0.8)
X <- model.matrix(~x)
Z <- cbind(z, 1-z)
eta <- X%*%beta + Z%*%b
pr <- 1/(1+exp(-eta))
y <- rbinom(n, 1, pr)
If we just wanted to fit a Logistic regression with no random parts, we could use the glm function:
glm(y~x, family="binomial")
glm(y~x, family="binomial")$coefficients
# (Intercept) x
# -0.2992785 2.1429825
Or constructing our own function of the log-likelihood
$$ \log\text{L}(\boldsymbol{\beta}) = \sum_{i=1}^n y_i \log \Lambda(\eta_i) + (1-y_i)\log(1-\Lambda(\eta_i)) $$
where $\Lambda(\eta_i)=\frac{1}{1+\exp(-\eta_i)}$ and $\eta_i=\boldsymbol{X_i' \beta}$
and use optim() to estimate the parameters that maximize it, as in the following example code:
ll.no.random <- function(theta,X,y){
beta <- theta[1:ncol(X)]
eta <- X%*%beta
p <- 1/(1+exp(-eta))
ll <- sum( y*log(p) + (1-y)*log(1-p) )
-ll
}
optim(c(0,1), ll.no.random, X=X, y=y)
optim(c(0,1), ll.no.random, X=X, y=y)$par
# -0.2992456 2.1427484
which of course gives the same estimates and maximizes the log-likelihood for the same value. For mixed effects, we would want something like
library(lme4)
glmer(y~x + (1|z), family="binomial")
But how can we do the same with our own function? Since the likelihood is
$$ L = \prod_{j=1}^J \int \text{Pr}(y_{1j},...,y_{n_jj}|\boldsymbol{x}, b_j) f(b_j)db_j $$
and the integral has no closed-form expression, we need to use numerical integration like Gaussian Quadrature. We can use the package statmod to get some quadratures, say 10
library(statmod)
gq <- gauss.quad(10)
w <- gq$weights
g <- gq$nodes
UPDATE: Using these quadrature locations $g_r$ and weights $w_r$ for the $r=1,...,R$ ($R=10$ here), we can approximate the integral over $b_j$ by a sum of the $R$ terms with $g_r$ substituted for $b_j$ and the whole term multiplied by the respective weights $w_r$. Thus, our likelihood function should be now
$$ L = \prod_{j=1}^J \sum_{r=1}^{R} Pr (y_{1j},...,y_{n_jj} | \boldsymbol{x}, b_j=g_r ) w_r $$
Also, we need to account for the variance of the random part, I read that this can be achieved by replacing the $b_j \sim N(0,\sigma_b^2)$ in our $\eta$ function with $\sigma_j \theta_j$ where $\theta_j \sim N(0,1)$, so in the likelihood function above we actually replace $\theta$'s with $g$'s and not $\beta$'s.
One computational issue I don't get is how to substitute the terms since the vectors won't be of the same length. But probably I don't understand that, because I'm missing something crucial here, or misunderstood how this method works.
