The most noticeable advantages of fitting mixed-effects models using the MixedModels package for Julia relative to, say, lme4 in R are speed and reliability.
For simple models fit to data sets of modest size, speed is not that much of an issue. If the model can be fit in less time than it takes to organize and plot the data to determine an appropriate model form then it is not too important if one implementation is faster than another.
However, in a simulation study that may involve tens or hundreds of thousands of model fits speed of fitting even simple models can be important. The purpose of this notebook is to show how to perform simulation studies in Julia quickly and efficiently. The emphasis will be on simulation of linear mixed-effects models but many of the techniques and software practices apply more generally.
It is almost trivial to state that the purpose of a simulation study is to generate a set of simulated values that we can use to assess whether a claim is reasonable. By the way, it is important to realize that a simulation study can only refute a claim by showing that the observed results are inconsistent with that claim. You cannot establish a result through a simulation study. You can only show that your results in a particular case are non inconsistent with the claim.
As with any other experiment it helps to determine up front what the scope of the study should be and what information from the study should be retained. The design space for the study can be very large, even for a simple case such as assessing the characteristics (Type I error and power) of a particular test, as described below. Also the amount of information generated can be very large and the expense of retaining a huge amount of information must be balanced with information loss by summarizing. There is a great temptation to perform too much data reduction, such as recording only the proportion of p-values falling below a threshold in a test. Retaining the simulated p-values themselves provides a richer source of information.
A basic rule in Julia is that, because of the way just-in-time (JIT) compilation works, any substantial computation should take place within a function. Thus the early phases of the simulation study are exercises in designing functions. In the later phases these functions should be wrapped in a package.
In Matuschek, Kliegl, Vasishth, Baayen and Bates (2017) we used a simulation study to explore "Balancing Type I Error and Power in Linear Mixed Models". The experimental design consisted of a single binary experimental factor, C, which varied within subjects (S) and within items (I). There are m items and n subjects.
Begin by loading the packages to be used. The MixedModels package is obviously needed to define and fit the models as is the DataFrames package for the data representation. The Distributions package provides the complementary cumulative distribution function (ccdf) to evaluate p-values from test statistics assuming a given reference distribution, and the Gadfly package is used for plotting.
In [1]:
using Distributions, DataFrames, Gadfly, Iterators, MixedModels
Define a mkdata function to generate a data frame according to the number of subjects, the number of items, and the coding of the levels of the covariates. (My coauthors in that paper overruled me and used a [-0.5, 0.5] encoding. Here I use the [-1, 1] encoding preferred by statisticians.)
In [2]:
function mkdata(m=10, n=30)
chars = vcat('A':'Z', 'a':'z')
tuples = collect(product(['N','Y'], ('a':'z')[1:m], chars[1:n]))
vecs = push!(Any[pool(getindex.(tuples, i)) for i in 1:3],
zeros(length(tuples)))
DataFrame(vecs, [:C, :I, :S, :Y])
end
Out[2]:
In [3]:
const dat = mkdata()
Out[3]:
In [4]:
typeof.(dat.columns) # check the column types
Out[4]:
The response, Y, is initialized to zeros so that it is available to construct LinearMixedModel objects. The const modifier in the assignment indicates that the type of dat will not change, although the actual data values may change (but in this case they won't). Using const, if appropriate, when assigning global variables is a good practice.
The lmm function creates, but does not fit, a LinearMixedModel. It is useful to have a model that has not yet been fit because it is the structure that can then be used to simulate! a response and refit! the model.
The formula for the linear mixed-effects model can be constructed with the @formula macro or using an explicit call to the Formula constructor. For example
In [5]:
srand(1234321)
refit!(simulate!(lmm(@formula(Y ~ 1 + (1|S) +(1|I)), dat),
β=[2000.], σ=300., θ=fill(1/3, 2)))
Out[5]:
or
In [6]:
srand(1234321)
const m0 = refit!(simulate!(lmm(Formula(:Y,:(1+(1|S)+(1|I))), dat),
β=[2000.], σ=300., θ=fill(1/3, 2)))
Out[6]:
We obtain the simulated response and fit the alternative model as
In [7]:
const y0 = model_response(m0)
const m1 = refit!(lmm(Formula(:Y, :(1+C+(1|S)+(1|I))), dat), y0)
Out[7]:
The likelihood ratio test statistic for H₀: β₂ = 0 versus Hₐ: β₂ ≠ 0 is the difference in the objective (negative twice the log-likelihood) values for the fitted models.
In [8]:
lrtval = objective(m0) - objective(m1)
Out[8]:
The p-value, assuming a Chisq(1) reference distribution, is
In [9]:
ccdf(Chisq(1), lrtval)
Out[9]:
Notice that the conclusions from the likelihoood ratio test are essentially the same as those from the Z test on β₂ in the coefficient summary table
Fixed-effects parameters:
Estimate Std.Error z value P(>|z|)
(Intercept) 1955.61 34.7294 56.31 <1e-99
C 5.81124 12.988 0.447432 0.6546for the alternative model m1. This is because the square of the z value, which is the ratio of the estimate of β₂ to its standard error, is essentially the same as lrtval
In [10]:
abs2(0.447432)
Out[10]:
and a Chisq(1) distribution is the square of a standard normal distribution.
Generally a hypothesis test based on fitting both the null model and the alternative model to the data (i.e. the likelihood ratio test) is more reliable than one based on fitting only the alternative model and inferring what the change to the null model fit will be (the Z-statistic). In this case they give essentially the same results.
Although it is tempting to continue with writing simulation functions, this is a good time to examine the responses that have been simulated. What do we expect them to look like?
To make it more convenient to plot the data we copy the contents of y0 vector to the :Y column of the data frame then look for a systematic shift
In [11]:
copy!(dat[:Y], y0)
plot(dat, x = :Y, Geom.density, Guide.xlabel("Simulated responses"))
Out[11]:
In [14]:
plot(dat, x = :Y, Geom.density, Guide.xlabel("Simulated responses"),
color = :C)
Out[14]:
There is no evidence of a systematic shift according to the level of the condition, C, which is what we would expect because we simulated the data from the null model.
Next consider the effect of the subjects and items
In [15]:
plot(dat, x = :S, y = :Y, Geom.boxplot, Guide.xlabel("Subject"),
Guide.ylabel("Simulated response"))
Out[15]:
In [16]:
plot(dat, x = :I, y = :Y, Geom.boxplot, Guide.xlabel("Item"),
Guide.ylabel("Simulated response"))
Out[16]:
In [17]:
plot(dat, x = :C, y = :Y, xgroup=:I, Geom.subplot_grid(Geom.boxplot),
Guide.ylabel("Simulated response"))
Out[17]: