Funções do R para ANOVA com 1 fator

---

Referências para leitura deste material:

• P. Morettin & W. Bussab. Estatística Básica. 6a ed. São Paulo: Saraiva, 2010.
• M. Magalhães & AC. Lima. Noções de Probabilidade e Estatística Básica. 5a ed. São Paulo: EDUSP, 2002.
• M. Campos. Introdução ao R. Vitória, 2018. Disponível em: https://maurocmcampos.github.io/ec/html/intror.html.
• M. Campos. Funções do R para estatística descritiva. Vitória, 2023. Disponível em: https://maurocmcampos.github.io/ec/html/edR.html.
• M. Campos. Funções do R para teste de hipóteses I. Vitória, 2023. Disponível em: https://maurocmcampos.github.io/ec/html/th1R.html.
• M. Campos. Funções do R para teste de hipóteses III. Vitória, 2023. Disponível em: https://maurocmcampos.github.io/ec/html/th3R.html.

---

1 Comparando duas populações normais independentes

1.1 O problema e os dados

Problema. (Morettin & Bussab (2010), Capítulo 15, Exemplo 15.1) Um psicólogo está investigando a relação entre o tempo que um indivíduo leva para reagir a um estímulo visual (\(y\)) e alguns fatores tais como: sexo (\(x_1\)), idade (\(x_2\)) e acuidade visual (\(x_3\), medida em porcentagem). Nesta seção, desejamos desenvolver um modelo descrevendo a relacão entre o tempo de reação a um estímulo visual e a variável sexo, que é um fator com dois níveis: sexo masculino (\(i=1\) e \(n_{1}=10\)) e sexo feminino (\(i=2\) e \(n_{2}=10\)). O problema fundamental em questão é o seguinte: será que o conhecimento do sexo de um indivíduo ajuda a melhorar a previsão do seu tempo de reação a um estímulo visual?

Dados. Os dados para este problema são como segue.

tempo <- c(96,92,106,100,98,104,110,101,116,106,109,100,112,105,118,108,113,112,127,117)
sexo <- factor(c("M","F","M","F","F","M","M","F","F","M","M","F","F","F","M","M","F","F","M","M"))
levels(sexo)=c("F","M")
# Para obter exatamente os mesmos resultados do livro-texto, 
# precisamos fixar o nivel M como o nivel de referencia
sexo <- relevel(sexo,ref="M") 
idade <- factor(rep(c("20","25","30","35","40"),c(4,4,4,4,4)))
levels=c("20","25","30","35","40")
acuidade <- c(90,100,80,90,100,90,80,90,70,90,90,80,90,80,70,90,90,90,60,80)
df1 <- data.frame(tempo,sexo,idade,acuidade)
str(df1)
## 'data.frame':    20 obs. of  4 variables:
##  $ tempo   : num  96 92 106 100 98 104 110 101 116 106 ...
##  $ sexo    : Factor w/ 2 levels "M","F": 1 2 1 2 2 1 1 2 2 1 ...
##  $ idade   : Factor w/ 5 levels "20","25","30",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ acuidade: num  90 100 80 90 100 90 80 90 70 90 ...
df1
##    tempo sexo idade acuidade
## 1     96    M    20       90
## 2     92    F    20      100
## 3    106    M    20       80
## 4    100    F    20       90
## 5     98    F    25      100
## 6    104    M    25       90
## 7    110    M    25       80
## 8    101    F    25       90
## 9    116    F    30       70
## 10   106    M    30       90
## 11   109    M    30       90
## 12   100    F    30       80
## 13   112    F    35       90
## 14   105    F    35       80
## 15   118    M    35       70
## 16   108    M    35       90
## 17   113    F    40       90
## 18   112    F    40       90
## 19   127    M    40       60
## 20   117    M    40       80
levels(df1$sexo)
## [1] "M" "F"
which(df1$sexo=="M")
##  [1]  1  3  6  7 10 11 15 16 19 20
which(df1$sexo=="F")
##  [1]  2  4  5  8  9 12 13 14 17 18

1.2 Resumo numérico e visualização dos dados

As funções aggregate() e boxplot() podem ser utilizadas para as tarefas de visualização dos dados e um resumo numérico.

head(df1)
##   tempo sexo idade acuidade
## 1    96    M    20       90
## 2    92    F    20      100
## 3   106    M    20       80
## 4   100    F    20       90
## 5    98    F    25      100
## 6   104    M    25       90
tail(df1)
##    tempo sexo idade acuidade
## 15   118    M    35       70
## 16   108    M    35       90
## 17   113    F    40       90
## 18   112    F    40       90
## 19   127    M    40       60
## 20   117    M    40       80
#
# Resumo numérico
#
aux1 <- aggregate(df1$tempo,by=list(df1$sexo),FUN=length)
aux2 <- aggregate(df1$tempo,by=list(df1$sexo),FUN=mean)
aux3 <- aggregate(df1$tempo,by=list(df1$sexo),FUN=var)
aux4 <- aggregate(df1$tempo,by=list(df1$sexo),FUN=sd)
tab.resumo <- cbind(aux1,aux2$x,aux3$x,aux4$x)
names(tab.resumo) <- c("Grupo","Tamanho","Media","Variancia","DP")
tab.resumo
##   Grupo Tamanho Media Variancia       DP
## 1     M      10 110.1  74.54444 8.633912
## 2     F      10 104.9  62.98889 7.936554
class(tab.resumo)
## [1] "data.frame"
dim(tab.resumo)
## [1] 2 5
#
# Visualizacao
#
boxplot(tempo~sexo,data=df1,xlab="Sexo",ylab="Tempo",col=c("lightblue","lightpink"))

dev.off()
## null device 
##           1

1.3 ANOVA com 1 fator de dois níveis

