Ajuste de modelos de regressão linear
Walmes Zeviani
Definições da sessão
##-----------------------------------------------------------------------------
## 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.1 plyr_1.8.1 alr3_2.0.5 car_2.0-20
## [5] reshape_0.8.5 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [9] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 MASS_7.3-33 methods_3.1.1
## [6] nnet_7.3-8 Rcpp_0.11.1 stringr_0.6.2 tools_3.1.1
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
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 88.6 <2e-16 ***
## Residuals 29 9824 339
## ---
## 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 88.6 <2e-16 ***
## Residuals 29 9824 339
## ---
## 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.5 664.3 687.9 711.5 758.8 821.9
##-----------------------------------------------------------------------------
## 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 25353 18.7 1.1e-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 9824
## 2 32 11340 -3 -1516 1.49 0.24
## 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.04 14.44 16.17 -29.04 11.61 -1.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 625.54 6.24 100.20 2.2e-06 ***
## A 964.47 86.76 11.12 0.0016 **
## I(A^2) -1362.07 198.34 -6.87 0.0063 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.5 on 3 degrees of freedom
## Multiple R-squared: 0.99, Adjusted R-squared: 0.983
## F-statistic: 147 on 2 and 3 DF, p-value: 0.00102
summary(m2)
##
## Call:
## lm(formula = Gain ~ A + I(A^2), data = turk0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.99 -16.54 2.19 12.79 36.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 625.54 5.23 119.7 < 2e-16 ***
## A 964.47 72.65 13.3 1.5e-14 ***
## I(A^2) -1362.07 166.08 -8.2 2.3e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.8 on 32 degrees of freedom
## Multiple R-squared: 0.929, Adjusted R-squared: 0.925
## F-statistic: 210 on 2 and 32 DF, p-value: <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))
