Comandos R para análises estatísticas

Capítulo 15: Inferência para Várias Populações

Exemplo 15.1

Para que o R saiba que as variáveis sexo e idade deverão ser consideradas como fatores, é preciso declará-las como fatores:

Tabela 15.1

tab15_1<-data.frame(
  individuo=1:20,
  tempo_Y=c(96,92,106,100,98,104,110,101,116,106,109,100,112,105,118,108,113,112,127,117),
  sexo_W=factor(c("H","M","H","M","M","H","H","M","M","H","H","M","M","M","H","H","M","M","H","H")),
  idade_X=factor(c(20,20,20,20,25,25,25,25,30,30,30,30,35,35,35,35,40,40,40,40)),
  acuidade_Z=c(90,100,80,90,100,90,80,90,70,90,90,80,90,80,70,90,90,90,60,80))
Tabela 15.1: Tempos de reação a um estímulo(Y) e acuidade visual(Z) de 20 indivíduos, segundo o sexo(W) e a idade(X)
individuo tempo_Y sexo_W idade_X acuidade_Z
1 96 H 20 90
2 92 M 20 100
3 106 H 20 80
4 100 M 20 90
5 98 M 25 100
6 104 H 25 90
7 110 H 25 80
8 101 M 25 90
9 116 M 30 70
10 106 H 30 90
11 109 H 30 90
12 100 M 30 80
13 112 M 35 90
14 105 M 35 80
15 118 H 35 70
16 108 H 35 90
17 113 M 40 90
18 112 M 40 90
19 127 H 40 60
20 117 H 40 80

Assim, quando acessamos cada uma delas, vemos que há a identificação dos níveis dos fatores de cada uma delas:

attach(tab15_1)
sexo_W
##  [1] H M H M M H H M M H H M M M H H M M H H
## Levels: H M
idade_X
##  [1] 20 20 20 20 25 25 25 25 30 30 30 30 35 35 35 35 40 40 40 40
## Levels: 20 25 30 35 40

Se quisermos calcular as estatísticas descritivas para cada subgrupo, masculino e feminino, podemos utilizar os subvetores de acordo com a variável sexo:

# Estatísticas descritivas para o grupo dos homens:
summary2(tempo_Y[sexo_W=="H"])
##              [,1]
## N        10.00000
## Min.     96.00000
## 1st Qu. 106.00000
## Median  108.00000
## Mean    110.00000
## 3rd Qu. 115.00000
## Max.    127.00000
## Tr Mean 109.75000
## Var      74.54444
## StDev     8.63391
## SE Mean   0.92919
# Estatísticas descritivas para o grupo das mulheres:
summary2(tempo_Y[sexo_W=="M"])
##              [,1]
## N        10.00000
## Min.     92.00000
## 1st Qu. 100.00000
## Median  103.00000
## Mean    105.00000
## 3rd Qu. 112.00000
## Max.    116.00000
## Tr Mean 105.12500
## Var      62.98889
## StDev     7.93655
## SE Mean   0.89087
# Uma maneira mais elegante é utilizar a função sapply:
## Estatísticas descrtivas do tempo de reacAo por sexo
summary_sexo<-sapply(levels(sexo_W),function(sex){summary2(tempo_Y[sexo_W==sex])})
summary_sexo
##               H         M
##  [1,]  10.00000  10.00000
##  [2,]  96.00000  92.00000
##  [3,] 106.00000 100.00000
##  [4,] 108.00000 103.00000
##  [5,] 110.00000 105.00000
##  [6,] 115.00000 112.00000
##  [7,] 127.00000 116.00000
##  [8,] 109.75000 105.12500
##  [9,]  74.54444  62.98889
## [10,]   8.63391   7.93655
## [11,]   0.92919   0.89087
## Estatísticas descrtivas do tempo de reacAo por idade
summary_idade<-sapply(levels(idade_X),function(age){summary2(tempo_Y[idade_X==age])})
summary_idade
##             20       25       30       35       40
##  [1,]   4.0000   4.0000   4.0000   4.0000   4.0000
##  [2,]  92.0000  98.0000 100.0000 105.0000 112.0000
##  [3,]  95.0000 100.0000 104.0000 107.0000 113.0000
##  [4,]  98.0000 102.0000 108.0000 110.0000 115.0000
##  [5,]  98.5000 103.0000 108.0000 111.0000 117.0000
##  [6,] 102.0000 106.0000 111.0000 114.0000 120.0000
##  [7,] 106.0000 110.0000 116.0000 118.0000 127.0000
##  [8,]  98.0000 102.5000 107.5000 110.0000 115.0000
##  [9,]  35.6667  26.2500  44.2500  31.5833  46.9167
## [10,]   5.9722   5.1235   6.6521   5.6199   6.8496
## [11,]   1.2219   1.1318   1.2896   1.1853   1.3086

