Parameter Estimation in Generalized Linear Mixed Models Let us assume a generalized linear mixed model with a binary dependent variable $y_{i, t } $ that is explained by a fixed effect matrix X and a simple random intercept for each individual $i$
$y_{i,t} = x_{i,t}'\beta +\alpha_i  + \epsilon_{i,t} \hspace{35pt} \alpha \sim N(0, \sigma_{\alpha})$
As far as I know, when we want to estimate the fixed and random effects via maximum likelihood, we have to integrate out the random effects over every individual so that the log likelihood becomes
$ LL(\beta) = \sum_{i=1}^N \log \left(\int_{-\infty}^{\infty} \prod_{t=1}^{T_i}F(y_{i,t}|\alpha_i, \beta) g(\alpha_i) d \alpha_i \right)$
However, if we integrate the random intercepts out, the $\alpha$ disappears and we can only solve for $\beta$. How do we estimate $\alpha$ then?
 A: When you have a binary outcome variable you typically use a link function to connect the probability of a positive response to the linear predictor that includes the random effects, i.e.,
\begin{equation}
\left \{ \begin{array}{l}
\log \displaystyle \frac{\pi_{it}}{1 - \pi_{it}} = x_{it}^\top \beta + z_{ij}^\top b_i,\\\\
\pi_{it} = \Pr(y_{it} = 1 \mid b_i),\\\\
b_i \sim \mathcal N(0, D),
\end{array}
\right.
\end{equation}
where $x_{it}$ is a design vector for the fixed effects coefficients $\beta$, and $z_{it}$ a design vector for the random effects $b_i$.
When you fit the model under maximum likelihood, the random effects are unobserved variables that are integrated out from the expression of the log-likelihood of the model, i.e.,
$$\ell(\theta) = \sum_{i = 1}^n \log \int \prod_t p(y_{it} \mid b_i; \beta) \; p(b_i; D) \; db_i,$$
where $\theta^\top = (\beta^\top, \mbox{vech}(D))$ denotes the parameters of the model. The integral does not have a closed-form solution, and is typically approximated using either the Laplace approximation or adaptive Gaussian quadrature (the latter is considered better).
The random effects are estimated in a second step using empirical Bayes methodology, that is, for $i = 1, \ldots, n$ they are estimated as the modes of the posterior distribution 
$$p(b_i \mid y_i) \propto p(y_i \mid b_i; \hat \beta) \; p(b_i; \hat D),$$
where $\hat \beta$ and $\hat D$ denote the maximum likelihood estimates of $\beta$ and $D$, respectively.
