Exercícios de Delineamento Inteiramente Casualizado

Prof. Fernando de Souza Bastos

26 de setembro de 2018

Exercício 4.1

Para comparar a produtividade de quatro variedades de milho, um agrônomo tomou vinte parcelas similares e distribuiu, inteiramente ao acaso, cada uma das 4 variedades em 5 parcelas experimentais. A partir dos dados experimentais fornecidos abaixo, é possível concluir que existe diferença significativa entre as variedades com relação a produtividade, utilizando o nível de significância de 5%?

\[\begin{array}{ccccc} \hline &&&Variedades&&\\ \hline & A & B & C & D \\ \hline & 25 & 31 & 22 &33 \\ & 26 & 25 & 26 &29 \\ & 20 & 28 & 28 &31 \\ & 23 & 27 & 25 &34 \\ & 21 & 24 & 29 &28 \\ \hline Totais& 115 & 135 & 130 &155 \\ \hline Médias& 23 & 27 & 26 &31 \\ \hline \end{array}\] 
library(readr)
if(!is.null(dev.list())) dev.off()
## null device 
##           1
da <- read_delim("https://raw.githubusercontent.com/maf261/maf261.github.io/master/Exercicios/Exercicios_noR/Exer_4_1.txt", 
                 "\t")
da
## # A tibble: 20 x 2
##    Var    Prod
##    <chr> <int>
##  1 A        25
##  2 A        26
##  3 A        20
##  4 A        23
##  5 A        21
##  6 B        31
##  7 B        25
##  8 B        28
##  9 B        27
## 10 B        24
## 11 C        22
## 12 C        26
## 13 C        28
## 14 C        25
## 15 C        29
## 16 D        33
## 17 D        29
## 18 D        31
## 19 D        34
## 20 D        28
str(da)
## Classes 'tbl_df', 'tbl' and 'data.frame':    20 obs. of  2 variables:
##  $ Var : chr  "A" "A" "A" "A" ...
##  $ Prod: int  25 26 20 23 21 31 25 28 27 24 ...
##  - attr(*, "spec")=List of 2
##   ..$ cols   :List of 2
##   .. ..$ Var : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   .. ..$ Prod: list()
##   .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
##   ..$ default: list()
##   .. ..- attr(*, "class")= chr  "collector_guess" "collector"
##   ..- attr(*, "class")= chr "col_spec"
da$Var <- as.factor(da$Var)
is.factor(da$Var)
## [1] TRUE
#Resumo total dos dados
summary(da)
##  Var        Prod      
##  A:5   Min.   :20.00  
##  B:5   1st Qu.:24.75  
##  C:5   Median :26.50  
##  D:5   Mean   :26.75  
##        3rd Qu.:29.00  
##        Max.   :34.00
#Média dos dados por variedade
res <- aggregate(Prod ~ Var, data = da, FUN = mean);
kable(res,format = "markdown",digits = 2, padding = 3,align = "c")
Var Prod
A 23
B 27
C 26
D 31

Sempre é interessante fazer uma análise descritiva dos dados para avaliar linearidade, independencia, normalidade e homocedásticidade.

O pacote lattice de Sarkar (2017) do R é valioso para tais análises:

library(lattice)

xyplot(Prod~Var, data=da, jitter.x=TRUE)

xyplot(Prod~reorder(Var,Prod), data=da, jitter.x=TRUE)

dotplot(Prod~Var, data=da, jitter.x=TRUE)

dotplot(Prod~reorder(Var,Prod), data=da, jitter.x=TRUE)

bwplot(Prod~Var, data=da, pch="|")

bwplot(Prod~reorder(Var,Prod), data=da, pch="|")

O gráfico dos resíduos versus valores ajustados (valores preditos) é uma das principais técnicas utilizadas para verificar as suposições dos resíduos. Além da detecção de heteroscedasticidade, esse gráfico pode indicar que não existe uma relação linear entra as variáveis explicativas com a variável resposta por meio de alguma tendência nos pontos. Por exemplo, se os pontos do gráfico formam uma parábola, é indicativo que termos de segundo grau sejam necessários.

Para o diagnóstico de heteroscedasticidade, tentamos encontrar alguma tendência no gráfico. Por isso, se os pontos estão aleatoriamente distribuídos em torno do 0, sem nenhum comportamento ou tendência, temos indícios de que a variância dos resíduos é homoscedástica. Já a presença de “funil” é um indicativo da presença de heteroscedasticidade.

