rm(list=ls())
par.ori <- par(no.readonly = TRUE)
##
## Regressão não paramétrica
##
## Exemplos do livro de Faraway
##

## Três conjuntos de dados para ilustrar o uso dos métodos,
## características, vantagens e desvantagens

## Uma função "complexa"
data(exa, package="faraway")
plot(y ~ x, exa, main="Example A")
lines(m ~ x, exa, lwd=2, col = "gray")

## Um dado "nulo" com pontos discrepantes para verificar sensibilidade a "outliers"
data(exb, package="faraway")
plot(y ~ x, exb, main="Example B")
lines(m ~ x, exb, lwd=2, col = "gray")

## Um conjunto "tradicional" na literatura que mostra dois grupos
plot(waiting ~ eruptions, faithful, main="Old Faithful")

## objetos para guardar predições dos modelos
exa.pred <- data.frame(x = seq(0, 1, length=201))
exb.pred <- data.frame(x = seq(0, 1, length=201))
faithful.pred <- data.frame(eruptions = seq(1,5.5,length=201))

##
## knn
##
## usando caret:knnreg() 
plot(y ~ x, data = exa, col=gray(0.75))
fit.knn5 <- caret::knnreg(y ~ x, data=exa, k=5)
exa.pred$knn5 <- predict(fit.knn5, newdata=exa.pred)
with(exa.pred, lines(x, knn5, col="red", lwd=2))

fit.knn10 <- caret::knnreg(y ~ x, data=exa, k=10)
exa.pred$knn10 <- predict(fit.knn10, newdata=exa.pred)
with(exa.pred, lines(x, knn10, col="blue", lwd=2))
legend("topright", legend=c("k = 5", "k = 10"), col=c("red","blue"), lwd=2, bty="n")

## usando FNN::knn.reg()
## Atenção: diferente de get.knnx() (abaixo), knn.reg() trata um "test"
## passado como VETOR simples (sem dim) como um único ponto de teste
## multidimensional, não como várias observações 1-D -- por isso "test"
## precisa ser um data.frame/matriz de uma coluna (exa.pred["x"]), nunca
## um vetor solto (exa.pred$x) como poderia parecer natural
plot(y ~ x, data = exa, col=gray(0.75))
fit.knn5.fnn <- with(exa, FNN::knn.reg(train = x, y = y, test = exa.pred["x"], k = 5))
exa.pred$knn5.fnn <- fit.knn5.fnn$pred
with(exa.pred, lines(x, knn5.fnn, col="red", lwd=2, lty=1))

fit.knn10.fnn <- with(exa, FNN::knn.reg(train = x, y = y, test = exa.pred["x"], k = 10))
exa.pred$knn10.fnn <- fit.knn10.fnn$pred
with(exa.pred, lines(x, knn10.fnn, col="blue", lwd=2, lty=1))
legend("topright", legend=c("k = 5", "k = 10"), col=c("red","blue"), lwd=2, bty="n")

## FNN::get.knnx(), que devolve os índices e as distâncias dos
## k vizinhos mais próximos de cada ponto de teste (knn.reg() não devolve
## essa informação); "data" e "query" precisam ter o mesmo número de
## colunas, por isso usamos apenas a coluna x de cada lado, não os
## data.frames inteiros
viz.knn5 <- with(exa, FNN::get.knnx(data = x, query = exa.pred$x, k = 5))
head(viz.knn5$nn.index)  ## índices (em exa) dos 5 vizinhos de cada linha de exa.pred
head(viz.knn5$nn.dist)   ## respectivas distâncias

##
## suavização por kernel (Naradaya-Watson)
##
## usando ksmooth() 
par(mfrow=c(1,3))
for(bw in c(0.1,0.5,2)){
    with(faithful,{
        plot(waiting ~ eruptions, col=gray(0.75))
        lines(ksmooth(eruptions,waiting, "normal", bw))
    })
    title(main = sprintf("bandwidth = %.1f", bw))
}
par(par.ori)

##
## usando sm::regression() 
## a função do pacote sm permite seleção do parâmetro de suavização por "cross-validation"
##
library(sm)
h.exa <- with(exa, h.select(x,y))
with(exa, sm.regression(x, y, h=h.exa))
title(sub = sprintf("h = %.3f (seleção por CV)", h.exa))
##with(exa, sm.regression(x, y, h=h.select(x,y), panel=TRUE))
h.exb <- with(exb, h.select(x,y))
with(exb, sm.regression(x, y, h=h.exb))
title(sub = sprintf("h = %.3f (seleção por CV)", h.exb))
##with(faithful, sm.regression(eruptions, waiting, h=h.select(eruptions,waiting)))

