естественный кубический сплайн регрессии с R

Question

я, кажется, есть проблема с splines::ns() в функции R.

Я создал простой фиктивную проблему

dat <- data.frame(t <- seq(0, 6, .01), 
        x <- rnorm(length(t), sd = 1), 
        y <- 5 + t - x^2 + rnorm(length(t), sd = .33)) 

lm(y ~ t + I(x^2), data = dat) 

library(splines) 
lm(y ~ t + ns(x, knots = c(0), Boundary.knots = c(-3, 3)), data = dat) 

В то время как первая модель работает отлично, второй не может чтобы правильно идентифицировать перехват. Что мне здесь не хватает?

Answer 1

Игнорирование нс() полностью вам не хватает двух вещей:

1) Комментарий выше, объясняющие, как определить dataframe:

t <- seq(0, 6, .01) 
x <- rnorm(length(t), sd = 1) 
y <- 5 + t - x^2 + rnorm(length(t), sd = .33) 
df <- data.frame(t, x, y)  
rm(t, x, y) 

2) То, как вы вызываете ваши модели:

lm(y ~ t + I(t^2), data=df) 
lm(y ~ splines::ns(t, knots = c(0), Boundary.knots = c(-3, 3)), data=df) 

Первая модель неправильно идентифицирует то, что, по вашему мнению, она делает.

Answer 2

Нет ничего плохого, потому что вы не подходите точно к одной и той же модели, и они даже не эквивалентны.

Чтобы объяснить различные результаты, которые вы видите, достаточно использовать более простой пример с одним ковариатором x. Мы генерируем данные из квадратичного многочлена: 5 + x + x^2, затем вставляем несколько моделей.

set.seed(0) 
x <- rnorm(500, mean = 1) ## `x` with non-zero mean 
y <- 5 + x + x * x + rnorm(500, sd = 0.5) 
library(splines) 

fit1 <- lm(y ~ x + I(x^2)) 
#(Intercept)   x  I(x^2) 
#  4.992  1.032  0.980 

fit2 <- lm(y ~ poly(x, degree = 2)) 
#(Intercept) poly(x, degree = 2)1 poly(x, degree = 2)2 
#  7.961    70.198    28.720 

fit3 <- lm(y ~ bs(x, degree = 2, df = 2)) 
#(Intercept) bs(x, degree = 2, df = 2)1 bs(x, degree = 2, df = 2)2 
#  6.583      -8.337      20.650 

fit4 <- lm(y ~ ns(x, df = 2)) 
#(Intercept) ns(x, df = 2)1 ns(x, df = 2)2 
#  5.523   10.737   21.265 

Первые 3 модели не то же самое, с точки зрения параметризации, но они эквивалентны: все они фитинг квадратичный полином с 3 степенью свободы. Для того, чтобы увидеть их эквивалентность, мы проверяем их подогнанные значения:

sum(abs(fit1$fitted - fit2$fitted)) 
# [1] 1.54543e-13 

sum(abs(fit1$fitted - fit3$fitted)) 
# [1] 2.691181e-13 

Чтобы увидеть разницу в параметризации, мы смотрим на матрице:

X1 <- model.matrix(~ x + I(x^2)) 
X2 <- model.matrix(~ poly(x, degree = 2)) 
X3 <- model.matrix(~ bs(x, degree = 2, df = 2)) 

par(mfrow = c(3,3), oma = rep.int(1,4), mar = c(4, 4, 0, 0)) 

plot(x, X1[, 1], cex = 0.2) 
plot(x, X1[, 2], cex = 0.2) 
plot(x, X1[, 3], cex = 0.2) 

plot(x, X2[, 1], cex = 0.2) 
plot(x, X2[, 2], cex = 0.2) 
plot(x, X2[, 3], cex = 0.2) 

plot(x, X3[, 1], cex = 0.2) 
plot(x, X3[, 2], cex = 0.2) 
plot(x, X3[, 3], cex = 0.2) 

Поскольку конструкция матрицы не являются то же (либо в форме, либо в масштабах), вы не получите тот же набор коэффициентов. В случае, если вы удивлены, давайте попробуем еще более простой пример:

x1 <- x - mean(x) 
test <- lm(y ~ x1 + I(x1^2)) 
#(Intercept)   x1  I(x1^2) 
#  7.003  2.991  0.980 

sum(abs(fit1$fitted - test$fitted)) 
# [1] 1.24345e-13 

Вот, я только что некоторые простые преобразования для x, то результат отличается (но все еще эквивалент).

4-я модель fit4, подходит для кубического многочлена с 3 степенями свободы, поэтому он не эквивалентен всем предыдущим моделям. Мы можем проверить установленные значения:

sum(abs(fit1$fitted - fit4$fitted)) 
# [1] 39.36563