Lecture 4: Linear Mixed-effects Models

1 Slides

2 Today’s Learning Goals

By the end of this lecture, you should be able to:

Identify the model assumptions in a linear Mixed-effects model.
Associate a term (or combination of terms) in a Mixed-effects model with the following quantities:
- Fixed effect estimates.
- Variances of the random effects.
- Regression coefficients for each group and population.
- Predictions on existing groups and a new group.
Fit a linear Mixed-effects model in R, and extract estimates of the above quantities.
Identify the consequences of fitting a fixed-effects linear regression model when there are groups, whether a slope parameter is pooled or fit separately per group.
Explain the difference between the distributional assumption on the random effects and the fixed effects estimates’ sampling distribution.

3 Linear Fixed-effects Model

So far, we have been working with regression models fitted with a training set of $n$ independent elements.

Let us start with different modelling techniques from the ones you learned in DSCI 561. Under a frequentist and classical Ordinary Least-squares (OLS) paradigm, given a set of $k$ regressors $X_{i,j}$ and a continuous response $Y_i$, we fit a model

\[Y_i = \beta_0 + \beta_1 X_{i,1} + \beta_2 X_{i,2} + \ldots + \beta_k X_{i,k} + \varepsilon_{i} \; \; \; \; \text{for} \; i = 1, \ldots, n.\]

Note parameters $\beta_0, \beta_1, \dots, \beta_k$ are fixed and constant for all the observed values $\left(x_{i,1}, \dots, x_{i,k}, y_i\right)$.

Definition of fixed effects

OLS parameters $\beta_0, \beta_1, \dots, \beta_k$ are called fixed effects in Regression Analysis. In an inferential framework, it is our interest to evaluate whether they are statistically significant for the response.

To illustrate today’s topic, let us introduce a suitable dataset.

3.1 Grunfeld’s Investment Dataset

Consider the following example: to study how gross investment depends on the firm’s value and capital stock, Grunfeld (1958) collected data from eleven different American companies over the years 1935-1954.

The Grunfeld’s Investment Dataset

The data frame Grunfeld, from package {AER}, contains 220 observations from a balanced panel of 11 sampled American firms from 1935 to 1954 (20 observations per firm). The dataset includes a continuous response investment subject to two explanatory variables, market_value and capital.

Definition of balanced panel

In statistical jargon, a panel refers to a dataset in which each individual (e.g., a firm) is observed within a timeframe. Furthermore, the term balanced indicates that we have the same number of observations per individual.

Firstly, we will load the data which has the following variables:

investment: the gross investment in millions of American dollars (additions to plant and equipment along with maintenance), a continuous response.
market_value: the firm’s market value in millions of American dollars, a continuous explanatory variable.
capital: stock of plant and equipment in millions of American dollars, a continuous explanatory variable.
firm: a nominal explanatory variable with eleven levels indicating the firm (General Motors, US Steel, General Electric, Chrysler, Atlantic Refining, IBM, Union Oil, Westinghouse, Goodyear, Diamond Match, and American Steel).
year: the year of the observation (it will not be used in our analysis, but take it into account for our next in-class question).

The code below only renames some columns from the original dataset.

We will start with an in-class question via iClicker.

Exercise 12

What class of data hierarchy do you observe in this dataset? Do you expect any class of correlation within the data points? Think about each hierarchical level as a possible sampling pool.

A. Yes, we have a data hierarchy with one level: firm. Still, there will not be a correlation among subsets of data points.

B. Yes, we have a data hierarchy with one level: firm. Hence, there will be a correlation among subsets of data points.

C. There is no data hierarchy at all. All observations in the training set are independent.

D. Yes, we have a data hierarchy with two levels: firm (level 1) and the corresponding yearly observations (level 2). Hence, there will be a correlation among subsets of data points.

Main Statistical Inquiries

We are interested in assessing the association of gross investment with market_value and capital in the population of American firms. Then, how can we fit a linear model to this data?

3.2 Exploratory Data Analysis

Let us plot the $n = 220$ data points of investment versus market_value and capital but facetted by firm and use geom_smooth() to fit sub-models by firm.

