I am using a linear mixed effects model with a natural spline function for age to describe the trajectory of an outcome y
(measured in grams) across time (age
in years) in a sample of individuals. My model includes random-effects for splines (with fewer degrees of freedom than the fixed-effects splines).
I would like to estimate the subject-specific (person-specific) age at peak velocity from this model.
A previous thread describes how to estimate the subject-specific age at peak velocity when the model contains a quadratic random-effects age term.
The code below provides the example dataset and shows how the model was fit, how the mean trajectory was plotted and how the mean velocity curve and mean age at peak velocity were calculated.
How can I estimate the age at peak velocity for each individual from this model?
# LOAD DATASET
library(RCurl)
dat <- read.csv(("https://raw.githubusercontent.com/aelhak/data/main/data.csv"))
# FIT MODEL
library(nlme)
library(splines)
model <- lme( # linear mixed-effects model
y ~ ns(age, df = 7), # FE: natural spline function for age with 6 knots
random = ~ ns(age, df = 3) | id, # RE: natural spline function for age with 2 knot
method = "ML", data = dat)
# PLOT MEAN FITTED TRAJECTORY CURVE
pred <- data.frame(age = seq(min(dat$age), max(dat$age), length = 100))
pred$pred <- predict(model, pred, level = 0)
library(tidyverse)
ggplot(data = pred, aes(x = age, y = pred)) + geom_line()
# PLOT MEAN VELOCITY CURVE
library(SplinesUtils)
spl_population <- RegSplineAsPiecePoly(model, "ns(age, df = 7)")
plot(spl_population, deriv = 1)
# ESTIMATE MEAN AGE AT PEAK VELOCITY AND PEAK VELOCITY
(apv_mean <- solve(spl_population, b = 0, deriv = 2))
(pv_mean <- predict(spl_population, apv_mean, deriv = 1))
(apv_pv_mean <- as.data.frame(cbind(apv_mean, pv_mean)))
(apv_pv_mean <- apv_pv_mean %>% top_n(1, pv_mean))
# apv_mean = 14.09897 years
# pv_mean = 317.1267 grams