O modelo para os dados. A variável \(y\) representa o tempo de reação a um estímulo visual e a variável sexo é um fator com dois níveis.

\[ y_{ij}=\mu_{i}+\varepsilon_{ij},\hspace{0.5cm}i=1,2;j=1,\ldots,n_{i} \]

onde \(\varepsilon_{ij}\) são erros aleatórios independentes e normalmente distribuídos com média comum zero e variância comum \(\sigma^{2}\). Este modelo possui três parâmetros, a saber: \(-\infty<\mu_{1},\mu_{2}<+\infty\) e \(\sigma^{2}>0\).

Estimação de parâmetros. Estimativas pontuais e intervalares para os parâmetros do modelo podem ser obtidas como segue.

#--------------------------------------
#
# Numero de niveis e tamanhos amostrais
#
I <- length(levels(df1$sexo)) # numero de niveis do fator
n1 <- length(which(df1$sexo=="M")) # tamanho da amostra da populacao 1
n2 <- length(which(df1$sexo=="F")) # tamanho da amostra da populacao 2
n <- n1+n2 # tamanho da amostral total
I
## [1] 2
c(n1,n2,n)
## [1] 10 10 20
tab.resumo
##   Grupo Tamanho Media Variancia       DP
## 1     M      10 110.1  74.54444 8.633912
## 2     F      10 104.9  62.98889 7.936554
#--------------------------------------
#
# Estimacao pontual para mu1, mu2 e sigma
#
m1.hat <- tab.resumo$Media[1]
m2.hat <- tab.resumo$Media[2]
sig1.hat2 <- tab.resumo$Variancia[1]
sig2.hat2 <- tab.resumo$Variancia[2]
sqr <- (n1-1)*sig1.hat2+(n2-1)*sig2.hat2 # soma de quadrados residual
sig.hat2 <- sqr/(n-I) # variancia amostral ponderada
sig.hat <- sqrt(sig.hat2) # DP amostral ponderado (residual standard error)
aux5 <- c("Média1","Média2","Variância dos erros","Desvio Padrão dos erros")
data.frame(Parâmetro=aux5,Estimativa=c(m1.hat,m2.hat,sig.hat2,sig.hat))
##                 Parâmetro Estimativa
## 1                  Média1 110.100000
## 2                  Média2 104.900000
## 3     Variância dos erros  68.766667
## 4 Desvio Padrão dos erros   8.292567
#
# Obs.: a estimativa pontual do desvio padrao dos erros eh tambem 
# conhecida como: erro padrao residual (residual standard error)
#
mi.hat <- c(m1.hat,m2.hat)
sigi.hat2 <- c(sig1.hat2,sig2.hat2)
mi.hat
## [1] 110.1 104.9
sigi.hat2
## [1] 74.54444 62.98889
sqr # soma de quadrados residual
## [1] 1237.8
n-I # graus de liberdade
## [1] 18
#--------------------------------------
#
# Estimacao intervalar para mu1 e mu2
#
t <- qt(p=0.05/2,df=n-I,lower.tail=FALSE)
ic.m1 <- c(m1.hat-t*sig.hat*sqrt(1/n1),m1.hat+t*sig.hat*sqrt(1/n1))
ic.m2 <- c(m2.hat-t*sig.hat*sqrt(1/n2),m2.hat+t*sig.hat*sqrt(1/n2))
ic.m1
## [1] 104.5907 115.6093
ic.m2
## [1]  99.39067 110.40933
#--------------------------------------
#
# Estimacao para dif=mu1-mu2
#
t <- qt(p=0.05/2,df=n-I,lower.tail=FALSE)
dif.hat <- m1.hat-m2.hat
ic.dif <- c(dif.hat-t*sig.hat*sqrt((1/n1)+(1/n2)),dif.hat+t*sig.hat*sqrt((1/n1)+(1/n2)))
ic.dif
## [1] -2.591372 12.991372
#--------------------------------------

Tabela ANOVA e o teste-F para comparação de médias. A função aov() pode ser usada para gerar o modelo. Consequentemente, a tabela ANOVA pode ser extraída e o teste de comparação de médias realizado. O teste de homocedasticidade é realizado com a função bartlett.test().

modelo <- aov(formula=tempo~sexo,data=df1)
names(modelo) # Extratores do modelo
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"
tab.anova <- summary(modelo)
tab.anova
##             Df Sum Sq Mean Sq F value Pr(>F)
## sexo         1  135.2  135.20   1.966  0.178
## Residuals   18 1237.8   68.77
bartlett.test(tempo~sexo,data=df1) # teste de homocedasticidade
## 
##  Bartlett test of homogeneity of variances
## 
## data:  tempo by sexo
## Bartlett's K-squared = 0.060404, df = 1, p-value = 0.8059
#
# Exemplos de extratores do modelo
#
modelo$coefficients
## (Intercept)       sexoF 
##       110.1        -5.2
modelo$fitted.values
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 110.1 104.9 110.1 104.9 104.9 110.1 110.1 104.9 104.9 110.1 110.1 104.9 104.9 
##    14    15    16    17    18    19    20 
## 104.9 110.1 110.1 104.9 104.9 110.1 110.1
modelo$residuals
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## -14.1 -12.9  -4.1  -4.9  -6.9  -6.1  -0.1  -3.9  11.1  -4.1  -1.1  -4.9   7.1 
##    14    15    16    17    18    19    20 
##   0.1   7.9  -2.1   8.1   7.1  16.9   6.9
modelo$df.residual
## [1] 18

Da tabela ANOVA vemos que a hipótese nula (\(H_{0}:\mu_{1}=\mu_{2}\)) não é rejeitada ao nível \(\alpha=5\%\) (o p-valor do teste-F é maior que o nível \(\alpha\) adotado). Isto significa que o fator sexo não afeta o tempo de reação de um indivíduo a um estímulo visual.

