rm(list=ls())
par.ori <- par(no.readonly = TRUE)
##
## Modelos segmentados
## Autor: PJ

## Conjunto de dados utilizado no "script"
xy01 <- read.csv("http://www.leg.ufpr.br/~paulojus/dados/xy01.csv")

## Inicializando objetos para armazenar modelos ajustados e predições destes modelos em um "grid" de valores
mods1 <- list()
nd1 <- data.frame(x = seq(0, 11, length=501))
## espaçamento da malha de pontos em x
(dx <- with(nd1, unique(diff(x)))[1])

## Visualização dos dados
plot(y ~ x, data = xy01, pch = 20)
##
## Modelos Globais (Polinomiais)
##
mods1$pol0 <- lm(y ~ 1, data = xy01)
mods1$pol0
nd1$ypol0 <- with(mods1, predict(pol0, newdata = nd1))
plot(y ~ x, data = xy01)
with(nd1, lines(ypol0 ~ x))

#mods1$pol1 <- lm(y ~ x, data = xy01)
mods1$pol1 <- lm(y ~ poly(x, 1), data = xy01)
mods1$pol1
nd1$ypol1 <- with(mods1, predict(pol1, newdata = nd1))
#plot(y ~ x, data = xy01)
with(nd1, lines(ypol1 ~ x))

#mods1$pol2 <- lm(y ~ x + I(x^2), data = xy01)
mods1$pol2 <- lm(y ~ poly(x, 2), data = xy01)
mods1$pol2
nd1$ypol2 <- with(mods1, predict(pol2, newdata = nd1))
#plot(y ~ x, data = xy01)
with(nd1, lines(ypol2 ~ x))

mods1$pol3 <- lm(y ~ poly(x, 3), data = xy01)
mods1$pol3
nd1$ypol3 <- with(mods1, predict(pol3, newdata = nd1))
#plot(y ~ x, data = xy01)
with(nd1, lines(ypol3 ~ x))

## ...

mods1$pol12 <- lm(y ~ poly(x, 12), data = xy01)
mods1$po12
nd1$ypol12 <- with(mods1, predict(pol12, newdata = nd1))
#plot(y ~ x, data = xy01)
with(nd1, lines(ypol12 ~ x))

## Ajustes e números de parâmetros dos modelos
sapply(mods1, logLik)
sapply(mods1, function(x) attr(logLik(x), "df"))

##
## Modelos segmentados
##

##
## 0. Constantes (grau 0)
##
xc <- c(1.7, 2.4, 4, 4.8, 5.7, 6.3, 8.5, 9.0, 9.7) 

plot(y ~ x, data = xy01)
abline(v = xc, lty = 2)

mods1$g0 <- lm(y ~ I(x > xc[1]) + I(x > xc[2]) + I(x > xc[3]) +
                   + I(x > xc[4]) + I(x > xc[5]) + I(x > xc[6]) +
                   + I(x > xc[7]) + I(x > xc[8]) + I(x > xc[9]),
               data = xy01)
mods1$g0
nd1$yg0 <- with(mods1, predict(g0, newdata = nd1))

plot(y ~ x, data = xy01, pch = 20)
with(nd1, lines(yg0 ~ x))
abline(v = xc, lty = 2)

rm(xc)

##
## 1. Lineares (grau 1)
##

## um conjunto menor de pontos de corte 
xc1 <- 1.8; xc2 <- 3.3; xc3 <- 5.3; xc4 <- 7.4

## 1.0 Sem restrição
plot(y ~ x, data = xy01)
abline(v = c(xc1, xc2, xc3, xc4), lty = 2)

mods1$g1.0 <- lm(y ~ x +
                     I(x > xc1) + I((x > xc1)*x) + 
                     I(x > xc2) + I((x > xc2)*x) + 
                     I(x > xc3) + I((x > xc3)*x) + 
                     I(x > xc4) + I((x > xc4)*x), 
               data = xy01)
mods1$g1.0
nd1$yg1.0 <- with(mods1, predict(g1.0, newdata = nd1))

plot(y ~ x, data = xy01, pch = 20)
with(subset(nd1, x <= xc1), lines(yg1.0 ~ x))
with(subset(nd1, x > xc1 & x <= xc2), lines(yg1.0 ~ x))
with(subset(nd1, x > xc2 & x <= xc3), lines(yg1.0 ~ x))
with(subset(nd1, x > xc3 & x <= xc4), lines(yg1.0 ~ x))
with(subset(nd1, x > xc4), lines(yg1.0 ~ x))