Tabela 15.4

Inicialmente, para calcularmos os resíduos, vamos construir dois vetores com as médias de cada grupo (por sexo e por idade) de acordo com as informações de cada indivíduo:

mean_idade<-rep(0,20)
for(age in levels(idade_X)){
  mean_idade[mean_idade==0]<-((idade_X==age)*mean(tempo_Y[idade_X==age]))[mean_idade==0]
}
mean_sexo<-rep(0,length(sexo_W))
for(sex in levels(sexo_W)){
  mean_sexo[mean_sexo==0]<-((sexo_W==sex)*mean(tempo_Y[sexo_W==sex]))[mean_sexo==0]
}
mean_idade
##  [1]  98.50  98.50  98.50  98.50 103.25 103.25 103.25 103.25 107.75 107.75
## [11] 107.75 107.75 110.75 110.75 110.75 110.75 117.25 117.25 117.25 117.25
mean_sexo
##  [1] 110.1 104.9 110.1 104.9 104.9 110.1 110.1 104.9 104.9 110.1 110.1
## [12] 104.9 104.9 104.9 110.1 110.1 104.9 104.9 110.1 110.1

Agora vamos calcular os resíduos e montar a Tabela 15.4:

e1<-tempo_Y-mean(tempo_Y)
e2<-tempo_Y-mean_sexo
e3<-tempo_Y-mean_idade

tab15_4<-cbind(tab15_1[,1:4],e1,e2,e3)
Tabela 15.4: Resíduos para vários modelos ajustados aos dados do Exemplo 15.1.
individuo tempo_Y sexo_W idade_X e1 e2 e3
1 96 H 20 -11.5000 -14.1000 -2.5000
2 92 M 20 -15.5000 -12.9000 -6.5000
3 106 H 20 -1.5000 -4.1000 7.5000
4 100 M 20 -7.5000 -4.9000 1.5000
5 98 M 25 -9.5000 -6.9000 -5.2500
6 104 H 25 -3.5000 -6.1000 0.7500
7 110 H 25 2.5000 -0.1000 6.7500
8 101 M 25 -6.5000 -3.9000 -2.2500
9 116 M 30 8.5000 11.1000 8.2500
10 106 H 30 -1.5000 -4.1000 -1.7500
11 109 H 30 1.5000 -1.1000 1.2500
12 100 M 30 -7.5000 -4.9000 -7.7500
13 112 M 35 4.5000 7.1000 1.2500
14 105 M 35 -2.5000 0.1000 -5.7500
15 118 H 35 10.5000 7.9000 7.2500
16 108 H 35 0.5000 -2.1000 -2.7500
17 113 M 40 5.5000 8.1000 -4.2500
18 112 M 40 4.5000 7.1000 -5.2500
19 127 H 40 19.5000 16.9000 9.7500
20 117 H 40 9.5000 6.9000 -0.2500
NA NA NA NA 8.5008 8.0714 5.3998

15.2.3 Intervalos de Confiança

Vamos construir a função IC_t para calcular IC para a distribuição t no caso de das subgrupos na amostra:

