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 |
Prof. Dr. Walmes M. Zeviani |
Embrapa Arroz e Feijão |
Lab. de Estatística e Geoinformação - LEG |
Departamento de Estatística - UFPR |
http://web.stanford.edu/class/stats191/notebooks/Diagnostics%20for%20multiple%20regression.pdf
##=============================================================================
## 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.
pkg <- c("lattice", "latticeExtra", "gridExtra", "car", "alr3", "asbio",
"plyr", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra car alr3 asbio
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr wzRfun
## TRUE TRUE
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Dados.
url <- "http://www.leg.ufpr.br/~walmes/data/
business_economics_dataset/EXAMPLES/REALESTA.DAT"
da <-
read.table(paste0(strwrap(url), collapse=""), header=FALSE)
da[,1] <- NULL
names(da) <- c("saleprice","landvalue","improvvalue","area")
str(da)
## 'data.frame': 20 obs. of 4 variables:
## $ saleprice : int 68900 48500 55500 62000 116500 45000 38000 83000 59000 47500 ...
## $ landvalue : int 5960 9000 9500 10000 18000 8500 8000 23000 8100 9000 ...
## $ improvvalue: int 44967 27860 31439 39592 72827 27317 29856 47752 39117 29349 ...
## $ area : int 1873 928 1126 1265 2214 912 899 1803 1204 1725 ...
summary(da)
## saleprice landvalue improvvalue area
## Min. : 22400 Min. : 1500 Min. : 5779 Min. : 899
## 1st Qu.: 40800 1st Qu.: 6965 1st Qu.:26351 1st Qu.:1054
## Median : 49250 Median : 8050 Median :31559 Median :1188
## Mean : 56660 Mean : 9213 Mean :35311 Mean :1383
## 3rd Qu.: 63725 3rd Qu.: 9625 3rd Qu.:41366 3rd Qu.:1744
## Max. :116500 Max. :23000 Max. :72827 Max. :2455
pairs(da)
scatterplotMatrix(da)
## Dividir para ficar mais tratável.
da <- da/1000
##-----------------------------------------------------------------------------
## Ver.
## p1 <- xyplot(saleprice~landvalue, data=da, type=c("p","smooth"))
## p2 <- xyplot(saleprice~improvvalue, data=da, type=c("p","smooth"))
## p3 <- xyplot(saleprice~area, data=da, type=c("p","smooth"))
## grid.arrange(p1, p2, p3)
xyplot.list(list(saleprice~improvvalue,
saleprice~landvalue,
saleprice~area),
data=da, x.same=FALSE, y.same=TRUE, layout=c(3,1),
xlab=list(c("Improvvalue","Landvalue","Area")),
ylab="Sales price (/1000)", type=c("p","smooth"))
##-----------------------------------------------------------------------------
## Ajuste.
m0 <- lm(saleprice~landvalue+improvvalue+area, data=da)
## class(m0)
methods(class="lm")
## [1] add1.lm* alias.lm* anova.lm*
## [4] Anova.lm* avPlot.lm* bootCase.lm*
## [7] Boot.lm* boxCox.lm* case.names.lm*
## [10] ceresPlot.lm* confidenceEllipse.lm* confint.lm
## [13] cooks.distance.lm* crPlot.lm* deltaMethod.lm*
## [16] deviance.lm* dfbeta.lm* dfbetaPlots.lm*
## [19] dfbetas.lm* dfbetasPlots.lm* drop1.lm*
## [22] dummy.coef.lm durbinWatsonTest.lm* effects.lm*
## [25] extractAIC.lm* family.lm* formula.lm*
## [28] hatvalues.lm* hccm.lm* infIndexPlot.lm*
## [31] influence.lm* influencePlot.lm* inverseResponsePlot.lm*
## [34] kappa.lm labels.lm* leveneTest.lm*
## [37] leveragePlot.lm* linearHypothesis.lm* logLik.lm*
## [40] mcPlot.lm* mmp.lm* model.frame.lm*
## [43] model.matrix.lm ncvTest.lm* nextBoot.lm*
## [46] nobs.lm* outlierTest.lm* plot.lm*
## [49] pod.lm* powerTransform.lm* predict.lm
## [52] print.lm* proj.lm* pureErrorAnova.lm*
## [55] qqPlot.lm* qr.lm* randomLinComb.lm*
## [58] residualPlot.lm* residualPlots.lm* residuals.lm
## [61] rstandard.lm* rstudent.lm* sigmaHat.lm*
## [64] simulate.lm* spreadLevelPlot.lm* summary.lm
## [67] variable.names.lm* vcov.lm*
##
## Non-visible functions are asterisked
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
par(mfrow=c(2,3)); plot(m0, which=1:6); layout(1)
## plot(m0, which=2)
residualPlots(m0)
## Test stat Pr(>|t|)
## landvalue 0.870 0.398
## improvvalue 1.742 0.102
## area 0.358 0.725
## Tukey test 1.652 0.099
## Teste para os efeitos.
anova(m0)
## Analysis of Variance Table
##
## Response: saleprice
## Df Sum Sq Mean Sq F value Pr(>F)
## landvalue 1 6102.2 6102.2 97.296 3.323e-08 ***
## improvvalue 1 2412.8 2412.8 38.470 1.267e-05 ***
## area 1 264.7 264.7 4.220 0.05666 .
## Residuals 16 1003.5 62.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
##
## Call:
## lm(formula = saleprice ~ landvalue + improvvalue + area, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.662 -2.155 1.027 2.722 15.976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4703 5.7463 0.256 0.80132
## landvalue 0.8145 0.5122 1.590 0.13137
## improvvalue 0.8204 0.2112 3.885 0.00131 **
## area 13.5286 6.5857 2.054 0.05666 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.919 on 16 degrees of freedom
## Multiple R-squared: 0.8974, Adjusted R-squared: 0.8782
## F-statistic: 46.66 on 3 and 16 DF, p-value: 3.9e-08
## Pouco efeito para o tamanho do imóvel? Não parece estrano?
## Quadros de somas de quadrados marginais.
drop1(m0, test="F")
## Single term deletions
##
## Model:
## saleprice ~ landvalue + improvvalue + area
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1003.5 86.310
## landvalue 1 158.58 1162.1 87.245 2.5285 0.131370
## improvvalue 1 946.60 1950.1 97.598 15.0929 0.001315 **
## area 1 264.67 1268.2 88.992 4.2200 0.056664 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(m0, type="III")
## Anova Table (Type III tests)
##
## Response: saleprice
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.11 1 0.0655 0.801316
## landvalue 158.58 1 2.5285 0.131370
## improvvalue 946.60 1 15.0929 0.001315 **
## area 264.67 1 4.2200 0.056664 .
## Residuals 1003.49 16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Ajuste do modelo, forma matricial, apenas para praticar.
## Matriz do modelo e vetor de observações.
X <- with(da, cbind(b0=1, b1=landvalue, b2=improvvalue, b3=area))
y <- with(da, cbind(saleprice))
head(X)
## b0 b1 b2 b3
## [1,] 1 5.96 44.967 1.873
## [2,] 1 9.00 27.860 0.928
## [3,] 1 9.50 31.439 1.126
## [4,] 1 10.00 39.592 1.265
## [5,] 1 18.00 72.827 2.214
## [6,] 1 8.50 27.317 0.912
## Produto de matrizes e inversa.
XlX <- crossprod(X)
Xly <- crossprod(X, y)
iXlX <- solve(XlX)
## Estimativas dos parâmetros.
b <- solve(XlX, Xly)
b
## saleprice
## b0 1.4702759
## b1 0.8144901
## b2 0.8204447
## b3 13.5286499
## Matrizes de projeção das observações (I) e do modelo (H).
I <- diag(nrow(y))
H <- X%*%solve(XlX, t(X))
## Função para calcular o traço de matriz.
trace <- function(matrix){
sum(diag(matrix))
}
## Traço das matrizes de projeção que são os graus de liberdade.
n <- trace(I)
p <- trace(H)
## Matriz de projeção da média (modelo nulo).
H1 <- (I*0+1)/n
## Soma de quadrados da regressão (corrigida) pra média e dos resíduos.
SQReg <- t(y)%*%(H-H1)%*%y
SQRes <- t(y)%*%(I-H)%*%y
## Quadrados médios.
QMReg <- SQReg/(p-1)
QMR <- SQRes/(n-p)
rQMR <- sqrt(QMR)
## Valor de F para ajuste do modelo.
F <- QMReg/QMR
pf(F, p-1, n-p, lower=FALSE)
## saleprice
## saleprice 3.899898e-08
## Valores preditos pelo modelo.
hy <- crossprod(H, y)
Os valores de alavancagem, de forma simples e direta, representam o peso que cada observação tem sobre o seu valor predito considerando o modelo construído. São os elementos da diagonal da matriz de projeção H. A soma desses valores é p, ou seja, é o traço da matriz. Sendo assim, em média, o peso de cada observação é p/n. Considera-se suspeitas as observações que possuem um valor de alavancagem duas vezes superior ao esperado.
\[ \begin{aligned} h_i &= H_{ii} \\ h &= \text{diag}(H) = \text{diag}(X(X^\top X)^{-1}X^\top) \end{aligned} \]
##-----------------------------------------------------------------------------
## Pontos de alavancagem.
hii <- diag(H)
plot(hii, type="h", ylim=c(0, 1.1*(max(hii))))
abline(h=2*p/n, lty=2, col=2)
## hatvalues() retorna os valores de alavancagem.
cbind(hatvalues(m0), hii)
## hii
## 1 0.29721478 0.29721478
## 2 0.14121680 0.14121680
## 3 0.08509103 0.08509103
## 4 0.08613511 0.08613511
## 5 0.36730340 0.36730340
## 6 0.13688377 0.13688377
## 7 0.14533094 0.14533094
## 8 0.52794088 0.52794088
## 9 0.10542731 0.10542731
## 10 0.20855291 0.20855291
## 11 0.17973851 0.17973851
## 12 0.09255776 0.09255776
## 13 0.38069000 0.38069000
## 14 0.09413692 0.09413692
## 15 0.12441213 0.12441213
## 16 0.22958162 0.22958162
## 17 0.18444624 0.18444624
## 18 0.09688042 0.09688042
## 19 0.21508782 0.21508782
## 20 0.30137165 0.30137165
Os resíduos crus são nada mais que a diferença entre valores observados e ajustados. Estão na escala da própria resposta. Têm valor esperado igual a zero. Embora seja uma suposição do modelo a independência condicional de y a x por meio do modelo, ao contrário do que se imagina, os resíduos crus não são independentes. O valor esperado é 0 mas a variância é \(\sigma^2(I-H)\). É por essa razão que não se recomenda aplicar testes de hipótese, como de normalidade, aos resíduos crus, por exemplo.
\[ \begin{aligned} \hat{e}_i &= y_i - \hat{y}_i\\ \hat{e} &= y - \hat{y}\\ \hat{e} &= y - X\hat{\beta} \end{aligned} \]
##-----------------------------------------------------------------------------
## Resíduos crus (ordinários).
e <- c(crossprod(I-H, y))
## e <- y-hy
## e <- y-crossprod(X, b)
plot(e); abline(h=0, lty=2)
## residuals() retorna os resíduos crus.
cbind(residuals(m0), e, da$saleprice-fitted(m0))
## e
## 1 0.34326507 0.34326507 0.34326507
## 2 4.28713670 4.28713670 4.28713670
## 3 5.26484739 5.26484739 5.26484739
## 4 2.78803439 2.78803439 2.78803439
## 5 10.66594524 10.66594524 10.66594524
## 6 1.85634162 1.85634162 1.85634162
## 7 -6.64364995 -6.64364995 -6.64364995
## 8 -0.77357950 -0.77357950 -0.77357950
## 9 2.55052449 2.55052449 2.55052449
## 10 -8.71683945 -8.71683945 -8.71683945
## 11 -14.48097731 -14.48097731 -14.48097731
## 12 -14.66237009 -14.66237009 -14.66237009
## 13 -1.97713293 -1.97713293 -1.97713293
## 14 2.69961720 2.69961720 2.69961720
## 15 -0.10013601 -0.10013601 -0.10013601
## 16 -2.68694921 -2.68694921 -2.68694921
## 17 15.97560222 15.97560222 15.97560222
## 18 -0.05204213 -0.05204213 -0.05204213
## 19 1.71028447 1.71028447 1.71028447
## 20 1.95207778 1.95207778 1.95207778
Os resíduos padronizados são resultado da padronização dos resíduos crus. Ao padronizar, ou seja, dividir pela correspondente variância \(\sigma^2(I-H)\), têm-se resíduos com variância unitária. Esses resíduos são chamados também de resíduos internamente studentizados.
\[ \hat{r}_i = \dfrac{\hat{e}_i}{s(\hat{e}_i)} = \dfrac{\hat{e}_i}{\hat{\sigma}\sqrt{1-h_{i}}} \]
##-----------------------------------------------------------------------------
## Resíduos padronizados ou internamente studentizados.
r <- e/sqrt(QMR*(1-hii))
plot(r); abline(h=0)
## rstandard() retorna os resíduos padronizados.
cbind(rstandard(m0), r)
## r
## 1 0.051703685 0.051703685
## 2 0.584155884 0.584155884
## 3 0.695024379 0.695024379
## 4 0.368264897 0.368264897
## 5 1.693186488 1.693186488
## 6 0.252305347 0.252305347
## 7 -0.907425428 -0.907425428
## 8 -0.142170649 -0.142170649
## 9 0.340506082 0.340506082
## 10 -1.237232437 -1.237232437
## 11 -2.018946947 -2.018946947
## 12 -1.943559660 -1.943559660
## 13 -0.317237864 -0.317237864
## 14 0.358157543 0.358157543
## 15 -0.013512746 -0.013512746
## 16 -0.386544354 -0.386544354
## 17 2.233747784 2.233747784
## 18 -0.006914895 -0.006914895
## 19 0.243759198 0.243759198
## 20 0.294901656 0.294901656
Os resíduos studentizados são considerados independentes pelo fato de serem resíduos decorrentes de procedimento leave-one-out. Para todos os efeitos, é como se o resíduo padronizado da observação i fosse calculado removendo-se o i-ésimo registro e ajustado o modelo. Não há necessidade de remover uma observação a cada vez para calcular tais resíduos pois tem-se fórmulas apropriadas para isso. Com isso, tem-se que \(y_i\) e \(y_{i(-i)}\) são independentes. Testes de influência consideram esses resíduos.
\[ \begin{aligned} \hat{t}_i &= \dfrac{\hat{e}_i}{s(\hat{e}_i)} = \dfrac{\hat{e}_i}{\hat{\sigma}_{-i}\sqrt{1-h_{i}}} = \hat{r}_i\left(\frac{n-p-1}{n-p-\hat{r}_i^2} \right)^{1/2}\\ \hat{\sigma}_{-i}^2 &= \dfrac{(n-p)\hat{\sigma}^2-\frac{\hat{e}_i^2}{1-h_{i}}}{(n-1)-p} \end{aligned} \]
##-----------------------------------------------------------------------------
## Resíduos studentizados.
t <- r*((n-p-1)/(n-p-r^2))^0.5
## rstudent() retorna os resíduos studentizados.
cbind(rstudent(m0), t)
## t
## 1 0.050066061 0.050066061
## 2 0.571736179 0.571736179
## 3 0.683349077 0.683349077
## 4 0.358091810 0.358091810
## 5 1.809532865 1.809532865
## 6 0.244781033 0.244781033
## 7 -0.902131048 -0.902131048
## 8 -0.137743171 -0.137743171
## 9 0.330894694 0.330894694
## 10 -1.259719431 -1.259719431
## 11 -2.264447340 -2.264447340
## 12 -2.153089760 -2.153089760
## 13 -0.308134853 -0.308134853
## 14 0.348183103 0.348183103
## 15 -0.013083735 -0.013083735
## 16 -0.376029863 -0.376029863
## 17 2.607226676 2.607226676
## 18 -0.006695329 -0.006695329
## 19 0.236458300 0.236458300
## 20 0.286316488 0.286316488
## Teste para outlier.
apropos("outlier")
## [1] "outlier.test" "outlierTest" "outlier.t.test"
outlierTest(m0)
##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 17 2.607227 0.019812 0.39624
## 2*c(t=1, bonferroni=n)*
## pt(max(abs(rstudent(m0))), df=n-p-1, lower.tail=FALSE)
Os \(n\) vetores de estimativas dos parâmetros considerando a remoção da \(i\)-ésima observação em cada vez são obtidas por
\[ \begin{aligned} \hat{beta}_{(-i)} = \hat{\beta}-\hat{e}_i\frac{(X^\top X)^{-1} x_{i}}{1-h_{i}}. \end{aligned} \]
A partir dessa medida, são obtidos os valores preditos, sem a \(i\)-ésima observação também, por
\[ \begin{aligned} \hat{y}_{(-i)} = x_i^\top\hat{\beta}_{(-i)}. \end{aligned} \]
##-----------------------------------------------------------------------------
## Estatísticas leave-one-out.
## Estimativas dos parâmetros sem a i-ésima observação.
bi <- 0*t(X)
for(i in 1:n){
## bi[,i] <- solve(crossprod(X[-i,]))%*%t(X[-i,])%*%y[-i,]
bi[,i] <- b-e[i]*crossprod(iXlX, X[i,])/(1-hii[i])
}
t(bi)
## b0 b1 b2 b3
## [1,] 1.5495414 0.8280208 0.8190054 13.400322
## [2,] 0.3553034 0.7589914 0.8170813 14.609679
## [3,] 0.6038282 0.7699289 0.8185515 14.292098
## [4,] 1.1775069 0.8085580 0.8082694 13.980312
## [5,] 5.2862179 0.8539770 0.6561456 14.091765
## [6,] 0.9837925 0.7937361 0.8186134 13.987550
## [7,] 3.2226525 0.8605318 0.8451833 11.604738
## [8,] 1.4400526 0.8843033 0.8099139 13.413586
## [9,] 1.1692898 0.8280933 0.8039131 13.974565
## [10,] 0.6624660 0.7969022 0.7273459 17.004255
## [11,] 4.0310386 0.6898487 0.9741052 9.223385
## [12,] 1.3149387 0.7403023 0.7627215 16.192478
## [13,] 0.5443722 0.8553346 0.8046348 14.444900
## [14,] 1.0614468 0.7940970 0.8341251 13.503081
## [15,] 1.4852686 0.8154461 0.8196796 13.535108
## [16,] 0.8715999 0.7514001 0.8434409 13.920655
## [17,] 0.6960190 1.0269219 0.9385098 8.951832
## [18,] 1.4802808 0.8142129 0.8204106 13.526217
## [19,] 1.2034635 0.8586896 0.8017954 13.824444
## [20,] 0.9792723 0.8281335 0.8485870 12.973374
## pairs(t(bi))
## Valores preditos sem a i-ésima observação.
## hyi <- diag(X%*%bi)
hyi <- numeric(n)
for(i in 1:n){
hyi[i] <- X[i,]%*%bi[,i]
}
plot(hyi~hy, asp=1); abline(a=0, b=1, lty=2)
## Quadrado médio sem a i-ésima observação.
QMRi <- ((n-p)*QMR-(e^2/(1-hii)))/(n-1-p)
##-----------------------------------------------------------------------------
## Resíduos studentizados ou externamente studentizados (ou forma de
## calcular).
s <- e/sqrt(QMRi*(1-hii))
cbind(rstudent(m0), s, t)
## s t
## 1 0.050066061 0.050066061 0.050066061
## 2 0.571736179 0.571736179 0.571736179
## 3 0.683349077 0.683349077 0.683349077
## 4 0.358091810 0.358091810 0.358091810
## 5 1.809532865 1.809532865 1.809532865
## 6 0.244781033 0.244781033 0.244781033
## 7 -0.902131048 -0.902131048 -0.902131048
## 8 -0.137743171 -0.137743171 -0.137743171
## 9 0.330894694 0.330894694 0.330894694
## 10 -1.259719431 -1.259719431 -1.259719431
## 11 -2.264447340 -2.264447340 -2.264447340
## 12 -2.153089760 -2.153089760 -2.153089760
## 13 -0.308134853 -0.308134853 -0.308134853
## 14 0.348183103 0.348183103 0.348183103
## 15 -0.013083735 -0.013083735 -0.013083735
## 16 -0.376029863 -0.376029863 -0.376029863
## 17 2.607226676 2.607226676 2.607226676
## 18 -0.006695329 -0.006695329 -0.006695329
## 19 0.236458300 0.236458300 0.236458300
## 20 0.286316488 0.286316488 0.286316488
## Teste t com correção de bonferroni para outlier.
pvals <- 2*pt(abs(s), df=n-1-p, lower.tail=FALSE)
p.adjust(pvals, method="bonferroni")
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [9] 1.0000000 1.0000000 0.7759051 0.9598583 1.0000000 1.0000000 1.0000000 1.0000000
## [17] 0.3962370 1.0000000 1.0000000 1.0000000
\[ \begin{aligned} D_i &= \frac{(\hat{\beta}_{(-i)}-\hat{\beta})^\top (X^\top X) (\hat{\beta}_{(-i)}-\hat{\beta})}{p\hat{\sigma}^2}\\ &= \dfrac{(\hat{y}-\hat{y}_{i(-i)})^\top (\hat{y}-\hat{y}_{i(-i)})}{p\hat{\sigma}^2} \\ &= \dfrac{1}{p}\cdot\dfrac{h_i}{(1-h_i)} \cdot\dfrac{\hat{e}_i^2}{\hat{\sigma}^2(1-h_i)}\\ &= \dfrac{1}{p}\cdot\dfrac{h_i}{(1-h_i)}\cdot r_i^2 \end{aligned} \]
##-----------------------------------------------------------------------------
## Distância de Cook.
D <- hii*(e^2)/(p*QMR*(1-hii)^2)
plot(D)
## cooks.distance() retorna a distância de Cook.
cbind(cooks.distance(m0), D)
## D
## 1 2.826382e-04 2.826382e-04
## 2 1.402815e-02 1.402815e-02
## 3 1.123171e-02 1.123171e-02
## 4 3.195647e-03 3.195647e-03
## 5 4.160821e-01 4.160821e-01
## 6 2.523920e-03 2.523920e-03
## 7 3.500435e-02 3.500435e-02
## 8 5.651306e-03 5.651306e-03
## 9 3.416074e-03 3.416074e-03
## 10 1.008410e-01 1.008410e-01
## 11 2.232948e-01 2.232948e-01
## 12 9.632292e-02 9.632292e-02
## 13 1.546584e-02 1.546584e-02
## 14 3.332619e-03 3.332619e-03
## 15 6.486198e-06 6.486198e-06
## 16 1.113138e-02 1.113138e-02
## 17 2.821145e-01 2.821145e-01
## 18 1.282336e-06 1.282336e-06
## 19 4.070585e-03 4.070585e-03
## 20 9.378872e-03 9.378872e-03
\[ dffits_i = \dfrac{\hat{y}_i-\hat{y}_{i(-i))}}{\hat{\sigma}_{-i}\sqrt{h_i}} = \hat{t}_i\left( \dfrac{h_i}{1-h_i} \right )^{1/2} \]
##-----------------------------------------------------------------------------
## DFfits.
## dffits <- (c(hy)-c(hyi))/sqrt(QMRi*hii)
dffits <- s*sqrt((hii/(1-hii)))
plot(dffits)
## dffits() retorna os valores de DFfits.
cbind(dffits(m0), dffits)
## dffits
## 1 0.032558720 0.032558720
## 2 0.231844656 0.231844656
## 3 0.208398965 0.208398965
## 4 0.109936900 0.109936900
## 5 1.378736269 1.378736269
## 6 0.097480803 0.097480803
## 7 -0.372005775 -0.372005775
## 8 -0.145668122 -0.145668122
## 9 0.113594823 0.113594823
## 10 -0.646652572 -0.646652572
## 11 -1.060001884 -1.060001884
## 12 -0.687636725 -0.687636725
## 13 -0.241586417 -0.241586417
## 14 0.112242260 0.112242260
## 15 -0.004931888 -0.004931888
## 16 -0.205270991 -0.205270991
## 17 1.239902094 1.239902094
## 18 -0.002192892 -0.002192892
## 19 0.123780417 0.123780417
## 20 0.188050485 0.188050485
\[ \begin{aligned} dbetas_i &= \dfrac{\hat{\beta}- \hat{\beta}_{(-i)}}{\hat{\sigma}_{(-i)} \sqrt{\text{diag}((X^\top X)^{-1})}} \end{aligned} \]
##-----------------------------------------------------------------------------
## DFbetas.
## dfbetas <- 0*t(X)
## for(i in 1:n){
## num <- e[i]*c(crossprod(iXlX, X[i,]))/(1-hii[i])
## den <- sqrt(QMRi[i]*diag(iXlX))
## dfbetas[,i] <- num/den
## }
## num <- sweep(-t(bi), 2, c(b), "+")
num <- t(sweep(iXlX%*%t(X), 2, e/(1-hii), "*"))
den <- outer(sqrt(QMRi), sqrt(diag(iXlX)), "*")
dfbetas <- num/den
## dfbetas() retorna os valores de DFbetas.
dfbetas(m0)
## (Intercept) landvalue improvvalue area
## 1 -0.013357206 -0.0255791757 0.0065994417 0.0188687060
## 2 0.189906987 0.1060459514 0.0155876421 -0.1606585314
## 3 0.148250023 0.0855351239 0.0088139472 -0.1139780947
## 4 0.049541486 0.0112613722 0.0560596793 -0.0666878921
## 5 -0.709697529 -0.0823869963 0.8314457891 -0.0913814452
## 6 0.082135170 0.0393095493 0.0084129640 -0.0676033687
## 7 -0.303176808 -0.0893623801 -0.1164585556 0.2904312440
## 8 0.005095804 -0.1320512114 0.0483122339 0.0169277796
## 9 0.050900401 -0.0258076893 0.0760704753 -0.0657985087
## 10 0.143133594 0.0349608993 0.4488526692 -0.5373439259
## 11 -0.499823291 0.2729255496 -0.8160871901 0.7332237670
## 12 0.029946754 0.1604506692 0.3027971856 -0.4480947782
## 13 0.156506181 -0.0774521758 0.0727144700 -0.1351353794
## 14 0.069164813 0.0387045021 -0.0629752868 0.0037743780
## 15 -0.002526252 -0.0018071770 0.0035078271 -0.0009495220
## 16 0.101350226 0.1198196523 -0.1059293160 -0.0579048016
## 17 0.157267714 -0.4840706407 -0.6525340542 0.8111621803
## 18 -0.001685814 0.0005240004 0.0001563683 0.0003576845
## 19 0.045041147 -0.0837058087 0.0856631047 -0.0435694922
## 20 0.082959029 -0.0258604588 -0.1293794865 0.0818610922
##-----------------------------------------------------------------------------
## Medidas de influência.
influencePlot(m0)
## StudRes Hat CookD
## 5 1.8095329 0.3673034 0.64504427
## 8 -0.1377432 0.5279409 0.07517517
## 17 2.6072267 0.1844462 0.53114456
dfbetasPlots(m0)
## Envelopes obtidos via simulação.
qqPlot(m0)
## qnorm(0.25)
## As medidas de influência.
im <- influence.measures(m0)
summary(im)
## Potentially influential observations of
## lm(formula = saleprice ~ landvalue + improvvalue + area, data = da) :
##
## dfb.1_ dfb.lndv dfb.impr dfb.area dffit cov.r cook.d hat
## 1 -0.01 -0.03 0.01 0.02 0.03 1.84_* 0.00 0.30
## 8 0.01 -0.13 0.05 0.02 -0.15 2.73_* 0.01 0.53
## 13 0.16 -0.08 0.07 -0.14 -0.24 2.04_* 0.02 0.38
## 20 0.08 -0.03 -0.13 0.08 0.19 1.81_* 0.01 0.30
m1 <- update(m0, data=da[-8,])
im <- influence.measures(m1)
summary(im)
## Potentially influential observations of
## lm(formula = saleprice ~ landvalue + improvvalue + area, data = da[-8, ]) :
##
## dfb.1_ dfb.lndv dfb.impr dfb.area dffit cov.r cook.d hat
## 1 -0.03 -0.06 0.02 0.04 0.07 2.09_* 0.00 0.37
## 13 0.27 -0.28 0.20 -0.19 -0.48 2.79_* 0.06 0.55
## 20 0.09 -0.04 -0.12 0.09 0.20 1.86_* 0.01 0.31
## TRUE indica influente.
str(im)
## List of 3
## $ infmat: num [1:19, 1:8] -0.0259 0.1776 0.1382 0.0453 -0.7261 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:19] "1" "2" "3" "4" ...
## .. ..$ : chr [1:8] "dfb.1_" "dfb.lndv" "dfb.impr" "dfb.area" ...
## $ is.inf: logi [1:19, 1:8] FALSE FALSE FALSE FALSE FALSE FALSE ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:19] "1" "2" "3" "4" ...
## .. ..$ : chr [1:8] "dfb.1_" "dfb.lndv" "dfb.impr" "dfb.area" ...
## $ call : language lm(formula = saleprice ~ landvalue + improvvalue + area, data = da[-8, ])
## - attr(*, "class")= chr "infl"
im$is.inf
## dfb.1_ dfb.lndv dfb.impr dfb.area dffit cov.r cook.d hat
## 1 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## 2 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 3 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 5 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 6 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 7 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 9 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 10 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 11 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 12 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 13 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## 14 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 15 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 16 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 17 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 18 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 19 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 20 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
##-----------------------------------------------------------------------------
## Fator de inflação da variância (VIF).
## Função vif da página da Professora Dra Sueli Giolo (DEST/UFPR).
## source("http://people.ufpr.br/~giolo/CE071/Exemplos/vif.R")
pairs(da)
## Correlação entre variáveis.
## colwise(is.numeric)(da)
cor(da)
## saleprice landvalue improvvalue area
## saleprice 1.0000000 0.7897767 0.9156607 0.8489815
## landvalue 0.7897767 1.0000000 0.7289003 0.6889384
## improvvalue 0.9156607 0.7289003 1.0000000 0.7881460
## area 0.8489815 0.6889384 0.7881460 1.0000000
## Covariância entre estimativas.
vcov(m0)
## (Intercept) landvalue improvvalue area
## (Intercept) 33.02024622 0.44219936 -0.05073683 -23.2527814
## landvalue 0.44219936 0.26236801 -0.04508018 -0.9162945
## improvvalue -0.05073683 -0.04508018 0.04459908 -0.8015259
## area -23.25278143 -0.91629452 -0.80152590 43.3711819
cov2cor(vcov(m0))
## (Intercept) landvalue improvvalue area
## (Intercept) 1.00000000 0.1502355 -0.04180905 -0.6144466
## landvalue 0.15023548 1.0000000 -0.41674200 -0.2716308
## improvvalue -0.04180905 -0.4167420 1.00000000 -0.5763071
## area -0.61444657 -0.2716308 -0.57630714 1.0000000
## Retorna um único valor.
## VIF(m0)
## Um valor para cada variável.
summary(m0)
##
## Call:
## lm(formula = saleprice ~ landvalue + improvvalue + area, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.662 -2.155 1.027 2.722 15.976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4703 5.7463 0.256 0.80132
## landvalue 0.8145 0.5122 1.590 0.13137
## improvvalue 0.8204 0.2112 3.885 0.00131 **
## area 13.5286 6.5857 2.054 0.05666 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.919 on 16 degrees of freedom
## Multiple R-squared: 0.8974, Adjusted R-squared: 0.8782
## F-statistic: 46.66 on 3 and 16 DF, p-value: 3.9e-08
vif(m0)
## landvalue improvvalue area
## 2.303501 3.194545 2.850020
## Calculando pelas correlações entre variáveis.
## X <- with(da, cbind(tar, nicotine, weight))
X <- model.matrix(m0)[,-1]
cX <- cor(X)
icX <- solve(cX)
vifs <- diag(icX); vifs
## landvalue improvvalue area
## 2.303501 3.194545 2.850020
vif(m0)
## landvalue improvvalue area
## 2.303501 3.194545 2.850020
## Regra empírica de decisão.
## vif_i ~= 1: variável i não envolvida em multicolinearidade;
## vif_i > 5: variável i envolvida em multicolinearidade, \beta_i com
## variância alta e fracamente estimado.
## Soluções para problemas de colinearidade:
## * Remover variáveis colineares. Remover as de menor significado,
## sujeitas a maior erro de medida ou mais caras/demoraras de
## adquirir.
## * Centralizar as variáveis.
## * Considerar regressão Ridge. Ver pacotes "ridge" e "genridge".
##-----------------------------------------------------------------------------
## Log verossimilhança.
## Quando o y~Gaussiano, após o ajuste pode-se calcular a
## log-verossimilhança por
n <- length(fitted(m0))
SQR <- deviance(m0); s2 <- SQR/n
ll <- -(n/2)*log(2*pi)-(n/2)*log(s2)-(n/2); ll
## [1] -67.53385
## logLik() retorna o valor de log-verossimilhança.
logLik(m0)
## 'log Lik.' -67.53385 (df=5)
## AIC (p é o total de parametros no modelo, \beta+\sigma^2).
p <- length(coef(m0))
pp <- p+1
-2*ll+2*pp
## [1] 145.0677
AIC(m0, k=2)
## [1] 145.0677
## BIC.
AIC(m0, k=log(n))
## [1] 150.0464
BIC(m0)
## [1] 150.0464
-2*ll+log(n)*pp
## [1] 150.0464
##-----------------------------------------------------------------------------
## R².
## Correlação entre observado e ajustado.
yobs <- m0$model[,1]
cor(yobs, fitted(m0))^2
## [1] 0.8974268
## Porção de soma de quadrados total atribuída ao modelo.
1-SQR/(sum(scale(y, scale=FALSE)^2))
## [1] 0.8974268
## R² ajustado para o número de parâmetros.
1-((n-1)/(n-p))*(1-summary(m0)$r.squared)
## [1] 0.8781943
summary(m0)$r.squared
## [1] 0.8974268
summary(m0)$adj.r.squared
## [1] 0.8781943
## O R² TAMBÉM É UMA MEDIDA RELATIVA DE AJUSTE! Por ser uma medida
## limitada em 0 e 1, acostumou-se a pensar que valores próximos de 1
## indicam um modelo adequadamente ajustado. Uma coisa não implica na
## na outra. Sempre deve ser feita uma análise de resíduos independente
## do valor do R². Ver conjunto de dados ascombe.
panel.regs <- function(x, y, ...){
panel.xyplot(x, y, ...)
m0 <- lm(y~x)
c0 <- c(coef(m0), summary(m0)$r.squared)
l0 <- sprintf("b0: %0.2f\n b1: %0.2f\n R²: %0.2f",
c0[1], c0[2], c0[3])
grid.text(x=0.95, y=0.05, label=l0, hjust=1, vjust=0)
}
xyplot.list(list("Linear"=y1~x1,
"Lack of fit"=y2~x2,
"Outlier"=y3~x3,
"Leverage"=y4~x4),
type=c("p","r"), panel=panel.regs,
xlab="x", ylab="y", layout=c(2,2), as.table=TRUE,
data=anscombe, x.same=TRUE, y.same=TRUE)
O PRESS (prediction residual error sum of squares) é uma medida ligada a capacidade preditiva do modelo dada por meio de estatísticas leave-one-out,
\[ \begin{aligned} \text{PRESS} = \sum (y_i-\hat{y}_{i(-i)})^2 = \sum \left(\frac{e_i}{1-h_i} \right)^2. \end{aligned} \]
##-----------------------------------------------------------------------------
## PRESS: predicted residual sum of squares.
## browseURL("http://en.wikipedia.org/wiki/PRESS_statistic")
sum((residuals(m0)/(1-hatvalues(m0)))^2)
## [1] 1549.209
press(m0)
## [1] 1549.209
## Um resumo com todas as medidas.
lm.select(list(m0))
## Model AIC AICc BIC Cp PRESS
## 1 saleprice ~ landvalue + improvvalue + area 145.0677 149.3534 150.0464 6 1549.209
##-----------------------------------------------------------------------------
str(swiss)
pairs(swiss)
splom(swiss, type=c("p","r"))
help(swiss, help_type="html")
## m0 <- lm(Fertility~Agriculture+Examination+Education+Catholic+Infant.Mortality, data=swiss)
m0 <- lm(Fertility~., data=swiss)
par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)
vif(m0)
m1 <- update(m0, .~.-Examination)
summary(m1)
anova(m0, m1)
##-----------------------------------------------------------------------------
prelink <- "http://www.leg.ufpr.br/~walmes/data/business_economics_dataset"
da <- read.table(paste(prelink, "/EXAMPLES/FTC.DAT", sep=""),
header=FALSE)
str(da)
names(da) <- c("tar", "nicotine", "weight", "co")
str(da)
## tar: conteúdo de alcatrão;
## nicotine: conteúdo de nicotina;
## weight: peso;
## co: monoxido de carbono;
## Os valores no data.frame são dos valores de alcatrão, nicotina e
## monoxido de carbono (mg) e peso (g) para uma amostra de 25 marcas de
## filtros testados. Deseja-se modelar o monoxido de carbono como função
## das demais variáveis.
m0 <- lm(co~tar+nicotine+weight, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)
residualPlots(m0)
im <- influence.measures(m0)
summary(im)
m1 <- update(m0, data=da[-3,])
par(mfrow=c(2,2)); plot(m1); layout(1)
summary(m1)
m2 <- update(m1, .~tar)
anova(m2, m1)
vif(m1)
summary(m2)
##-----------------------------------------------------------------------------
## Dados de gasto com alimentação (foodcons) em função da receita da
## família (income) e do número de membros (size).
##
## Ajuste um modelo com foodcons~income+size. Faça uma análise
## completa do modelo e verifique se há necessidade de modificações. Se
## sim, proceda e justifique.
url <- "http://www.leg.ufpr.br/~walmes/data/
business_economics_dataset/EXERCISE/DCFOOD.DAT"
da <- read.table(paste0(strwrap(url), collapse=""), header=FALSE)
names(da) <- c("house","foodcons","income","size")
str(da)
##-----------------------------------------------------------------------------
## Engenheiros mediram o comprimento (length, cm) e o peso (weight, g) e
## o nível de DDT (ddt, ppm) para 144 peixes capturados. Além do mais, a
## distância de captura rio acima (mile), o rio (river) e a espécie do
## peixe (species) também foram registrados.
##
## Ajuste um modelo com ddt~mile+length+weight. Faça uma análise
## completa do modelo e verifique se há necessidade de modificações. Se
## sim, proceda e justifique.
url <- "http://www.leg.ufpr.br/~walmes/data/
business_economics_dataset/EXERCISE/DDT.DAT"
da <- read.table(paste0(strwrap(url), collapse=""), header=FALSE)
names(da) <- c("river","mile","species","length","weight","ddt")
str(da)
## http://westfall.ba.ttu.edu/isqs5349/Rdata/turtles.txt
##-----------------------------------------------------------------------------
## Informações da sessão.
Sys.time()
## [1] "2014-12-11 20:54:12 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] tcltk grid stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.4 plyr_1.8.1 asbio_1.1-1 alr3_2.0.5
## [5] car_2.0-22 gridExtra_0.9.1 latticeExtra_0.6-26 RColorBrewer_1.0-5
## [9] lattice_0.20-29 rmarkdown_0.3.3 knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] deSolve_1.11 digest_0.6.4 doBy_4.5-12 evaluate_0.5.5
## [5] formatR_1.0 htmltools_0.2.6 MASS_7.3-35 Matrix_1.1-4
## [9] methods_3.1.2 multcomp_1.3-7 multcompView_0.1-5 mvtnorm_1.0-1
## [13] nnet_7.3-8 pixmap_0.4-11 plotrix_3.5-10 Rcpp_0.11.3
## [17] sandwich_2.3-2 scatterplot3d_0.3-35 splines_3.1.2 stringr_0.6.2
## [21] survival_2.37-7 TH.data_1.0-5 tkrplot_0.0-23 tools_3.1.2
## [25] yaml_2.1.13 zoo_1.7-11