Para obtenção dos resíduos padronizados fazemos uma mudança de escala dos resíduos, tend uma linha horizontal temos que a variabilidade é constante

Pontos ao redor da linha significam normalidade.

Mesma avaliação do plot 1, mais utilizado quando utilizamos regressão linear

m0 <- lm(Prod~Var, data=da)
par(mfrow=c(2,2))
plot(m0)

layout(1)

\[ \begin{align} H_{0}&:\mu_{1}=\mu_{2}=\mu_{3}=\mu_{4}\\ H_{1}&: \textrm{Não}\ H_{0}\ (\textrm{Pelo menos uma média difere das demais}) \end{align} \]

alpha <- 0.05
#Variedades: A, B, C e D
t1 <- c(25,26,20,23,21)
t2 <- c(31,25,28,27,24)
t3 <- c(22,26,28,25,29)
t4 <- c(33,29,31,34,28)
tab <- matrix(c(t1,t2,t3,t4),5,4)
(n <- length(tab))
## [1] 20
(I <- dim(tab)[2])
## [1] 4
(J <- dim(tab)[1])
## [1] 5
(bart1 <- mean(t1))
## [1] 23
(bart2 <- mean(t2))
## [1] 27
(bart3 <- mean(t3))
## [1] 26
(bart4 <- mean(t4))
## [1] 31
(Tott1 <- sum(t1))
## [1] 115
(Tott2 <- sum(t2))
## [1] 135
(Tott3 <- sum(t3))
## [1] 130
(Tott4 <- sum(t4))
## [1] 155
(mu <- mean(tab))
## [1] 26.75
####################################
difTot <- matrix(NA,dim(tab)[1],dim(tab)[2])
for (j in 1:dim(tab)[2]) {
  for (i in 1:dim(tab)[1]) {
    difTot[i,j] <- ((tab[i,j]-mu)^2)
  }
  
  }
(SQTot <- sum(difTot))
## [1] 275.75
####################################
difTrat <- matrix(NA,1,I)
for (i in 1:I) {
  difTrat[i] <- ((mean(tab[,i])-mu)^2)
}
(SQTrat <- 5*sum(difTrat))
## [1] 163.75
####################################
difRes <- matrix(NA,dim(tab)[1],dim(tab)[2])
for (i in 1:dim(tab)[2]) {
  for (j in 1:dim(tab)[1]) {
    difRes[j,i] <- (tab[j,i]-mean(tab[,i]))^2
     }
}
(SQRes <- sum(difRes))
## [1] 112
####################################
#Observem que:

#Soma de quadrado totais = Soma de quadrados dos resíduos + Soma de quadrados de tratamentos
FV <- cbind("Trat","Res","Total")
GL <- cbind(I-1,I*(J-1),I*J-1)
SQ <- cbind(SQTrat,SQRes,SQTot)
QM <- cbind(SQTrat/GL[1],SQRes/GL[2],SQTot/GL[3])
(Fcal <- QM[1,1]/QM[1,2])
## [1] 7.797619
(Ftab <- qf(alpha,GL[1],GL[2],lower.tail = FALSE))
## [1] 3.238872
(pvalor <- pf(Fcal,GL[1],GL[2],lower.tail = FALSE))
## [1] 0.001975594
RR <- "Rejeita-se H0 ao nível alpha de significancia"
RNR <- "Não rejeita-se H0 ao nível alpha de significancia"
dif <- "Pelo menos uma média difere das demais"
ndif <- "Todas as médias são iguais"
ifelse(Fcal>Ftab,RR,RNR)
## [1] "Rejeita-se H0 ao nível alpha de significancia"
ifelse(Fcal>Ftab,dif,ndif)
## [1] "Pelo menos uma média difere das demais"
m0 <- lm(Prod~Var, data=da)
anova(m0)
## Analysis of Variance Table
## 
## Response: Prod
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## Var        3 163.75  54.583  7.7976 0.001976 **
## Residuals 16 112.00   7.000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(ExpDes.pt)
dic <- dic(da$Var,
  da$Prod,
  quali = TRUE,
  mcomp = "tukey",
  nl = FALSE,
  hvar = 'bartlett',
  sigT = 0.05,
  sigF = 0.05)
