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


Medidas de diagnóstico e influência

http://web.stanford.edu/class/stats191/notebooks/Diagnostics%20for%20multiple%20regression.pdf


Diagnóstico visual

##=============================================================================
## 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 de forma matricial

##-----------------------------------------------------------------------------
## 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)

Cálculo das medias de diagnóstico e influência

Pontos de alavancagem

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

Resíduos crus (ordinários)

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

Resíduos padronizados

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

Resíduos studentizados

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)

Estatísticas baseadas em leave-one-out

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

Distância de Cook

\[ \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

\[ 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

DFbetas

\[ \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

##-----------------------------------------------------------------------------
## 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

Medida de colinearidade

##-----------------------------------------------------------------------------
## 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".

Medidas de ajuste

Log-verossimilhança, AIC e BIC

  • A log-verossimilhança para distribuição gaussiana é dada por \[ \begin{aligned} ll &= -\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\hat{\sigma}^2)-\frac{1}{2\hat{\sigma}^2}\sum (y_i-\hat{y}_i)^2\\ &= -\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\hat{\sigma}^2)-\frac{n}{2}, \end{aligned} \] em que \(\hat{\sigma}^2 = (1/n)\sum (y_i-\hat{y}_i)^2\).
  • O critério de informação de Akaike, AIC, é dado por \[ AIC = -2\cdot ll +2\cdot p, \] em que \(p\) é o total de parâmetros estimados no modelo, ou seja, incluí aqueles parâmetros que representam dispersão/covariância também
  • O BIC é dado por \[ BIC = -2\cdot ll +\log(n)\cdot p. \]
##-----------------------------------------------------------------------------
## 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

Coeficientes de determinação

  • O coeficiente de determinação é dados por \[ \begin{aligned} R^2 &= \text{cor}(y, \hat{y})^2\\ &= 1-\frac{\sum (y_i-\hat{y}_i)^2}{\sum (y_i-\bar{y}_i)^2}\\ &= 1-\frac{SQR}{SQTc}, \end{aligned} \] em que \(SQR\) e \(SQTc\) representam a soma de quadrados dos resíduos e a soma de quadrados total corrigida para a média (modelo com apenas intercepto).
  • O coeficiente de determinação ajustado é dado por \[ \begin{aligned} R_{adj}^2 = 1-\frac{n-1}{n-p}(1-R^2). \end{aligned} \]
##-----------------------------------------------------------------------------
## 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)

PRESS

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

Praticar

##-----------------------------------------------------------------------------

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