boxcox.R

library(MASS)
library(car)

## Loading required package: carData

source('https://dl.dropbox.com/s/bhy3tk9dalv2pyt/tnormal.R',
       encoding = 'UTF-8',
       echo = F)

source('https://dl.dropbox.com/s/ms9sn1ge4a3k2co/normal.R',
       encoding = 'UTF-8',
       echo = F)

av.b <- aov(Days ~ Eth*Sex*Age*Lrn, data = quine)

t.normal(resid(av.b))

## Loading required package: nortest

## 
##  Média amostral = 1.690524e-16 
##  Desvio padrão amostral = 12.7334

## Warning in ks.test(x, "pnorm", xb, sx): ties should not be present for the
## Kolmogorov-Smirnov test

##                                                Estatística         p
## One-sample Kolmogorov-Smirnov test                 0.08832 0.2048393
## Lilliefors (Kolmogorov-Smirnov) normality test     0.08832 0.0072630
## Cramer-von Mises normality test                    0.28178 0.0005288
## Shapiro-Wilk normality test                        0.96122 0.0003913
## Shapiro-Francia normality test                     0.95915 0.0004865
## Anderson-Darling normality test                    1.63276 0.0003267

par(mfrow=c(2,2)) 
plot(av.b)

## Warning: not plotting observations with leverage one:
##   32, 60, 105, 136

## Warning: not plotting observations with leverage one:
##   32, 60, 105, 136

par(mfrow=c(1,1)) 
qqp(resid(av.b))

## [1]  59 104

box.tr <- boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda =  seq(-1, 1, 1/10))

(lambda <- box.tr$x[which(box.tr$y == max(box.tr$y))])

## [1] 0.2121212

av.c <- aov((Days+1)^(lambda) ~ Eth*Sex*Age*Lrn, data = quine)

t.normal(resid(av.c))

## 
##  Média amostral = -8.272982e-18 
##  Desvio padrão amostral = 0.2791861

## Warning in ks.test(x, "pnorm", xb, sx): ties should not be present for the
## Kolmogorov-Smirnov test

##                                                Estatística       p
## One-sample Kolmogorov-Smirnov test                 0.06849 0.49995
## Lilliefors (Kolmogorov-Smirnov) normality test     0.06849 0.09043
## Cramer-von Mises normality test                    0.06054 0.36989
## Shapiro-Wilk normality test                        0.98908 0.31201
## Shapiro-Francia normality test                     0.98992 0.32121
## Anderson-Darling normality test                    0.37696 0.40579

par(mfrow=c(2,2)) 
plot(av.c)

## Warning: not plotting observations with leverage one:
##   32, 60, 105, 136

## Warning: not plotting observations with leverage one:
##   32, 60, 105, 136

par(mfrow=c(1,1)) 
qqp(resid(av.c))

## [1] 98 61