Finding an optimal mixed effects model

After processing your MRI data you may be wondering to yourself, what is the most appropriate statistical model to explain the variance in my data?

Here I will go through the steps necessary to compare models based on best fit of the data. For example, if you have longitudinal data and are unsure how to fit the effect of time, you may want to compare a linear to a quadratic fit (old example). Alternatively, you could compare two models with varying complexity, as in the first example. The preferred method implements the AIC as described first, but the older method is included at the bottom as well.

AIC

Here we compare one model that interacts timepoint (age) and treatment including sex as a covariate with a second model that interacts age, sex, and treatment. While this example is written for mincLmer, it applies mostly the same for vertexLmer or anatLmer.

You will need to load the following libraries in R:

> library(RMINC)
> library(lme4)

First define a function to calculate AIC:

AIC_summary <-
  function(mmod)
    c(RMINC:::fixef_summary(mmod), AIC = extractAIC(mmod)[2])

Then, run both of the models you are interested in comparing.

mlm <-
  mincLmer(jacobians ~ Treatment * age + Sex + (1|ID)
          , data = data
          , REML=FALSE
          , summary_type = AIC_summary)
mlm2 <-
  mincLmer(jacobians ~ Treatment * age * Sex + (1|ID)
          , data = data
          , REML=FALSE
          , summary_type = AIC_summary)

When they have run completely, write the AIC columns to minc files thus:

mincWriteVolume(mlm, “model1_AIC.mnc”, column = ‘AIC’)
mincWriteVolume(mlm2, “model2_AIC.mnc”, column = ‘AIC’)

Then, open a terminal and load minc-toolkit

In your terminal, type:

minccalc -expression 'A[0]<A[1]?1:2' model1_AIC.mnc model2_AIC.mnc winnermap.mnc

In this expression, A[0] is the first .mnc file you enter (model1_AIC.mnc). A[1] is the second (model2_AIC.mnc). This expression will compare A[0] to A[1] at each voxel and create an output mnc (winnermap.mnc) showing at each voxel which model is better. For AIC, a smaller value indicates a better model. If A[0] is less (better) than A[1], that voxel receives a value of 1. If not, that voxel receives a value of 2. Then, we can Display winnermap.

Display winnermap.mnc

Areas that are black indicate model1 is a better fit. Areas that are white indicate model2 is a better fit.

In some cases it may be obvious which model is better, but in this example both maps account for a fairly equal amount of the image. Two things can help you decide which model to end up using. First: pay attention to areas that may be subject to motion artefacts, for example dorsal, posterior areas. Second: check if one model is more symmetrical. In this case, I chose model 1.

Running AIC on linear models

If you are simply working with a linear model, and not a linear mixed effects model, you can do the following: First run your two models and store them as model1 and model 2:

model1 = mincLm(rel_jac_overall ~ new_groups*sex,filtered.data, mask=votedmask)
model2 = mincLm(rel_jac_overall ~ new_groups+sex,filtered.data, mask=votedmask)

Then you can run compare_models to compare the two models using AIC:

model_comparison = compare_models(model1, model2, metric=AICc) The output looks as follows, where the first column is the AIC for the first model (for each voxel), and the second column corresponds to the second model:

model_comparison
model_comparison matrix, showing  6  of  24749075  rows, and  2  of  2 columns
print more by supplying n and width to the print function
showing models:
showing models:
1: rel_jac_overall ~ new_groups * sex
2: rel_jac_overall ~ new_groups + sex
         [,1]     [,2]
[1,] 6.189573 8.285714
[2,] 6.189573 8.285714
[3,] 6.189573 8.285714
[4,] 6.189573 8.285714
[5,] 6.189573 8.285714
[6,] 6.189573 8.285714

If you want the AIC map for each model you can store model_comparison[,1] or model_comparison[,2] as minc files: mincWriteVolume(buffer=model_comp[,1], output.filename = "AIC_model1.mnc", like.filename = "your_nlin3_average.mnc")

Alternatively, if you want a map as the one described for the lmer method, you can run the following line:

winner = apply(compare_models(model1,model2), 1, which.min)

As above, areas in black will be best explained by model 1, and areas in white are best explained by model 2.

Old Comparison method (still valid)

You will need to load the following libraries in R:

> library(RMINC)
> library(splines)
> library(lme4)

Once you've set up your data frames, you will run your first, most simple model. Here I am modelling a linear age term using nonlinear splines. When comparing models with different fixed effects, as we will in this example, you should set REML=FALSE. You would use REML=TRUE if you are comparing models with different random effects, or when performing your actual statistical analysis.

mod_lin = mincLmer(abs_jac ~ ns(age,1) + (1|mouse_id), mask=mask, REML=FALSE)

Next, you will estimate your degrees of freedom, and perform a False Discovery Rate (FDR) correction on your model: fdr_mod_lin = mincFDR(mincLmerEstimateDF(mod_lin))

Next, you will run your more complex model. In this case, it I am fitting a quadratic age term followed by the degrees of freedom estimation and FDR correction:

mod_quad = mincLmer(abs_jac ~ ns(age,2) + (1|mouse_id), mask=mask, REML=FALSE) fdr_mod_quad = mincFDR(mincLmerEstimateDF(mod_quad))

Now you are ready to compare models. We will use the log likelihood ratio test. The simpler model goes first, the more complex model goes second:

linvsquad = mincFDR(mincLogLikRatio(mod_lin,mod_quad))

Here is an example result from this comparison. In this case, the quadratic model is significantly better at explaining the variance in my data at q=0.01:

Degrees of Freedom: 3 
FDR Thresholds:
       mod_quad
0.01 11.955936
0.05  8.221689
0.1   6.567381
0.15  5.577953
0.2   4.862967

Next, you have to visualize where in the brain your new model best fits your data. To do so, you will output your model comparison results as a minc file: mincWriteVolume(buffer = linvsquad, column = "mod_quad", output.filename = "model_comp_lin_vs_quad.mnc", like.filename = "second_level-nlin-3.mnc")

Finally, you will use Display to overlay your model_comp_lin_vs_quad.mnc file over your nlin-3.mnc average.

Here I have set my q-value threshold to 0.05. Areas in white/color are areas in which the more complex model (i.e. quadratic) better fits the data.