## ------------------------------------------------------------------------
## Quadro da analise de variancia
## ------------------------------------------------------------------------
##            GL     SQ     QM     Fc     Pr>Fc
## Tratamento  3 163.75 54.583 7.7976 0.0019756
## Residuo    16 112.00  7.000                 
## Total      19 275.75                        
## ------------------------------------------------------------------------
## CV = 9.89 %
## 
## ------------------------------------------------------------------------
## Teste de normalidade dos residuos 
## Valor-p:  0.2358736 
## De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Teste de homogeneidade de variancia 
## valor-p:  0.9981235 
## De acordo com o teste de bartlett a 5% de significancia, as variancias podem ser consideradas homogeneas.
## ------------------------------------------------------------------------
## 
## Teste de Tukey
## ------------------------------------------------------------------------
## Grupos Tratamentos Medias
## a     D   31 
## ab    B   27 
##  b    C   26 
##  b    A   23 
## ------------------------------------------------------------------------
#Gráficos de Resíduos do Pacote ExpDes.pt
plotres(dic)

Exercício 4.2

Um treinador de corrida rústica, objetivando melhorar o desempenho de seus atletas, testou três novas técnicas de preparação. Para tanto trabalhou com um grupo de 15 atletas completamente homogêneos para as características essenciais. A designação das técnicas de preparação aos atletas foi feita totalmente ao acaso e de tal forma que o número de atletas avaliados em cada uma das técnicas fosse o mesmo. Os resultados obtidos, após um determinado período de tempo de aprendizado da técnica pelos atletas, foram os seguintes (minutos / 25 Km):

\[\begin{array}{ccccc} \hline &&\textrm{Técnicas de Preparação}&&\\ \hline \textrm{Repetições} & 1 & 2 & 3 \\ \hline 1 & 130 & 125 & 135 \\ 2 & 129 & 131 & 129 \\ 3 & 128 & 130 & 131 \\ 4 & 126 & 129 & 128 \\ 5 & 130 & 127 & 130 \\ \hline Totais& 643 & 642 & 653 \\ \hline \textrm{Médias}& 128.6 & 128.4 & 130.6 \\ \hline \end{array}\] 
library(readr)
da <- read_delim("https://raw.githubusercontent.com/maf261/maf261.github.io/master/Exercicios/Exercicios_noR/Exer_4_2.txt", 
                 "\t")
da
## # A tibble: 15 x 2
##    Tec   Tempo
##    <chr> <int>
##  1 X1      130
##  2 X1      129
##  3 X1      128
##  4 X1      126
##  5 X1      130
##  6 X2      125
##  7 X2      131
##  8 X2      130
##  9 X2      129
## 10 X2      127
## 11 X3      135
## 12 X3      129
## 13 X3      131
## 14 X3      128
## 15 X3      130
str(da)
## Classes 'tbl_df', 'tbl' and 'data.frame':    15 obs. of  2 variables:
##  $ Tec  : chr  "X1" "X1" "X1" "X1" ...
##  $ Tempo: int  130 129 128 126 130 125 131 130 129 127 ...
##  - attr(*, "spec")=List of 2
##   ..$ cols   :List of 2
##   .. ..$ Tec  : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   .. ..$ Tempo: list()
##   .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
##   ..$ default: list()
##   .. ..- attr(*, "class")= chr  "collector_guess" "collector"
##   ..- attr(*, "class")= chr "col_spec"
da$Tec <- as.factor(da$Tec)
is.factor(da$Tec)
## [1] TRUE
#Resumo total dos dados
summary(da)
##  Tec        Tempo      
##  X1:5   Min.   :125.0  
##  X2:5   1st Qu.:128.0  
##  X3:5   Median :129.0  
##         Mean   :129.2  
##         3rd Qu.:130.0  
##         Max.   :135.0
#Média dos dados por variedade
res <- aggregate(Tempo ~ Tec, data = da, FUN = mean);
kable(res,format = "markdown",digits = 2, padding = 3,align = "c")
Tec Tempo
X1 128.6
X2 128.4
X3 130.6

Usando o pacote lattice de Sarkar (2017):

library(lattice)

xyplot(Tempo~Tec, data=da, jitter.x=TRUE)

xyplot(Tempo~reorder(Tec,Tempo), data=da, jitter.x=TRUE)

dotplot(Tempo~Tec, data=da, jitter.x=TRUE)

