What is the most reliable way to visualize the mixed model?

Question

this may be a beginner’s question, but any help is appreciated! First, how to interpret the interaction between continuous variables? Such as WORD_LENGTH_SCALE:Freq_SCALE and WORD_LENGTH_SCALE:Freq_SCALE:Freq_2nd?

Here are my codes (Freq_2nd is not applicable to a word length condition. That is, when the word length = 1, Freq_2nd is NA. Thus, i wrote is as (WORD_LENGTH_SCALE:Freq_2nd)):

fm = lmer (response  ~ WORD_LENGTH_SCALE * Freq_SCALE * 
    (WORD_LENGTH_SCALE:Freq_2nd) + (1| subjects) + ( 1 | trail_no),  
    data = data_01, REML = TRUE)

and the results:

 Fixed effects:
                                       Estimate Std. Error t value
(Intercept)                            4.205600   0.015334 274.272
WORD_LENGTH_SCALE                     -0.025104   0.015568  -1.613
Freq_SCALE                            -0.004509   0.001953  -2.309
WORD_LENGTH_SCALE:Freq_SCALE          -0.019313   0.008061  -2.396
WORD_LENGTH_SCALE:Freq_2nd             0.004601   0.002513   1.831
WORD_LENGTH_SCALE:Freq_SCALE:Freq_2nd  0.003309   0.001370   2.416

Correlation of Fixed Effects:
                           (Intr) WORD_LENGTH_SCALE F_SCAL WORD_LENGTH_SCALE:Fr_SCALE WORD_LENGTH_SCALE:F_2
WORD_LENGTH_SCALE          -0.029                                                                          
Freq_SCALE                  0.145  0.093                                                                   
WORD_LENGTH_SCALE:Fr_SCALE -0.015  0.849             0.051                                                 
WORD_LENGTH_SCALE:F_2      -0.024 -0.958            -0.295 -0.825                                          
WORD_LENGTH_SCALE:F_SCALE: -0.033 -0.785            -0.284 -0.957                      0.839  

Second, I used the ggeffects, sjPlot and prediction function to plot graphs, but they seem different from the model. In the model, the Freq_SCALE and WORD_LENGTH_SCALE are negative, while in the plot, the Freq_SCALE is negative but the WORD_LENGTH_SCALE is positive. Why is this the case?

Third, what is the most reliable way to visualize the mixed model? sjPlot and the prediction function give me different plots. So I am not sure what to do ...

sjPlot

prediction:

Answer 1

You will have to apply your knowledge of the subject matter to decide the most reliable way to proceed. First some general considerations, then some warnings about your application.

In general, with interactions among three predictors there is no single association of any predictor with outcome. The association of one of them with outcome depends on the values of the other two.

For example, the single-predictor coefficient for WORD_LENGTH_SCALE of -0.025104 only holds when the values of the other 2 predictors are 0. Some particular values of Freq_SCALE and Freq_2nd were necessarily assumed to produce the first plot of response versus WORD_LENGTH_SCALE. Presumably, the assumed values of those other predictors were such that the interactions of WORD_LENGTH_SCALE with the other predictors led to a net positive slope.

I don't use sjPlot and the details of your prediction method aren't completely clear, but what you have to do is independent of those details. To understand the plots that you show, you need to read the help pages to find out the particular assumptions those methods are making about all of the predictors. Your plots show at most values of the predictor along the horizontal axis and one other predictor. You need to know values of all 3 predictors to see if the predictions make sense.

For comparison, use your coefficients to write the point estimates of response from your model directly in terms of WORD_LENGTH_SCALE (ws), Freq_SCALE (fs), and Freq_2nd (f2):

response <- function(ws,fs,f2) {4.205600 - (0.025104 * ws) - (0.004509 * fs) - (0.019313 * ws * fs) +  (0.004601 * ws * f2) +  (0.003309 * ws * fs * f2)}

and plug in appropriate numbers.

The most useful way to present results of models with interactions is to illustrate with intelligently chosen combinations of predictor values that are expected in practice. Don't depend on hidden assumptions that software packages make about predictor values that you don't directly specify. Use your knowledge of the subject matter to make sets of explicit choices about all 3 predictors that will best illustrate your results to your audience. Software should be able to make predictions based on your own choices of values.

Now for the warnings. First, if WORD_LENGTH_SCALE (or any other predictor) is derived from only a small set of integer values, you probably don't want to show plots that imply a continuous relationship with outcome. (You might even consider whether modeling as a single continuous predictor makes sense.) Separate plots and standard errors for individual values would be preferable.

Second, you say:

Freq_2nd is not applicable to a word length condition. That is, when the word length = 1, Freq_2nd is NA. Thus, I wrote is as (WORD_LENGTH_SCALE:Freq_2nd).

I'm not sure that your model handles this situation properly. In general, software will remove any data rows that have NA values in any variable included in the model before fitting. You might need to tweak your model and your data coding accordingly.