Properly plotting an lm or lmer model predicted curve in R with ggplot

A challenge when running lm and lmer models in R is how does one properly visualize the "significant" effects found in a model when there are multiple covariates also included in the model.

The quick way to attempt this is ggplot(data = data, aes...) + geom_stat(method="lm") etc. however this method cannot handle covariates, which means any model line you produce will not be representative of the linear model you actually ran.

The effects package offers a solution to this, which I will show in the code below.

The way that this visualization is achieved is by "simulating" or "predicting" the model effects for hypothetical subjects using the Effect function. For any predictor not specified in the Effect options to stratify, the prediction will be made for the mean of the population. This means predictors you don't intend to graph should be left out of the Effect call.

In the below example, we run a model for total GM volume of subject brains, with a poly age effect a sex effect, a covariate of total brain volume, and random effects of field strength and subject. When generating an Effects output, we stratify for sex, and generate an age prediction every 1 year (the default sampling is too sparse). The resulting prediction curve will be for the mean brain volume of the population subject, since it is not included in the Effect stratification call, similarly for random effects.

The Effect output contains the confidence interval for the prediction, so we can plot proper bounds.

library(lme4)
library(ggplot2)
library(effects)
library(lmerTest)

#Load data and such here ##

# Here gmvolume could by ananatomy volume from MAGET, a voxel value generated from mincGetVoxel or mincGetWorldVoxel from DBM data
# Or a vertex value from vertexTable
gmvolume_model = lmer(gmvolume ~ poly(age, 2) + sex + brainvolume + (1|DICOM_field_strength) + (1|subject), data = data)
# Note that, for the effects call below, ONLY THE PREDICTORS YOU WISH TO GRAPH should be called. If you include more, R will chunk the data by other predictors and you will not get a single line, but many fractured ones.
fitdata_gmvolume = as.data.frame(Effect(c("sex", "age"),gmvolume_model, xlevels=list(age=seq(6,53,1))))

ggplot(data = data, aes(y=gmvolume,x=age, color=sex, group=sex, fill=sex)) +
  geom_point() +
  geom_line(data = fitdata_gmvolume, aes(y=fit)) +
  geom_ribbon(data = fitdata_gmvolume, aes(y=fit, ymin=lower, ymax=upper), alpha=0.4) +
  xlab("Age") +
  ylab(bquote('Gray Matter Volume'~(mm^3))) + 
  labs(fill = "Sex", group="Sex", color="Sex")

While this example is done for a lmer call, the lm method is exactly the same.

THE CODE BELOW IS DEPRECIATED AND NOT VALIDATED

Plotting a voxel's trajectory after running mincLMER If you want to use the effects package to plot the trajectory of the model in one specific voxel, this is what you need to do:

1. After you run your mincLMER and FDR correction, you can save the stats of interest and visualize them with Display. If you don't know how to save your stats go here, if you don't know how to visualize them with Display go here.

2. In Display, locate the voxel of interest that you would like to plot and write down the world coordinates. The world coordinates (Xw, Yw, Zw) appear at the bottom left side of the Display interface when hovering over the voxel.

3. You'll need to run an LMER (not mincLMER) in the specified voxel, therefore, you need to load your data exactly as you loaded your data previously, to run the mincLMER. In my case I subsetted my data, so I would need to make sure I subset the data in the same way and indicate the same baselines.

Example of code:

#model I ran
modns <- mincLmer(resampled_log_full_det ~ Genotype*ns(Age_days,2) + Sex + (1|ID), 
                  data = jacdata, 
                  mask = mask,
                  parallel = c("sge", 200),
                  REML = TRUE)

#fdr correction
fdrmodns  <- mincFDR(mincLmerEstimateDF(model = modns), mask = mask)

#save stats of interest
mincWriteVolume(modns,"GenotypeModNs.mnc", column = 'tvalue-Genotype3Tg:ns(Age_days, 2)2')

After I selected the voxel of interest, which in this case is located in the Left Olfactory Bulb, I would enter coordinates:

#coordinates
Xw <- -1.0
Yw <- 6.1
Zw <- 2.0

Now, I will need to run again the LMER in the specified voxel, so first I would create a variable containing the voxel information from my data (my data is 'jacdata$resampled_log_full_det'). Then I unlogged the values.

#get voxel data
voxel <- mincGetWorldVoxel(jacdata$resampled_log_full_det, Xw, Yw, Zw)
#unlogged values
voxel_unlogged <- 10^(voxel)

Then I printed the voxel coordinates because I can use them in the future to create slice figures:

#print voxel coordinates and save coordinates of the voxel "attr(,"voxelCoord")"
print(voxel_unlogged)

This is what I saved from printing the voxel coordinates (I usually save this within my script to easily find the information in case I need it later):

#dimension 1, 2, 3
#[1]  62 143  53    

now I run the LMER in the voxel of interest, see how my model is exactly the same I used with mincLMER:

model_voxel = lmer(voxel_unlogged ~ Genotype*ns(Age_days,2) + Sex + (1|ID), data = jacdata)

Now I can easily use the Effect function as instructed before:

#stratify for Genotype, and generate an age prediction every 1 year
fitdata_jacob= as.data.frame(Effect(c("Genotype", "Age_days"),model_voxel, xlevels=list(Age_days=seq(50,200,1))))

png("genotypeModNsLeftOlfactoryBulbFittedWithEffects.png", units= "in", width = 4, height =3, res = 300)  #this are indications to save my image as png

ggplot(data = jacdata, aes(x=Age_days, y=voxel_unlogged, color=Genotype, group=Genotype)) +
  geom_point(alpha=0.8, size = 1) +
  geom_line(data = fitdata_jacob, aes(y=fit),size=1, alpha=0.8) +
  geom_ribbon(data = fitdata_jacob, aes(y=fit, ymin=lower, ymax=upper, color=NULL), alpha=0.1) +
  scale_color_brewer(palette = "Set2")+
  labs(x = "Age in days", y= "Jacobian value", title = "Peak Voxel in Olfactory Bulb-Left")+
  theme(plot.title = element_text(hjust = 0.5))

dev.off()  # This function closes the specified png plot I indicated above

This is how the figure looks like:

The following code can be used to plot a single slice. This slice is the same slice I was looking at when selecting the voxel I previously plotted.

#Single Slice from the previous figure, see voxel coordinates for slice number
#dimension 1 is Sagital view, dimension 2 is Coronal view, dimension 3 is a horizontal view

png("genotypeSliceLeftOlfBulb.png", units= "in", width = 5, height =3, res = 300) 

sliceSeries(nrow = 1, slices = 143, dimension = 2) %>%
  anatomy(avol1, low = 700, high = 1400) %>%
  overlay(mincArray(GenotypeModNs, "tvalue-Genotype3Tg:ns(Age_days, 2)2"),    
  low=2.74, high=3.56, symmetric = T) %>%       # t score range; symmetric = T means show positive and negative t values
  legend("t-statistics 5% to 1%") %>%
  contourSliceIndicator(avol1, c(700,1400)) %>% 
  draw()

dev.off()

This is how the figure looks like: