Dear Russel and Henrik,

thank you very much for both your responses.
@ Henrik: I indeed was not aware of the Effects package, so it took me some time to follow up your suggestion. Indeed it worked and I have solved my issue. I can’t share my data for reasons of privacy, so below, I am providing a solution using the Hox (2002) popularity data based on your comments.

@ Russel: you are right, it makes much more sense for the use with the Effect package, NOT to explicate the interaction as a variabele, but as: var1:var2 .
With regard to your second comment: Thanks, it is indeed required to set the variable as factor (if not numeric), but please note that my understanding of the documentation is that you cannot do it as you proposed: There is a warning in the documentation of the Effects package not to use as.factor(var1), but to explicitly set the var to factor before running the model.
Also thank you for your note on the lsmeans package, I still need to follow up on that, but it sounds very promising.

I will now try and see whether I can get the interaction effects plotted directly from the predict function in the lme4 package, as suggested by Henrik. One reason to still pursue this, is that the Effects package is not flexible in changing labels since, according to the author, the plots were never actually ment to be used for publication purposes.

Finally, below my solution based on Henrik and Russel’s comments, using a reproducible example:

Best, 

Richard


# File: glmer_interaction_example01.R
# Date: July 13, 2014
# Purpose: to show how to visualize interactions from a glmer model
# Last change: -

# load packages
library(lme4)
library(effects)

#### load data ####
# load data from Joop Hox text book incorporated on Jeroen Ooms github
# see: http://jeroenooms.github.io/multilevel/ and 
# http://joophox.net/mlbook2/MLbook.htm
# The data are about popularity of pupils in classrooms and how this depends
# on characteristics such as extraversity
# for full description see: http://joophox.net/mlbook2/DataStories.pdf

pop2 <- read.csv("http://jeroenooms.github.io/multilevel/data/chapter2/popular2.csv")
summary(pop2)

#### creating a dummy var for whether you are popular or not ####
pop2$popd <- pop2$popteach > mean(pop2$popteach)

# setup model
m0 <- glmer(popd ~ (1 | class), 
            data = pop2, family= "binomial") # null model
#summary(m0)
m1 <- glmer(popd ~ sex + Cextrav + Ctexp + (1 | class),
            data = pop2, family= "binomial") # model with fixed effects
#summary(m1)
m2 <- glmer(popd ~ sex + Cextrav + Ctexp + Cextrav:Ctexp + (1 | class),
            data = pop2, family= "binomial") # model with interaction)
summary(m2)

#### estimating effects with effects package ####
m2fx <- allEffects(m2)

# plot the interaction
iact1 <- plot(m2fx, 'Cextrav:Ctexp',
                   multiline = TRUE, 
                   rescale.axis = FALSE, # have equal steps on y-axis
                   x.var="Cextrav", # set var for x-axis
                   xlab = "Extraversion Rating (centered)",
                   ylab ="Probability of being popular",
                   ylim = c(0, 1)) # set y-axis range between 0-1
iact1
# so basically:
# - extraverts are more likely to be popular
# - as teacher experience grows, the gap between levels of extraversion and the
# probability of being popular decreases: to put it boldly: more experienced
# teachers attach less value to levels of extraversion when rating their 
# students' popularity

# let's repeat this, but now have an interaction with a dummy variable
# create the dummy for extraversion
summary(pop2$extrav)
pop2$extrvd <- pop2$extrav > mean(pop2$extrav)

# run model with extraversion dummy
m3 <- glmer(popd ~ sex + extrvd + Ctexp + extrvd:Ctexp + (1 | class),
            data = pop2, family= "binomial") # model with interaction)
# eigen value warning, let's ignore it for the moment
summary(m3)

# estimating effects
m3fx <- allEffects(m3)
# results in error, since extrvd is no factor or numeric
pop2$extrvd <- as.factor(pop2$extrvd)
m3fx <- allEffects(m3)

# plot the interaction
iact2 <- plot(m3fx, 'extrvd:Ctexp',
              multiline = TRUE, 
              rescale.axis = FALSE, # have equal steps on y-axis
              x.var="extrvd", # set var for x-axis
              xlab = "Being extravert (TRUE) or not (FALSE)",
              ylab ="Probability of being popular",
              ylim = c(0, 1)) # set y-axis range between 0-1
iact2
# again the plot shows, that on average those who are extraverts are more
# likely to be popular than those who are not
# but the fact that the slope of the lines decreases (become more horizontal)
# with teacher experience, indicates that being extravert or not matters less
# for the probability that a teacher rates a pupil as popular

# EOF



On 12 Jul 2014, at 17:53, Lenth, Russell V <russell-lenth@uiowa.edu> wrote:

> I meant to chime-in on this a couple of days ago. It seems to me you are going to have trouble getting what you want with the model you have fitted. Instead, fit a model that has the interaction explicitly so that the effects package (or whatever you use) can tell you want to look at an interaction. That is, fit a model of the form
>  
>       model2 <- glmer(dep.var ~ var1 + var2 + ? + time + var1:time + (1 | reff), …)
>  
> If var1 and var2 are both numeric predictors, then I think the effects package probably offers what you want. However, I think you have var1 as a factor (you say a dichotomous variable). If you coded this as a dummy variable (0 or 1), then again you’ll have better luck changing it to a factor, e.g.
>  
>         model3<- glmer(dep.var ~ factor(var1)*time + var2 + ? + (1 | reff), …)
>  
> In that case, the lsmeans package might prove handy, as it has an lstrends function that can estimate and compare the time trends for the two levels of var1:
>  
>     library(lsmeans)
>     lstrends(model3, “var1”, “time”)
>     pairs(.Last.value)  ### For comparing the slopes
>  
> Hope this helps.
>  
> Russ
>  
> Russell V. Lenth  -  Professor Emeritus
> Department of Statistics and Actuarial Science  
> The University of Iowa  -  Iowa City, IA 52242  USA  
> Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017
>  
>  
> Message: 1
> Date: Thu, 10 Jul 2014 11:46:06 +0100
> From: Richard Zijdeman <richard.zijdeman@gmail.com>
> To: r-sig-mixed-models@r-project.org
> Subject: [R-sig-ME] visualising interaction with predict.merMod
> Message-ID: <56910434-600C-4ECB-B339-F3079F3DC193@gmail.com>
> Content-Type: text/plain; charset=windows-1252
>  
> Dear all,
>  
> I?ve fitted a model with the lme4 (version 1.1.7) package, using the glmer command. The only random component in the model is the intercept. The model contains an interaction with time and a dichotomous variable and I would like to visualise this interaction (the model actually contains more interactions and variables, but I don?t think this matters):
>  
> model1 <- glmer(dep.var ~ var1 + var2 + ? + time + var1_time + (1 | reff), data = data, family = binomial("logit?)) # where var1_time is the interaction
>  
> To visualise the interaction, I think I would need to predict new values from the model. To do so, I have used the actual data, but I don?t think that?s right:
>  
> pred.values <- predict(model1, type = "response")
>  
> I think I should create a new dataset, in which all variables are held at constant value (e.g. the mean), except for, in this case the main effects of the interaction and the interaction.
>  
> My questions are:
> - could someone please comment on my reflection (do I indeed need to create a new dataset)?
> - If I indeed would need to create such a new dataset, with average values, how would I do that (my model actually contains a lot of variables)?
>  
> Thanks in advance,
>  
> Richard


	[[alternative HTML version deleted]]

