Get some data

Heights of Boys in Oxford. Data comes from the {nlme} package. This data frame contains the following columns:

data(Oxboys, package = "nlme")

quick plot

Height obviously appears to be associated with age.

library(ggplot2)
ggplot(Oxboys) +
  aes(age, height) +
  geom_line(aes(group = Subject)) + 
  geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

fit a non-linear model

Model height as a function of age. A non-linear model is probably not necessary here, but we’ll do it anyway.

library(lme4)
m <- lmer(height ~ poly(age, 2) + (1 | Subject), data = Oxboys)
summary(m, corr = F)
## Linear mixed model fit by REML ['lmerMod']
## Formula: height ~ poly(age, 2) + (1 | Subject)
##    Data: Oxboys
## 
## REML criterion at convergence: 922.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.58149 -0.63386 -0.02158  0.57060  2.26937 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 65.569   8.097   
##  Residual              1.641   1.281   
## Number of obs: 234, groups:  Subject, 26
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    149.519      1.590  94.022
## poly(age, 2)1   64.526      1.281  50.371
## poly(age, 2)2    4.194      1.281   3.274

create a simple effect plot

library(ggeffects)
plot(ggpredict(m, terms = "age [all]"))

and now do post-hoc comparison

Compare height at standardized age of 0.5 to -0.5. Notice we need to specify the age values to compare.

library(emmeans)
emm <- emmeans(m, specs = "age", at = list(age = c(0.5, -0.5)))
contrast(emm, method = "pairwise")
##  contrast           estimate   SE  df t.ratio p.value
##  age0.5 - (age-0.5)     6.52 0.13 206  50.293  <.0001
## 
## Degrees-of-freedom method: kenward-roger