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?"