Especificação, suposições e estimação
Especificação do modelo
- De momentos (média e variância) \[
\begin{aligned}
Y|X &= X\beta+\epsilon, \\
\text{E}(\epsilon) &= 0, \\
\text{V}(\epsilon) &= \sigma^2.
\end{aligned}
\]
- Especificação completa (distribuição de probabilidades) \[
\begin{aligned}
Y|X &\sim \text{Normal}(\mu, \sigma^2), \\
\mu &= X\beta.
\end{aligned}
\]
##=============================================================================
## Aplicação de modelos de regressão linear e
## não linear em ciências agrárias
##
## 09 à 11 de Dezembro de 2014 - Goiânia/GO
## Embrapa Arroz e Feijão
##
## Prof. Dr. Walmes M. Zeviani
## LEG/DEST/UFPR
##=============================================================================
##-----------------------------------------------------------------------------
## Definições da sessão.
require(rpanel)
##-----------------------------------------------------------------------------
## Modelo de regressão linear simples.
rls <- function(panel){
with(panel, {
curve(b0+b1*x, xlim[1], xlim[2], ylim=ylim, lwd=2,
ylab=ylab, xlab=xlab, main=main)
ncur <- seq(xlim[1], xlim[2], length=nc)
ncur <- ncur[-length(ncur)]
lapply(ncur,
function(x){
xx <- seq(ylim[1],ylim[2], l=100)
mu <- b0+b1*x
s <- sqrt(s2)
y <- dnorm(xx, mu, s)/scale
ind <- y>0.0001
if(curvas){
lines(y[ind]+x, xx[ind], col="green4")
segments(x, mu, x+dnorm(mu, mu, s)/scale, mu,
col="green4", lty=2)
segments(x, mu-4*s, x, mu+4*s,
col="gray70", lty=2)
}
})
if(simula){
x <- rep(ncur, nr)
y <- rnorm(x, mean=b0+b1*x, sd=s2)
points(x, y)
if(ajusta){
m0 <- lm(y~x)
abline(m0, col=2)
}
}
})
panel
}
panel <- rp.control(xlim=c(0,10), ylim=c(0,10), nc=7, scale=1,
xlab="x", ylab="[Y|x]",
main="Regressão linear simples")
rp.slider(panel, b0, 0, 10, initval=5, showvalue=TRUE,
action=rls, res=0.05, title="b0")
rp.slider(panel, b1, -4, 4, initval=0, showvalue=TRUE,
action=rls, res=0.05, title="b1")
rp.slider(panel, s2, 0, 2, initval=0.1, showvalue=TRUE,
action=rls, res=0.01, title="s²")
rp.doublebutton(panel, variable=nr, initval=2, range=c(1, NA), step=1,
action=rls, showvalue=TRUE,
title="Repetições")
rp.doublebutton(panel, variable=nc, initval=7, range=c(2, NA), step=1,
action=rls, showvalue=TRUE,
title="Valores de x")
rp.checkbox(panel, curvas, labels=c("Modelo teórico"),
initval=FALSE, action=rls)
rp.checkbox(panel, simula, labels=c("Simular do modelo"),
initval=FALSE, action=rls)
rp.checkbox(panel, ajusta, labels=c("Ajustar o modelo"),
initval=FALSE, action=rls)
x11()
rp.do(panel, action=rls)
Estimação
Mínimos quadrados (ordinários)
- Princípio geométrico, projeção de espaços vetoriais.
- Objetivo de Minimizar soma dos desvios ao quadrado \[
\underset{\beta \in \mathbb{R}}{\arg \min}\, S(\beta) = \sum_{i=1}^{n}
(y_{i}-x_{i}^\top\beta)^2 = (y-X\beta)^{\top}(y-X\beta)
\]
- Equações normais \[
X^\top X \beta = X^\top y
\]
- Solução \[
\hat{\beta} = (X^\top X)^{-1} X^\top y
\]
Baseada na função de verossimilhança
- Princípio mais geral.
- Objetivo de Maximizar a verossilhança (log-verossimilhança). Se considerar o modelo tradicional de regressão (distribuição Gaussiana, independência, variância homogênea) \[
\begin{aligned}
\underset{\beta \in \mathbb{R},\, \sigma >0}{\arg \max}\, L(\beta,
\sigma^2) &= \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}}\cdot
\exp\left\{\frac{(y_{i}-x_{i}^\top\beta)^2}{2\sigma^2}\right\}\\
\underset{\beta \in \mathbb{R},\, \sigma >0}{\arg \max}\, \log
L(\beta, \sigma^2) &= -\frac{n}{2}
\log(2\pi\sigma^2)-\frac{1}{2\sigma^2} \sum_{i=1}^{n}
(y_{i}-x_{i}^\top\beta)^2
\end{aligned}
\]
- Solução (é a mesma em \(\beta\)) \[
\hat{\beta} = (X^\top X)^{-1} X^\top y
\]
Mínimos quadrados generalizados
- Permite acomodar casos onde a suposição de independendia ou variância constante não são verificadas.
Demonstrando a estimação
##-----------------------------------------------------------------------------
## Dados.
## Comprimento da diagonal de aparelhos celulares (cm).
x <- c(12.5, 13.2, 14.4, 11.3, 11.2, 11.2, 11.4, 11.0, 10.6, 12.4, 12.4)
## Peso dos aparelhos (g).
y <- c(136.7, 144.8, 152.6, 120.9, 113.6, 126.9, 91.3, 107.9, 103.6,
116.9, 117.3)
plot(y~x,
xlab="Comprimento diagonal (cm)",
ylab="Peso do aparelho celular (g)")
grid()
##-----------------------------------------------------------------------------
## Modelo: E(y|x) = b0+b1*x.
y <- matrix(y, ncol=1)
X <- cbind(b0=1, b1=x)
## De forma literal.
XlX <- t(X)%*%X
Xly <- t(X)%*%y
solve(XlX)%*%Xly
## [,1]
## b0 -33.22079
## b1 12.90219
## Equivalente porém numericamente superior.
solve(crossprod(X), crossprod(X, y))
## [,1]
## b0 -33.22079
## b1 12.90219
##-----------------------------------------------------------------------------
## Usando a função lm().
m0 <- lm(y~1+x)
m0 <- lm(y~x)
coef(m0)
## (Intercept) x
## -33.22079 12.90219
summary(m0)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.5642 -5.1349 0.0575 8.0189 15.6162
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.221 37.878 -0.877 0.40326
## x 12.902 3.153 4.092 0.00271 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.3 on 9 degrees of freedom
## Multiple R-squared: 0.6504, Adjusted R-squared: 0.6115
## F-statistic: 16.74 on 1 and 9 DF, p-value: 0.00271
##-----------------------------------------------------------------------------
## Observados e ajustados.
betas <- coef(m0)
plot(y~x,
xlab="Comprimento diagonal (cm)",
ylab="Peso do aparelho celular (g)")
grid()
abline(a=betas[1], b=betas[2], col=2)
##-----------------------------------------------------------------------------
## Visualização.
rp.regression(x=x, y=y)
Métodos numéricos considerados na estimação
##-----------------------------------------------------------------------------
## Usando decomposição QR de X.
QR <- qr(X, LAPACK=TRUE)
Q <- qr.Q(QR, complete=TRUE)
## Q é ortonormal.
round(t(Q)%*%Q, 3)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 1 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0 0
## [7,] 0 0 0 0 0 0 1 0 0 0 0
## [8,] 0 0 0 0 0 0 0 1 0 0 0
## [9,] 0 0 0 0 0 0 0 0 1 0 0
## [10,] 0 0 0 0 0 0 0 0 0 1 0
## [11,] 0 0 0 0 0 0 0 0 0 0 1
## R é triangular superior.
R <- qr.R(QR, complete=TRUE); R
## b1 b0
## [1,] -39.84043 -3.3031771
## [2,] 0.00000 0.2983644
## [3,] 0.00000 0.0000000
## [4,] 0.00000 0.0000000
## [5,] 0.00000 0.0000000
## [6,] 0.00000 0.0000000
## [7,] 0.00000 0.0000000
## [8,] 0.00000 0.0000000
## [9,] 0.00000 0.0000000
## [10,] 0.00000 0.0000000
## [11,] 0.00000 0.0000000
## Estimação.
z <- t(Q)%*%y; z
## [,1]
## [1,] -404.294816
## [2,] -9.911903
## [3,] -5.000805
## [4,] 7.483397
## [5,] 1.608694
## [6,] 14.908694
## [7,] -23.541900
## [8,] -1.240712
## [9,] 0.160475
## [10,] -12.194868
## [11,] -11.794868
betas <- backsolve(R, z); betas
## [,1]
## [1,] 12.90219
## [2,] -33.22079
##-----------------------------------------------------------------------------
## Usando decomposição de Cholesky.
XlX <- t(X)%*%X
XlX <- crossprod(X); XlX
## b0 b1
## b0 11.0 131.60
## b1 131.6 1587.26
## U é a triangular superior (upper), L é a inferior (lower).
U <- chol(XlX)
L <- t(U)
## LU = XlX.
L%*%U
## b0 b1
## b0 11.0 131.60
## b1 131.6 1587.26
## Estimação.
d <- t(X)%*%y
z <- forwardsolve(L, d); z
## [,1]
## [1,] 401.76387
## [2,] 46.24218
betas <- backsolve(U, z); betas
## [,1]
## [1,] -33.22079
## [2,] 12.90219
str(m0)
## List of 12
## $ coefficients : Named num [1:2] -33.2 12.9
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
## $ residuals : Named num [1:11] 8.6434 7.7118 0.0292 8.326 2.3162 ...
## ..- attr(*, "names")= chr [1:11] "1" "2" "3" "4" ...
## $ effects : Named num [1:11] -401.76 -46.24 -4.78 7.28 1.39 ...
## ..- attr(*, "names")= chr [1:11] "(Intercept)" "x" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:11] 128 137 153 113 111 ...
## ..- attr(*, "names")= chr [1:11] "1" "2" "3" "4" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:11, 1:2] -3.317 0.302 0.302 0.302 0.302 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:11] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "x"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## ..$ qraux: num [1:2] 1.3 1.31
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 9
## $ xlevels : Named list()
## $ call : language lm(formula = y ~ x)
## $ terms :Classes 'terms', 'formula' length 3 y ~ x
## .. ..- attr(*, "variables")= language list(y, x)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "x"
## .. .. .. ..$ : chr "x"
## .. ..- attr(*, "term.labels")= chr "x"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, x)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "nmatrix.1" "numeric"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
## $ model :'data.frame': 11 obs. of 2 variables:
## ..$ y: num [1:11, 1] 137 145 153 121 114 ...
## ..$ x: num [1:11] 12.5 13.2 14.4 11.3 11.2 11.2 11.4 11 10.6 12.4 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 y ~ x
## .. .. ..- attr(*, "variables")= language list(y, x)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "x"
## .. .. .. .. ..$ : chr "x"
## .. .. ..- attr(*, "term.labels")= chr "x"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, x)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "nmatrix.1" "numeric"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
## - attr(*, "class")= chr "lm"
##-----------------------------------------------------------------------------
## Informações da sessão.
Sys.time()
## [1] "2014-12-11 20:52:15 BRST"
sessionInfo()
## R version 3.1.2 (2014-10-31)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] rmarkdown_0.3.3 knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 stringr_0.6.2
## [6] tools_3.1.2 yaml_2.1.13