IC_t<-function(x_barra1,x_barra2,s1,s2,n1,n2,gama=0.95){
  left_tail<-(1-gama)/2
  gamma_adj<-left_tail+gama # ajuste para somar a cauda inferior da distribuição somente para encontrar o quantil gama usando a função dist. acumulada
  Se=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2)
  llim1<-x_barra1-qt(gamma_adj,df = (n1+n2-2))*sqrt(Se)/sqrt(n1) # calcula o limite inferior
  ulim1<-x_barra1+qt(gamma_adj,df = (n1+n2-2))*sqrt(Se)/sqrt(n1) # calcula o limite inferior
  llim2<-x_barra2-qt(gamma_adj,df = (n1+n2-2))*sqrt(Se)/sqrt(n2) # calcula o limite inferior
  ulim2<-x_barra2+qt(gamma_adj,df = (n1+n2-2))*sqrt(Se)/sqrt(n2) # calcula o limite inferior

    #writeLines(paste("IC_1(",gama,"):[",round(llim1,3),";",round(ulim1,3),"]\n",sep=""))
    #writeLines(paste("IC_2(",gama,"):[",round(llim2,3),";",round(ulim2,3),"]\n",sep=""))
  return(list(ic_calc1=c(inf.lim=llim1,upper.lim=ulim1),ic_calc2=c(inf.lim=llim2,upper.lim=ulim2)))
}

IC_t(mean(tempo_Y[sexo_W=="H"]),mean(tempo_Y[sexo_W=="M"]),var(tempo_Y[sexo_W=="H"]),var(tempo_Y[sexo_W=="M"]),10,10,gama=0.95)
## $ic_calc1
##   inf.lim upper.lim 
##    104.59    115.61 
## 
## $ic_calc2
##   inf.lim upper.lim 
##    99.391   110.409

Para o IC de diferença de médias, utilizaremos a mesma função IC_diff_medias que construímos no capítulo anterior:

IC_diff_medias(mean(tempo_Y[sexo_W=="M"]),mean(tempo_Y[sexo_W=="H"]),
var(tempo_Y[sexo_W=="M"]),var(tempo_Y[sexo_W=="H"]),10,10,gama=0.95)
## [1] "IC(0.95):]-2.5914;12.9914["

15.2.4 Tabela de Análise de Variância

Para construirmos uma tabela de ANOVA, devemos referenciá-la a um modelo linear, que será calculado com a função lm. Assim, a tabela de ANOVA para o modelo 15.4 que contém o fator sexo é:

fit15_4_sexo<-anova(lm(tempo_Y~sexo_W))
Tabela 15.6: Tabela ANOVA para o Exemplo 15.1.
FV gl SQ QM F
sexo_W 1 135.2 135.200 1.9661 0.17788
Residuals 18 1237.8 68.767 NA NA

Uma outra função que também calcula a tabela de ANOVA é a função aov. Nela é possível também aplicar o método plot e verificar a análise de resíduos do modelo:

par(mfrow=c(2,2))
plot(aov(lm(tempo_Y~sexo_W)))

Figura 15.1.A: Análise de resíduos para o Modelo de Estímulo em função do sexo.

Exemplo 15.4

qf(0.95,1,18)
## [1] 4.4139
if(fit15_4_sexo$`F value`[1]>qf(0.95,1,18)){
  print("Rejeito H_0")
} else print( "Aceito H0")
## [1] "Aceito H0"

Figura 15.2

O boxplot pode ser contruído com a notação y~x1+...+xk, em que y representa a variável de interesse e x1+...+xk as variáveis agrupadoras. Assim:

boxplot(tempo_Y~idade_X, col="lightblue3")

Figura 15.2: Box plots para a variável Y(estímulo) para cada nível de idade.

Podemos também construir tal figura utilizando o pacote ggplot2:

library(ggplot2)
p <- ggplot(tab15_1,aes(idade_X, tempo_Y))
p + geom_boxplot(stat = "boxplot", fill = "brown4",
             position = "dodge", colour = "#858585", 
             outlier.shape = 16, outlier.size = 1, notch = FALSE, 
             notchwidth = 0.5) +
  xlab("Idade") + ylab("Estímulo")+
             theme(panel.background = element_rect(fill = "azure2"))

Figura 15.2: Box plots para a variável Y(estímulo) para cada nível de idade utilizando o pacote ggplot2.

Para calcular a ANOVA para o fator idade, teremos:

