Intuitive explanation of "integrate out random effect" We are trying to figure out an intuitive reasoning behind integrate out the unobserved random effect. The specific formula is:
$$f\left(y_i|x_i;\beta, \sigma_c\right)=\int_{-\infty}^{+\infty}\left(\prod_{i=1}^{T}f_t(y_{it}|x_{it};c,\beta)\right)(1/\sigma_c)\phi(c/\sigma_c)~\mathrm dc. $$
I think we get most of the idea: you give a random effect to a person to model the unobserved effect such as IQ or his health. But this random effect is unobserved, so you integrate it out.
But here comes my question which makes me doubt my understanding: do you calculate the likelihood for all the possible values the random effect can take?
Why do you then say that you have integrated it out? Does it mean that you let the random effect take all possible values and just pick the most probably $\beta$'s, at whatever value of the random effect?
 A: I'll offer an explanation based on marginal vs joint probability distributions. First, let's introduce some notation.
Let $Y_1,\ldots,Y_n$ be the random sample, with $Y_i$ the $i$th (random) observation vector,  with $Y_i = (Y_{i1},\ldots,Y_{iT})$; suppose also that we have some observed covariates $x_i$. Thus the observable data are $(Y_i,x_i)$, for $i=1,\ldots,n$.
Typically, but not always, the values in $Y_i$ are repeated measures taken in time or space, pertaining to the same individual $i$. Thus, towards building the likelihood function, while we may assume independence between $Y_i,Y_j$, we cannot assume independence across $Y_i$.
Here is where the random effects come in. We assume that for each of the individuals we have an unobserved random effect $c_i$; where unobserved means that we cannot measure it and, random means that it is a random variable. Since all individuals are independent, we assume that $c_i$ are also independent and that $c_i\sim N(0,\sigma_c^2)$; without loss of generality, we can assume $\sigma_c^2$ is known.
Therefore the complete data for individual $i$ is $(Y_i, x_i,c_i)$ and thus the joint probability distribution is
$$f(y_i, c_i;x_i,\beta) = f(y_i|c_i;x_i,\beta_i)\phi(c_i).$$
If we were given a value of $c_i$, we could use this probability distribution to estimate $\beta$. It's clear that no one will ever give us such a $c_i$ so we have to come out with a better strategy.
The common$^1$ approach is work unconditionally or marginally with respect to $c_i$. That is, the idea is to take the likelihood
$$\bar{f}(y_i;x_i,\beta) = \int f(y_i|c_i;x_i,\beta_i)\phi(c_i)\,dc_{i}.$$
This would guarantee that the estimator for $\beta$ will not depend on the particular $c_i$ although it will depend on the distribution of $c_i$, thus we are taking care of the dependence across $Y_i$.
To sum up, the random effects are just uninteresting random data that we cannot measure. To deal with this missingness issue, the most reasonable choice is to eliminate, i.e. integrate, them from the joint distribution.
$^1$ There is another frequentist approach called Hierarchical Generalized Linear Models in which, essentially, $c_i$ is treated as a parameter and is estimated jointly with $\beta$, using the so-called $h$-likelihood.
A: 
Do you calculate the likelihood for all the possible values the random effect can take?


Does it mean that you let the random effect take all possible values and just pick the most probably 

It’s kind of both. This is a common technique in statistics that’s often not explained explicitly: you “integrate out” a random variable that you don’t know or don’t observe, so you can get the distribution of the random variable you’re interested in.
I’ll take a step back to basic probability theory to draw a parallel to the expected value. Remember the definition of expected value for discrete random variables is
$$\mathbb{E}(X) = \sum_{x \in \mathcal{X}} x P(X = x).$$
To get the average value (expectation) of a random variable, you take a weighted average of every possible value of X. The weight given to each $x$ is how likely each event is to occur, e.g. $P(X=x)$. Remembering that integrals act like sums over infinite sets, the “weighted average” of a continuous variable is
$$\mathbb{E}(X) = \int_{-\infty}^{\infty} x f(x)dx.$$
Like $P(X=x)$ did in the discrete case, the $f(x)$ term tells us how much weight to assign to each value $x$ when we “add” everything together with the integral.
It’s the same idea in your example. We want to know the distribution $f(y | x; \beta, \sigma_c)$, but we don’t observe $c$. So, we “integrate out” $y$’s dependence on $c$:
$$
\int_{-\infty}^{+\infty}f(y|x;c,\beta) (1/\sigma_c)\phi(c/\sigma_c)~\mathrm dc. $$
Notice that only $c$ varies in the integral. $y, x$, and $\beta$ are fixed. We’re plugging in every possible value of $c$ (using an integral over the real numbers), weighing each value by its likelihood ($\phi(c/\sigma_c))$, and averaging over everything. That gives us the most plausible distribution of $y | x; \beta$ after we’ve removed (“integrated out”) everything we know about $c$.
