Trabalho 1 - Extensões de Modelos de Regressão - CE092 - 2017/2

Professor Paulo Justiniano Ribeiro Jr.

Bruna Wundervald

2017-09-04


Introdução

library(lattice)
library(latticeExtra)


xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       ylab = "Y - Renda",
       xlab = "X - Tempo de estudo",
       pch = 16)
Figura 1: Gráfico de dispersão dos dados: estudo da relação entre X e Y

Figura 1: Gráfico de dispersão dos dados: estudo da relação entre X e Y

Regressão Linear Simples

Onde \(\beta_0\) pode ser interpretado como o intercepto da equação e \(\beta_1\) como o quanto muda \(y\) quando o valor de \(x\) é alterado em uma unidade. O \(\epsilon_i\) representa os valores dos resíduos, isto é, tudo aquilo que não é explicado pela parte determinística da equação. Temos que \(\epsilon_i \sim N(0, \sigma^2)\).

aic <- c()
m0 <- lm(y ~ x, data = da)
summary(m0)

Call:
lm(formula = y ~ x, data = da)

Residuals:
     Min       1Q   Median       3Q      Max 
-1083.13  -277.41   -34.15   241.46  1238.17 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1931.105    462.739   4.173 0.000125 ***
x             47.162      8.616   5.474 1.58e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 487.1 on 48 degrees of freedom
Multiple R-squared:  0.3843,    Adjusted R-squared:  0.3715 
F-statistic: 29.96 on 1 and 48 DF,  p-value: 1.579e-06
(aic[1] <- AIC(m0))
[1] 764.7086
xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       ylab = "Y - Renda",
       xlab = "X - Tempo de estudo",
       pch = 16,
       main = paste0("Regressão linear simples - AIC = ", 
                     round(aic[1], 3)),
       panel = function(x, y, ...) {
    panel.xyplot(x, y, ...)
    panel.abline(lm(y~x), col = 2, lwd = 3)
  })
Figura 2: Reta de regressão linear simples sobreposta aos dados

Figura 2: Reta de regressão linear simples sobreposta aos dados

Os resultados acima mostram que teríamos uma reta na qual: \[ y_i = 1931.105 + 47.162 x_i + \epsilon_i,\] Cuja demonstração está no gráfico anterior.

Regressão Por Partes

\[ \begin{align} y_i = \beta_{0} + \beta_{1} x_i , & \quad x \leq c \\ y_i = \beta_{0} + \beta_{1} x_i + \beta_{2} (x_i - c), & \quad x > c \end{align} \]

Onde \(c\) é o valor do “ponto de corte”, ou seja, é o ponto aonde a equação de estimação é alterada.

#-------- Encontrando pontos de corte
g <- c()

c.est <- for (i in 1:20) {
  da$v <- with(da, ifelse(x <= (40+i), 0, x - (40+i)))
  m1 <- lm(y ~ x + v, data = da)
  g[i] <- AIC(m1)
}
which.min(g) # Corresponde a 40 + 10 = 50
[1] 10
da$v <- with(da, ifelse(x <= 53, 0, x - 53))
m1 <- lm(y ~ x + v, data = da)
(aic[2] <- AIC(m1))
[1] 763.7651
da$p <- with(da, ifelse(x <= 50, 0, x - 50))
m2 <- lm(y ~ x + p, data = da)
(aic[3] <- AIC(m2))
[1] 762.4402
#--------------

xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       ylab = "Y - Renda",
       xlab = "X - Tempo de estudo",
       pch = 16,
       main = paste0("Regressão Segmentada - AIC = ", 
                     round(aic[2], 3), " - \nPonto de corte em 53"),
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         le <- lm(y ~ x, da,
                  subset = x <= 53)
         llines(da$x[da$x <= 53],
                predict(le), col = 2, lwd = 3)
         ld <- lm(y ~ x + v, da,
                  subset = x > 53)
         llines(da$x[da$x > 53],
          predict(ld), col = 3, lwd = 3)
})
Figura 3: Modelo de regressão segmentada com ponto de corte em 53

