Transformação de dados

Exercícios

1. O artigo “Origin of precambrian iron formations” (Econ. Geology, 1964:1025-1057) relata os seguintes dados sobre o total de Fe de quatro tipos de formação de ferro (t1=carbonato, t2=silicato, t3=magnetita, t4=hematita).

t1	t2	t3	t4
20.5	26.3	29.5	36.5
28.1	24.0	34.0	44.2
27.8	26.2	27.5	34.1
27.0	20.2	29.4	30.3
28.0	23.7	27.9	31.4
25.2	34.0	26.2	33.1
25.3	17.1	29.9	34.1
27.1	26.8	29.5	32.9
20.5	23.7	30.0	36.3
31.3	24.9	35.6	25.5

Após uma análise de resíduos, foi constatado desvio de normalidade, violando deste modo um dos pressupostos da anova. Logo ,uma transformação se faz necessário.

1. Utilize as trasformações raiz quadrada, logarítmica e boxcox e apresente os dados.
Resposta

Para facilitar a vida do usuário, irei fazer uma função que faz a transformação boxcox “quase” que automaticamente. Em alguns casos pode ser necessário alterar o intervalo de lambda para que o lambda estimado seja o mais preciso possível. Logo, aconselho o uso da função acompanhado do gráfico para que não haja erro. Segue a função:

library(MASS)
bx <- function(x, data, lambda = seq(-2,2,1/10),plotit=TRUE,...){ #um objeto da classe lm, aov
  if(!require(MASS)){
   install.packages('MASS')
  }

 aux <- MASS::boxcox(object = x,data = data, lambda = lambda,plotit=plotit,...)
 aux2 <- do.call('cbind',aux)
 aux3 <- which(aux2 == max(aux2[,2]), arr.ind=TRUE)
 aux4 <- aux2[aux3[1,1],]
 lambda <- aux4[1]

 if(class(x) == 'formula'){
 
    oldvar <- data[[all.vars(x)[1]]] 
 
 } else {

   oldvar <- x$model$res 
 
 }
 newvar <- (oldvar^lambda - 1)/lambda
 res <- list(resp = newvar,
             lambda = lambda)
 return(res)
}

Fazendo as transformações tem-se:

t1 = c(20.5,28.1,27.8,27,28,25.2,25.3,27.1,20.5,31.3)
t2 = c(26.3,24,26.2,20.2,23.7,34,17.1,26.8,23.7,24.9)
t3 = c(29.5,34,27.5,29.4,27.9,26.2,29.9,29.5,30,35.6)
t4 = c(36.5,44.2,34.1,30.3,31.4,33.1,34.1,32.9,36.3,25.5)

dados1 <- data.frame(Tipos = rep(c('t1','t2','t3','t4'),rep(10,4)),
                    resp  = c(t1,t2,t3,t4)) 

yraiz <- sqrt(dados1$resp)
ylog  <- log(dados1$resp)

# No caso da função boxcox, devemos informar o modelo estatístico.
ybox <- bx(x = resp ~ Tipos, data = dados1, plotit=F)[[1]]

# Colocando as novas variáveis na base de dados
dados1$yraiz <- round(yraiz,4)
dados1$ylog  <- round(ylog,4)
dados1$ybox  <- round(ybox,4)

# Logo tem-se os seguintes dados:
dados1

##    Tipos resp  yraiz   ylog    ybox
## 1     t1 20.5 4.5277 3.0204  8.5404
## 2     t1 28.1 5.3009 3.3358 10.6664
## 3     t1 27.8 5.2726 3.3250 10.5873
## 4     t1 27.0 5.1962 3.2958 10.3745
## 5     t1 28.0 5.2915 3.3322 10.6401
## 6     t1 25.2 5.0200 3.2268  9.8862
## 7     t1 25.3 5.0299 3.2308  9.9137
## 8     t1 27.1 5.2058 3.2995 10.4012
## 9     t1 20.5 4.5277 3.0204  8.5404
## 10    t1 31.3 5.5946 3.4436 11.4909
## 11    t2 26.3 5.1284 3.2696 10.1862
## 12    t2 24.0 4.8990 3.1781  9.5529
## 13    t2 26.2 5.1186 3.2658 10.1591
## 14    t2 20.2 4.4944 3.0057  8.4505
## 15    t2 23.7 4.8683 3.1655  9.4685
## 16    t2 34.0 5.8310 3.5264 12.1606
## 17    t2 17.1 4.1352 2.8391  7.4881
## 18    t2 26.8 5.1769 3.2884 10.3209
## 19    t2 23.7 4.8683 3.1655  9.4685
## 20    t2 24.9 4.9900 3.2149  9.8035
## 21    t3 29.5 5.4314 3.3844 11.0315
## 22    t3 34.0 5.8310 3.5264 12.1606
## 23    t3 27.5 5.2440 3.3142 10.5078
## 24    t3 29.4 5.4222 3.3810 11.0057
## 25    t3 27.9 5.2820 3.3286 10.6137
## 26    t3 26.2 5.1186 3.2658 10.1591
## 27    t3 29.9 5.4681 3.3979 11.1346
## 28    t3 29.5 5.4314 3.3844 11.0315
## 29    t3 30.0 5.4772 3.4012 11.1602
## 30    t3 35.6 5.9666 3.5723 12.5474
## 31    t4 36.5 6.0415 3.5973 12.7619
## 32    t4 44.2 6.6483 3.7887 14.5179
## 33    t4 34.1 5.8395 3.5293 12.1850
## 34    t4 30.3 5.5045 3.4111 11.2370
## 35    t4 31.4 5.6036 3.4468 11.5161
## 36    t4 33.1 5.7533 3.4995 11.9398
## 37    t4 34.1 5.8395 3.5293 12.1850
## 38    t4 32.9 5.7359 3.4935 11.8904
## 39    t4 36.3 6.0249 3.5918 12.7145
## 40    t4 25.5 5.0498 3.2387  9.9685

2. Extraia os resíduos da anova considerando a variável na escola original e com as transformações executadas. Faça o gráfico qqplot para cada resíduos colocando-os lado a lado com o intuito de comparação. Todas as transformações resolveram o problema de normalidade?
Resposta

Precisamos fazer uma anova para cada variável transformada. Em seguida vamos extrair os resíduos e fazer os gráficos.

modorig <- aov(resp ~ Tipos, data=dados1)
modraiz <- aov(yraiz ~ Tipos, data=dados1)
modlog  <- aov(ylog ~ Tipos, data=dados1) 
modbox <- aov(ybox ~ Tipos, data=dados1) 

errosorig <- residuals(modorig)
errosraiz <- residuals(modraiz) 
erroslog  <- residuals(modlog) 
errosbox  <- residuals(modbox)

par(mfrow=c(2,2))
qqnorm(errosorig,main='Escala original')
qqline(errosorig)

qqnorm(errosraiz,main='Escala raiz')
qqline(errosraiz)
 
qqnorm(erroslog, main='Escala logarítmica')
qqline(erroslog)
  
qqnorm(errosbox, main='Escala BoxCox')
qqline(errosbox)



Como podemos observar nos gráficos, parece que as transformações não resolveram a normalidade. Vamos apresentar o teste de Lilliefors para todas as variáveis (original e transformada) para que possamos julgar melhor. Segue os p-valores:

##    Original        Raiz         Log         Box 
## 0.028903432 0.013253690 0.004051292 0.016000014

Os pvalores atestam a análise gráfica, ou seja, nenhuma transformação resolveu o problema da normalidade. Neste caso recomenda-se um outro tipo de análise que não faz uso de tal pressuposto como a análise não-paramétrica.

2. Quatro serviços de entrega durante a noite são testados para “fragilidade”, enviando ítens frágeis. As taxas de quebra observadas são dadas abaixo:

A	B	C	D
17	7	11	5
20	11	9	4
15	15	5	3
21	10	12	7
28	10	6	6

A análise de resíduos mostrou que as variâncias são heterocedásticas. Diante disto responda.

1. Utilize as trasformações raiz quadrada, logarítmica e boxcox e apresente os dados.
Resposta

Fazendo as transformações tem-se:

A <- c(17,20,15,21,28)
B <- c(7,11,15,10,10)
C <- c(11,9,5,12,6)
D <- c(5,4,3,7,6)

dados22 <- data.frame(Tipos=rep(LETTERS[1:4],rep(5,4)),
                      resp  =c(A,B,C,D)) 

yraiz <- sqrt(dados22$resp)
ylog  <- log(dados22$resp)

# No caso da função boxcox, devemos informar o modelo estatístico.
ybox <- bx(x = resp ~ Tipos, data = dados22, plotit=F)[[1]]

# Colocando as novas variáveis na base de dados
dados22$yraiz <- round(yraiz,4)
dados22$ylog  <- round(ylog,4)
dados22$ybox  <- round(ybox,4)