Estimação de parâmetros. Extratores do objeto modelo também podem ser usados para gerar estimativas pontuais e intervalos de confiança para os parâmetros do modelo.

#--------------------------------------
#
# Estimacao via extratores do objeto modelo
#
coefficients(modelo) 
## (Intercept)       sexoF 
##       110.1        -5.2
confint(modelo)
##                 2.5 %     97.5 %
## (Intercept) 104.59067 115.609332
## sexoF       -12.99137   2.591372
#--------------------------------------
#
# Numero de niveis e tamanhos amostrais
#
I <- length(levels(df1$sexo)) # numero de niveis do fator
n1 <- length(which(df1$sexo=="M")) # tamanho da amostra da populacao 1
n2 <- length(which(df1$sexo=="F")) # tamanho da amostra da populacao 2
n <- n1+n2 # tamanho da amostral total
I
## [1] 2
c(n1,n2,n)
## [1] 10 10 20
#--------------------------------------
#
# Estimacao pontual para mu1, mu2 e sigma
#
mu1.hat <- as.numeric(coefficients(modelo)[1])
mu2.hat <- sum(coefficients(modelo)[c(1,2)])
sq.residual <- sum(residuals(modelo)^(2)) # soma de quadrados residual
sigma.hat2 <- sq.residual/df.residual(modelo)
sigma.hat <- sqrt(sigma.hat2)
aux5 <- c("Média1","Média2","Variância dos erros","Desvio Padrão dos erros")
data.frame(Parâmetro=aux5,Estimativa=c(mu1.hat,mu2.hat,sigma.hat2,sigma.hat))
##                 Parâmetro Estimativa
## 1                  Média1 110.100000
## 2                  Média2 104.900000
## 3     Variância dos erros  68.766667
## 4 Desvio Padrão dos erros   8.292567
mui.hat <- c(mu1.hat,mu2.hat)
sigmai.hat2 <- c(sum(residuals(modelo)[which(df1$sexo=="M")]^(2))/(n1-1),
                 sum(residuals(modelo)[which(df1$sexo=="F")]^(2))/(n2-1))
mui.hat
## [1] 110.1 104.9
sigmai.hat2
## [1] 74.54444 62.98889
sq.residual # soma de quadrados residual
## [1] 1237.8
deviance(modelo) # funcao alternativa para a soma de quadrados residual
## [1] 1237.8
df.residual(modelo) # graus de liberdade
## [1] 18
#--------------------------------------
#
# Estimacao intervalar para mu1 e mu2
#
ic.mu1 <- confint(modelo)[1,]
ic.mu2 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[2])
ic.mu1
##    2.5 %   97.5 % 
## 104.5907 115.6093
ic.mu2
##     2.5 %    97.5 % 
##  99.39067 110.40933
#--------------------------------------
#
# Estimacao para d=mu1-mu2
#
d.hat <- mu1.hat-mu2.hat
ic.d <- -confint(modelo)[2,2:1]
names(ic.d) <- names(ic.d)[2:1]
ic.d
##     2.5 %    97.5 % 
## -2.591372 12.991372
#--------------------------------------

Análise dos resíduos. ANOVA pressupõe que os erros \(\varepsilon_{ij}\) são independentes e normalmente distribuídos com média comum zero e variância comum \(\sigma^{2}\). Essas suposições devem ser verificadas através da análise dos resíduos.

n
## [1] 20
sq.residual
## [1] 1237.8
sqt <- sum((df1$tempo-mean(df1$tempo))^(2))
sqt
## [1] 1373
r2 <- 1-(sq.residual/sqt) # coeficiente de explicacao do modelo
r2
## [1] 0.0984705
par(mfrow=c(2,2))
boxplot(modelo$residuals~modelo$model[,2],xlab="Sexo",ylab=expression(r),
        names=c("M","F"),col=c("lightblue","lightpink"))
plot(modelo$residuals~modelo$fitted.values,xlab="Valores ajustados",ylab=expression(r))
qqnorm(modelo$residuals/sigma.hat,xlab="Quantis teóricos",ylab=expression(u==r/hat(sigma)),main="")
qqline(modelo$residuals/sigma.hat,col="red")
plot(modelo$fitted.values,sqrt(abs(modelo$residuals/sigma.hat)),
     xlab="Valores ajustados",ylab=expression(sqrt(abs(u))))

dev.off()
## null device 
##           1

1.4 Equivalência com os resultados do teste-t

No caso de duas populações normais independentes, ANOVA com 1 fator de dois níveis é equivalente ao teste-t. Portanto, a função t.test() também pode ser usada para realizar o teste \(H_{0}:\mu_{1}=\mu_{2}\).

x <- df1[which(df1$sexo=="M"),1] # Tempos de reacao dos individuos do sexo masculino
y <- df1[which(df1$sexo=="F"),1] # Tempos de reacao dos individuos do sexo feminino
cd1 <- list(Masculino=x,Feminino=y)
cd1
## $Masculino
##  [1]  96 106 104 110 106 109 118 108 127 117
## 
## $Feminino
##  [1]  92 100  98 101 116 100 112 105 113 112
out <- t.test(x,y,alternative="two.sided",mu=0,paired=FALSE,var.equal=TRUE,conf.level=0.95)
names(out)
##  [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
##  [6] "null.value"  "stderr"      "alternative" "method"      "data.name"
out
## 
##  Two Sample t-test
## 
## data:  x and y
## t = 1.4022, df = 18, p-value = 0.1779
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.591372 12.991372
## sample estimates:
## mean of x mean of y 
##     110.1     104.9
out$statistic^2
##        t 
## 1.966069
#
# Note que a estatistica t elevada ao quadrado é igual a estatistica F na tabela ANOVA
#