dotplot(Tempo~reorder(Tec,Tempo), data=da, jitter.x=TRUE)

bwplot(Tempo~Tec, data=da, pch="|")

bwplot(Tempo~reorder(Tec,Tempo), data=da, pch="|")

m0 <- lm(Tempo~Tec, data=da)
par(mfrow=c(2,2))
plot(m0)

layout(1)

\[ \begin{align} H_{0}&:\mu_{1}=\mu_{2}=\mu_{3}\\ H_{1}&: \textrm{Não}\ H_{0}\ (\textrm{Pelo menos uma média difere das demais}) \end{align} \] - Usando a função interna do R

m0 <- lm(Tempo~Tec, data=da)
anova(m0)
## Analysis of Variance Table
## 
## Response: Tempo
##           Df Sum Sq Mean Sq F value Pr(>F)
## Tec        2   14.8     7.4  1.3962  0.285
## Residuals 12   63.6     5.3
library(ExpDes.pt)
dic <- dic(da$Tec,
  da$Tempo,
  quali = TRUE,
  mcomp = "tukey",
  nl = FALSE,
  hvar = 'bartlett',
  sigT = 0.01,
  sigF = 0.01)
## ------------------------------------------------------------------------
## Quadro da analise de variancia
## ------------------------------------------------------------------------
##            GL   SQ  QM     Fc Pr>Fc
## Tratamento  2 14.8 7.4 1.3962 0.285
## Residuo    12 63.6 5.3             
## Total      14 78.4                 
## ------------------------------------------------------------------------
## CV = 1.78 %
## 
## ------------------------------------------------------------------------
## Teste de normalidade dos residuos 
## Valor-p:  0.9262885 
## De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Teste de homogeneidade de variancia 
## valor-p:  0.663661 
## De acordo com o teste de bartlett a 5% de significancia, as variancias podem ser consideradas homogeneas.
## ------------------------------------------------------------------------
## 
## De acordo com o teste F, as medias nao podem ser consideradas diferentes.
## ------------------------------------------------------------------------
##   Niveis Medias
## 1     X1  128.6
## 2     X2  128.4
## 3     X3  130.6
## ------------------------------------------------------------------------
plotres(dic)

Exercício do livro do Banzato (Exer 3.2.1)

Agora vamos usar o pacote Labest data produzido pela UFPR (2016)

# library(devtools)
# 
# # Do repositório de divulgação no GitHub.
# install_github(repo = "labestData",
#                username = "pet-estatistica",
#                ref = "master", build_vignettes = TRUE)

library(labestData)
#labestDataView()
str(BanzattoQd3.2.1)
## 'data.frame':    30 obs. of  3 variables:
##  $ trat   : Factor w/ 5 levels "Azinfos etilico",..: 5 1 3 4 2 5 1 3 4 2 ...
##  $ rept   : int  1 1 1 1 1 2 2 2 2 2 ...
##  $ pulgoes: int  2370 1282 562 173 193 1687 1527 321 127 71 ...
library(lattice)

data(BanzattoQd3.2.1)

aggregate(pulgoes ~ trat,  data = BanzattoQd3.2.1,
          FUN = function(x) { c(mean = mean(x), var = var(x)) })
##                   trat pulgoes.mean  pulgoes.var
## 1      Azinfos etilico   1075.00000  75558.80000
## 2        Diazinon 60CE     91.33333   2791.86667
## 3 Supracid 40CE dose 1    527.16667  40126.16667
## 4 Supracid 40CE dose 2    156.33333   1502.26667
## 5           Testemunha   2477.00000 233749.60000
ban <- BanzattoQd3.2.1

xyplot(pulgoes ~ trat, data = ban,
       xlab = "Tratamentos",
       ylab = "Número de pulgões 36h após pulverização")

xyplot(pulgoes ~ reorder(trat,pulgoes), data = ban,
       xlab = "Tratamentos",
       ylab = "Número de pulgões 36h após pulverização")

m0 <- lm(pulgoes~trat,data = ban)
plot(m0)

#Suponha que esteja tudo bem!
library(ExpDes.pt)
dic <- dic(ban$trat,
  ban$pulgoes,
  quali = TRUE,
  mcomp = "tukey",
  nl = FALSE,
  hvar = 'bartlett',
  sigT = 0.05,
  sigF = 0.05)
