受限立方样条拐点求解失败求助（R语言RMS包）

报错原因 #

uniroot(f, interval=c(20,35))报错f(lower) * f(upper) > 0，核心问题有两个：

• 你传入的f是模型的预测因变量函数，而非找拐点需要的二阶导数函数，目标方向错误。
• 即使是找极值点（一阶导数为0），当前区间[20,35]内函数值符号未发生变化，uniroot无法找到有效根。

正确寻找拐点的步骤 #

拐点是曲线二阶导数为0的点，对应rcs(BMI)拟合曲线的曲率变化点，需按以下操作执行：

1. 先可视化拟合曲线，确认拐点存在性及合理区间 #

先画出BMI与BMD_ALL的拟合曲线及二阶导数曲线，观察曲率变化，确定uniroot的有效区间：

# 生成覆盖数据实际范围的BMI序列
bmi_seq <- seq(min(test$BMI, na.rm=TRUE), max(test$BMI, na.rm=TRUE), length.out=100)
# 固定其他协变量（示例为gender取参考水平、age取均值，可根据需求调整）
new_data <- data.frame(BMI=bmi_seq, gender=factor(levels(test$gender)[1]), age=mean(test$age, na.rm=TRUE))
# 预测二阶导数
preds_second <- predict(test_model, newdata=new_data, deriv=2)
# 绘制二阶导数曲线，查看是否穿越0线
plot(bmi_seq, preds_second, type="l", xlab="BMI", ylab="二阶导数")
abline(h=0, col="red", lty=2)

通过此图可观察到二阶导数在哪个区间内穿过0，以此作为uniroot的输入区间。

2. 定义二阶导数函数，用uniroot找根 #

方法一：利用rms内置的导数计算

直接用predict的deriv参数构造二阶导数函数：

# 固定其他协变量的二阶导数函数
f_second_deriv <- function(bmi) {
  predict(test_model, 
          newdata=data.frame(BMI=bmi, gender=factor(levels(test$gender)[1]), age=mean(test$age, na.rm=TRUE)), 
          deriv=2)
}

# 替换为你从可视化中得到的有效区间
uniroot_out <- uniroot(f = f_second_deriv, interval = c(25, 30))
uniroot_out

方法二：用Deriv包自动求导

若需要更灵活的求导方式，可使用Deriv包：

install.packages("Deriv")
library(Deriv)

# 定义固定协变量的预测函数
f_pred <- function(bmi) {
  predict(test_model, 
          newdata=data.frame(BMI=bmi, gender=factor(levels(test$gender)[1]), age=mean(test$age, na.rm=TRUE)))
}

# 求二阶导数
f_second_deriv <- Deriv(Deriv(f_pred))

# 传入有效区间找根
uniroot_out <- uniroot(f = f_second_deriv, interval = c(25, 30))
uniroot_out

3. 关键注意事项 #

• 必须固定其他协变量（gender、age）的取值，拐点是控制其他变量时BMI与BMD_ALL关系的曲率变化点，不同协变量取值可能对应不同拐点。
• 若可视化后发现二阶导数始终大于0或小于0，说明拟合曲线无拐点，此时uniroot必然报错，无需强行寻找。
• 确保interval区间内二阶导数两端符号相反（一正一负），否则需调整区间范围。

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