Below we demonstrate how to fit linear mixed models using the Ruby gem mixed_models. We also show how to access various parameters estimated by the linear mixed model, and we use them to assess the goodness of fit of the resulting model.
We will fit the model with the user friendly method LMM#from_formula, which mimics the behaviour of the function lmer from the R package lme4. The formula language used by mixed_models is in fact a subset of the lme4 formula interface. It allows to fit all of the same models, but does not allow for cetain shortcuts, namely the symbols *, || and /. However, for all of these shortcuts longer alternative formulations exist in all cases (e.g. you need to use the equivalent formulation a + b + a:b instead of a*b). An expanation of the lme4 formula interface can be found in the lme4 vignette. Additionally, there is a multitude of posts on stackexchange (such as this) discussing the formula language of lme4.
Example data
Linear mixed model
Model attributes
Fitted values and residuals
Assessing the quality of the model fit
The data set, which is simulated, contains two numeric variables Age and Aggression, and two categorical variables Location and Species. These data are available for 100 (human and alien) individuals.
The data is supplied to LMM#from_formula as a Daru::DataFrame (from the excellent Ruby gem daru). We load the data, and display the first 10 lines with:
In [1]:
require 'daru'
alien_species = Daru::DataFrame.from_csv '../examples/data/alien_species.csv'
# mixed_models expects that all variable names in the data frame are ruby Symbols:
alien_species.vectors = Daru::Index.new(alien_species.vectors.map { |v| v.to_sym })
alien_species.head
Out[1]:
We model the Aggression level of an individual as a linear function of the Age (Aggression decreases with Age), with a different constant added for each Species (i.e. each species has a different base level of aggression). Moreover, we assume that there is a random fluctuation in Aggression due to the Location that an individual is at. Additionally, there is a random fluctuation in how Age affects Aggression at each different Location.
Thus, the Aggression level of an individual of Species $spcs$ who is at the Location $lctn$ can be expressed as: $$Aggression = \beta_{0} + \gamma_{spcs} + Age \cdot \beta_{1} + b_{lctn,0} + Age \cdot b_{lctn,1} + \epsilon,$$ where $\epsilon$ is a random residual, and the random vector $(b_{lctn,0}, b_{lctn,1})^T$ follows a multivariate normal distribution (the same distribution but different realizations of the random vector for each Location). That is, we have a linear mixed model with fixed effects $\beta_{0}, \beta_{1}, \gamma_{Dalek}, \gamma_{Ood}, \dots$, and random effects $b_{Asylum,0}, b_{Asylum,1}, b_{Earth,0},\dots$.
We fit this model in mixed_models using a syntax familiar from the R package lme4 (as described above), and display the estimated fixed and random effects coefficients:
In [2]:
require 'mixed_models'
model_fit = LMM.from_formula(formula: "Aggression ~ Age + Species + (Age | Location)",
data: alien_species)
puts "Fixed effects:"
puts model_fit.fix_ef
puts "Random effects:"
puts model_fit.ran_ef
A summary of some important information about the fixed and random effects of the model can be conveniently displayed with the methods LMM#fix_ef_summary and LMM#ran_ef_summary.
LMM#fix_ef_summary contains the fixed effects coefficient estimates, the standard deviations of the estimates, as well as the corresponding z scores and Wald Z p-values testing the significance of each fixed effects term.
In [3]:
model_fit.fix_ef_summary
Out[3]:
LMM#ran_ef_summary summarizes the correlation structure of the random effects terms, that is, the correlation matrix of the vector $(b_{lctn,0}, b_{lctn,1})^T$ in the present data analysis.
In [4]:
model_fit.ran_ef_summary
Out[4]:
In particular, we see that the standard deviation of the constant term corresponding to the variable Location is rather big (104.26...), and therefore Location must contribute hugely to the variance in Aggression. The variability of the effect of Age with respect to Location, however, seems rather small (0.15567...) in comparison.
Apart from the fixed and random effects coefficients (seen above), we can access many attributes of the fitted model. Among others:
fix_ef_names and ran_ef_names are Arrays of names of the fixed and random effects.
reml is an indicator whether the profiled REML criterion or the profiled deviance function was optimized by the model fitting algorithm.
formula returns the R-like formula used to fit the model as a String.
model_data, optimization_result and dev_fun store the various model matrices in an LMMData object, the results of the utilized optimization algorithm, and the corresponding objective function as a Proc.
sigma2 is the residual variance (unless weights was specified in the model fit).
sigma_mat is the covariance matrix of the multivariate normal random effects vector.
We can look at some of these parameters for our example model:
In [5]:
puts "REML criterion used: \t#{model_fit.reml}"
puts "Residual variance: \t#{model_fit.sigma2}"
puts "Formula: \t" + model_fit.formula
puts "Variance of the intercept due to 'location' (i.e. variance of b0): \t#{model_fit.sigma_mat[0,0]}"
puts "Variance of the effect of 'age' due to 'location' (i.e. variance of b1): \t#{model_fit.sigma_mat[1,1]}"
puts "Covariance of b0 and b1: \t#{model_fit.sigma_mat[0,1]}"
Some further convenience methods are (apart from #fix_ef_summary and #ran_ef_summary as seen above):
#sigma returns the square root of sigma2.
#theta returns the optimal solution of the minimization of the deviance function or the REML criterion (whichever was used to fit the model).
#deviance returns the value of the deviance function or the REML criterion at the optimal solution.
#ran_ef_cov is a version of #ran_ef_summary that returns the covariance structure of the random effects rather than the correlation structure.
In [6]:
puts "Residual standard deviation: \t#{model_fit.sigma}"
puts "REML criterion: \t#{model_fit.deviance}"
In [7]:
puts "Fitted values at the population level:"
model_fit.fitted(with_ran_ef: false)
Out[7]:
In [8]:
puts "Model residuals:"
model_fit.residuals
Out[8]:
We can assess the goodness of the model fit (to some extent) by plotting the residuals agains the fitted values, and checking for unexpected patterns. We use the gem gnuplotrb for plotting.
In [9]:
require 'gnuplotrb'
include GnuplotRB
x, y = model_fit.fitted, model_fit.residuals
fitted_vs_residuals = Plot.new([[x,y], with: 'points', pointtype: 6, notitle: true],
xlabel: 'Fitted', ylabel: 'Residuals')
Out[9]:
We see that the residuals look more or less like noise, which is good.
We can further analyze the validity of the linear mixed model somewhat, by checking if the residuals appear to be approximately normally distributed.
In [10]:
bin_width = (y.max - y.min)/10.0
bins = (y.min..y.max).step(bin_width).to_a
rel_freq = Array.new(bins.length-1){0.0}
y.each do |r|
0.upto(bins.length-2) do |i|
if r >= bins[i] && r < bins[i+1] then
rel_freq[i] += 1.0/y.length
end
end
end
bins_center = bins[0...-1].map { |b| b + bin_width/2.0 }
residuals_hist = Plot.new([[bins_center, rel_freq], with: 'boxes', notitle: true],
style: 'fill solid 0.5')
Out[10]:
The histogram does not appear to be too different from a bell shaped curve, although it might be slightly skewed to the right.
We can further explore the validity of the normality assumption by looking at the Q-Q plot of the residuals.
In [11]:
require 'distribution'
observed = model_fit.residuals.sort
n = observed.length
theoretical = (1..n).to_a.map { |t| Distribution::Normal.p_value(t.to_f/n.to_f) * model_fit.sigma}
qq_plot = Plot.new([[theoretical, observed], with: 'points', pointtype: 6, notitle: true],
['x', with: 'lines', notitle: true],
xlabel: 'Normal theoretical quantiles', ylabel: 'Observed quantiles',
title: 'Q-Q plot of the residuals')
Out[11]:
The straight line in the above plot is simply the diagonal. We see that the observed quantiles aggree with the theoretical values fairly well, as expected from a "good" model.