fit15_4_idade<-anova(lm(tempo_Y~idade_X))
par(mfrow=c(2,2))
plot(aov(lm(tempo_Y~idade_X)))

Figura 15.1.B: Análise de resíduos para o Modelo de Estímulo em função do Idade.

Tabela 15.7: Tabela ANOVA para o Exemplo 15.1.
FV gl SQ QM F
idade_X 4 819 204.750 5.5438 0.00605
Residuals 15 554 36.933 NA NA 
plot(as.vector(idade_X),e3,col="darkred",pch=20,xaxt="n")
axis(1,at=c(19.6,24.6,29.6,34.6,39.6),labels=c("20 anos","25 anos","30 anos","35 anos","40 anos"))
abline(v=c(19.6,24.6,29.6,34.6,39.6), col="darkgrey", lty=2)

Figura 15.3: Resíduos do modelo \(y_{ij}=\mu_i+e_{ij}\) para o fator idade.

qf(0.95,4,15)
## [1] 3.0556
if(fit15_4_idade$`F value`[1]>qf(0.95,4,15)){
  print("Rejeito H_0")
} else print( "Aceito H0")
## [1] "Rejeito H_0"

Exemplo 15.5

medias_idade<-sapply(levels(idade_X),function(age){mean(tempo_Y[idade_X==age])})
medias_idade
##     20     25     30     35     40 
##  98.50 103.25 107.75 110.75 117.25

Para calcular a diferença entre as idades conforme calculado na Tabela 15.8, utilize a função diff

diff(medias_idade)
##   25   30   35   40 
## 4.75 4.50 3.00 6.50
print(kable(cbind(c("média","diferença"),rbind(medias_idade,c(0,diff(medias_idade)))),caption="**Tabela 15.8**:Médias e diferenças de médias para os diversos grupos de idade para o Exemplo 15.1.",col.names = c("","20 anos","25 anos","30 anos","35 anos","40 anos"),row.names = FALSE))
Tabela 15.8:Médias e diferenças de médias para os diversos grupos de idade para o Exemplo 15.1.
20 anos 25 anos 30 anos 35 anos 40 anos
média 98.5 103.25 107.75 110.75 117.25
diferença 0 4.75 4.5 3 6.5

Exemplo 15.6

Novamente, vamos construir uma função para construir os ICs:

IC_t_diff<-function(x_barra1,x_barra2,SQ,n,n1,n2,gama=0.95){
  left_tail<-(1-gama)/2
  gamma_adj<-left_tail+gama # ajuste para somar a cauda inferior da distribuição somente para encontrar o quantil gama usando a função dist. acumulada
  Se=SQ/(n-2)
  llim1<-(x_barra1-x_barra2)-qt(gamma_adj,df = (n-2))*sqrt(Se)*sqrt(1/n1+1/n2) # calcula o limite inferior
  ulim1<-(x_barra1-x_barra2)+qt(gamma_adj,df = (n-2))*sqrt(Se)*sqrt(1/n1+1/n2) # calcula o limite inferior
  #writeLines(paste("t_c:",qt(gamma_adj,df = (n-2))))
  #writeLines(paste("IC_1(",gama,"):[",round(llim1,3),";",round(ulim1,3),"]\n",sep=""))
  #writeLines(paste("IC_2(",gama,"):[",round(llim2,3),";",round(ulim2,3),"]\n",sep=""))
  return(ic_calc=c(inf.lim=llim1,upper.lim=ulim1))
}

Calculando agora todos os ICs de diferenças entre grupos de idades, temos:

v<-levels(idade_X)
SQ<-fit15_4_idade$`Sum Sq`[1]  #fit15_4_idade ;e o modelo ajustado e $`Sum Sq`[1] é a SQ para a idade.
ics_diff_idades<-matrix(NA,10,2)
cont=1
lab_ics<-character(0)
vec_diff_idades<-numeric()
for( i in 1:4)
  for(j in (i+1):5){
    ics_diff_idades[cont,]<-IC_t_diff(mean(tempo_Y[idade_X==v[i]]),mean(tempo_Y[idade_X==v[j]]),SQ,20,4,4,gama = 0.995)
    vec_diff_idades[cont]<-mean(tempo_Y[idade_X==v[i]])-mean(tempo_Y[idade_X==v[j]])
    lab_ics[cont]<-paste(v[i],v[j],sep="_")
    cont=cont+1
  }
