I have a data structure that requires a mixed effects model - specifically, I have 3 crossed (rather than nested) factors with random intercepts. The number of levels for each of the 3 random factors is:

level 1: 25

level 2: 29

level 3: 58

I am trying to determine how linear is the relationship between a variable X and a dependent variable Y. The first thing I thought of is a scatterplot. However, because of the data structure, I am not sure that I can just plot X v. Y. How to address? Should I plot separate scatterplots by each level? For example for level 3, show the relationship for each of the 58 members of level 3? How to do that?

In the mixed effects model, I will then include X and X^2 in the model, and determine if X^2 is statistically significant, and then if it is, include a cubed term, and so on. When assessing non-linearity, is there anything else I should do?

Thank you!

Note: this is not a longitudinal dataset.

  • 1
    $\begingroup$ At what level does X vary? Does it vary within each of the levels of each your grouping variables? $\endgroup$ – Ben Bolker Nov 11 '18 at 0:17
  • $\begingroup$ Thank you. Yes, it varies within each of the levels of each of the grouping variables. $\endgroup$ – Laura Nov 11 '18 at 1:38

Graphically, I would suggest three spaghetti plots - scatterplots with values connected in order of X, with lines grouped by your three grouping variables. For example, in R:

data("sleepstudy", package="lme4") 
lattice::xyplot(Reaction~Days, groups=Subject,sleepstudy,type="b")
## or
ggplot(sleepstudy, aes(x=Days, y=Reaction, group=Subject)) +

I personally prefer this to putting each group in its own facet (Reaction~Days|Subject in lattice or facet_wrap(~Subject) in ggplot2), but tastes differ.

Since your levels are crossed, it's probably not practical to try to graphically decompose the effects of grouping any further.

As for the modeling: why are you assuming that only the intercept varies among groups? If X varies within groups, you should probably include the full term in your random-effects model, e.g.

Y ~ poly(X,2) + (poly(X,2)|L1var) + (poly(X,2)|L2var) + poly(X,2)|L3var)

although the practicality of this will depend on the size of your data set (this model fits 3 fixed + 3*(3*(3+1)/2) = 18 random-effects parameters ...)

  • $\begingroup$ Thank you!!! I will think more about the model set up. $\endgroup$ – Laura Nov 11 '18 at 14:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.