Important

Only for plotting, we will transform both $x$ and $y$-axes on the logarithmic scale in base 10 (trans = "log10"). This allows us to compare those firms under dissimilar market values, capital, and gross investments.

Exercise 13

What do you observe in the above plots?

3.3 OLS Modelling Framework

Before digging into today’s new modelling approach, let us try a classical approach via OLS for comparative purposes.

Important

Always keep in mind the main statistical inquiry when taking any given modelling approach.

3.3.1 Regression Alternatives

Based on what you have seen via OLS in DSCI 561, there might be four possible approaches:

Take the average for each firm, and fit an OLS regression on the averages. This is not an ideal approach, since we would lose valuable data points if we summarize them using a summary statistic of central tendency such as the average. Recall the frequentist paradigm: the more data points, the better! Therefore, we will not use this approach in this lecture.
We could ignore firm, and fit an OLS regression with the other columns in the training set. Nonetheless, this is not an ideal approach since we will not be properly defining the corresponding systematic component of our population of interest: American firms.
Allow different intercepts for each firm. This might be an approach worth taking using OLS modelling techniques.
Allow a different slope and intercept for each firm (i.e., an interaction model!). This might be an approach worth taking using OLS modelling techniques.

3.3.2 OLS Ignoring Firm

Let us start with this basic OLS model to warm up our modelling skills regarding setting up equations. Suppose we ignore factor firm. Then, we will estimate an OLS regression with investment as a response to market_value and capital as regressors.

The regression equation for the $i$th sampled observation will be:

\[\texttt{investment}_{i} = \beta_0 + \beta_1 \texttt{marketValue}_{i} + \beta_2\texttt{capital}_{i} + \varepsilon_{i} \; \; \; \; \; \text{for} \; i = 1, \ldots, 220.\]

We use lm() and glance() to get the corresponding outputs.

With $\alpha = 0.05$, we have evidence to state see that market_value and capital are statistically associated to investment since the $p\text{-values} < .001$.

3.3.3 OLS Regression with Varying Intercept

Now, let us estimate another OLS regression model with investment as a response to market_value and capital as regressors but with varying intercepts by each firm.

We will do this with the lm() function by adding - 1 on the right-hand side of the argument formula. This - 1 will allow the baseline firm to have its intercept (i.e., replacing the usual (Intercept) in column estimate with this specific baseline company’s intercept). In this case, General Motors is the baseline company (as it appears on the left-had side of the levels() output).

By checking the adj.r.squared, we see that model_varying_intercept has a larger value (0.959) than ordinary_model (0.816) (i.e., the first fitted model without firm as a regressor). This indicates that a model with estimated intercepts by firm fits the data better than a model without taking firm into account (at least by looking at the metrics!).

Note that model_varying_intercept is equivalent to just fitting the OLS model using formula = investment ~ market_value + capital + firm. However, in this below primary_OLS, General Motors estimated intercept would be called (Intercept) whereas the remaining ten firms will have their corresponding estimated intercepts as (Intercept) + firmCompany as in the above model_varying_intercept. Moreover, the regression parameter estimates for market_value and capital (along with their corresponding standard errors, test statistics, and $p\text{-values}$) stay the same.

Important

Even though the intercept estimates by firm can be easily computed using column estimate from primary_OLS as previously mentioned; standard errors, test statistics, and $p\text{-values}$ will be different given the model parametrization in the regression equation.

Going back to model_varying_intercept and ordinary_model, we can test if there is a gain in considering a varying intercept versus fixed intercept. Hence, we will make a formal $F$-test to check whether the model_varying_intercept fits the data better than the ordinary_model.

We obtain a $p\text{-value} < .001$. Thus, with $\alpha = 0.05$, we have evidence to conclude that model_varying_intercept fits the data better than the ordinary_model.

However, this costs us one extra degree of freedom per firm except for the baseline. Therefore, we lose another 10 degrees of freedom (column DF in the anova() output).

Important

In this case, losing 10 degrees of freedom is not a big deal with 220 data points. Nonetheless, in other cases, when data is scarce or the model’s complexity demands a large amount of regression parameters, this could be an issue regarding inferential inquiries since we need them to perform the corresponding hypothesis testing.