## 1.1 Com restrição de continuidade na função
plot(y ~ x, data = xy01)
abline(v = c(xc1, xc2, xc3, xc4), lty = 2)

mods1$g1.1 <- lm(y ~ x + I((x > xc1)*(x-xc1)) + 
                   I((x > xc2)*(x-xc2)) +
                   I((x > xc3)*(x-xc3)) +
                   I((x > xc4)*(x-xc4)),
               data = xy01)
mods1$g1.1
nd1$yg1.1 <- with(mods1, predict(g1.1, newdata = nd1))

plot(y ~ x, data = xy01)
with(nd1, lines(yg1.1 ~ x))

## Derivadas
nd1 <- transform(nd1, yg1.0.d1 = c(NA,diff(yg1.0)/dx))
nd1 <- transform(nd1, yg1.1.d1 = c(NA,diff(yg1.1)/dx))

## Ajustes e 1as derivadas
par(mfcol=c(2,2), mar=c(3,3,0.5, 0.5), mgp=c(1.7,0.7, 0))
plot(y ~ x, data = xy01)
with(subset(nd1, x <= xc1), lines(yg1.0 ~ x))
with(subset(nd1, x > xc1 & x <= xc2), lines(yg1.0 ~ x))
with(subset(nd1, x > xc2 & x <= xc3), lines(yg1.0 ~ x))
with(subset(nd1, x > xc3 & x <= xc4), lines(yg1.0 ~ x))
with(subset(nd1, x > xc4), lines(yg1.0 ~ x))
with(nd1, plot(yg1.0.d1 ~ x, cex=0.25, pch=19, ylim = c(-40, 40)))
##
plot(y ~ x, data = xy01, pch = 20)
with(nd1, lines(yg1.1 ~ x))
with(nd1, plot(yg1.1.d1 ~ x, cex=0.25, pch=19, ylim = c(-40, 40)))

par(par.ori)
## dev.off()

## ============================================================
## 1.2 A restrição de continuidade como hipótese linear testável
## ============================================================
## Em cada nó, o modelo SEM restrição (g1.0) tem dois parâmetros livres:
##   c_k (salto no intercepto)  e  d_k (mudança na inclinação)
## de forma que a reta à esquerda e a reta à direita de x = xc_k podem
## assumir valores DIFERENTES em xc_k (descontinuidade).
##
## Impor continuidade em x = xc_k exige que os dois trechos coincidam
## nesse ponto:
##   (b0 + b1*xc_k) == (b0 + c_k) + (b1 + d_k)*xc_k
##   0 == c_k + d_k*xc_k   =>   c_k == -d_k*xc_k
##
## Substituindo essa restrição nos termos de g1.0:
##   c_k*I(x>xc_k) + d_k*I(x>xc_k)*x
##     = -d_k*xc_k*I(x>xc_k) + d_k*I(x>xc_k)*x
##     =  d_k*I(x>xc_k)*(x - xc_k)
## que é exatamente a base truncada (x-xc_k)+ usada em g1.1.
##
## Ou seja: g1.1 é g1.0 com 4 restrições lineares impostas (uma por nó),
## reduzindo 10 para 6 parâmetros. Como o espaço de g1.1 está contido no
## espaço de g1.0 (mesma combinação linear das colunas), os modelos são
## ANINHADOS e a restrição pode ser testada com um teste F:
with(mods1, anova(g1.1, g1.0))

## conferindo "na mão" que a diferença de graus de liberdade é 4:
length(coef(mods1$g1.0)) - length(coef(mods1$g1.1))

## ============================================================
## 1.3 Conferindo com bs()
## ============================================================
require(splines)          ## usado aqui e também mais adiante no arquivo
mods1$g1.bs <- lm(y ~ bs(x, degree = 1, knots = c(xc1, xc2, xc3, xc4)), data = xy01)
logLik(mods1$g1.1)
logLik(mods1$g1.bs)               ## mesmo ajuste (log-verossimilhanças iguais)
max(abs(fitted(mods1$g1.1) - fitted(mods1$g1.bs)))  ## ~ 0: mesmos valores ajustados

