Блог

Я понимаю, что rcs() (из rms пакета) использует усеченную степенную основу для представления естественных (ограниченных) кубических сплайнов. В качестве альтернативы я мог бы использовать ns() (из splines пакета), который использует базис B-сплайна.

Однако я заметил, что прогнозы обучения и тестирования могут сильно отличаться (особенно при x экстраполяции). Я пытаюсь понять различия между rcs() и ns() и могу ли я использовать функции взаимозаменяемо.

Поддельные нелинейные данные.

 library(tidyverse)
library(splines)
library(rms)

set.seed(100)

xx <- rnorm(1000)
yy <- 10   5*xx - 0.5*xx^2 - 2*xx^3   rnorm(1000, 0, 4)
df <- data.frame(x=xx, y=yy)
  

Подгоните одну модель к ns другой rcs с одинаковыми узлами.

 ns_mod <- lm(y ~ ns(x, knots=c(-2, 0, 2)), data=df)

ddist <- datadist(df)
options("datadist" = "ddist")

trunc_power_mod <- ols(y ~ rcs(x, knots=c(-2, 0, 2)), data=df)
  

Изучите их соответствия (MSE).

 mean(ns_mod$residuals^2)
mean(trunc_power_mod$residuals^2)

df$pred_ns <- ns_mod$fitted.values
df$pred_trunc_power <- trunc_power_mod$fitted.values

df_melt <- df %>% 
  gather(key="model", value="predictions", -x, -y)

ggplot(df_melt, aes(x=x, y=y))  
  geom_point(alpha=0.1)  
  geom_line(aes(x=x, y=predictions, group=model, linetype=model))
  

Создайте набор тестовых данных и постройте прогнозы между двумя моделями.

 newdata <- data.frame(x=seq(-10, 10, 0.1))

pred_ns_new <- predict(ns_mod, newdata=newdata)
pred_trunc_new <- predict(trunc_power_mod, newdata=newdata)

newdata$pred_ns_new <- pred_ns_new
newdata$pred_trunc_new <- pred_trunc_new

newdata_melted <- newdata %>% 
  gather(key="model", value="predictions", -x)

ggplot(newdata_melted, aes(x=x, y=predictions, group=model, linetype=model))  
  geom_line()
  

Существует довольно простое объяснение: knots это не аргумент rcs() . Он хочет, чтобы узлы были указаны с помощью параметра parms . Другая проблема заключается в том, что knots параметр to ns() не указывает «граничные узлы», которые по умолчанию range(x) . Итак, чтобы получить одинаковые прогнозы, вам нужно

 trunc_power_mod <- ols(y ~ rcs(x, parms=c(min(x), -2, 0, 2, max(x))), data=df)