Exercise 14

What is the sample’s regression equation for model_varying_intercept?

\[ \begin{align*} \texttt{investment}_{i,j} &= \beta_{0} + \beta_1 \texttt{marketValue}_{i,j} + \beta_2\texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \quad \text{for} \; i = 1, \ldots, 20 \; \; \text{and} \; \; j = 1, \ldots, 11. \end{align*} \]

\[ \begin{align*} \texttt{investment}_{i} &= \beta_{0} + \beta_1 \texttt{marketValue}_{i} + \beta_2\texttt{capital}_{i} + \varepsilon_{i} \\ & \qquad \qquad \quad \text{for} \; i = 1, \ldots, 220. \end{align*} \]

\[ \begin{align*} \texttt{investment}_{i,j} &= \beta_{0,j} + \beta_1 \texttt{marketValue}_{i,j} + \beta_2\texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \quad \text{for} \; i = 1, \ldots, 20 \; \; \text{and} \; \; j = 1, \ldots, 11. \end{align*} \]

3.3.4 OLS Regression for Each Firm

We can make the model more complex with two interactions (market_value * firm and capital * firm). This will estimate a linear regression by firm with its own slopes.

In this case, we are fitting eleven linear regressions, each with 20 points.

Important

We have plenty of data points in this case to estimate each parameter with the required degrees of freedom. Nonetheless, this might not hold for other small datasets. Also, we must consider other modelling scenarios in which there may be more than two regressors of interest (not only market_value and capital) with different natures (discrete and continuous), which may require a large number of degrees of freedom.

Exercise 15

What is the sample’s regression equation for model_by_firm?

\[ \begin{align*} \texttt{investment}_{i,j} &= \beta_{0,j} + \beta_{1,j} \texttt{marketValue}_{i,j} + \beta_{2,j} \texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \quad \text{for} \; i = 1, \ldots, 20 \; \; \text{and} \; \; j = 1, \ldots, 11. \end{align*} \]

\[ \begin{align*} \texttt{investment}_{j} &= \beta_{0} + \beta_1 \texttt{marketValue}_{j} + \beta_2\texttt{capital}_{j} + \varepsilon_{j} \\ & \qquad \qquad \quad \text{for} \; j = 1, \ldots, 11. \end{align*} \]

\[ \begin{align*} \texttt{investment}_{i,j} &= \beta_{0,i} + \beta_{1,i} \texttt{marketValue}_{i,j} + \beta_{2,i} \texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \quad \text{for} \; i = 1, \ldots, 20 \; \; \text{and} \; \; j = 1, \ldots, 11. \end{align*} \]

We went as far as fitting a separate linear regression for each firm.

How to interpret the coefficients in this model_by_firm?

Each regression coefficient is associated with a firm. For example, $\texttt{firmUS Steel:capital} = 0.02$ means that the variable capital has a slope of $0.02 + 0.37 = 0.39$ for US Steel. We can double-check this by estimating an individual linear regression for US Steel:

Exercise 16

Now, after all these OLS modelling approaches, are we in line with the main modelling objective with all these linear regression models?

4 Linear Mixed-effects Model

The fact that we have a sample of 11 companies and our inquiry aims to make inference on the population of American companies, along with the correlated data structure in this panel, paves the way to linear Mixed-effects modelling.

Important

In linear Mixed-effects modelling, the $n$ rows in our training set will not be independent anymore. Let us view it this way: possibly, your data subsets of elements share a correlation structure (e.g., observations over time on a given response and regressors for a specific subject in a panel). This is a fundamental idea in Mixed-effects regression.

Let us take a step back and think about a population of companies. For instance, all American companies. Grunfeld did not collect data on all the American companies but sampled 11 companies from this population. The author was interested in assessing whether market_value and capital were related to investment and by how much.

Let us assume that the $j$th sampled firm has its own intercept $b_{0,j}$ and the overall fixed intercept is $\beta_0$ for all American companies. Therefore, for the $j$th firm, we define the following mixed intercept:

\[ \beta_{0,j} = \beta_0 + b_{0,j}. \]

The intercept $b_{0,j}$ is specifically for the $j$th firm that was sampled. It will change due to chance since it is linked to the $j$th sampled firm which would make it a random effect. This is the deviation of the $j$th firm from the overall fixed intercept $\beta_0$.

Definition of random effect

A random effect in a regression model is a term that captures systematic differences across units within a group (e.g., repeated measurements within a single patient, students within a school, etc.) by allowing each group to have its own intercept and/or coefficients, treated as draws from a common distribution.

The regression paradigm of estimating a fixed unknown intercept $\beta_0$ will change now. Moreover, the intercept $\beta_{0,j}$ is what we call a mixed effect:

\[ \begin{align*} \texttt{investment}_{i,j} &= \overbrace{\beta_{0,j}}^{\text{Mixed Effect}} + \beta_1 \texttt{marketValue}_{i,j} + \beta_2\texttt{capital}_{i,j} + \varepsilon_{i,j} \\ &= (\beta_0 + b_{0,j}) + \beta_1 \texttt{marketValue}_{i,j} + \beta_2\texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \qquad \; \; \; \; \text{for} \; i = 1, \ldots, n_j \; \; \text{and} \; \; j = 1, \ldots, 11. \end{align*} \tag{1}\]

Important

Note that $n_j$ in Equation 1 is making the model even more flexible by allowing different numbers of observations $n_j$ in each $j$th firm.

Now,

\[b_{0,j}\sim \mathcal{N}(0, \sigma_0^2)\]

is called a random effect and we assume it is independent of the error component

\[\varepsilon_{i,j}\sim \mathcal{N}(0, \sigma^2).\]

The observations for the same firm (group) share the same random effect making a correlation structure.

The variance of the $i$th response for the $j$th firm will be

\[ \text{Var}(\texttt{investment}_{i,j}) = \text{Var}(b_{0,j}) + \text{Var}(\varepsilon_{i,j}) = \sigma_0^2 + \sigma^2. \]

For the $k$th and $l$th responses, within the $j$th firm, the correlation is given by:

\[ \text{Corr}(\texttt{investment}_{k,j}, \texttt{investment}_{l,j}) = \frac{\sigma^2_0}{\sigma_0^2 + \sigma^2}. \]

We could even go further and model random slopes, along with the existing fixed ones, as follows:

\[ \begin{align*} \texttt{investment}_{i,j} &= \overbrace{\beta_{0,j}}^{\text{Mixed Effect}} + \overbrace{\beta_{1,j}}^{\text{Mixed Effect}} \times \texttt{marketValue}_{i,j} + \\ & \qquad \overbrace{\beta_{2,j}}^{\text{Mixed Effect}} \times \texttt{capital}_{i,j} + \varepsilon_{i,j} \\ &= (\beta_0 + b_{0,j}) + (\beta_1 + b_{1,j}) \times \texttt{marketValue}_{i,j} + \\ & \qquad (\beta_2 + b_{2,j}) \times \texttt{capital}_{i,j} + \varepsilon_{i,j} \\ & \qquad \qquad \qquad \text{for} \; i = 1, \ldots, n_j \; \; \text{and} \; \; j = 1, \ldots, 11; \end{align*} \tag{2}\]

with $(b_{0,j}, b_{1,j}, b_{2,j})^{T} \sim \mathcal{N}_3(\mathbf{0}, \mathbf{D})$, where $\mathbf{0} = (0, 0, 0)^T$ and $\mathbf{D}$ is a generic covariance matrix.

Important

Note that the random effects $b_{0,j}$, $b_{1,j}$, and $b_{2,j}$ in Equation 2 have a multivariate Gaussian (or Normal) distribution of $d = 3$. You can review further details about this distribution in this appendix.

4.1 What is this Generic Covariance Matrix $\mathbf{D}$?

This is a standard form in linear Mixed-effects modelling. Hence, this matrix becomes:

\[ \begin{equation*} \mathbf{D} = \begin{bmatrix} \sigma_{0}^2 & \rho_{01} \sigma_{0} \sigma_{1} & \rho_{02} \sigma_{0} \sigma_{2}\\ \rho_{01} \sigma_{0} \sigma_{1} & \sigma_{1}^2 & \rho_{12} \sigma_{1} \sigma_{2}\\ \rho_{02} \sigma_{0} \sigma_{2} & \rho_{12} \sigma_{1} \sigma_{2} & \sigma_{2}^2 \end{bmatrix} = \begin{bmatrix} \sigma_{0}^2 & \sigma_{0, 1} & \sigma_{0, 2} \\ \sigma_{0, 1} & \sigma_{1}^2 & \sigma_{1,2}\\ \sigma_{0, 2} & \sigma_{1, 2} & \sigma_{2}^2 \end{bmatrix}, \end{equation*} \]

where $\sigma_{0}^2$, $\sigma_{1}^2$, and $\sigma_{2}^2$ are the variances of $b_{0,j}$, $b_{1,j}$, and $b_{2,j}$ respectively.

Moreover $\rho_{uv} \in [0,1]$ is the correlation between the $u$th and the $v$th random effects. We can reexpress the covariances as $\sigma_{u, v}$.

Important

$\rho_{uv}$ indicates a Pearson correlation.
While the random effects are assumed to follow a mulivariate Normal distribution, this is different from the sampling distribution of the estimates of the fixed effects.
The multivariate Normal distribution explains the variability of random regression intercepts and coefficients. The spread does not change when we collect more data.
The sampling distribution explains the uncertainty in the fixed regression estimates and gets narrower as we collect more data.

4.2 Model Fitting, Inference, and Coefficient Interpretation

Let us estimate the regression model with a mixed intercept only (mixed_intercept_model) via the function lmer() from package lme4. Note that (1 | firm) allows the model to have a random intercept by firm.

Recall the regression equation for mixed_intercept_model:

where

\[b_{0,j}\sim \mathcal{N}(0, \sigma_0^2)\]

and

\[\varepsilon_{i,j}\sim \mathcal{N}(0, \sigma^2).\]

The section Random effects will show the estimated standard deviations $\hat{\sigma}_0$ and $\hat{\sigma}$ as 82.10 and 50.27, respectively. Estimated Fixed effects $\hat{\beta}_0$, $\hat{\beta}_1$, and $\hat{\beta}_2$ are -54.0318, 0.1094, and 0.3082 respectively.

Now, let us estimate the Mixed-effects regression model with mixed intercept and slopes (full_mixed_model). Note that (market_value + capital | firm) allows the model to have a random intercept and slopes by firm.

Recall the regression equation for full_mixed_model:

with $(b_{0,j}, b_{1,j}, b_{2,j})^{T} \sim \mathcal{N}_3(\mathbf{0}, \mathbf{D})$, where $\mathbf{0} = (0, 0, 0)^T$ and $\mathbf{D}$ is a generic covariance matrix.

The section Random effects will show the estimated standard deviations $\hat{\sigma}_0$, $\hat{\sigma}_1$, $\hat{\sigma}_2$, and $\hat{\sigma}$ as 15.15612, 0.05235, 0.12291, and 40.77647 respectively. Estimated Fixed effects $\hat{\beta}_0$, $\hat{\beta}_1$, and $\hat{\beta}_2$ are -7.79756, 0.06118, and 0.22694 respectively.

Let us proceed with inference using mixed_intercept_model. We now assess whether the fixed effects are statistically associated with investment in each model via summary(). We will use the package lmerTest along with function summary().

We do the same with full_mixed_model.

We can see that market_value and capital are significant with $\alpha = 0.05$ in both models. Moreover, the regression coefficients’ interpretation for the fixed effects will be on the effect these regressors have on the population investment mean of the American companies.

We can obtain the estimated coefficients by firm along with the intercepts for both models via coef().

Column (Intercept) is the sum $\hat{\beta}_0 + \hat{b}_{0,j}$. Note the estimated regression coefficients for market_value and capital are the same since mixed_intercept_model only has $\beta_1$ and $\beta_2$ as its general modelling setup.