Obviamente que os resultados do teste-t (via função t.test()) concordam com os resultados do teste-F na tabela ANOVA (via função aov()).

2 Comparando \(I\) populações normais independentes

2.1 O problema e os dados

Problema. (Morettin & Bussab (2010), Capítulo 15, Exemplo 15.1) Considere novamente os dados do estudo apresentado seção anterior onde um psicólogo está investigando a relação entre o tempo que um indivíduo leva para reagir a um estímulo visual (\(y\)) e alguns fatores tais como: sexo (\(x_1\)), idade (\(x_2\)) e acuidade visual (\(x_3\)). Na presente seção, desejamos desenvolver um modelo descrevendo a relacão entre o tempo de reação a um estímulo visual e a variável idade, que neste estudo, será considerada como um fator com cinco níveis: indivíduos com 20 anos (\(i=1\) e \(n_{1}=4\)), indivíduos com 25 anos (\(i=2\) e \(n_{2}=4\)), indivíduos com 30 anos (\(i=3\) e \(n_{3}=4\)), indivíduos com 35 anos (\(i=4\) e \(n_{4}=4\)) e, finalmente, indivíduos com 40 anos (\(i=5\) e \(n_{5}=4\)). Novamente, o problema fundamental em questão é o seguinte: será que o conhecimento da idade de um indivíduo ajuda a melhorar a previsão do seu tempo de reação a um estímulo visual?

Dados. Os dados para o problema em questão são como segue.

str(df1)
## 'data.frame':    20 obs. of  4 variables:
##  $ tempo   : num  96 92 106 100 98 104 110 101 116 106 ...
##  $ sexo    : Factor w/ 2 levels "M","F": 1 2 1 2 2 1 1 2 2 1 ...
##  $ idade   : Factor w/ 5 levels "20","25","30",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ acuidade: num  90 100 80 90 100 90 80 90 70 90 ...
df1
##    tempo sexo idade acuidade
## 1     96    M    20       90
## 2     92    F    20      100
## 3    106    M    20       80
## 4    100    F    20       90
## 5     98    F    25      100
## 6    104    M    25       90
## 7    110    M    25       80
## 8    101    F    25       90
## 9    116    F    30       70
## 10   106    M    30       90
## 11   109    M    30       90
## 12   100    F    30       80
## 13   112    F    35       90
## 14   105    F    35       80
## 15   118    M    35       70
## 16   108    M    35       90
## 17   113    F    40       90
## 18   112    F    40       90
## 19   127    M    40       60
## 20   117    M    40       80
levels(df1$idade)
## [1] "20" "25" "30" "35" "40"
y1 <- df1[which(df1$idade=="20"),1] # Tempos de reacao dos individuos com 20 anos
y2 <- df1[which(df1$idade=="25"),1] # Tempos de reacao dos individuos com 25 anos
y3 <- df1[which(df1$idade=="30"),1] # Tempos de reacao dos individuos com 30 anos
y4 <- df1[which(df1$idade=="35"),1] # Tempos de reacao dos individuos com 35 anos
y5 <- df1[which(df1$idade=="40"),1] # Tempos de reacao dos individuos com 40 anos
cd2 <- list(I20=y1,I25=y2,I30=y3,I35=y4,I40=y5)
cd2
## $I20
## [1]  96  92 106 100
## 
## $I25
## [1]  98 104 110 101
## 
## $I30
## [1] 116 106 109 100
## 
## $I35
## [1] 112 105 118 108
## 
## $I40
## [1] 113 112 127 117

2.2 Resumo numérico e visualização dos dados

#
# Resumo numérico
#
aux1 <- aggregate(df1$tempo,by=list(df1$idade),FUN=length)
aux2 <- aggregate(df1$tempo,by=list(df1$idade),FUN=mean)
aux3 <- aggregate(df1$tempo,by=list(df1$idade),FUN=var)
aux4 <- aggregate(df1$tempo,by=list(df1$idade),FUN=sd)
tab.resumo <- cbind(aux1,aux2$x,aux3$x,aux4$x)
names(tab.resumo) <- c("Grupo","Tamanho","Media","Variancia","DP")
tab.resumo
##   Grupo Tamanho  Media Variancia       DP
## 1    20       4  98.50  35.66667 5.972158
## 2    25       4 103.25  26.25000 5.123475
## 3    30       4 107.75  44.25000 6.652067
## 4    35       4 110.75  31.58333 5.619905
## 5    40       4 117.25  46.91667 6.849574
#
# Visualizacao
#
boxplot(df1$tempo~df1$idade,xlab="Idade",ylab="Tempo",col=c("tan","tan1","tan2","tan3","tan4"))

dev.off()
## null device 
##           1

2.3 ANOVA com 1 fator de \(I\) níveis (\(I\geq 2\))

2.3.1 O modelo

A variável \(y\) representa o tempo de reação a um estímulo visual e a variável idade é um fator com cinco níveis (\(I=5\)).

\[ y_{ij}=\mu_{i}+\varepsilon_{ij},\hspace{0.5cm}i=1,\ldots,I;j=1,\ldots,n_{i}. \]

onde \(\varepsilon_{ij}\) são erros aleatórios independentes, normalmente distribuídos com média comum zero e variância comum \(\sigma^{2}\). Neste modelo, temos agora seis parâmetros, a saber: \(-\infty<\mu_{1},\mu_{2},\mu_{3},\mu_{4},\mu_{5}<+\infty\) e \(\sigma^{2}>0\).

2.3.2 A tabela ANOVA e o teste-F