## ============================================================
## 1.4 Função genérica: spline linear com K nós arbitrários
## ============================================================
## Generaliza g1.0/g1.1 para qualquer vetor de nós, com (restrict=TRUE)
## ou sem (restrict=FALSE) a restrição de continuidade.
## Os nós entram como constantes literais na fórmula (via I(...)), e não
## como colunas pré-computadas, para que predict(newdata=...) funcione
## normalmente (só precisa de x em newdata, como nos modelos g1.0/g1.1).
linspline.f <- function(knots, data, restrict = TRUE) {
    if (restrict) {
        termos <- sprintf("I(ifelse(x > %.10g, x - %.10g, 0))", knots, knots)
    }
    else {
        termos <- c(sprintf("I(as.numeric(x > %.10g))", knots),
                    sprintf("I(as.numeric(x > %.10g) * x)", knots))
    }
    form <- as.formula(paste("y ~ x +", paste(termos, collapse = " + ")))
    lm(form, data = data)
}

## conferindo que reproduz exatamente g1.0 e g1.1 originais
knots4 <- c(xc1, xc2, xc3, xc4)
g1.0b <- linspline.f(knots4, xy01, restrict = FALSE)
g1.1b <- linspline.f(knots4, xy01, restrict = TRUE)
logLik(mods1$g1.0); logLik(g1.0b)
logLik(mods1$g1.1); logLik(g1.1b)

## ============================================================
## 1.5 Nós desconhecidos: estimando a posição por otimização numérica
## ============================================================
## Generaliza a busca em grade usada em 01partes.R e 03partes.R
## (1 nó, por grade em Ks) para K nós simultâneos, minimizando a
## soma de quadrados dos resíduos via optim().
rss.knots <- function(knots, data) {
    if (any(knots <= min(data$x)) || any(knots >= max(data$x))) return(1e10)
    if (any(diff(sort(knots)) <= 0)) return(1e10)  ## evita nós fora de ordem
    fit <- linspline.f(sort(knots), data, restrict = TRUE)
    sum(residuals(fit)^2)
}

start <- as.numeric(quantile(xy01$x, probs = c(0.2, 0.4, 0.6, 0.8)))
opt <- optim(start, rss.knots, data = xy01, method = "Nelder-Mead")
(knots.est <- sort(opt$par))

g1.1.est <- linspline.f(knots.est, xy01, restrict = TRUE)
nd1$yg1.1.est <- predict(g1.1.est, newdata = nd1)

par(par.ori)
plot(y ~ x, data = xy01, pch = 20)
abline(v = knots.est, lty = 2, col = "red")
with(nd1, lines(yg1.1.est ~ x, col = "red"))
with(nd1, lines(yg1.1 ~ x, col = 4, lty = 2))
legend("topright", c("nós estimados", "nós fixos (1.8, 3.3, 5.3, 7.4)"),
       col = c(2, 4), lty = c(1, 2))

## comparando ajuste: nós "à mão" vs. nós estimados
AIC(mods1$g1.1, g1.1.est)
## Obs 3: a otimização é sensível ao ponto de partida (`start`) e pode
## convergir a um ótimo local pior que os nós escolhidos manualmente -
## é o problema clássico de "free-knot splines". Vale tentar outros
## pontos de partida e comparar.

rm(xc1, xc2, xc3, xc4, knots4)
##
## 2. Quadraticos (grau 2)
##

## Modelos com um ponto de corte
xc <- 5

## 2.0 Sem restrição
mods1$g2.0 <- lm(y ~ x + I(x^2) +
                     I((x > xc)) +
                     I((x > xc)*(x-xc)) + 
                   I((x > xc)*(x^2-xc^2)),
               data = xy01)
mods1$g2.0
nd1$yg2.0 <- with(mods1, predict(g2.0, newdata = nd1))

plot(y ~ x, data = xy01, pch = 20)
#with(nd1, lines(yg2.0 ~ x))
with(subset(nd1, x <= xc), lines(yg2.0 ~ x))
with(subset(nd1, x >  xc), lines(yg2.0 ~ x))

## 2.1 Uma restrição (continuidade da função)
mods1$g2.1 <- lm(y ~ x + I(x^2) +
                     I((x > xc)*(x-xc)) + 
                     I((x > xc)*(x^2-xc^2)),
               data = xy01)
