如何从线性混合自然样条模型获取个体峰值速度及对应年龄？

你遇到的异常结果（负值pv、不合理的apv）主要是因为错误地将随机二次项直接叠加到固定效应样条的分段多项式系数上——你的模型结构是固定效应为自然样条，随机效应为年龄的二次多项式，个体的完整拟合曲线应该是这两部分的和，而非修改样条本身的系数。下面是修正后的完整流程：

核心思路 #

每个个体的拟合函数为：
$$f_i(\text{age}) = \text{固定效应样条预测值} + \text{随机二次项预测值}$$
我们需要找到该函数的一阶导数（即速度）的最大值：

• 峰值速度年龄（apv）：使一阶导数取得最大值的年龄
• 峰值速度（pv）：该最大值对应的导数数值

具体实现代码 #

1. 拟合模型（保留你的原始代码） #

library(nlme)
library(splines)
library(tidyverse)
library(numDeriv) # 用于数值求导，或手动计算导数

# 拟合模型
nspline_model <- lme(y ~ ns(age, df = 3), data = dat, random = ~ poly(age, 2) | id)

# 提取固定效应和随机效应
fixef_coef <- fixef(nspline_model)
ranef_coef <- ranef(nspline_model)

2. 定义个体水平的拟合函数及其导数 #

我们手动构造函数，结合固定效应样条和个体随机二次项：

# 定义单个个体的拟合函数
individual_fit <- function(age, id) {
  # 固定效应部分：ns(age,3)的预测值，复用原模型的样条参数
  ns_design <- ns(age, df = 3, 
                  knots = attr(ns(dat$age, df=3), "knots"), 
                  Boundary.knots = attr(ns(dat$age, df=3), "Boundary.knots"))
  fixed_pred <- ns_design %*% fixef_coef
  
  # 随机效应部分：poly(age,2)的预测值，保持与模型一致的正交多项式
  poly_design <- poly(age, degree = 2, raw = FALSE)
  random_pred <- poly_design %*% t(ranef_coef[id, ])
  
  # 总预测值
  as.numeric(fixed_pred + random_pred)
}

# 定义一阶导数函数（即速度）
individual_velocity <- function(age, id) {
  # 用数值求导计算速度，稳定且易实现
  grad(func = individual_fit, x = age, id = id)
}

3. 为每个个体估算apv和pv #

我们用optimize()函数在数据的实际年龄范围内寻找速度的最大值，避免无意义的外推：

# 获取数据中的年龄范围，限制优化区间
age_range <- range(dat$age, na.rm = TRUE)

# 批量处理所有个体
apv_pv_results <- map_dfr(unique(dat$id), function(id) {
  # 寻找速度的最大值点
  opt_result <- optimize(f = individual_velocity, 
                         interval = age_range, 
                         id = id, 
                         maximum = TRUE)
  
  # 整理结果为数据框
  tibble(
    id = id,
    apv = opt_result$maximum,
    pv = opt_result$objective
  )
})

# 查看结果
head(apv_pv_results)

4. 关键注意事项 #

• 正交多项式一致性：模型中poly(age,2)默认使用正交多项式，构造随机效应预测时要保持raw=FALSE，避免系数不匹配。
• 年龄区间限制：必须将优化范围限定在数据的实际年龄区间内，外推到无数据区域会导致不合理的结果。
• 解析导数优化（可选）：如果追求更高效率，可以手动计算解析导数：
  • 用SplinesUtils提取固定效应样条的分段多项式导数
  • 正交二次多项式的一阶导数是一次多项式，可直接通过系数组合计算

为什么你的原始代码出错？ #

你尝试修改固定效应样条的分段多项式系数，但随机效应是独立的二次项，并非样条结构的一部分。正确的个体曲线是**样条（固定） + 二次多项式（随机）**的叠加，而非修改样条本身的结构。这种叠加后的曲线导数是样条导数加一次函数，需要整体寻找最大值，而非单独处理样条部分。

内容的提问来源于stack exchange，提问作者aelhak