modelo <- aov(formula=tempo~idade,data=df1)
names(modelo) # Extratores do modelo
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"
tab.anova <- summary(modelo)
tab.anova
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## idade        4    819  204.75   5.544 0.00605 **
## Residuals   15    554   36.93                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bartlett.test(tempo~idade,data=df1) # teste de homocedasticidade
## 
##  Bartlett test of homogeneity of variances
## 
## data:  tempo by idade
## Bartlett's K-squared = 0.29867, df = 4, p-value = 0.9899

A tabela ANOVA mostra que a hipótese nula (\(H_{0}:\mu_{1}=\ldots=\mu_{5}\)) deve ser rejeitada ao nível \(\alpha=5\%\) (o p-valor do teste-F é menor que o nível \(\alpha\) adotado). Isto significa que há evidências de que o fator idade afeta o tempo de reação de um indivíduo a um estímulo visual, de tal forma que os tempos médios de reação para os diversos grupos de idade não são todos iguais.

2.3.3 Estimação pontual e intervalar

#--------------------------------------
coefficients(modelo) 
## (Intercept)     idade25     idade30     idade35     idade40 
##       98.50        4.75        9.25       12.25       18.75
confint(modelo)
##                   2.5 %    97.5 %
## (Intercept) 92.02329205 104.97671
## idade25     -4.40944822  13.90945
## idade30      0.09055178  18.40945
## idade35      3.09055178  21.40945
## idade40      9.59055178  27.90945
#--------------------------------------
#
# Numero de niveis e tamanhos amostrais
#
I <- length(levels(df1$idade)) # numero de niveis do fator
n1 <- length(which(df1$idade=="20")) # tamanho da amostra da populacao 1
n2 <- length(which(df1$idade=="25")) # tamanho da amostra da populacao 2
n3 <- length(which(df1$idade=="30")) # tamanho da amostra da populacao 3
n4 <- length(which(df1$idade=="35")) # tamanho da amostra da populacao 4
n5 <- length(which(df1$idade=="40")) # tamanho da amostra da populacao 5
n <- n1+n2+n3+n4+n5 # tamanho da amostral total
I
## [1] 5
c(n1,n2,n3,n4,n5,n)
## [1]  4  4  4  4  4 20
#--------------------------------------
#
# Estimacao pontual das médias e de sigma
#
mu1.hat <- as.numeric(coefficients(modelo)[1])
mu2.hat <- sum(coefficients(modelo)[c(1,2)])
mu3.hat <- sum(coefficients(modelo)[c(1,3)])
mu4.hat <- sum(coefficients(modelo)[c(1,4)])
mu5.hat <- sum(coefficients(modelo)[c(1,5)])
mui.hat <- c(mu1.hat,mu2.hat,mu3.hat,mu4.hat,mu5.hat)
sq.residual <- sum(residuals(modelo)^(2)) # soma de quadrados residual
sigma.hat2 <- sq.residual/df.residual(modelo)
sigma.hat <- sqrt(sigma.hat2)
aux5 <- c("Média1","Média2","Média3","Média4","Média5",
          "Variância dos erros","Desvio Padrão dos erros")
data.frame(Parâmetro=aux5,Estimativa=c(mui.hat,sigma.hat2,sigma.hat))
##                 Parâmetro Estimativa
## 1                  Média1   98.50000
## 2                  Média2  103.25000
## 3                  Média3  107.75000
## 4                  Média4  110.75000
## 5                  Média5  117.25000
## 6     Variância dos erros   36.93333
## 7 Desvio Padrão dos erros    6.07728
sigmai.hat2 <- c(sum(residuals(modelo)[which(df1$idade=="20")]^(2))/(n1-1),
                 sum(residuals(modelo)[which(df1$idade=="25")]^(2))/(n2-1),
                 sum(residuals(modelo)[which(df1$idade=="30")]^(2))/(n3-1),
                 sum(residuals(modelo)[which(df1$idade=="35")]^(2))/(n4-1),
                 sum(residuals(modelo)[which(df1$idade=="40")]^(2))/(n5-1))
mui.hat
## [1]  98.50 103.25 107.75 110.75 117.25
sigmai.hat2
## [1] 35.66667 26.25000 44.25000 31.58333 46.91667
sq.residual # soma de quadrados residual
## [1] 554
deviance(modelo) # funcao alternativa para a soma de quadrados residual
## [1] 554
df.residual(modelo) # graus de liberdade
## [1] 15
#--------------------------------------
#
# Estimacao intervalar para as médias
#
ic.mu1 <- confint(modelo)[1,]
ic.mu2 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[2])
ic.mu3 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[3])
ic.mu4 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[4])
ic.mu5 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[5])
ic.mu1
##     2.5 %    97.5 % 
##  92.02329 104.97671
ic.mu2
##     2.5 %    97.5 % 
##  96.77329 109.72671
ic.mu3
##    2.5 %   97.5 % 
## 101.2733 114.2267
ic.mu4
##    2.5 %   97.5 % 
## 104.2733 117.2267
ic.mu5
##    2.5 %   97.5 % 
## 110.7733 123.7267
#--------------------------------------

2.3.4 Comparações múltiplas

Podemos usar o método de Tukey para realizar comparações múltiplas entre as médias dos diferentes níveis do fator.