mods1$g2.1
nd1$yg2.1 <- with(mods1, predict(g2.1, newdata = nd1))

plot(y ~ x, data = xy01, xlim = c(1,11))
with(subset(nd1, x <= xc), lines(yg2.1 ~ x))
with(subset(nd1, x > xc), lines(yg2.1 ~ x))

## 2.2 Duas restrições (continuidade da função e de sua 1a derivada)
mods1$g2.2 <- lm(y ~ x + I(x^2) +
                     I((x > xc)*(x-xc)^2),
                 data = xy01)
mods1$g2.2
nd1$yg2.2 <- with(mods1, predict(g2.2, newdata = nd1))

plot(y ~ x, data = xy01)
with(subset(nd1, x <= xc), lines(yg2.2 ~ x))
with(subset(nd1, x >= xc), lines(yg2.2 ~ x))

#rm(xc)

## 2.2.2 Quadratico (com as duas restrições) e com dois pontos de corte
xc1 <- 4; xc2 <- 6

mods1$g2.2.2 <- lm(y ~ x + I(x^2) +
                     I((x > xc1)*(x-xc1)^2) +
                     I((x > xc2)*(x-xc2)^2),
                 data = xy01)
mods1$g2.2.2
nd1$yg2.2.2 <- with(mods1, predict(g2.2.2, newdata = nd1))

plot(y ~ x, data = xy01)
with(subset(nd1, x <= xc1), lines(yg2.2.2 ~ x))
with(subset(nd1, x > xc1 & x <= xc2), lines(yg2.2.2 ~ x))
with(subset(nd1, x >= xc2), lines(yg2.2.2 ~ x))

#rm(xc1, xc2)

##
## Examinando derivadas dos modelos quadráticos ajustados
##
nd1 <- transform(nd1, yg2.0.d1 = c(NA,diff(yg2.0)/dx), yg2.0.d2 = c(NA, NA,diff(yg2.0, differ = 2)/dx))
nd1 <- transform(nd1, yg2.1.d1 = c(NA,diff(yg2.1)/dx), yg2.1.d2 = c(NA, NA,diff(yg2.1, differ = 2)/dx))
nd1 <- transform(nd1, yg2.2.d1 = c(NA,diff(yg2.2)/dx), yg2.2.d2 = c(NA, NA,diff(yg2.2, differ = 2)/dx))
nd1 <- transform(nd1, yg2.2.2.d1 = c(NA,diff(yg2.2.2)/dx), yg2.2.2.d2 = c(NA, NA,diff(yg2.2.2, differ = 2)/dx))

## Colocando tudo junto (funções quadraticas ajustadas e suas derivadas)
par(mfcol=c(3,4), mar=c(3,3,0.5, 0.5), mgp=c(1.7,0.7, 0))
plot(y ~ x, data = xy01, pch = 20, xlim = c(0,11))
with(subset(nd1, x <= xc), lines(yg2.0 ~ x))
with(subset(nd1, x >  xc), lines(yg2.0 ~ x))
plot(yg2.0.d1 ~ x, type = "n",  data = nd1, ylim = c(-30,50)); abline(h=0, lty=3)
lines(yg2.0.d1 ~ x,  data = subset(nd1, x <= xc))
lines(yg2.0.d1 ~ x,  data = subset(nd1, x > xc+dx))
plot(yg2.0.d2 ~ x, data = nd1, ylim = c(-0.5, 0.5), pch=19, cex=0.25)
##
plot(y ~ x, data = xy01, pch = 20, xlim = c(0,11))
with(subset(nd1, x <= xc), lines(yg2.1 ~ x))
with(subset(nd1, x >= xc), lines(yg2.1 ~ x))
plot(yg2.1.d1 ~ x, type = "n",  data = nd1, ylim = c(-30,50)); abline(h=0, lty=3)
lines(yg2.1.d1 ~ x,  data = subset(nd1, x <= xc))
lines(yg2.1.d1 ~ x,  data = subset(nd1, x > xc+dx))
plot(yg2.1.d2 ~ x, data = nd1, ylim = c(-0.5, 0.5), pch=19, cex=0.25)
##
plot(y ~ x, data = xy01, pch = 20, xlim = c(0,11))
with(subset(nd1, x <= xc), lines(yg2.2 ~ x))
with(subset(nd1, x >= xc), lines(yg2.2 ~ x))
plot(yg2.2.d1 ~ x, type = "l",  data = nd1, ylim = c(-30,50)); abline(h=0, lty=3)
plot(yg2.2.d2 ~ x, data = nd1, ylim = c(-0.5, 0.5), pch=19, cex=0.25)
##
plot(y ~ x, data = xy01, pch = 20, xlim = c(0,11))
with(subset(nd1, x <= xc1), lines(yg2.2.2 ~ x))
with(subset(nd1, x > xc1 & x <= xc2), lines(yg2.2.2 ~ x))
with(subset(nd1, x >= xc2), lines(yg2.2.2 ~ x))
plot(yg2.2.2.d1 ~ x, type = "l",  data = nd1, ylim = c(-30,50)); abline(h=0, lty=3)
plot(yg2.2.2.d2 ~ x, data = nd1, ylim = c(-0.5, 0.5), pch=19, cex=0.25)