rownames(ics_diff_idades)<-lab_ics
ics_diff_idades
##          [,1]    [,2]
## 20_25 -19.997 10.4967
## 20_30 -24.497  5.9967
## 20_35 -27.497  2.9967
## 20_40 -33.997 -3.5033
## 25_30 -19.747 10.7467
## 25_35 -22.747  7.7467
## 25_40 -29.247  1.2467
## 30_35 -18.247 12.2467
## 30_40 -24.747  5.7467
## 35_40 -21.747  8.7467
### Grafico com os ICs
plot(1:10,vec_diff_idades, ylim=c(min(ics_diff_idades)-1,max(ics_diff_idades)+1), type="p",xaxt="n",xlab="pares",ylab="ICs") 
axis(1,at=1:10,labels=lab_ics) # labels no eixo x
for (i in 1:10) {
  if(ics_diff_idades[i,1]>0 || ics_diff_idades[i,2]<0) {color="red"} else color="darkblue"
  lines(x =c(i,i), y=t(ics_diff_idades)[,i], col=color) # Plotando os ICs
}

Figura 15.6.A: Intervalos de confiança simultâneos para as diferenças de médias da variável Y(estímulo) por idade.

Exemplo 15.7

O teste de Bartlett é calculado no R a partir da função bartlett.test:

bartlett.test(x, g)
# x: Vetor com a variável de interesse
# g: vetor contendo a qual grupo cada elemento de x pertence
bartlett.test(x=tempo_Y, g=idade_X)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  tempo_Y and idade_X
## Bartlett's K-squared = 0.299, df = 4, p-value = 0.99

Quadro 15.1

O Quadro 15.1 é formado pela saída da análise de variância, que já calculamos anteriormente (fit15_4_idade) e os ICs de cada um dos grupos de idades:

fit15_4_idade
## Analysis of Variance Table
## 
## Response: tempo_Y
##           Df Sum Sq Mean Sq F value Pr(>F)
## idade_X    4    819   204.8    5.54 0.0061
## Residuals 15    554    36.9
ics_idades<-matrix(NA,5,2)
cont=1
for(age in levels(idade_X)){
  ics_idades[cont,]<-IC(n = 4, mean(tempo_Y[idade_X==age]),sd(tempo_Y[idade_X==age]),gama=0.95)
  cont=cont+1
}
## IC(0.95):[92.647;104.353]
## 
## IC(0.95):[98.229;108.271]
## 
## IC(0.95):[101.231;114.269]
## 
## IC(0.95):[105.243;116.257]
## 
## IC(0.95):[110.538;123.962]
print(kable(cbind(medias_idade,ics_idades),caption="Média e Intervalos de confiança por grupo de idade",col.names = c("Média","Limite Inferior", "Limite Superior")))
## 
## 
## Table: Média e Intervalos de confiança por grupo de idade
## 
##        Média   Limite Inferior   Limite Superior
## ---  -------  ----------------  ----------------
## 20     98.50            92.647            104.35
## 25    103.25            98.229            108.27
## 30    107.75           101.231            114.27
## 35    110.75           105.243            116.26
## 40    117.25           110.538            123.96
### gráfico dos ICs
plot(c(20,25,30,35,40),medias_idade, ylim=c(min(ics_idades)-1,max(ics_idades)+1), type="p",xlab="Idade",ylab="ICs") 
for (i in 1:5) lines(x =c(v[i],v[i]), y=t(ics_idades)[,i]) # Plotando os ICs

Figura 15.6.A: ICs para médias dos grupos de idades.

Figuras 15.4 e 15.5

boxplot(e3~as.vector(idade_X),col="lightblue3",pch=20)

Figura 15.4: Box plot para os resíduos por nível do fator idade.

boxplot(e3,col = "lightblue3")

Figura 15.5: Box plot para os resíduos de todas as idades.

Capítulo Anterior | Introdução | Próximo Capítulo