Ajuste de modelo de regressão linear
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "reshape", "car", "alr3",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra reshape car alr3 plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
source("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()
## R version 3.1.1 (2014-07-10)
## 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] wzRfun_0.3 plyr_1.8.1 alr3_2.0.5 car_2.0-21
## [5] reshape_0.8.5 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [9] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] MASS_7.3-34 methods_3.1.1 nnet_7.3-8 Rcpp_0.11.3 stringr_0.6.2
## [11] tools_3.1.1 yaml_2.1.13
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
Distância para parada em função da velocidade
##-----------------------------------------------------------------------------
## Explorar os dados.
## Estrutura do Arquivo.
str(cars)
## 'data.frame': 50 obs. of 2 variables:
## $ speed: num 4 4 7 7 8 9 10 10 10 11 ...
## $ dist : num 2 10 4 22 16 10 18 26 34 17 ...
## Visualização.
xyplot(dist~speed, data=cars, type=c("p","smooth"))
##-----------------------------------------------------------------------------
## Estimação por mínimos quadrados via solução matricial.
X <- cbind(b0=1, b1=cars$speed)
head(X)
## b0 b1
## [1,] 1 4
## [2,] 1 4
## [3,] 1 7
## [4,] 1 7
## [5,] 1 8
## [6,] 1 9
y <- cbind(cars$dist)
## Estimativa dos parâmetros.
solve(t(X)%*%X)%*%t(X)%*%y
## [,1]
## b0 -17.579095
## b1 3.932409
##-----------------------------------------------------------------------------
## Usando a lm().
## y ~ Normal(modelo linear, sigma²)
## modelo linear: b0+b1*x.
m0 <- lm(dist~speed, data=cars)
## m0 <- lm(dist~0+speed, data=cars)
c0 <- coef(m0)
c0
## (Intercept) speed
## -17.579095 3.932409
## Sobrepondo ajuste às observações.
xyplot(dist~speed, data=cars, xlim=c(0,NA), ylim=c(c0[1],NA))+
layer(panel.abline(a=c0[1], b=c0[2], col=2))
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Falta de ajuste. Extender para um modelo que permita curvatura da
## função.
##-----------------------------------------------------------------------------
## Ajustar polinômio de segundo grau.
## modelo linear: b0+b1*x+b2*x^2.
m1 <- lm(dist~speed+I(speed^2), data=cars)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Mostra uma possível relação média~variância.
## Testa o abandono do termo extra por meio da mudança na soma de
## quadrados.
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: dist ~ speed + I(speed^2)
## Model 2: dist ~ speed
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 10825
## 2 48 11354 -1 -528.81 2.296 0.1364
summary(m1)
##
## Call:
## lm(formula = dist ~ speed + I(speed^2), data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.720 -9.184 -3.188 4.628 45.152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.47014 14.81716 0.167 0.868
## speed 0.91329 2.03422 0.449 0.656
## I(speed^2) 0.09996 0.06597 1.515 0.136
##
## Residual standard error: 15.18 on 47 degrees of freedom
## Multiple R-squared: 0.6673, Adjusted R-squared: 0.6532
## F-statistic: 47.14 on 2 and 47 DF, p-value: 5.852e-12
## Verifica se cabe uma tranformação.
MASS::boxcox(m1); abline(v=0.5, col=2)
## Usar lambda=0.5 como valor para transformar.
xyplot(sqrt(dist)~speed, data=cars)
##-----------------------------------------------------------------------------
## Com a variável transformada.
## Modelo quadrático com transformação na resposta.
m2 <- lm(sqrt(dist)~speed+I(speed^2), data=cars)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m2); layout(1)
## Teste para as estimativas, h0: beta==0.
summary(m2)
##
## Call:
## lm(formula = sqrt(dist) ~ speed + I(speed^2), data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0734 -0.7260 -0.1833 0.6369 3.1159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.990337 1.086712 0.911 0.367
## speed 0.365587 0.149193 2.450 0.018 *
## I(speed^2) -0.001429 0.004838 -0.295 0.769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.113 on 47 degrees of freedom
## Multiple R-squared: 0.71, Adjusted R-squared: 0.6976
## F-statistic: 57.52 on 2 and 47 DF, p-value: 2.334e-13
##-----------------------------------------------------------------------------
## Volta pro modelo mais simples.
## Modelo que abandona b2, ou seja, b2==0.
m3 <- lm(sqrt(dist)~speed, data=cars)
summary(m3)
##
## Call:
## lm(formula = sqrt(dist) ~ speed, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0684 -0.6983 -0.1799 0.5909 3.1534
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.27705 0.48444 2.636 0.0113 *
## speed 0.32241 0.02978 10.825 1.77e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.102 on 48 degrees of freedom
## Multiple R-squared: 0.7094, Adjusted R-squared: 0.7034
## F-statistic: 117.2 on 1 and 48 DF, p-value: 1.773e-14
## Teste para o abandono de b2 (o mesmo que o t do summary).
anova(m3, m2)
## Analysis of Variance Table
##
## Model 1: sqrt(dist) ~ speed
## Model 2: sqrt(dist) ~ speed + I(speed^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 48 58.334
## 2 47 58.226 1 0.10814 0.0873 0.7689
## Diagnóstico para o modelo final.
par(mfrow=c(2,2)); plot(m3); layout(1)
##-----------------------------------------------------------------------------
## Representa os resultados.
## Gráfico do modelo final.
c3 <- coef(m3)
xyplot(sqrt(dist)~speed, data=cars)+
layer(panel.abline(a=c3[1], b=c3[2], col=2))
## Na escala natural.
xyplot(dist~speed, data=cars, xlim=c(0,NA), ylim=c(0,NA))+
layer(panel.curve((c3[1]+c3[2]*x)^2, col=2))
##-----------------------------------------------------------------------------
## Com bandas de confiança.
pred <- data.frame(speed=seq(0, 30, length.out=30))
pred <- cbind(pred, predict(m3, newdata=pred, interval="confidence"))
str(pred)
## 'data.frame': 30 obs. of 4 variables:
## $ speed: num 0 1.03 2.07 3.1 4.14 ...
## $ fit : num 1.28 1.61 1.94 2.28 2.61 ...
## $ lwr : num 0.303 0.695 1.086 1.477 1.867 ...
## $ upr : num 2.25 2.53 2.8 3.08 3.35 ...
xyplot(sqrt(dist)~speed, data=cars)+
as.layer(xyplot(fit~speed, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))
##-----------------------------------------------------------------------------
## Na escala original.
i <- c("fit","lwr","upr")
pred[,i] <- pred[,i]^2
xyplot(dist~speed, data=cars, xlim=c(0,NA), ylim=c(0,NA))+
as.layer(xyplot(fit~speed, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))
Ganho de peso em perus em função da metionina na dieta
##-----------------------------------------------------------------------------
## Ganho de peso de perus em função de metionina na dieta.
str(turk0)
## 'data.frame': 35 obs. of 2 variables:
## $ A : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Gain: int 644 631 661 624 633 610 615 605 608 599 ...
## help(turk0, help_type="html")
## A
## Amount of methionine supplement (percent of diet)
## Gain
## Pen weight increase (g)
## Diagrama de dispersão com linha de tendência.
xyplot(Gain~A, data=turk0, type=c("p","smooth"),
xlab="Metionina (% da dieta)",
ylab="Ganho de peso (g)")
##-----------------------------------------------------------------------------
## Ajuste do modelo.
## Ajuste do modelo de regressão linear simples: b0+b1*x+b2*x^2;
m0 <- lm(Gain~A, turk0)
## Observados vs ajustados.
xyplot(Gain~A, data=turk0,
xlab="Metionina (% da dieta)",
ylab="Ganho de peso (g)")+
layer(panel.abline(m0))
## Resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Indiscutível falta de ajuste. Considerar um modelo que permita
## curvatura na função.
##-----------------------------------------------------------------------------
## Ajuste do modelo saturado (considerando A como fator). Considerar A
## como fator ocupa o mesmo espaço vetorial que ajustar um polinômio de
## grau k-1, em que k é o número de níveis. A equação da reta é o
## polinômio de grau um.
turk0$Afat <- factor(turk0$A)
m1 <- lm(Gain~poly(A, nlevels(Afat)-1), turk0)
anova(m1)
## Analysis of Variance Table
##
## Response: Gain
## Df Sum Sq Mean Sq F value Pr(>F)
## poly(A, nlevels(Afat) - 1) 5 150041 30008.2 88.587 < 2.2e-16 ***
## Residuals 29 9824 338.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m1 <- lm(Gain~Afat, turk0)
anova(m1)
## Analysis of Variance Table
##
## Response: Gain
## Df Sum Sq Mean Sq F value Pr(>F)
## Afat 5 150041 30008.2 88.587 < 2.2e-16 ***
## Residuals 29 9824 338.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Note que embora os modelos tenham especificações diferentes o quadro
## de anova é o mesmo. Isso porque ambos exploram o mesmo espaço
## vetorial de k-1 dimensões, apenas os vetores que formam a base desse
## espaço é que são diferentes. O conjunto de vetores que formam a base
## são as colunas da matriz do modelo, X.
## Resíduos.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Resíduos ok.
##-----------------------------------------------------------------------------
## Médias amostrais e os valores preditos.
## Médias amostrais.
aggregate(Gain~A, data=turk0, mean)
## A Gain
## 1 0.00 623.0
## 2 0.04 668.4
## 3 0.10 715.6
## 4 0.16 732.0
## 5 0.28 794.0
## 6 0.44 785.4
## Valores preditos pelo modelo saturado.
predict(m1, newdata=list(Afat=levels(turk0$Afat)))
## 1 2 3 4 5 6
## 623.0 668.4 715.6 732.0 794.0 785.4
## Valores preditos pelo modelo m0
predict(m0, newdata=list(A=unique(turk0$A)))
## 1 2 3 4 5 6
## 648.4899 664.2535 687.8988 711.5441 758.8348 821.8890
##-----------------------------------------------------------------------------
## Teste de falta de ajuste para o modelo m0.
anova(m0, m1)
## Analysis of Variance Table
##
## Model 1: Gain ~ A
## Model 2: Gain ~ Afat
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 33 35176
## 2 29 9824 4 25352 18.711 1.062e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Não passou na falta de ajuste conforme já se havia antecipado.
##-----------------------------------------------------------------------------
## Ajuste do modelo de segundo grau: b0+b1*x+b2*x^2;
m2 <- lm(Gain~A+I(A^2), turk0)
## m2 <- lm(Gain~poly(A, 2), turk0)
## Resíduos.
par(mfrow=c(2,2)); plot(m2); layout(1)
##-----------------------------------------------------------------------------
## Teste de falta de ajuste para o modelo m2.
anova(m1, m2)
## Analysis of Variance Table
##
## Model 1: Gain ~ Afat
## Model 2: Gain ~ A + I(A^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 29 9823.6
## 2 32 11339.9 -3 -1516.2 1.492 0.2374
## Model 1: Gain ~ Afat
## Model 2: Gain ~ A + I(A^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 29 9823.6
## 2 32 11339.9 -3 -1516.2 1.492 0.2374
## Não apresentou falta de ajuste. Então é um modelo que pode ser mantido.
##-----------------------------------------------------------------------------
## Gráfico dos valores preditos com bandas de confiança.
## range(turk0$A)
pred2 <- data.frame(A=seq(0, 0.45, by=0.025))
a <- predict(m2, newdata=pred2, interval="confidence")
pred2 <- cbind(pred2, a)
str(pred2)
## 'data.frame': 19 obs. of 4 variables:
## $ A : num 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 ...
## $ fit: num 626 649 670 690 708 ...
## $ lwr: num 615 640 663 682 700 ...
## $ upr: num 636 657 678 698 717 ...
xyplot(Gain~A, data=turk0,
xlab="Metionina (% da dieta)",
ylab="Ganho de peso (g)")+
as.layer(xyplot(fit~A, data=pred2, type="l",
ly=pred2$lwr, uy=pred2$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))
## Sobrepondo as estimativas de médias do modelo m1.
pred1 <- data.frame(Afat=levels(turk0$Afat))
a <- predict(m1, newdata=pred1, interval="confidence")
pred1 <- cbind(pred1, a)
pred1$A <- as.numeric(as.character(pred1$Afat))
## pred1
xyplot(Gain~A, data=turk0,
xlab="Metionina (% da dieta)",
ylab="Ganho de peso (g)")+
as.layer(xyplot(fit~A, data=pred2, type="l",
ly=pred2$lwr, uy=pred2$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))+
as.layer(xyplot(fit~(A+0.005), data=pred1, pch=19,
ly=pred1$lwr, uy=pred1$upr, cty="bars",
prepanel=prepanel.cbH, panel=panel.cbH))
##-----------------------------------------------------------------------------
## Uso de pesos no ajuste do modelo de regressão.
## Imagine que não tem os dados completos, apenas as médias e
## frequências.
tu <- ddply(turk0, .(A), summarise,
mGain=mean(Gain), nGain=length(Gain))
tu
## A mGain nGain
## 1 0.00 623.0 10
## 2 0.04 668.4 5
## 3 0.10 715.6 5
## 4 0.16 732.0 5
## 5 0.28 794.0 5
## 6 0.44 785.4 5
m3 <- lm(mGain~A+I(A^2), data=tu, weights=nGain)
summary(m3)
##
## Call:
## lm(formula = mGain ~ A + I(A^2), data = tu, weights = nGain)
##
## Weighted Residuals:
## 1 2 3 4 5 6
## -8.037 14.442 16.171 -29.042 11.611 -1.816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 625.542 6.243 100.201 2.19e-06 ***
## A 964.471 86.763 11.116 0.00156 **
## I(A^2) -1362.069 198.338 -6.867 0.00632 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.48 on 3 degrees of freedom
## Multiple R-squared: 0.9899, Adjusted R-squared: 0.9832
## F-statistic: 146.9 on 2 and 3 DF, p-value: 0.001016
summary(m2)
##
## Call:
## lm(formula = Gain ~ A + I(A^2), data = turk0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.988 -16.542 2.193 12.788 36.059
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 625.542 5.227 119.665 < 2e-16 ***
## A 964.471 72.651 13.275 1.46e-14 ***
## I(A^2) -1362.069 166.077 -8.201 2.28e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.82 on 32 degrees of freedom
## Multiple R-squared: 0.9291, Adjusted R-squared: 0.9246
## F-statistic: 209.6 on 2 and 32 DF, p-value: < 2.2e-16
## Observe as estimativas dos parâmetros. Por que o valor pontual é
## igual nos dois modelos? Por que o erro padrão é maior com o uso das
## médias? Por que o R² é maior com o uso das médias? O que iria
## acontecer se não fossem usados os pesos?
##-----------------------------------------------------------------------------
## Sobrepondo os valores preditos com bandas de confiança obtido com o
## ajuste com as médias.
pred3 <- data.frame(A=seq(0, 0.45, by=0.025))
a <- predict(m3, newdata=pred3, interval="confidence")
pred3 <- cbind(pred3, a)
## pred3
xyplot(Gain~A, data=turk0,
xlab="Metionina (% da dieta)",
ylab="Ganho de peso (g)")+
as.layer(xyplot(fit~A, data=pred2, type="l",
ly=pred2$lwr, uy=pred2$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))+
as.layer(xyplot(fit~A, data=pred3, type="l",
ly=pred3$lwr, uy=pred3$upr, cty="bands",
prepanel=prepanel.cbH, panel=panel.cbH))
