Linear models in R estimate parameters that are considered fixed or non-random and are called fixed-effects. In contrast, random-effect parameters 
assume data share a common error distribution, and can produce different estimates when there are small amounts of data or outliers. 
Models with both fixed and random-effects are mixed-effect models or linear mixed-effect regression.

The lme4 package fits mixed-effect models (models with both fixed- and random-effects) with lmer(), which uses a formula similar to lm().
But, random-effect intercepts use special syntax:

lmer(y ~ x + (1 | random-effect), data = my_data)

lmer() function requires the model to include a random-effect, otherwise the model gives you an error. Here, you will fit a lm() and a lmer(), 
and then graphically compare the fitted models using a subset of the data. We provide this code because of the advanced data wrangling, 
which is required because random-effects are usually not plotted (ggplot2 also does not include nice plot options for mixed-effect models). 
In this plot, notice how the dashed lines from random-effect slopes compare to the solid lines from the fixed-effect slopes.

Note: broom.mixed is required because the broom package does not support lme4.

# Build a liner model including class as fixed-effect model
lm_out <- lm(mathgain ~ classid + mathkind, data = student_data)

# Build a mixed-effect model including class id as a random-effect
lmer_out <- lmer(mathgain ~ mathkind + (1 | classid), data = student_data)

# Extract out the slope estimate for mathkind
tidy(lm_out) %>%
    filter(term == "mathkind")
    
tidy(lmer_out) %>%
    filter(term == "mathkind")
    
    
# Re-run the models to load their outputs
lm_out <- lm(mathgain ~ classid + mathkind, data = student_data)
lmer_out <- lmer(mathgain ~ mathkind + (1 | classid), data = student_data)

# Add the predictions to the original data
student_data_subset <-
    student_data %>%
    mutate(lm_predict = predict(lm_out),
           lmer_predict = predict(lmer_out)) %>%
    filter(schoolid == "1")

# Plot the predicted values
ggplot(student_data_subset,
       aes(x = mathkind, y = mathgain, color = classid)) +
    geom_point() +
    geom_line(aes(x = mathkind, y = lm_predict)) +
    geom_line(aes(x = mathkind, y = lmer_predict), linetype = 'dashed') +
    xlab("Kindergarten math score") +
    ylab("Math gain later in school") +
    theme_bw() +
    scale_color_manual("Class ID", values = c("red", "blue"))