The coefficient summary changes in full_mixed_model given that we also include random effects for market_value and capital, as shown below. Columns market_value and capital are the sums $\hat{\beta}_1 + \hat{b}_{1,j}$ and $\hat{\beta}_2 + \hat{b}_{2,j}$, respectively. Column (Intercept) is the sum $\hat{\beta}_0 + \hat{b}_{0,j}$.

Exercise 17

Note the standard errors for the estimated slopes in market_value and capital behave in a really particular way when comparing the OLS model_varying_intercept and the Mixed-effects full_mixed_model.

Therefore, what are the advantages of a Mixed-effects model over an OLS model with fixed-effects only?

4.3 Estimation

Mixed-effects models are still fit using likelihood ideas, but there is more to estimate than in OLS. In OLS, we estimate the fixed-effect intercept and coefficients and a single overall variance $\sigma^2$. In a Mixed-effects model, we also estimate additional variance components: how much groups/subjects vary through the random effects (and sometimes how random effects move together), plus the aforementioned overall variance $\sigma^2$.

Restricted maximum likelihood (REML) is an estimation approach for Mixed-effects regression. The basic idea is that estimating the fixed effects “uses up” some information in the data, and plain MLE can then underestimate the random-effect variances. REML adjusts for that by focusing the likelihood on the part of the data that reflects variation left over after accounting for the fixed effects.

In practice, REML proceeds in two steps:

Estimate the variance components (random-effect variances/covariances and the overall variance $\sigma^2$).
Estimate the fixed-effects given those variance estimates (essentially a weighted regression that accounts for within-group correlation).

Note

If you would like a more technical treatment of REML (including how it differs from MLE and how the restricted likelihood is constructed), see the documentation for the lme4 package.

4.4 Prediction

We can make two classes of predictions with Mixed-effects models:

To predict on an existing group, we find that group’s regression coefficients (and therefore model function) by summing the fixed effects and (if present) the random effects, then use that model function to make predictions.
To predict on a new group (using a mean prediction), we use the fixed effects as the regression coefficients (because the random effects are assumed to have a mean equal to zero) and use that model function to make predictions.

For predictions on an existing group in our training set, via full_mixed_model, we have:

If we wanted to predict the investment for General Motors with a market_value of USD $\$2,000$ million and capital of USD $\$1,000$ million, then our answer would be USD $\$537.4$ million. This prediction uses $\hat{\beta}_0$, $\hat{\beta}_1$, $\hat{\beta}_2$, $\hat{b}_{0, j}$, $\hat{b}_{1, j}$, and $\hat{b}_{2, j}$ from Equation 2.

Whereas that for predictions on American companies in general, via full_mixed_model, we have:

If we wanted to predict the MEAN investment for American companies with a market_value of USD $\$2,000$ million and capital of USD $\$1,000$ million, then our answer would be USD $\$341.51$ million. This prediction only uses $\hat{\beta}_0$, $\hat{\beta}_1$, and $\hat{\beta}_2$ from Equation 2.

5 Wrapping Up on Mixed-effects Modelling

In many different cases, when there is a correlation structure in our observations, OLS models are not suitable for our inferential or predictive inquiries.
Therefore, linear Mixed-effects models are suitable for correlated observations. Nonetheless, the model’s complexity will also be in function of our specific inquiries.
We can even extend the Mixed-effects approach to generalized linear models (GLMs)!

1 Slides

2 Today’s Learning Goals

3 Linear Fixed-effects Model

3.1 Grunfeld’s Investment Dataset

3.2 Exploratory Data Analysis

3.3 OLS Modelling Framework

3.3.1 Regression Alternatives

3.3.2 OLS Ignoring Firm

3.3.3 OLS Regression with Varying Intercept

3.3.4 OLS Regression for Each Firm

4 Linear Mixed-effects Model

4.1 What is this Generic Covariance Matrix \(\mathbf{D}\)?

4.2 Model Fitting, Inference, and Coefficient Interpretation

4.3 Estimation

4.4 Prediction

5 Wrapping Up on Mixed-effects Modelling