##
## splines suavizadores
##
with(exa,{
    plot(y ~ x, col=gray(0.75))
    lines(x,m)  # modelo "verdadeiro"
})
fit.ssp.a <- with(exa, smooth.spline(x,y))
exa.pred$ssp.a <- with(exa.pred, predict(fit.ssp.a, x)$y)
with(exa.pred, lines(x, ssp.a, lty=2))
legend("topright", legend=c("m(x) verdadeira", "smoothing spline"), lty=c(1,2), bty="n")

with(exb,{
    plot(y ~ x, col=gray(0.75))
    lines(x,m)  # modelo "verdadeiro"
})
fit.ssp.b <- with(exb, smooth.spline(x,y))
exb.pred$ssp.b <- with(exb.pred, predict(fit.ssp.b, x)$y)
with(exb.pred, lines(x, ssp.b, lty=2))
legend("topright", legend=c("m(x) verdadeira", "smoothing spline"), lty=c(1,2), bty="n")

plot(waiting ~ eruptions, col=gray(0.75), data = faithful)
fit.ssp.f <-   with(faithful, smooth.spline(eruptions, waiting))
faithful.pred$ssp.f <- with(faithful.pred, predict(fit.ssp.f, eruptions)$y)
with(faithful.pred, lines(eruptions, ssp.f, lty=2))

##
## regressão segmentada
##
## funções base para regressão segmentada (9 "knots")
rhs <- function(x,c) ifelse(x>c,x-c,0)
curve(rhs(x, 0.5), 0, 1)  ## função base para 1 nó
(knots <- 1:9/10)
dm  <- with(exa, outer(x, knots, rhs))
with(exa, matplot(x, dm, type="l", col=1, xlab="x", ylab=""))
lmod <- with(exa, lm(y ~ dm))
plot(y ~ x, exa, col=gray(0.75))
with(exa, lines(x, predict(lmod), lwd=2))

newknots <- c(0,0.5,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95)
dmn  <- with(exa, outer(x, newknots, rhs))
lmod <- with(exa, lm(y ~ dmn))
#plot(y ~ x, exa, col=gray(0.75))
with(exa, lines(x, predict(lmod), col="red", lwd=2, lty=2))
legend("topright", legend=c("9 nós uniformes", "10 nós escolhidos"),
       col=c("black","red"), lty=c(1,2), lwd=2, bty="n")

## agora usando pacote: spline::bs()
## splines (segmentada com polinômios de grau 1, contínua)
library(splines)
## funções base (12 nós) grau = 1
matplot(bs(seq(0, 1, length=1000), df=12, degree=1), type="l", ylab="", col=1)

fit.sp1 <- lm(y ~ bs(x, 12, degree = 1, Boundary.knots = c(0,1)), exa)
exa.pred$sp1 <- predict(fit.sp1, newdata = exa.pred)
plot(y ~ x, exa, col=gray(0.75))
lines(m ~ x, exa, col = "gray")       ## verdadeiro
with(exa.pred, lines(x, sp1, lty=2))  ## splines lineares
legend("topright", legend=c("m(x) verdadeira", "spline linear (bs, 12 nós)"),
       col=c("gray","black"), lty=c(1,2), bty="n")

## funções base (12 nós) grau = 3
matplot(bs(seq(0,1,length=1000), df=12), type="l", ylab="",col=1)

fit.sp2 <- lm(y ~ bs(x, 12, degree = 3, Boundary.knots = c(0,1)), exa)
exa.pred$sp2 <- predict(fit.sp2, newdata = exa.pred)
plot(y ~ x, exa, col=gray(0.75))
lines(m ~ x, exa, col = "gray")       ## verdadeiro
with(exa.pred, lines(x, sp2, lty=2))  ## splines cúbicos

## natural splines
fit.sp3 <- lm(y ~ ns(x, df = 12, Boundary.knots = c(0,1)), data = exa)
exa.pred$sp3 <- predict(fit.sp3, newdata = exa.pred)
##plot(y ~ x, exa, col=gray(0.75))
##lines(m ~ x, exa, col = "gray")       ## verdadeiro
with(exa.pred, lines(x, sp3, lty=3))   ## natural splines
## legenda cobre as três curvas do mesmo gráfico (bs cúbico + ns, acima)
legend("topright", legend=c("m(x) verdadeira", "B-spline cúbico (bs, 12 nós)", "natural spline (ns, 12 nós)"),
       col=c("gray","black","black"), lty=c(1,2,3), bty="n")

##
## polinômios locais (loess ou lowess)
##
with(exa,{
    plot(y ~ x, col=gray(0.75))
    lines(m ~ x, col = "gray")
    })
