Кубические сплайны с серией для манекенов: добавление списка манекенов к регрессии - PullRequest
0 голосов
/ 19 сентября 2019

У меня есть кривая с возмущением в середине (при t = 9), которое вызывает спад (t <9) и подъем (t> 9) в моих данных.Я хотел бы подогнать кубические сплайны и добавить фиктивные переменные в области возмущения, чтобы исключить эти значения при прогнозировании.

У меня есть две проблемы, которые могут быть одинаковыми:

  1. Я хотел бы автоматически добавить серию манекенов в регрессию.Обычно я использовал бы фактор-переменные, но пакет сплайнов не допускает факторы.Поэтому я явно создаю серию манекенов (dat5, dat6, dat7, ...), которые я хотел бы легко добавить к регрессии.Обратите внимание, что это игрушечный пример, и фактическое количество манекенов будет меняться в зависимости от модели.

  2. Я хотел бы предсказать основную кривую БЕЗ возмущения.Для этого я сначала вычисляю прогноз для полной модели и вычитаю значения для фиктивных коэффициентов.(см. ближе к концу кода).На данный момент это действительно громоздко, есть ли хороший способ автоматизировать это?

y <-c(185, 160, 145, 127, 117, 74, 76, 78, 101, 115, 
                  120, 70, 64, 65, 64, 63, 62, 60, 59, 58)
t <- seq(1,20,1)
z <- ifelse(t>=5 & t<13, t, NA)
dat = data.frame(y=y, t=t, z=z)

# dummies
library(dummies)
dat <- cbind(dat, dummy(dat$z, sep = ""))

# splines:
library(splines)
fit <- lm(y ~ ns(t, knots =c( 4, 10, 16)) + dat5 +dat6 + dat7 +dat8 
+ dat9 + dat10 + dat11, data = dat)
summary(fit)


# save coefficients:
d5<-fit[["coefficients"]][["dat5"]]
d6<-fit[["coefficients"]][["dat6"]]
d7<-fit[["coefficients"]][["dat7"]]
d8<-fit[["coefficients"]][["dat8"]]
d9<-fit[["coefficients"]][["dat9"]]
d10<-fit[["coefficients"]][["dat10"]]
d11<-fit[["coefficients"]][["dat11"]]

# predict:
dat$pred <- predict(fit, dat) - d5*dat$dat5 - d6*dat$dat6 - d7*dat$dat7 - d8*dat$dat8 - 
                                   d9*dat$dat9 - d10*dat$dat10 - d11*dat$dat11

library(ggplot2)
ggplot(dat, aes(x=t, y=y)) +
  geom_line() +
  geom_line(aes(x=t, y=pred), color='blue')

Так выглядит окончательный прогноз.

enter image description here

...