0

I am fitting mixed model in R. My model is:

model <- lmer(provision_rate~breeding_type+nestling_age+time+sex:nestling_age)+(1|nest)+(1|individual), data = provision)

Sex is the sex of parent. I want to graph the relation between the interaction term (sex:nestling_age) and provision_rate. I don't know how to add the predicted line according to model with other parameters set to their mean values to scatter figure in R?

laylarenee
  • 3,276
  • 7
  • 32
  • 40
Yu jin
  • 11
  • 1
  • what have you tried? Also, can you provide a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)? Dummy data is okay. – alexwhitworth Sep 14 '15 at 16:56

1 Answers1

0

One approach to this problem is to use the data manipulation tools in the merTools package in order to create a data.frame that simulates the relationship you want to explore. Here is an example using the VerbAgg data from the lme4 package, but you should be able to easily adopt it to your data and model.

data(VerbAgg)
fmVA <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 +
       (1|id) + (1|item), family = binomial, 
       data = VerbAgg)

We prep the data using the draw function in merTools. Here we draw the average observation from the model frame. We then wiggle the data by expanding the dataframe to include the same observation repeated but with different values of the variable specified by the var parameter. Here, we expand the dataset to all values of btype, situ, and Anger subsequently.

# Select the average case
newData <- draw(fmVA, type = "average")
newData <- wiggle(newData, var = "btype", values = unique(VerbAgg$btype))
newData <- wiggle(newData, var = "situ", values = unique(VerbAgg$situ))
newData <- wiggle(newData, var = "Anger", values = unique(VerbAgg$Anger))

head(newData, 10)

   r2 Anger Gender btype  situ id        item
1   N    20      F curse other  5 S3WantCurse
2   N    20      F scold other  5 S3WantCurse
3   N    20      F shout other  5 S3WantCurse
4   N    20      F curse  self  5 S3WantCurse
5   N    20      F scold  self  5 S3WantCurse
6   N    20      F shout  self  5 S3WantCurse
7   N    11      F curse other  5 S3WantCurse
8   N    11      F scold other  5 S3WantCurse
9   N    11      F shout other  5 S3WantCurse
10  N    11      F curse  self  5 S3WantCurse

In the next step we simply pass this new dataset to predictInterval in order to generate predictions for these counterfactuals and we bind these predictions to the original data.

plotdf <- predictInterval(fmVA, newdata = newData, type = "probability", 
        stat = "median", n.sims = 1000)
plotdf <- cbind(plotdf, newData)

Then we plot the predicted values against the continuous variable, Anger, and facet and group on the two categorical variables situ and btype respectively.

ggplot(plotdf, aes(y = fit, x = Anger, color = btype, group = btype)) + 
  geom_point() + geom_smooth(aes(color = btype), method = "lm") + 
  facet_wrap(~situ) + theme_bw() +
  labs(y = "Predicted Probability")

Which produces the following plot of fitted values:

enter image description here

jknowles
  • 479
  • 3
  • 13