f1a <- with(exa, loess(y ~ x)) 
exa.pred$lo1 <- predict(f1a, newdata = exa.pred)
with(exa.pred, lines(x, lo1, lty=2))

f2a <- with(exa, loess(y ~ x, span=0.22))
exa.pred$lo2 <- predict(f2a, newdata = exa.pred)
with(exa.pred, lines(x, lo2, lty=3))
legend("topright", legend=c("m(x) verdadeira", "loess (span padrão)", "loess (span = 0.22)"),
       col=c("gray","black","black"), lty=c(1,2,3), bty="n")

with(exb, {
    plot(y ~ x, col=gray(0.75))
    lines(m ~ x, col = "gray")
})
f1b <- loess(y ~ x, data = exb)
exb.pred$lo <- predict(f1b, newdata = exb.pred)
with(exb.pred, lines(x, lo, lty=2))
legend("topright", legend=c("m(x) verdadeira", "loess (span padrão)"),
       col=c("gray","black"), lty=c(1,2), bty="n")

with(faithful, plot(waiting ~ eruptions, col=gray(0.75)))
fit.lo <- loess(waiting ~ eruptions, data = faithful)
faithful.pred$lo <- predict(fit.lo, newdata = faithful.pred)
with(faithful.pred, lines(eruptions, lo, lty=2))

##
## modelos aditivos generalizados (GAM), univariado
##
require(mgcv)
with(exa,{
    plot(y ~ x, col=gray(0.75))
    lines(m ~ x, col = "gray")
})
fit.gam.a <- gam(y ~ s(x), data = exa)
exa.pred$gam <- predict(fit.gam.a, newdata = exa.pred)
with(exa.pred, lines(x, gam, lty=2))
legend("topright", legend=c("m(x) verdadeira", "gam: s(x)"),
       col=c("gray","black"), lty=c(1,2), bty="n")

with(exb,{
    plot(y ~ x, col=gray(0.75))
    lines(m ~ x, col = "gray")
})
fit.gam.b <- gam(y ~ s(x), data = exb)
exb.pred$gam <- predict(fit.gam.b, newdata = exb.pred)
with(exb.pred, lines(x, gam, lty=2))
legend("topright", legend=c("m(x) verdadeira", "gam: s(x)"),
       col=c("gray","black"), lty=c(1,2), bty="n")

with(faithful, plot(waiting ~ eruptions, col=gray(0.75)))
fit.gam.f <- gam(waiting ~ s(eruptions), data = faithful)
faithful.pred$gam <- predict(fit.gam.f, newdata = faithful.pred)
with(faithful.pred, lines(eruptions, gam, lty=2))


##
## métodos "embutidos" em opções do ggplot e adicionando bandas de incerteza
##
require(ggplot2)
ggplot(exa, aes(x=x,y=y)) +
    geom_point(alpha=0.25) +
    geom_smooth(method="loess", span=0.22) +
    geom_line(aes(x=x,y=m),linetype=2)

ggplot(exa, aes(x=x,y=y)) +
    geom_point(alpha=0.25) +
    geom_smooth(method="gam", formula=y ~ s(x, k=20)) +
    geom_line(aes(x=x,y=m),linetype=2)

## Ajuste em duas dimensões
data(savings, package="faraway")
x <- with(savings, cbind(pop15, ddpi))
y <- with(savings, sr)
sm.regression(x,y,h=c(1,1),xlab="pop15",ylab="growth",zlab="savings rate")
sm.regression(x,y,h=c(5,5),xlab="pop15",ylab="growth",zlab="savings rate")

sm.regression(x,y,h=c(5,5),xlab="pop15",ylab="growth",zlab="savings rate",
              panel=TRUE, display="persp", theta=-35, ticktype="detailed")
sm.regression(x,y,h=c(5,5),xlab="pop15",ylab="growth",zlab="savings rate",
              panel=TRUE, display="rgl", theta=-35, ticktype="detailed")


## Modelos aditivos generalizados (combina suavizadores em s()), bivariado
library(mgcv)
data(savings, package="faraway")
amod <- gam(sr ~ s(pop15,ddpi), data=savings)
vis.gam(amod, col="gray", ticktype="detailed",theta=-35)
lomod <- loess(sr ~ pop15 + ddpi, data=savings)
xg <- seq(21,48,len=20)
yg <- seq(0,17,len=20)
zg <- expand.grid(pop15=xg,ddpi=yg)
persp(xg, yg, predict(lomod, zg), theta=-35, ticktype="detailed", xlab="pop15", ylab="growth", zlab="savings rate", col="gray")