tab.resumo
##   Grupo Tamanho  Media Variancia       DP
## 1    20       4  98.50  35.66667 5.972158
## 2    25       4 103.25  26.25000 5.123475
## 3    30       4 107.75  44.25000 6.652067
## 4    35       4 110.75  31.58333 5.619905
## 5    40       4 117.25  46.91667 6.849574
c(I,n,n1,n2,n3,n4,n5)
## [1]  5 20  4  4  4  4  4
sig.hat2 <- sum((tab.resumo$Tamanho-1)*tab.resumo$Variancia)/(n-I)
prob.c <- 0.95 # the family-wise probability of coverage
t <- qtukey(p=prob.c,nmeans=I,df=n-I);t
## [1] 4.366985
difmu <- NULL
hsd <- NULL # Tukey's honestly significant difference
ic.inf <- NULL
ic.sup <- NULL
for(i2 in 1:4){
  for(i1 in (i2+1):5){
    difmu.aux <- tab.resumo$Media[i1]-tab.resumo$Media[i2]
    hsd.aux <- t*sqrt((sig.hat2/2)*((1/tab.resumo$Tamanho[i1])+(1/tab.resumo$Tamanho[i2])))  
    difmu <- c(difmu,difmu.aux)
    hsd <- c(hsd,hsd.aux)
    ic.inf <- c(ic.inf,difmu.aux-hsd.aux)
    ic.sup <- c(ic.sup,difmu.aux+hsd.aux)
  }
}
#
# ICs para as diferencas entre medias
#
aux <- c("25-20","30-20","35-20","40-20","30-25",
         "35-25","40-25","35-30","40-30","40-35")
data.frame(aux,difmu,hsd,ic.inf,ic.sup)
##      aux difmu      hsd      ic.inf   ic.sup
## 1  25-20  4.75 13.26969  -8.5196946 18.01969
## 2  30-20  9.25 13.26969  -4.0196946 22.51969
## 3  35-20 12.25 13.26969  -1.0196946 25.51969
## 4  40-20 18.75 13.26969   5.4803054 32.01969
## 5  30-25  4.50 13.26969  -8.7696946 17.76969
## 6  35-25  7.50 13.26969  -5.7696946 20.76969
## 7  40-25 14.00 13.26969   0.7303054 27.26969
## 8  35-30  3.00 13.26969 -10.2696946 16.26969
## 9  40-30  9.50 13.26969  -3.7696946 22.76969
## 10 40-35  6.50 13.26969  -6.7696946 19.76969

Por método de Tukey, apenas duas diferenças são significativas, a saber: entre os niveis 5 (40 anos) e 1 (20 anos) e entre os niveis 5 (40 anos) e 2 (25 anos). De fato, existe uma função pronta para realizar comparações múltiplas pelo método de Tukey, a saber: TukeyHSD(). Utilize a função help() para entender melhor essa função e seus argumentos. Note que os resultados obtidos via função TukeyHSD() são exatamente os mesmos que foram obtidos acima, com a vantagem que a saída da função TukeyHSD() pode ser visualizada facilmente usando a função plot. A desvantagem é que os rótulos do gráfico estão escritos em Inglês e não podem ser modificados (pelo menos, não facilmente).

out <- TukeyHSD(modelo,conf.level=prob.c)
out
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = tempo ~ idade, data = df1)
## 
## $idade
##        diff         lwr      upr     p adj
## 25-20  4.75  -8.5196946 18.01969 0.8012355
## 30-20  9.25  -4.0196946 22.51969 0.2495379
## 35-20 12.25  -1.0196946 25.51969 0.0773018
## 40-20 18.75   5.4803054 32.01969 0.0043044
## 30-25  4.50  -8.7696946 17.76969 0.8296822
## 35-25  7.50  -5.7696946 20.76969 0.4379674
## 40-25 14.00   0.7303054 27.26969 0.0363399
## 35-30  3.00 -10.2696946 16.26969 0.9537673
## 40-30  9.50  -3.7696946 22.76969 0.2282694
## 40-35  6.50  -6.7696946 19.76969 0.5703582
plot(out,las=2)

dev.off()
## null device 
##           1

Alternativamente, podemos realizar comparações múltiplas usando a função pairwise.t.test(). Novamente, utilize a função help() para entender melhor essa função e seus argumentos.

x <- df1$tempo # resposta
f <- df1$idade # fator
pairwise.t.test(x,f,p.adjust.method="bonferroni",pool.sd=T,paired=F,alternative="two.sided")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  x and f 
## 
##    20     25     30     35    
## 25 1.0000 -      -      -     
## 30 0.4804 1.0000 -      -     
## 35 0.1215 1.0000 1.0000 -     
## 40 0.0056 0.0530 0.4301 1.0000
## 
## P value adjustment method: bonferroni

Os resultados obtidos via função pairwise.t.test() concordam quase que completamente com os resultados obtidos via função TukeyHSD.

2.3.5 Análise dos resíduos

n
## [1] 20
sq.residual <- sum(residuals(modelo)^(2)) # soma de quadrados residual
sq.residual
## [1] 554
sqt <- sum((df1$tempo-mean(df1$tempo))^(2))
sqt
## [1] 1373
r2 <- 1-(sq.residual/sqt) # coeficiente de explicacao do modelo
r2
## [1] 0.596504
par(mfrow=c(2,2))
boxplot(modelo$residuals~modelo$model[,2],xlab="Idade",ylab=expression(r),
        names=c("20","25","30","35","40"),col=c("tan","tan1","tan2","tan3","tan4"))
plot(modelo$residuals~modelo$fitted.values,xlab="Valores ajustados",ylab=expression(r))
qqnorm(modelo$residuals/sigma.hat,xlab="Quantis teóricos",ylab=expression(u==r/hat(sigma)),main="")
qqline(modelo$residuals/sigma.hat,col="red")
plot(modelo$fitted.values,sqrt(abs(modelo$residuals/sigma.hat)),
     xlab="Valores ajustados",ylab=expression(sqrt(abs(u))))

dev.off()
## null device 
##           1

2.4 ANOVA via função lm()

2.4.1 A função stack() para estruturar os dados

Suponha que o conjunto de dados para realizar uma ANOVA esteja estruturado como uma lista conforme o conjunto cd2 no script abaixo. Este é exatamente o mesmo conjunto de dados descrito no início da presente Seção.