# Logo tem-se os seguintes dados:
dados22

##    Tipos resp  yraiz   ylog   ybox
## 1      A   17 4.1231 2.8332 4.4652
## 2      A   20 4.4721 2.9957 4.8549
## 3      A   15 3.8730 2.7081 4.1778
## 4      A   21 4.5826 3.0445 4.9756
## 5      A   28 5.2915 3.3322 5.7245
## 6      B    7 2.6458 1.9459 2.6426
## 7      B   11 3.3166 2.3979 3.5105
## 8      B   15 3.8730 2.7081 4.1778
## 9      B   10 3.1623 2.3026 3.3175
## 10     B   10 3.1623 2.3026 3.3175
## 11     C   11 3.3166 2.3979 3.5105
## 12     C    9 3.0000 2.1972 3.1106
## 13     C    5 2.2361 1.6094 2.0689
## 14     C   12 3.4641 2.4849 3.6915
## 15     C    6 2.4495 1.7918 2.3726
## 16     D    5 2.2361 1.6094 2.0689
## 17     D    4 2.0000 1.3863 1.7191
## 18     D    3 1.7321 1.0986 1.3013
## 19     D    7 2.6458 1.9459 2.6426
## 20     D    6 2.4495 1.7918 2.3726

2. Extraia os resíduos da anova considerando a variável na escola original e com as transformações executadas. Faça o gráfico qqplot para cada resíduo colocando-os lado a lado com o intuito de comparação. As transformações mantiveram a normalidade dos resíduos?
Resposta

Precisamos fazer uma anova para cada variável transformada. Em seguida vamos extrair os resíduos e fazer os gráficos.

modorig <- aov(resp ~ Tipos, data=dados22)
modraiz <- aov(yraiz ~ Tipos, data=dados22)
modlog  <- aov(ylog ~ Tipos, data=dados22) 
modbox <- aov(ybox ~ Tipos, data=dados22) 

errosorig <- residuals(modorig)
errosraiz <- residuals(modraiz) 
erroslog  <- residuals(modlog) 
errosbox  <- residuals(modbox)

par(mfrow=c(2,2))
qqnorm(errosorig,main='Escala original')
qqline(errosorig)

qqnorm(errosraiz,main='Escala raiz')
qqline(errosraiz)
 
qqnorm(erroslog, main='Escala logarítmica')
qqline(erroslog)
  
qqnorm(errosbox, main='Escala BoxCox')
qqline(errosbox)



Como podemos observar nos gráficos, as transformações não interferiram na normalidade dos resíduos. Vamos apresentar o teste de Lilliefors para todas as variáveis (original e transformada) para ratificar o que foi apresentado nos gráficos. Segue os p-valores:

##  Original      Raiz       Log       Box 
## 0.6023316 0.8947168 0.7614876 0.8223737

Os pvalores atestam a análise gráfica.

3. Apresente a análise gráfica referente a homocedasticidade tanto do resíduo na escala original quanto das variáveis transformadas. A transformação resolveu o problema da heterocedasticidade? Apresente o teste F-máximo de Hartley.
Resposta

Os gráficos são:

par(mfrow=c(2,2))
plot(errosorig ~ fitted(modorig),main='Escala original')
abline(h=0,col='red')
plot(errosraiz ~ fitted(modraiz),main='Escala raiz')
abline(h=0,col='red')
plot(erroslog ~ fitted(modlog),main='Escala logarítmica')
abline(h=0,col='red')
plot(errosbox ~ fitted(modbox),main='Escala BoxCox')
abline(h=0,col='red')



Todas as transformações resolveram o problema da heterocedasticidade. Segue o teste F-máximo de Hartley para confirmar a análise gráfica.

varior <- with(dados22,tapply(resp,Tipos,var))
varira <- with(dados22,tapply(yraiz,Tipos,var)) 
varilo <- with(dados22,tapply(ylog,Tipos,var)) 
varibo <- with(dados22,tapply(ybox,Tipos,var))  

fmaxor <- max(varior)/min(varior)
fmaxra <- max(varira)/min(varira) 
fmaxlo <- max(varilo)/min(varilo) 
fmaxbo <- max(varibo)/min(varibo) 

pf(fmaxor,4,4,lower.tail=F)

## [1] 0.02379042

pf(fmaxra,4,4,lower.tail=F)

## [1] 0.2273168

pf(fmaxlo,4,4,lower.tail=F)

## [1] 0.1886888

pf(fmaxbo,4,4,lower.tail=F)

## [1] 0.294552