Figura 3: Modelo de regressão segmentada com ponto de corte em 53

xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       ylab = "Y - Renda",
       xlab = "X - Tempo de estudo",
       pch = 16,
       main = paste0("Regressão Segmentada - AIC = ", 
                     round(aic[3], 3), " - \n Ponto de corte em 50"),
       panel = function(x, y, ...){
         panel.xyplot(x, y, ...)
         le <- lm(y ~ x, da,
                  subset = x <= 50)
         llines(da$x[da$x <= 50],
                predict(le), col = 2, lwd = 3)
         ld <- lm(y ~ x + p, da,
                  subset = x > 50)
         llines(da$x[da$x > 50],
          predict(ld), col = 3, lwd = 3)
})
Figura 4: Modelo de regressão segmentada com ponto de corte em 50

Figura 4: Modelo de regressão segmentada com ponto de corte em 50

Suavização por Splines

Isto porque o primeiro termo, de perda, estimula \(f\) a se ajustar aos dados, enquanto o termo de penalidade justamente penaliza a variabilidade em \(f\).

p1 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise",
             ylab = "Y - Renda",
             xlab = "",
             pch = 16,
             main = "15 Graus de Liberdade")+
  layer(panel.lines(smooth.spline(da$x, da$y, df = 15),
                    lwd = 2,
                    col = 2))

p2 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise", 
             xlab = "X - Tempo de estudo",
             ylab = "",
             pch = 16,
             main = "10 Graus de Liberdade")+
  layer(panel.lines(smooth.spline(da$x, da$y, df = 10),
                    lwd = 2,
                    col = 3))

p3 <- xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       pch = 16,
       ylab = "",
       xlab = "",
       main = "5 Graus de Liberdade")+
  layer(panel.lines(smooth.spline(da$x, da$y, df = 5),
                    lwd = 2,
                    col = 2))

print(p1, position = c(0, 0, 1/3, 1), more = TRUE)
print(p2, position = c(1/3, 0, 2/3, 1), more = TRUE)
print(p3, position = c(2/3, 0, 1, 1))
Figura 5: Suavização por Splines com variação dos graus de liberdade

Figura 5: Suavização por Splines com variação dos graus de liberdade

Suavização por Kernel

\[ \hat{f}_{\lambda} (x) = \frac{1}{n \lambda} \sum_{j = 1}^{n} k \left( \frac{x - x_{j}}{\lambda} \right) Y_{j} = \frac{1}{n} \sum_{j = 1}^{n} w_{j} Y_{j}\]

Onde \(w_{j} = k \left( \frac{x - x_{j}}{\lambda} \right) / \lambda\), e \(K\) representa o Kernel, com: \[\int K(u)d(u) = 1.\]

p1 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise",
             xlab = "",
             ylab = "Y - Renda",
             pch = 16,
             main = expression(lambda~"= 1"))+
  layer(panel.lines(ksmooth(da$x, da$y, bandwidth = 1, 
                            kernel = "normal"),
                    lwd = 2,
                    col = 3))

p2 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise", 
             xlab = "X - Tempo de Estudo",
             ylab = "",
             pch = 16,
             main = expression(lambda~"= 4"))+
  layer(panel.lines(ksmooth(da$x, da$y, bandwidth = 4, 
                            kernel = "normal"),
                    lwd = 2,
                    col = 2))

p3 <- xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       pch = 16,
       ylab = "",
       xlab = "",
       main = expression(lambda~"= 7"))+
  layer(panel.lines(ksmooth(da$x, da$y, bandwidth = 7, 
                            kernel = "normal"),
                    lwd = 2,
                    col = 3))

print(p1, position = c(0, 0, 1/3, 1), more = TRUE)
print(p2, position = c(1/3, 0, 2/3, 1), more = TRUE)
print(p3, position = c(2/3, 0, 1, 1))
Figura 6: Suavização por Kernel com variação do parâmetro de suavização

Figura 6: Suavização por Kernel com variação do parâmetro de suavização

Modelo Aditivo Generalizado