par(par.ori)
##
## 3. Cúbicos
##
## um ponto de corte
xc <- 5

## 3.0 Sem restrições
mods1$g3.0 <- lm(y ~ x + I(x^2) + I(x^3) +
                     I(x > xc) + 
                     I((x > xc)*x) + 
                     I((x > xc)*(x^2)) +
                     I((x > xc)*(x^3)),
                 data = xy01)
mods1$g3.0

nd1$yg3.0 <- with(mods1, predict(g3.0, newdata = nd1))

plot(y ~ x, pch = 20, data = xy01)
with(subset(nd1, x <= xc), lines(yg3.0 ~ x))
with(subset(nd1, x >= xc), lines(yg3.0 ~ x))

## 3.1 Uma restrição de continuidade na função
mods1$g3.1 <- lm(y ~ x + I(x^2) + I(x^3) +
                     I((x > xc)*(x - xc)) + 
                     I((x > xc)*(x^2 - xc^2)) + 
                     I((x > xc)*(x^3 - xc^3)),
                 data = xy01)
mods1$g3.1
nd1$yg3.1 <- with(mods1, predict(g3.1, newdata = nd1))

plot(y ~ x, pch = 20, data = xy01)
with(nd1, lines(yg3.1 ~ x))

## 3.2 Duas restrições: continuidade na função e na 1a derivada
mods1$g3.2 <- lm(y ~ x + I(x^2) + I(x^3) +
                     I((x > xc)*((x - xc)^2)) + 
                     I((x > xc)*(x^3 - 3*(xc^2)*x + 2*xc^3)),
                 data = xy01)
mods1$g3.2
nd1$yg3.2 <- with(mods1, predict(g3.2, newdata = nd1))

plot(y ~ x, data = xy01)
with(nd1, lines(yg3.2 ~ x))


## três restrições (continuidade da função e 1a e 2a derivada). Um ponto de corte
mods1$g3.3 <- lm(y ~ x + I(x^2) + I(x^3) +
                     I((x > xc)*(x-xc)^3),
                 data = xy01)
mods1$g3.3

nd1$yg3.3 <- with(mods1, predict(g3.3, newdata = nd1))

plot(y ~ x, data = xy01)
with(nd1, lines(yg3.3 ~ x))
#with(subset(nd1, x <= xc), lines(yg3.3 ~ x))
#with(subset(nd1, x >= xc), lines(yg3.3 ~ x))

## três restrições/nó e dois pontos de corte
xc1 <- 4; xc2 <- 7

mods1$g3.3.2 <- lm(y ~ x + I(x^2) + I(x^3) +
                     I((x > xc1)*(x-xc1)^3) +
                     I((x > xc2)*(x-xc2)^3),
                 data = xy01)
mods1$g3.3.2

nd1$yg3.3.2 <- with(mods1, predict(g3.3.2, newdata = nd1))

plot(y ~ x, data = xy01)
with(nd1, lines(yg3.3.2 ~ x))

## Colocando ajustes juntos
plot(y ~ x, data = xy01)
with(subset(nd1, x <= xc), lines(yg3.0 ~ x))
with(subset(nd1, x > xc), lines(yg3.0 ~ x))
#with(nd1, lines(yg3.0 ~ x))
with(nd1, lines(yg3.1 ~ x, col = 2))
with(nd1, lines(yg3.2 ~ x, col = 3))
with(nd1, lines(yg3.3 ~ x, col = 4))
with(nd1, lines(yg3.3.2 ~ x, col = 5))