#
# Tempos de reacao de individuos a um estimulo visual por idade
#
cd2
## $I20
## [1]  96  92 106 100
## 
## $I25
## [1]  98 104 110 101
## 
## $I30
## [1] 116 106 109 100
## 
## $I35
## [1] 112 105 118 108
## 
## $I40
## [1] 113 112 127 117
class(cd2)
## [1] "list"

Podemos usar a função stack() para reestruturar este conjunto de dados como uma dataframe.

d <- stack(cd2)
names(d)
## [1] "values" "ind"
d
##    values ind
## 1      96 I20
## 2      92 I20
## 3     106 I20
## 4     100 I20
## 5      98 I25
## 6     104 I25
## 7     110 I25
## 8     101 I25
## 9     116 I30
## 10    106 I30
## 11    109 I30
## 12    100 I30
## 13    112 I35
## 14    105 I35
## 15    118 I35
## 16    108 I35
## 17    113 I40
## 18    112 I40
## 19    127 I40
## 20    117 I40
class(d)
## [1] "data.frame"
str(d)
## 'data.frame':    20 obs. of  2 variables:
##  $ values: num  96 92 106 100 98 104 110 101 116 106 ...
##  $ ind   : Factor w/ 5 levels "I20","I25","I30",..: 1 1 1 1 2 2 2 2 3 3 ...
boxplot(values~ind,data=d,xlab="Idade",ylab="Tempo",names=c("20","25","30","35","40"),
        col=c("tan","tan1","tan2","tan3","tan4"))

dev.off()
## null device 
##           1

2.4.2 ANOVA com lm()

Estimação, tabela ANOVA e o teste-F.

head(d)
##   values ind
## 1     96 I20
## 2     92 I20
## 3    106 I20
## 4    100 I20
## 5     98 I25
## 6    104 I25
tail(d)
##    values ind
## 15    118 I35
## 16    108 I35
## 17    113 I40
## 18    112 I40
## 19    127 I40
## 20    117 I40
modelo <- lm(values~ind,data=d)
summary(modelo)
## 
## Call:
## lm(formula = values ~ ind, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.750 -4.500 -1.000  2.812  9.750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   98.500      3.039  32.416 2.65e-15 ***
## indI25         4.750      4.297   1.105 0.286423    
## indI30         9.250      4.297   2.153 0.048045 *  
## indI35        12.250      4.297   2.851 0.012151 *  
## indI40        18.750      4.297   4.363 0.000556 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.077 on 15 degrees of freedom
## Multiple R-squared:  0.5965, Adjusted R-squared:  0.4889 
## F-statistic: 5.544 on 4 and 15 DF,  p-value: 0.006055
anova(modelo)
## Analysis of Variance Table
## 
## Response: values
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## ind        4    819 204.750  5.5438 0.006055 **
## Residuals 15    554  36.933                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bartlett.test(values~ind,data=d) # teste de homocedasticidade
## 
##  Bartlett test of homogeneity of variances
## 
## data:  values by ind
## Bartlett's K-squared = 0.29867, df = 4, p-value = 0.9899

Estimação pontual para as medias e o coeficiente de explicação do modelo.

#
# Estimacao pontual para as medias
#
coefficients(modelo)
## (Intercept)      indI25      indI30      indI35      indI40 
##       98.50        4.75        9.25       12.25       18.75
mu1.hat <- as.numeric(coefficients(modelo)[1])
mu2.hat <- sum(coefficients(modelo)[c(1,2)])
mu3.hat <- sum(coefficients(modelo)[c(1,3)])
mu4.hat <- sum(coefficients(modelo)[c(1,4)])
mu5.hat <- sum(coefficients(modelo)[c(1,5)])
mui.hat <- c(mu1.hat,mu2.hat,mu3.hat,mu4.hat,mu5.hat);mui.hat
## [1]  98.50 103.25 107.75 110.75 117.25
#
# Estimacao pontual para a variancia do erro e 
# o coeficiente de explicacao do modelo (R^2)
#
y <- d$values
tss <- sum((y-mean(y))^2)
rss <- deviance(modelo) # soma de quadrados residual
sigma2hat <- rss/df.residual(modelo)
sigmahat <- sqrt(sigma2hat)
sigma2hat 
## [1] 36.93333
sigmahat # erro padrao residual (residual standard error)
## [1] 6.07728
R2 <- 1-(rss/tss);R2 # coeficiente de explicacao do modelo
## [1] 0.596504

Estimação intervalar para as medias.

confint(modelo)
##                   2.5 %    97.5 %
## (Intercept) 92.02329205 104.97671
## indI25      -4.40944822  13.90945
## indI30       0.09055178  18.40945
## indI35       3.09055178  21.40945
## indI40       9.59055178  27.90945
ic.mu1 <- confint(modelo)[1,]
ic.mu2 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[2])
ic.mu3 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[3])
ic.mu4 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[4])
ic.mu5 <- confint(modelo)[1,]+as.numeric(coefficients(modelo)[5])
ic.mu1
##     2.5 %    97.5 % 
##  92.02329 104.97671
ic.mu2
##     2.5 %    97.5 % 
##  96.77329 109.72671
ic.mu3
##    2.5 %   97.5 % 
## 101.2733 114.2267
ic.mu4
##    2.5 %   97.5 % 
## 104.2733 117.2267
ic.mu5
##    2.5 %   97.5 % 
## 110.7733 123.7267

Comparações múltiplas pelo método de Tukey.

modelo.aov <- aov(values~ind,data=d)
out <- TukeyHSD(modelo.aov)
out
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = values ~ ind, data = d)
## 
## $ind
##          diff         lwr      upr     p adj
## I25-I20  4.75  -8.5196946 18.01969 0.8012355
## I30-I20  9.25  -4.0196946 22.51969 0.2495379
## I35-I20 12.25  -1.0196946 25.51969 0.0773018
## I40-I20 18.75   5.4803054 32.01969 0.0043044
## I30-I25  4.50  -8.7696946 17.76969 0.8296822
## I35-I25  7.50  -5.7696946 20.76969 0.4379674
## I40-I25 14.00   0.7303054 27.26969 0.0363399
## I35-I30  3.00 -10.2696946 16.26969 0.9537673
## I40-I30  9.50  -3.7696946 22.76969 0.2282694
## I40-I35  6.50  -6.7696946 19.76969 0.5703582
plot(out,las=2)