#------ Instalação do pacote através do GitHub
# library(devtools)
# install_github("cran/mgcv")
#------

mg.gam1 <- mgcv::gam(y ~ log(x), data = da)
aic[4] <- AIC(mg.gam1)

mg.gam2 <- mgcv::gam(y ~ s(x), data = da)
aic[5] <- AIC(mg.gam2)

p1 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise", 
             pch = 16,
             ylab = "Y - Renda",
             xlab = "X - Anos de Estudo",
             main = paste0("GAM com log em x - AIC = ", round(aic[4], 3)),
             panel = function(x, y, ...) {
               panel.xyplot(x, y, ...)
               
               le <- mgcv::gam(da$y ~ log(da$x))
               p <- data.frame(predict(le))
               pp <- sort(p$predict.le., decreasing = FALSE)
               x <- sort(da$x)
               llines(x, pp, col = 3, lwd = 3)
             })

p2 <- xyplot(y ~ x, data = da, type = c("p", "g"),
             col = "turquoise", 
             pch = 16,
             ylab = "Y - Renda",
             xlab = "X - Anos de Estudo",
             main = paste0("GAM com spline em x - AIC = ", round(aic[5],
                                                                 3)),
             
             panel = function(x, y, ...) {
               panel.xyplot(x, y, ...)
               
               le <- mgcv::gam(da$y ~ s(da$x))
               p <- data.frame(predict(le))
               pp <- sort(p$predict.le., decreasing = FALSE)
               x <- sort(da$x)
               llines(x, pp, col = 2, lwd = 3)
             })

print(p1, position = c(0, 0, 0.5, 1), more = TRUE)
print(p2, position = c(0.5, 0, 1, 1), more = TRUE)
Figura 7: Modelos Lineares Generalizados: variações não lineares na resposta

Figura 7: Modelos Lineares Generalizados: variações não lineares na resposta

Resultados e Conclusão

xyplot(y ~ x, data = da, type = c("p", "g"),
       col = "turquoise", 
       pch = 16,
       ylab = "Y - Renda",
       xlab = "X - Tempo de estudo",      
       main = "Curvas Ajustadas",

       key=list(space="bottom",
                lines= list(col = c(2, 3, "yellow", "violetred",
                                    "orange", 4),
                            lwd=6),
                text=list(c("Reg. Linear",
                            "Reg. Por Partes",
                            expression("Suav. por Kernel - "~lambda~"= 7"),
                            "Suav. por Splines - 5 graus de liberdade",
                            "GAM com log em x",
                            "GAM com spline em x")), 
                columns = 2
       ),
       panel = function(x, y, ...) {
         
         panel.xyplot(x, y, ...)
         panel.abline(lm(y~x), col = 2, lwd = 3)
         
    le <- lm(y ~ x, da,
             subset = x <= 50)
    llines(da$x[da$x <= 50],
           predict(le), col = 3, lwd = 3)
    ld <- lm(y ~ x + p, da,
             subset = x > 50)
    llines(da$x[da$x > 50],
           predict(ld), col = 3, lwd = 3)
    #-------
    gam1 <- mgcv::gam(da$y ~ s(da$x))
    p <- data.frame(predict(gam1))
    pp <- sort(p$predict.gam1, decreasing = FALSE)
    x <- sort(da$x)
    llines(x, pp, col = "orange", lwd = 3)
    
    #-------
    gam2 <- mgcv::gam(da$y ~ log(da$x))
    p <- data.frame(predict(gam2))
    pp <- sort(p$predict.gam2, decreasing = FALSE)
    x <- sort(da$x)
    llines(x, pp, col = 4, lwd = 3)
    
  })+
  layer(panel.lines(ksmooth(da$x, da$y, bandwidth = 7, 
                            kernel = "normal"),
                    lwd = 3,
                    col = "yellow"))+
    layer(panel.lines(smooth.spline(da$x, da$y, df = 5),
                    lwd = 3,
                    col = "violetred"))
Figura 8: Comparação de ajustes por diferentes técnicas

Figura 8: Comparação de ajustes por diferentes técnicas

Referências Bibliográficas