Skip to content

LMER: Visualizing three-way interaction

5 messages · Klemens Knöferle, Evan Palmer-Young, Thierry Onkelinx +2 more

#
Hi all,

I'm trying to visualize a three-way interaction from a rather complex
linear mixed model in R (lmer function from the lme4 package; the model has
a complex random-effects structure). The interaction consists of two
continuous variables and one categorical variable (two experimental
conditions).

So far, I have graphed the interaction via two 3D-surface plots using
visreg2d from the visreg package. But my reviewers found these plots
confusing and asked for a different illustration, such as conditional
coefficient plots (i.e., plots of the strength of coefficient 1 as
coefficient 2 increases).

I've tried to find a package that allows me to create these kind of plots,
but failed. The existing packages only allow coefficient plots for two-way
interactions (for instance the interplot package;
https://cran.r-project.org/web/packages/interplot/vignettes/interplot-vignette.html).
That means I only get a conditional coefficient plot of the two-way
interaction, collapsed across both levels of the categorical variable.

Is there a package for my case? If not, I probably have to manually extract
fitted values from my model (e.g., using broom) and somehow plot them in
ggplot2. But I don't really know how to do this, whether or not to take
into account random effects (and how), etc. Any ideas would be much
appreciated...

Klemens Kn?ferle, Ph.D.
Associate Professor - Department of Marketing
BI Norwegian Business School
Visiting address: Nydalsveien 37, 0484 Oslo
#
Dear Klemens,
I also have been asked to remove plot3d plots that I thought were
stupendous, but I guess not everybody likes them!
How about two two-way interaction plots for the regressions of continuous
variable 1 vs 2, with separate panels for the levels of the categorical
variable?
You can get fitted model values and CI's with the lsmeans function, using
the "at" argument to specify covariate values.
If the plot is too bland, you can sprinkle the raw data on top as an extra
layer of points.
Best wishes, Evan


On Wed, Apr 12, 2017 at 5:41 AM, Klemens Kn?ferle <knoeferle at gmail.com>
wrote:

  
    
#
Dear Klemens,

Here is my ?0.02

library(ggplot2)
dataset <- expand.grid(
  x = seq(0, 1, length = 41),
  y = seq(0, 1, length = 41),
  z =factor( c("A", "B"))
)
dataset$fit <- with(dataset,
  ifelse(
    z == "A",
     x - 2 * x^2 + 0.5 * x * y + 3 * y - y ^ 2,
    -x - 1 * x^2 - 2 * x * y - 3 * y + y ^ 2
  )
)
ggplot(dataset, aes(x = x, y = y, fill = fit)) +
  geom_raster() +
  facet_wrap(~ z) +
  scale_fill_gradient2()
ggplot(dataset, aes(x = x, y = y, z = fit)) +
  geom_contour(aes(colour = ..level..), binwidth = 0.25) +
  facet_wrap(~ z) +
  scale_colour_gradient2()
ggplot(dataset, aes(x = x, y = y)) +
  geom_raster(aes(fill = fit)) +
  geom_contour(aes(z = fit), binwidth = 0.25) +
  facet_wrap(~ z) +
  scale_fill_gradient(low = "black", high = "white")

Best regards,

Thierry
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given body
of data. ~ John Tukey


2017-04-12 19:06 GMT+02:00 Evan Palmer-Young <ecp52 at cornell.edu>:
#
Hello Klemens,

Have a look at the effects package and this document:

https://www.jstatsoft.org/article/view/v008i15/effect-displays-revised.pdf

Using the Owles data (two continuous variables and one categorical, but
fit to a bernoulli response, which doesn't matter for your problem,
because it results also in a continuous variable for the fit), you can
cut your data in sections:

library(effects)

data("Cowles")
cowles.mod <- glm(volunteer ~ sex * neuroticism * extraversion,
                  data=Cowles, family=binomial)

# standard effects plot via lattice:
plot(effect("sex*neuroticism*extraversion", cowles.mod,
            xlevels = list(neuroticism = 0:24,
		           extraversion = seq(0, 24, 6))),
     multiline = TRUE, ylab = "Probability(Volunteer)")

# easily to get a dataframe to do it yourself with ggplot2:
ne.effect <- data.frame(effect("sex*neuroticism*extraversion",
                               cowles.mod, 		
			       xlevels=list(neuroticism = seq(0, 24, 6),
			       extraversion=0:24)))

library(ggplot2)
ne.effect$neuroticismF <- factor(ne.effect$neuroticism,
				labels = paste0("extraversion: ",
				unique(ne.effect$neuroticism)))

ggplot(ne.effect, aes(x = extraversion, y = fit, group = sex)) +
theme_bw() + ylim(0,1) +
    geom_line(aes(linetype = sex), size = 1.2) +
    geom_ribbon(aes(ymin = lower, ymax = upper, linetype = sex),
    color = "black", fill = "grey", alpha = 0.1) +
    facet_wrap(~ neuroticismF)

Extracting data with different levels for continuous variables does also
work with objects from lmer if you want to only concentrate on the fixed
effects estimates. And this attempt at plotting three variables is the
most standard one for the every day reviewer from my opinion...

Good luck!



Am 12.04.2017 um 20:54 schrieb Thierry Onkelinx:
--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Rainer Schoppik
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING
#
Dear Klemens,


You might want to take a look at the effects package, by John Fox and
collaborators. It allows to visualize different types of interactions,
including conditional plots, with a lot of control over what and how, deals
with fits from lmer, and can also plot partial residuals, which I find
extremely useful for diagnostic purposes.


Best,

R.
On Wed, 12-04-2017, at 11:41:00, Klemens Kn?ferle <knoeferle at gmail.com> wrote: