Загрузим данные Boston, где y-crim (уровень преступности по городу), x-nox(концентрация оксида азота).

my.seed <- 1
data(Boston)
attach(Boston)

Рассмотрим наблюдения, при этом выясним, сколько наблюдений показателя crim превышает 30.

gp <- ggplot(data = Boston, aes(x = nox, y = crim))
gp <- gp + geom_point() + geom_abline(slope = 0, intercept = 30, col = 'red')
gp

У нас имеется группа наблюдений, которые превысили черту. Построим локальную регрессию. Зададим величину окон s=0.8, так как график практически полностью лег в интревальные границы. ДЛя начала была построена модель fit, кубический полином.

#полином 3ой степени
fit <- lm(crim ~ poly(nox, 3), data = Boston)
round(coef(summary(fit)), 2)
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)       3.61       0.32   11.24        0
## poly(nox, 3)1    81.37       7.23   11.25        0
## poly(nox, 3)2   -28.83       7.23   -3.99        0
## poly(nox, 3)3   -60.36       7.23   -8.34        0
# границы изменения переменной nox
noxlims <- range(nox)
# значения nox, для которых делаем прогноз 
nox.grid <- seq(from = noxlims[1], to = noxlims[2], length.out = 506)
# рассчитать прогнозы и их стандартные ошибки
preds <- predict(fit, newdata = list(nox = nox.grid), se = T)

# границы доверительного интервала для уровня пеступности
se.bands <- cbind(lower.bound = preds$fit - 2*preds$se.fit,
                  upper.bound = preds$fit + 2*preds$se.fit)
# смотрим результат
round(head(se.bands), 2)
##   lower.bound upper.bound
## 1        0.40        5.06
## 2        0.34        4.93
## 3        0.28        4.79
## 4        0.22        4.66
## 5        0.17        4.53
## 6        0.11        4.41
par(mfrow=c(1,1))
plot(nox, crim, xlim = noxlims, cex = 0.5, col = 'darkgrey')

title('Локальная регрессия')

# подгоняем модель c окном 0.8
fit1 <- loess(crim ~ nox, span = 0.8, data = Boston)


# рисум модели
lines(nox.grid, predict(fit1, data.frame(nox = nox.grid)),
      col = 'red', lwd = 2)

# доверительные интервалы прогноза
matlines(x = nox.grid, y = se.bands, lwd = 1, col = 'black', lty = 3)
# легенда
legend('topright', 
       c('s = 0.8'),
       col = c('red'), lty = 1, lwd = 2, cex = 0.8)

Построим модель с непрерывным откликом GAM на локальной регрессии.

gam.lo <- gam(crim ~  lo(nox, span = 0.7) , 
              data = Boston)
par(mfrow = c(1, 1))
plot(gam.lo, se = T, col = 'green')

summary(gam.lo)
## 
## Call: gam(formula = crim ~ lo(nox, span = 0.7), data = Boston)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8194 -1.8851 -0.2366 -0.1355 79.5854 
## 
## (Dispersion Parameter for gaussian family taken to be 53.8853)
## 
##     Null Deviance: 37363.22 on 505 degrees of freedom
## Residual Deviance: 27063.41 on 502.2413 degrees of freedom
## AIC: 3459.06 
## 
## Number of Local Scoring Iterations: 2 
## 
## Anova for Parametric Effects
##                         Df  Sum Sq Mean Sq F value    Pr(>F)    
## lo(nox, span = 0.7)   1.00  6621.4  6621.4  122.88 < 2.2e-16 ***
## Residuals           502.24 27063.4    53.9                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                     Npar Df Npar F     Pr(F)    
## (Intercept)                                     
## lo(nox, span = 0.7)     1.8 38.815 7.661e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1