## Derivadas
nd1 <- transform(nd1, yg3.0.d1 = c(NA,diff(yg3.0)/dx), yg3.0.d2 = c(NA, NA,diff(yg3.0, difference = 2)/dx))
nd1 <- transform(nd1, yg3.1.d1 = c(NA,diff(yg3.1)/dx), yg3.1.d2 = c(NA, NA,diff(yg3.1, difference = 2)/dx))
nd1 <- transform(nd1, yg3.2.d1 = c(NA,diff(yg3.2)/dx), yg3.2.d2 = c(NA, NA,diff(yg3.2, difference = 2)/dx))
nd1 <- transform(nd1, yg3.3.d1 = c(NA,diff(yg3.3)/dx), yg3.3.d2 = c(NA, NA,diff(yg3.3, difference = 2)/dx))
nd1 <- transform(nd1, yg3.3.2.d1 = c(NA,diff(yg3.3.2)/dx), yg3.3.2.d2 = c(NA, NA,diff(yg3.3.2, difference = 2)/dx))

## Colocando tudo junto (funções cúbicas ajustadas e suas derivadas)
par(mfcol=c(3,5), mar=c(3,3,0.5, 0.5), mgp=c(1.7,0.7, 0))
plot(y ~ x, data = xy01, xlim = c(0,11))
with(subset(nd1, x <= xc), lines(yg3.0 ~ x))
with(subset(nd1, x > xc), lines(yg3.0 ~ x))
with(nd1, plot(yg3.0.d1 ~ x, ylim = c(-40, 20), pch=19, cex=0.25))
with(nd1, plot(yg3.0.d2 ~ x, ylim = c(-2, 2), pch=19, cex=0.25))
##
plot(y ~ x, pch = 20, data = xy01, xlim = c(0,11))
with(nd1, lines(yg3.1 ~ x))
with(nd1, plot(yg3.1.d1 ~ x, pch=19, cex=0.25))
with(nd1, plot(yg3.1.d2 ~ x, ylim = c(-1, 1), pch=19, cex=0.25))
##
plot(y ~ x, pch = 20, data = xy01, xlim = c(0,11))
with(nd1, lines(yg3.2 ~ x))
with(nd1, plot(yg3.2.d1 ~ x, pch=19, cex=0.25))
with(nd1, plot(yg3.2.d2 ~ x, pch=19, cex=0.25))
##
plot(y ~ x, pch = 20, data = xy01, xlim = c(0,11))
with(nd1, lines(yg3.3 ~ x))
with(nd1, plot(yg3.3.d1 ~ x, pch=19, cex=0.25))
with(nd1, plot(yg3.3.d2 ~ x, pch=19, cex=0.25))
##
plot(y ~ x, pch = 20, data = xy01, xlim = c(0,11))
with(nd1, lines(yg3.3.2 ~ x))
with(nd1, plot(yg3.3.2.d1 ~ x, pch=19, cex=0.25))
with(nd1, plot(yg3.3.2.d2 ~ x, pch=19, cex=0.25))

par(par.ori)
##
## Todos são "splines" !
##
require(splines)

#par(mfcol=c(2,2), mar=c(3,3,0.5, 0.5), mgp=c(1.7,0.7, 0))
## = polinomio cubico global
mods1$sp3 <- lm(y ~ bs(x, degree = 3), data = xy01)
mods1$sp3
nd1$ysp3 <- with(mods1, predict(sp3, newdata = nd1))
plot(y ~ x, pch = 20, data = xy01)
with(nd1, lines(ysp3 ~ x))
with(nd1, lines(ypol3 ~ x), col = 4)

## linear por partes (= spline grau 1)
mods1$sp1 <- lm(y ~ bs(x, knots = c(1.8, 3.3, 5.3, 7.4), degree = 1), data = xy01)
mods1$sp1
nd1$ysp1 <- with(mods1, predict(sp1, newdata = nd1))
plot(y ~ x, pch = 20, data = xy01)
with(nd1, lines(ysp1 ~ x))

