B Путаница со сплайнами

## simulate some data - from mgcv::magic
set.seed(1)
n <- 400
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
y <- f + rnorm(n, 0, sd = 2)

## load the splines package - comes with R
require(splines)

Вы используете функцию bs() в формуле для lm, когда вам нужны оценки OLS. bs предоставляет базовые функции, заданные узлами, степенью многочлена и т. д.

mod <- lm(y ~ bs(x, knots = seq(0.1, 0.9, by = 0.1)))

Вы можете относиться к этому так же, как к линейной модели.

> anova(mod)
Analysis of Variance Table

Response: y
                                        Df Sum Sq Mean Sq F value    Pr(>F)    
bs(x, knots = seq(0.1, 0.9, by = 0.1))  12 2997.5 249.792  65.477 < 2.2e-16 ***
Residuals                              387 1476.4   3.815                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Несколько советов по размещению узлов. bs имеет аргумент Boundary.knots, по умолчанию Boundary.knots = range(x), поэтому, когда я указал выше аргумент knots, я не включил граничные узлы.

Прочтите ?bs для получения дополнительной информации.

Создание графика подобранного сплайна

В комментариях я обсуждаю, как нарисовать подогнанный сплайн. Один из вариантов — упорядочить данные в терминах ковариаты. Это хорошо работает для одной ковариаты, но не обязательно для двух или более ковариат. Еще одна проблема заключается в том, что вы можете оценить подобранный сплайн только при наблюдаемых значениях x — это нормально, если вы произвели плотную выборку ковариаты, но в противном случае сплайн может выглядеть странно с длинными линейными участками.

Более общее решение состоит в том, чтобы использовать predict для создания прогнозов модели для новых значений ковариаты или ковариат. В приведенном ниже коде я показываю, как это сделать для приведенной выше модели, прогнозируя 100 равномерно распределенных значений в диапазоне x.

pdat <- data.frame(x = seq(min(x), max(x), length = 100))
## predict for new `x`
pdat <- transform(pdat, yhat = predict(mod, newdata = pdat))

## now plot
ylim <- range(pdat$y, y) ## not needed, but may be if plotting CIs too
plot(y ~ x)
lines(yhat ~ x, data = pdat, lwd = 2, col = "red")

Это производит

Основываясь на примере в ответе, более простым способом построения подогнанного сплайна было бы использование пакета effects.

## simulate some data - from mgcv::magic
set.seed(1)
n <- 400
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
y <- f + rnorm(n, 0, sd = 2)

## load the splines package - comes with R
require(splines)
require(car)
require(effects)

## estimate model
mod <- lm(y ~ bs(x, knots = seq(0.1, 0.9, by = 0.1)))

Затем вы можете использовать Anova из car:

> Anova(mod)
Anova Table (Type II tests)

Response: y
                                       Sum Sq  Df F value    Pr(>F)    
bs(x, knots = seq(0.1, 0.9, by = 0.1)) 2997.5  12  65.477 < 2.2e-16 ***
Residuals                              1476.4 387                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

И вы можете легко построить подогнанный сплайн, используя пакет effects.

plot(allEffects(mod))

Что выведет это:

Смотрите также:

• Представление предельных эффектов логита с фиксированными эффектами