Solved – how to plot interactions with continuous and categorical predictors in mixed models

ggplot2interactionmixed modelmultilevel-analysisr

I'm very unsure how to plot mixed-level data consisting of a mixture of categorical and continuous predictors, so any help would be appreciated.

This is the data

subject <- factor(rep(c(1,2,3,4,5,6),each=12))
dep <- c(0.3763244126,0.2185001692,0.4191841742,0.9812978664,0.7429273683,0.6254715486,0.6200958213,0.4693300191,0.1779032899,0.0035873980,0.8821949826,0.4818012617,0.0008013437,0.6280700732,0.7126500814,0.4984349359,0.2457996449,0.3085733312,0.5903398243,0.3704800352,0.8215325437,0.0445236221,0.1849731791,0.3670945817,0.0022268933,0.1630332691,0.9734050406,0.2638539758,0.8550054496,0.9413964085,0.4548943471,0.0440815873,0.5222098769,0.6553600784,0.6853486744,0.0571945074,0.0923124240,0.6976544929,0.9257440316,0.5658967043,0.0636543627,0.1038574059,0.0662497468,0.9165439918,0.5200087291,0.9528015053,0.5347318368,0.1848373057,0.9948602219,0.9633110918,0.1482162909,0.9000614029,0.0898618386,0.7975253051,0.8334557347,0.8629821099,0.0001795699,0.2488384889,0.6382902598,0.1103540359,0.2199716354,0.2737281912,0.5694398067,0.7940423761,0.4906677457,0.5191186895,0.4770589883,0.2823238128,0.2458788699,0.6363522802,0.0306954833,0.6979198116)

f1 <- factor(rep(c("Female","Male","Female","Male","Male","Female"), each=12))
f2 <- c(0.098788608,0.934606288,0.145045152,0.841969882,0.498234471,0.562897249,0.359740488,0.082046687,0.183987342,0.082418820,0.173424633,0.799291329,0.041450568,0.686708743,0.352092230,0.823550310,0.650857094,0.331705763,0.659111451,0.745187314,0.066165065,0.870759966,0.154977488,0.031703774,0.065251788,0.707452073,0.564604314,0.224798417,0.656363138,0.047954841,0.500513114,0.923316812,0.706629266,0.561530974,0.670860932,0.414969178,0.709973062,0.452946384,0.187624344,0.278656351,0.562138433,0.655193272,0.014868182,0.518697012,0.414113229,0.273464316,0.844080831,0.962636550,0.952739605,0.728627219,0.761122951,0.309977150,0.755239042,0.208627128,0.481429897,0.376021223,0.753871400,0.164842337,0.921081061,0.859677311,0.600462073,0.119193708,0.276722102,0.854752641,0.962710853,0.956277061,0.228313179,0.920405764,0.001594131,0.104930433,0.241548888,0.549643015)
f3 <- factor(rep(c("day1","day2","day3","day4"),each=3, times=6))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2, f3=f3)

m <- lmer(dep ~ f1*f2*f3 + (1|sub), data=data)

how can I plot the 3-way interaction f1*f2*f3?

I tried using ggplot2 like this

ggplot(aes(x=f2,y=dep,color=f1),data=data) + geom_smooth(method="lm") + facet_grid(".~f3")

but I don't think this is taking into account the fact that this is a mixed-model design, is it?

Best Answer

You can specify that you'd like a line for each interaction of f1 and f3 using the group aesthetic and interaction function in R:

ggplot(data, aes(x=f2, y=dep, group = interaction(f1, f3))) +
    geom_smooth(method="lm")

This will create a plot with (in your case 2 * 4 = 8 lines). Then you can use features other than color to distinguish between the factor levels. In your example, you facet, which is OK, but makes it a bit hard to compare the lines. You can also use shape (for points) or line weight / type. Here's an example that uses color for the day and line type for gender:

ggplot(aes(x=f2, y=dep, group = interaction(f1, f3))) + 
    geom_smooth(method="lm", alpha = 0.1, aes(lty=f1, color=f3))

Note that I added alpha = 0.1 to make it a bit easier on the eyes.

Related Solutions

Solved – Interpreting 5-way Mixed Model ANOVA

Generally, you should start from the highest order interactions. You are probably aware that it is usually not sensible to interpret a main effect A when that effect is also involved in an interaction A:B. This is because the interaction tells you that the effect of A actually depends on the level of B, rendering any simple main effect interpretation of A impossible. In the same way, if you have factors A, B, C, then A:B should not be interpreted if A:B:C is significant.

Thus, when you have a 5-way interaction, none of the lower-order interactions can be sensibly interpreted. Therefore, if I understand you correctly and you have interpreted your lower order interactions, you should probably not continue along those lines.

Rather, what you can do is to split up your data set and continue to analyze factor levels of your data set separately. Which of the factors you use to split up the dataset is arbitrary, but often it is very useful to split up the data for each variable and assess what you see. In your example, you might start with sex, and calculate an ANOVA for males, and another one for females (each ANOVA contains the 4 remaining factors). Just as well, you could split up the data according to ethnicity (one ANOVA for Asian, one for Caucasian). You could also split up by one of the within-subject factors.

I will assume that you have decided to split the data by sex (just to continue with the example here). Then, assume that for males, you get a 4-way interaction. You would then go on to split up the male data by one of the remaining variables (say, ethnicity). You would then calculate ANOVAs for male Asians (over the remaining 3 factors), and for male Caucasians.

Importantly, if you get only a lower-order interaction, then you are only "allowed" to analyze these further. This is because the other factors did not show significant differences. Thus, if your males ANOVA gives you only a 2-way interaction, then you would average over the other factors and calculate only an ANOVA over the 2 interacting factors (and, because we are in the male part of the ANOVAs, this would be for the males alone).

For the females, everything may look different, and so the decision which follow-up ANOVAs to calculate is separate for this group. So, what you did for males should be done for females in the same way ONLY if you got the same interactions.

Thus, you will potentially have a lot of ANOVAs, and it might not be easy to decide which ones to report. You should report 1 complete line down from the hightest interaction to the last effects (possibly t-tests to compare only 1 of your factors at the end). You should not usually report several lines (e.g., one starting the split-up by sex, then another one starting by ethnicity). However, you must report a complete line, and cannot simply choose to report only some of the ANOVAs of that line. So, you report one complete analysis, not more, not less. Which way to go in terms of splitting up / follow-up ANOVA is a subjective decision (unless you have clear hypotheses you can follow), and might depend on which results can be understood best etc.

Solved – Comparing slopes in mixed-effect model

One short thing: Day + Hour + Day*Hour can and should be replaced by Day*Hour (or Day + Hour + Day:Hour), because * allready includes the main effects see the formula help page. Second, conditioning on ID doesn't make sense, because you only have one observation per ID. You should condition on genotype and treatment. As you wrote down a random intercepts model (without random slopes) I stick to this in my proposed model:

lmer(DV ~ 1 + Day*Hour + (1|treatment) + (1|genotype), data)

Best Answer

Related Solutions

Solved – Interpreting 5-way Mixed Model ANOVA

Solved – Comparing slopes in mixed-effect model

Related Question