## quadratico com 2 pontos de corte
mods1$sp2 <- lm(y ~ bs(x, knots = c(4, 7), degree = 2), data = xy01)
mods1$sp2
nd1$ysp2 <- with(mods1, predict(sp2, newdata = nd1))
plot(y ~ x, pch = 20, data = xy01)
with(nd1, lines(ysp2 ~ x), col = 4)
##
xc1 <- 4; xc2 <- 7
with(subset(nd1, x < xc1), lines(yg2.2.2 ~ x), col = 4)
with(subset(nd1, x > xc1 & x <= xc2), lines(yg2.2.2 ~ x), col = 4)
with(subset(nd1, x >= xc2), lines(yg2.2.2 ~ x), col = 4)

## cubico com 2 pontos de corte
mods1$sp3.2 <- lm(y ~ bs(x, knots = c(4, 7), degree = 3), data = xy01)
mods1$sp3.2
nd1$ysp3.2 <- with(mods1, predict(sp3.2, newdata = nd1))
plot(y ~ x, pch = 20, data = xy01)
with(nd1, lines(ysp3.2 ~ x))
##
xc1 <- 4; xc2 <- 7
with(subset(nd1, x < xc1), lines(yg3.2 ~ x), col = 4)
with(subset(nd1, x > xc1 & x <= xc2), lines(yg3.2 ~ x), col = 4)
with(subset(nd1, x >= xc2), lines(yg3.2 ~ x), col = 4)

rm(xc1, xc2)

##
## Splines e intervalos de predição
##
require(splines)
xy01 <- read.csv("xy01.csv")
modelos <- list()
nd <- data.frame(x = seq(0, 11, length=501))

## Polinômio global
#modelos$pol12 <- lm(y ~ poly(x, 12), data = xy01)
modelos$pol12 <- lm(y ~ bs(x, degree = 12), data = xy01)
nd$ypol12 <- with(modelos, predict(pol12, newdata = nd, interval = "prediction"))
plot(y ~ x, data = xy01, xlim = c(0,12), ylim = c(-5,45))
with(nd, matlines(x, ypol12, col = 1, lty=c(1,2,2)))

## Spline linear
#xc1 <- 1.8; xc2 <- 3.3; xc3 <- 5.3; xc4 <- 7.4
#modelos$g1 <- lm(y ~ x + I((x > xc1)*(x-xc1)) + 
#                       I((x > xc2)*(x-xc2)) +
#                       I((x > xc3)*(x-xc3)) +
#                       I((x > xc4)*(x-xc4)),
#                   data = xy01)
modelos$g1 <- lm(y ~ bs(x, degree = 1, knots = c(1.8, 3.3, 5.3, 7.4)), data = xy01)
nd$yg1 <- with(modelos, predict(g1, newdata = nd, interval = "prediction"))
plot(y ~ x, data = xy01, xlim = c(0,12), ylim = c(-5,45))
with(nd, matlines(x, yg1, col = 1, lty=c(1,2,2)))

## Spline quadrático
##     2 nós
#xc1 <- 4; xc2 <- 6
#modelos$g2k2 <- lm(y ~ x + I(x^2) +
#                     I((x > xc1)*(x-xc1)^2) +
#                     I((x > xc2)*(x-xc2)^2),
#                 data = xy01)
modelos$g2k2 <- lm(y ~ bs(x, degree = 2, knots = c(4, 6)), data = xy01)
nd$yg2k2 <- with(modelos, predict(g2k2, newdata = nd, interval = "prediction"))
#xc1 <- 2.5; xc2 <- 5; xc3 <- 7.5
#modelos$g2k3 <- lm(y ~ x + I(x^2) +
#                     I((x > xc1)*(x-xc1)^2) +
#                     I((x > xc2)*(x-xc2)^2) +
#                     I((x > xc3)*(x-xc3)^2),
#                 data = xy01)
modelos$g2k3 <- lm(y ~ bs(x, degree = 2, knots = c(2.5, 5, 7.5)), data = xy01)
nd$yg2k3 <- with(modelos, predict(g2k3, newdata = nd, interval = "prediction"))
##     3 nós
plot(y ~ x, data = xy01, xlim = c(0,12), ylim = c(-5,45))
with(nd, matlines(x, yg2k2, col = 1, lty=c(1,2,2)))
with(nd, matlines(x, yg2k3, col = 4, lty=c(1,2,2)))
legend("topright", c("2 nós", "3 nós"), col = c(1,4), lty=1)

