Equivalent to ANOVA for non-linear (quadratic) relationships

Question

I am trying to perform the equivalent of a repeated-measures ANOVA using data that have a non-linear relationship. The independent variable position runs from -20 to +20, and we have measured a dependent variable performance in 17 subjects. What I'd like to know is the extent to which position affects performance. I'm not interested in whether there is a specific position at which all subjects have a better performance; in fact I have reason to suspect that the optimal position for each subject is idiosyncratic.

Assumptions

• The relationship between position and performance is quadratic.
• Each subject's performance is at its lowest at a specific position, and gets better either side of this (i.e. we will always see a 'U' shape; never '∩'). (The quadratic equation is always positive.)
• The position at which any given subject's performance is at it's lowest is idiosyncratic, and in some cases, our range of positions (from -20 to +20) may not include this position (i.e. the inflection point may be outside the sampled positions). (The x-intercept varies between subjects.)
• Some subjects' performance is better than others, even at their lowest position. (The y-intercept varies between subjects.)
• The rate at which performance improves either side of this idiosyncratic lowest position is also idiosyncratic. (The quadratic term of the equation varies between subjects.)

To illustrate, here are some examples:

Each curve above can each be described by a quadratic equation, and, with the exception that they are all U-shaped, they differ in every other conceivable way.

Given this, I need to come up with a statistical test to prove (or disprove) that changing position has an effect on performance.

The problems

I see two problems. First, even if the relationship between position and performance were linear, the idiosyncratic nature of the position at which performance is lowest means that even an ANOVA wouldn't be appropriate. i.e., although performance is our dependent variable, we are not interested in whether there is a between-subjects effect (i.e. a position where performance is worst), but rather we are looking to know the extent to which position affects performance, regardless of whether the direction of the effect is different between individuals. In other words, some subjects may have better performance at lower positions, and others at high positions. It's not the absolute position we're interested in, but the relative differences within any given subjects. One possible solution to this would be to run a different test for each of our subjects (there are 17 of them). However, this increases the chance of a Type I error.

The second issue is of course the fact that the relationship is non-linear. So I need to find a statistical method that can model quadratic relationships. I'm working in R. Two possibilities that I'm aware of are to use a non-linear mixed model, or a generalised additive model with a quadratic link function (whatever that is called in real terms). However, my knowledge of both is poor, and I may be barking up the wrong tree entirely.

Answer 1

What you are looking for is a regression model that can handle within-subjects observations. An ANOVA is not appropriate here since you have a continuous predictor. You are looking for a multilevel model (also called mixed-models or hierarchical models).

What you can do is specify a model where your dependent variable is a function of an intercept, coefficient for a main effect, and then a coefficient for a quadratic relationship. You can let the intercept and coefficients vary by individuals. That is, the coefficients then are predicted by other variables themselves. In the present case, you are just letting them vary (i.e., only an intercept and error will be predicting these coefficients). I highly recommend reading Chapter 2 from Joop Hox's Multilevel Modeling, which can be found here.

Here is how you would do it in R:

First, let's generate some data that somewhat resemble the data that you have. I am assuming that you have 17 participants, and each participant is measured at 30 different occasions (I call these "trials").

set.seed(1839)
n1 <- 30 # number of observations per person
n2 <- 17 # number of people in the study
dat <- data.frame(expand.grid(id = 1:n2, trial = 1:n1)) # make matrix of p ids and trial ids

# generate random slopes and intercepts, one for each person
b0 <- rnorm(n2, 0, 5)
b1 <- rnorm(n2, 0, 5)
b2 <- runif(n2, 0, 1)

# initialize empty columns for position and performance
dat$position <- 0
dat$performance <- 0

# generate data
for (i in 1:n2) {
  x <- runif(n1, -20, 20) # make position
  dat$position[dat$id == i] <- x # assign to data
  # create performance, with a little bit of error
  dat$performance[dat$id == i] <- b0[i] + b1[i] * x + b2[i] * x ^ 2 + rnorm(n1, 0, 10)
}

Let's look at the head and tail of the data:

> rbind(head(dat), tail(dat))
    id trial    position performance
1    1     1 -15.9075065    8.259001
2    2     1  11.0343883   56.287865
3    3     1   7.5321399   62.576056
4    4     1  11.7032195  157.332692
5    5     1   1.5328204   -8.821568
6    6     1   8.1800251   52.460136
505 12    30  -8.5917467   -7.888843
506 13    30 -15.1133658   90.140682
507 14    30  -1.5184124  -20.156739
508 15    30  17.1828493  128.683005
509 16    30  15.2308364  149.121488
510 17    30  -0.8624032   18.709641

You can see that this is in "long" format (or, what it is called in the data science world, "tidy" format), because every measurement is its own row. Note that we have an id variable that tells us what participant the data came from. We will use this in multilevel modeling to tell us what to group the observations by.

We can then graph the quadratic relationship, separating each line by id. This looks somewhat like what you drew in your post:

library(ggplot2)
ggplot(dat, aes(x = position, y = performance, group = factor(id), color = factor(id))) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, size = .7)

It is OK if your data do not look exactly like this, as the model will hold regardless.

Now, you are ready to run a multilevel model. I suggest reading this great post, Using R and lme/lmer to fit different two- and three-level longitudinal models, by Kristoffer Magnusson, which covers an introduction on how to use the lme4 and lmerTest R packages. Don't let the "longitudinal" title fool you—you can use these same models to nest within person that are cross-sectional.

# run mixed model with quadratic effects
library(lme4)
library(lmerTest)
model <- lmer(performance ~ position + I(position ^ 2) + # fixed effects
                (1 + position + I(position ^ 2) | id), # random effects
              data = dat)

performance ~ position + I(position ^ 2) tells lme4 that you want to estimate fixed effects for position and position ^ 2 (i.e., the quadratic term). This means that it will give you the average coefficient for each of these in the summary of the model.

Next, the (1 + position + I(position ^ 2) | id) tells lme4 that you want to estimate the same effects (including an intercept, which is 1) as random effects that are grouped by (using |) the id. The random effects means that you are letting these coefficients differ by person. The summary object will give you the average coefficients, but you can extract the coefficients for each person, too.

Let's look at the summary:

> summary(model)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to degrees of
  freedom [lmerMod]
Formula: performance ~ position + I(position^2) + (1 + position + I(position^2) |      id)
   Data: dat

REML criterion at convergence: 3975.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8371 -0.6241 -0.0525  0.6124  2.8863 

Random effects:
 Groups   Name          Variance Std.Dev. Corr       
 id       (Intercept)   20.49317 4.5269              
          position      20.79349 4.5600   -0.06      
          I(position^2)  0.06644 0.2578    0.19 -0.17
 Residual               87.47573 9.3528              
Number of obs: 510, groups:  id, 17

Fixed effects:
              Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)    0.92948    1.27400 16.09000   0.730    0.476    
position       1.60834    1.10658 16.00100   1.453    0.165    
I(position^2)  0.36147    0.06262 16.00300   5.772 2.85e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) positn
position    -0.049       
I(positn^2)  0.142 -0.165

The Fixed effects: part is where the hypothesis tests for a study are usually done. We can see that the quadratic relationship between position and performance is = 0.36147. This is significant, as the p-value is 2.85e-05.

What is cool, though, is that you can look at each own person's coefficient! Under the Random effects: part, you can see that the variance of this quadratic effect is 0.06644. This means that the variance around the mean effect of .36 is about .07. We can look at every person's individual coefficient for the random slope:

> coef(model)$id
   (Intercept)    position I(position^2)
1    5.9534256  6.49409694    0.45250260
2   -0.9840823  2.41624005    0.23704314
3    3.1271227  7.17657290    0.03244345
4   -5.0344391  6.75137341    0.63935341
5   -3.8292241 -3.26936324    0.02005057
6    0.8090050 -1.87472471    0.83077897
7   -6.7428209  0.92931149    0.05555921
8    1.1568852  3.25817546    0.49851941
9    2.4236226 -1.89218156    0.37627493
10  -1.6531994  8.23137550    0.07024075
11   1.4485008 -5.25734473    0.32639520
12   5.5198297 -0.07633506    0.10475536
13   0.8987842 -2.18005992    0.27164383
14  -0.0215908  9.12063733    0.32520260
15   3.7644442 -3.77068022    0.64417213
16   0.6251266 -0.69074713    0.71894030
17   8.3397560  1.97541719    0.54108187

The participant with the id of 5 had the smallest quadratic effect (.02). You can actually see that in the plot above, as the line for id of 5 is the closest to being a striaght, linear line.

One of the nice parts about multilevel modeling is that you could, in turn, predict this variance in coefficients, if you are interested. You could now turn to questions like: "What kind of people show the quadratic effect, and which do not?"