dev.off()
## null device 
##           1

Análise dos resíduos.

par(mfrow=c(2,2))
boxplot(modelo$residuals~modelo$model[,2],xlab="Idade",ylab=expression(r),
        names=c("20","25","30","35","40"),col=c("tan","tan1","tan2","tan3","tan4"))
plot(modelo$residuals~modelo$fitted.values,xlab="Valores ajustados",ylab=expression(r))
qqnorm(modelo$residuals/sigma.hat,xlab="Quantis teóricos",ylab=expression(u==r/hat(sigma)),main="")
qqline(modelo$residuals/sigma.hat,col="red")
plot(modelo$fitted.values,sqrt(abs(modelo$residuals/sigma.hat)),
     xlab="Valores ajustados",ylab=expression(sqrt(abs(u))))

dev.off()
## null device 
##           1

3 Uma abordagem não-paramétrico

3.1 O problema e os dados

Uma atração natural nos Estados Unidos é o gêiser Old Faithful, que fica localizado no Parque Nacional de Yellowstone. O conjunto de dados no script abaixo apresenta intervalos de tempo (em minutos) entre erupções sucessivas para 4 anos diferentes (1951, 1985, 1995 e 1996). Uma análise descritiva dos dados sugere que essas amostras se originam de populações com médias diferentes. Uma questão fundamental agora é a seguinte: será que essas diferenças amostrais são estatisticamente significativas?

y1 <- c(74,60,74,42,74,52,65,68,62,66,62,60)
y2 <- c(89,90,60,65,82,84,54,85,58,79,57,88)
y3 <- c(86,86,62,104,62,95,79,62,94,79,86,85)
y4 <- c(88,86,85,89,83,85,91,68,91,56,89,94)
cd3 <- list("1951"=y1,"1985"=y2,"1995"=y3,"1996"=y4)
d <- stack(cd3)
str(d)
## 'data.frame':    48 obs. of  2 variables:
##  $ values: num  74 60 74 42 74 52 65 68 62 66 ...
##  $ ind   : Factor w/ 4 levels "1951","1985",..: 1 1 1 1 1 1 1 1 1 1 ...
aux1 <- aggregate(d$values,by=list(d$ind),FUN=length)
aux2 <- aggregate(d$values,by=list(d$ind),FUN=mean)
aux3 <- aggregate(d$values,by=list(d$ind),FUN=var)
aux4 <- aggregate(d$values,by=list(d$ind),FUN=sd)
tab.resumo <- cbind(aux1,aux2$x,aux3$x,aux4$x)
names(tab.resumo) <- c("Grupo","Tamanho","Media","Variancia","DP")
tab.resumo
##   Grupo Tamanho    Media Variancia        DP
## 1  1951      12 63.25000  89.29545  9.449627
## 2  1985      12 74.25000 200.75000 14.168627
## 3  1995      12 81.66667 188.24242 13.720147
## 4  1996      12 83.75000 119.11364 10.913919
boxplot(values~ind,data=d,xlab="Ano",ylab="Intervalo de tempo",col=c("red","blue","green","yellow"))

dev.off()
## null device 
##           1

A análise descritiva sugere que a hipótese de homocedasticidade ou variância constante pode não ser válida (note que a amplitude das caixas são diferentes). Porém, o mais importante, é que a análise descritiva também sugere que a hipótese de normalidade pode não ser válida (note a falta de simetria das caixas e principalmente presença de pontos exteriores, destoantes das demais observações). Esses achados da análise descritiva sugerem que o teste não-paramétrico de Kruskal-Wallis pode, neste caso, ser uma abordagem mais apropriada para efetuar a comparação dessas populações independentes.

3.2 O teste de Kruskal-Wallis para populações independentes

d
##    values  ind
## 1      74 1951
## 2      60 1951
## 3      74 1951
## 4      42 1951
## 5      74 1951
## 6      52 1951
## 7      65 1951
## 8      68 1951
## 9      62 1951
## 10     66 1951
## 11     62 1951
## 12     60 1951
## 13     89 1985
## 14     90 1985
## 15     60 1985
## 16     65 1985
## 17     82 1985
## 18     84 1985
## 19     54 1985
## 20     85 1985
## 21     58 1985
## 22     79 1985
## 23     57 1985
## 24     88 1985
## 25     86 1995
## 26     86 1995
## 27     62 1995
## 28    104 1995
## 29     62 1995
## 30     95 1995
## 31     79 1995
## 32     62 1995
## 33     94 1995
## 34     79 1995
## 35     86 1995
## 36     85 1995
## 37     88 1996
## 38     86 1996
## 39     85 1996
## 40     89 1996
## 41     83 1996
## 42     85 1996
## 43     91 1996
## 44     68 1996
## 45     91 1996
## 46     56 1996
## 47     89 1996
## 48     94 1996
modelo.kw <- kruskal.test(values~ind,data=d)
print(modelo.kw)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  values by ind
## Kruskal-Wallis chi-squared = 14.479, df = 3, p-value = 0.00232
names(modelo.kw)
## [1] "statistic" "parameter" "p.value"   "method"    "data.name"

Com base no resultado do teste Kruskal-Wallis, rejeitamos a hipótese de populações com médias idênticas. Ou seja, existem evidências que os tempos médios entre erupções sucessivas não são todos iguais para os diferentes anos investigados.