## Spline cúbico (base spline)
##     2 nós
#xc1 <- 4; xc2 <- 7
#modelos$g3k2 <- lm(y ~ x + I(x^2) + I(x^3) +
#                     I((x > xc1)*(x-xc1)^3) +
#                     I((x > xc2)*(x-xc2)^3),
#                 data = xy01)
modelos$g3k2 <- lm(y ~ bs(x, knots = c(4, 7)), data = xy01)
nd$yg3k2 <- with(modelos, predict(g3k2, newdata = nd, interval = "prediction"))
##     3 nós
#xc1 <- 2.5; xc2 <- 5; xc3 <- 7.5
#modelos$g3k3 <- lm(y ~ x + I(x^2) + I(x^3) +
#                     I((x > xc1)*(x-xc1)^3) +
#                     I((x > xc2)*(x-xc2)^3) + 
#                     I((x > xc3)*(x-xc3)^3),
#                 data = xy01)
modelos$g3k3 <- lm(y ~ bs(x, knots = c(2.5, 5, 7.5)), data = xy01)
nd$yg3k3 <- with(modelos, predict(g3k3, newdata = nd, interval = "prediction"))
plot(y ~ x, data = xy01, xlim = c(0,12), ylim = c(-5,45))
with(nd, matlines(x, yg3k2, col = 1, lty=c(1,2,2)))
with(nd, matlines(x, yg3k3, col = 4, lty=c(1,2,2)))
legend("topright", c("2 nós", "3 nós"), col = c(1,4), lty=1)

## Spline natural cúbico - ns()
modelos$g3k3n <- lm(y ~ ns(x, knots = c(2.5, 5, 7.5)), data = xy01)
nd$yg3k3n <- with(modelos, predict(g3k3n, newdata = nd, interval = "prediction"))
plot(y ~ x, data = xy01, xlim = c(0,12), ylim = c(-5,45))
with(nd, matlines(x, yg3k3, col = 1, lty=c(1,2,2)))
with(nd, matlines(x, yg3k3n, col = "darkorange", lty=c(1,2,2)))
abline(v = c(2.5, 5, 7.5), lty = 3, col = "darkgreen")
legend("topright", c("bs 3 nós", "ns 3 nós"), col = c("black", "darkorange"), lty=1)


## Funções "base" 
with(nd, matplot(x, cbind(I((x > 1.8)*(x-1.8)),
                          I((x > 3.3)*(x-3.3)),
                          I((x > 5.3)*(x-5.3)),
                          I((x > 7.4)*(x-7.4))),
                 type = "l", lty = 2))
with(nd, matplot(x, cbind(I((x > 4)*(x-4)^2),
                          I((x > 6)*(x-6)^2)),
                 type = "l", lty = 2))
with(nd, matplot(x, cbind(I((x > 2.5)*(x-2.5)^2),
                          I((x > 5)*(x-5)^2),
                          I((x > 7.5)*(x-7.5)^2)),
                 type = "l", lty = 2))
with(nd, matplot(x, cbind(I((x > 4)*(x-4)^3),
                          I((x > 7)*(x-7)^3)),
                 type = "l", lty = 2))
with(nd, matplot(x, cbind(I((x > 2.5)*(x-2.5)^3),
                          I((x > 5)*(x-5)^3),
                          I((x > 7.5)*(x-7.5)^3)),
                 type = "l", lty = 2))


## Funções "base" "equivalentes" em bs() e ns()
with(nd, matplot(x, bs(x, degree = 12), type = "l"))
with(nd, matplot(x, bs(x, degree = 1, knots = c(1.8, 3.3, 5.3, 7.4)), type = "l"))
with(nd, matplot(x, bs(x, degree = 2, knots = c(4, 6)), type = "l"))
with(nd, matplot(x, bs(x, degree = 2, knots = c(2.5, 5, 7.5)), type = "l"))
with(nd, matplot(x, bs(x, knots = c(4, 7)), type = "l"))
with(nd, matplot(x, bs(x, knots = c(2.5, 5, 7.5)), type = "l"))
with(nd, matplot(x, ns(x, knots = c(2.5, 5, 7.5)), type = "l"))