## ------------------------------------------------------------------------
## Quadro da analise de variancia
## ------------------------------------------------------------------------
##            GL       SQ      QM    Fc      Pr>Fc
## Tratamento  4 23145259 5786315 81.79 5.5054e-14
## Residuo    25  1768644   70746                 
## Total      29 24913903                         
## ------------------------------------------------------------------------
## CV = 30.74 %
## 
## ------------------------------------------------------------------------
## Teste de normalidade dos residuos 
## Valor-p:  0.006742331 
## ATENCAO: a 5% de significancia, os residuos nao podem ser considerados normais!
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Teste de homogeneidade de variancia 
## valor-p:  0.9616506 
## De acordo com o teste de bartlett a 5% de significancia, as variancias podem ser consideradas homogeneas.
## ------------------------------------------------------------------------
## 
## Teste de Tukey
## ------------------------------------------------------------------------
## Grupos Tratamentos Medias
## a     Testemunha      2477 
##  b    Azinfos etilico     1075 
##   c   Supracid 40CE dose 1    527.1667 
##   c   Supracid 40CE dose 2    156.3333 
##   c   Diazinon 60CE   91.33333 
## ------------------------------------------------------------------------
#Avalie os resultados
plotres(dic)

\[ y_{ij}=\Biggl\{ \begin{array}{c} \dfrac{y^{\lambda}-1}{\lambda},\ \textrm{se}\ \lambda\neq 0\\ \log{(y)},\ \textrm{se}\ \lambda=0 \end{array} \] O valor de \(\lambda\) será o valor que maxima a função log de verossimilhança (medida de compatibilidade) do modelo m0. Não precisa ser o valor exato, escolha um ponto aproximado entre as linhas mais externas abaixo:

Escolha um número que seja próximo, o ideal é que seja um número do tipo \(1/a,\) onde \(a \in \mathbb{Z}.\)

MASS::boxcox(m0)

MASS::boxcox(m0)
abline(v=0.2,col="red")
abline(v=0.25,col="blue")

m1 <- lm(pulgoes^0.2~trat,data = ban)
plot(m1)

#Suponha que esteja tudo bem!
library(ExpDes.pt)
dic <- dic(ban$trat,
  ban$pulgoes^0.2,
  quali = TRUE,
  mcomp = "tukey",
  nl = FALSE,
  hvar = 'bartlett',
  sigT = 0.05,
  sigF = 0.05)
## ------------------------------------------------------------------------
## Quadro da analise de variancia
## ------------------------------------------------------------------------
##            GL      SQ     QM     Fc      Pr>Fc
## Tratamento  4 21.6288 5.4072 118.82 6.9784e-16
## Residuo    25  1.1377 0.0455                  
## Total      29 22.7665                         
## ------------------------------------------------------------------------
## CV = 6.13 %
## 
## ------------------------------------------------------------------------
## Teste de normalidade dos residuos 
## Valor-p:  0.6206532 
## De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
## ------------------------------------------------------------------------
## 
## ------------------------------------------------------------------------
## Teste de homogeneidade de variancia 
## valor-p:  0.9824495 
## De acordo com o teste de bartlett a 5% de significancia, as variancias podem ser consideradas homogeneas.
## ------------------------------------------------------------------------
## 
## Teste de Tukey
## ------------------------------------------------------------------------
## Grupos Tratamentos Medias
## a     Testemunha      4.759716 
##  b    Azinfos etilico     4.022682 
##   c   Supracid 40CE dose 1    3.4685 
##    d      Supracid 40CE dose 2    2.736617 
##    d      Diazinon 60CE   2.421748 
## ------------------------------------------------------------------------
#Avalie os resultados
plotres(dic)

Referências

Ferreira, Eric Batista, Portya Piscitelli Cavalcanti, and Denismar Alves Nogueira. 2018. ExpDes.pt: Pacote Experimental Designs (Portuguese). https://CRAN.R-project.org/package=ExpDes.pt.

Sarkar, Deepayan. 2017. “Lattice: Trellis Graphics for R.” março. https://cran.r-project.org/web/packages/lattice/index.html.

UFPR, PET Estatística. 2016. LabestData: Biblioteca de Dados Para Ensino de Estatística. http://gitlab.c3sl.ufpr.br/pet-estatistica/labestData, https://github.com/pet-estatistica/labestData.