Comandos R para análises estatísticas

Capítulo 13: Inferência para duas populações

Exemplo 13.2

maqA<-c(145,127,136,142,141,137)
maqB<-c(143,128,132,138,142,132)
s2_a<-var(maqA)
s2_b<-var(maqB)

Desejamos testar:

  • \(H_0:\sigma_A^2=\sigma_B^2=\sigma^2\)
  • \(H_1:\sigma_A^2\ne\sigma_B^2\)

Utilizamos o comando var.test para comparar duas variâncias amostrais:

teste_13_2<-var.test(x = maqA, y = maqB, conf.level = 0.9, alternative = "two.sided")
teste_13_2
## 
##  F test to compare two variances
## 
## data:  maqA and maqB
## F = 1.08, num df = 5, denom df = 5, p-value = 0.93
## alternative hypothesis: true ratio of variances is not equal to 1
## 90 percent confidence interval:
##  0.21425 5.46474
## sample estimates:
## ratio of variances 
##             1.0821

A opção “two.sided” define a hipótese alternativa \(H_0:\sigma_A^2\ne\sigma_B^2\). Se quiséssemos as alternativas \(H_0:\sigma_A^2 > \sigma_B^2\). ou \(H_0:\sigma_A^2<\sigma_B^2\), utilizaríamos “greater” ou “less”, respectivamente. O valor teste_13_2$statistic nos dá o valor da estatística do teste, neste caso, 1.08206. Como podemos ver, temos 5 graus e liberdade (df) no numerador e 5 no denominador, resultando num valor-p de 0.93315 (variável teste_13_2$p.value). Como este teste avalia a razão \(\sigma_A^2 /\ \sigma_B^2\), rejeitamos a hipótese nula para valores distantes de 1, neste caso, razões menores do que 0.21425 ou maiores do que 5.46474. Estes valores estão armazenados na variável teste_13_2$conf.int.

Para testar igualdade de variâncias ou de médias, porém, não há, no R, funções específicas para tais testes. Desta forma, para os Exemplos 13.3 a 13.5, se desejamos testar um valor específico da estatística W, precisamos escrever a função do teste.

Exemplo 13.3

Vamos escrever uma função para realizar o teste de variâncias de acordo com as estatísticas de teste apresentadas no livro:

teste_var<- function(s1,s2,n1,n2,alpha){
  #s1: Variancia amostral da primeira amostra
  #s2: Variancia amostral da primeira amostra
  #n1: tamanho da primeira amostra
  #n2: tamanho da seunda amostra
  #alpha: Nivel de significância
  
  f_c<-qf(p = alpha/2,n1,n2) # Calculando o quantil da distribruicao F
  print("Teste para igualdade de variâncias")
  print("H0:Sigma_A/Sigma_B=1")
  if(f_c*s1/s2>1 || 1/f_c*s1/s2<1)
    print("Rejeito H0: Variâncias iguais")
  else print("Aceito H0: Variâncias iguais")
print(paste("Intervalo de aceitação: (",round(f_c*s1/s2,5)," ; ",round(1/f_c*s1/s2,5),")",sep=""))
}
teste_var(s1=85,s2=8,n1 = 6,n2=6,alpha = 0.1)
## [1] "Teste para igualdade de variâncias"
## [1] "H0:Sigma_A/Sigma_B=1"
## [1] "Rejeito H0: Variâncias iguais"
## [1] "Intervalo de aceitação: (2.48024 ; 45.51607)"

Para os Exemplos 13.4 3 13.5, utilizaremos a seguinte função para realizar os testes de médias através das estatísticas amostrais:

teste_t<-function(x_barra1,x_barra2,s1,s2,n1,n2,alpha=0.05,H_1="A!=B",var_iguais=TRUE){
  
  #x_barra1,x_barra2 : médias amostrais das amostras 1 e 2 respectivamente
  #s1,s2             : variancias amostrais das amostras 1 e 2 respectivamente
  #n1,n2             : Tamanhos das amostras 1 e 2 respectivamente
  #alpha             : Nivel de Significância
  #H_1               : Hipótese Alternativa
  #var_iguais        : Se TRUE calcula as estatisticas supondo variâncias iguais e desconhecidas.
  #                    Caso contrário, supõe variâncias desiguais e desconhecidas
  
  if(var_iguais==TRUE){ # Caso com variâncias iguais e desconhecidas
  
      s_p=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2)
      t_obs=(x_barra2-x_barra1)/sqrt(s_p*(1/n1+1/n2)) # Calculando estatística do teste
      print("Teste para igualdade de médias: Variâncias iguais e desconhecidas")
    
      if (H_1=="A!=B"){
        t_c=qt(p = 1-alpha/2, df = n1 + n2 - 2) # Calculando o valor para a RC
        if(t_obs < -t_c || t_obs > t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (",round(-t_c,4)," ;",
                    round(t_c,4),")",sep=""))
        }
      else if (H_1=="A<B"){
        t_c=qt(p = 1-alpha, df = n1 + n2 - 2) # Calculando o valor para a RC
        if(t_obs > t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (-infinito ;", round(t_c,4),")",sep=""))
      }
      else if (H_1=="A>B"){
        t_c=qt(p = alpha, df = n1 + n2 - 2) # Calculando o valor para a RC
        if(t_obs < t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (", round(t_c,4),";+infinito)",sep=""))
      }
      else {
        print("'H_1'deve ser 'A!=B', 'A>B' ou 'A<B'")
        return()
      }
      print(paste("Estatística do teste:" , round(t_obs,4)))
  }
  else{ # Caso com variâncias desiguais e desconhecidas
      t_obs=(x_barra2-x_barra1)/sqrt(s1/n1+s2/n2) # Calculando estatística do teste
      gl=round((s1/n1+s2/n2)^2/((s1/n1)^2/(n1-1)+(s2/n2)^2/(n2-1)),0)
      print("Teste para igualdade de médias: Variâncias desiguais e desconhecidas")
    
      if (H_1=="A!=B"){
        t_c=qt(p = 1-alpha/2, df = gl) # Calculando o valor para a RC
        if(t_obs < -t_c || t_obs > t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (",round(-t_c,4)," ;",
                    round(t_c,4),")",sep=""))
        }
      else if (H_1=="A<B"){
        t_c=qt(p = 1-alpha, df = gl) # Calculando o valor para a RC
        if(t_obs > t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (-infinito ;", round(t_c,4),")",sep=""))
      }
      else if (H_1=="A>B"){
        t_c=qt(p = alpha, df = gl) # Calculando o valor para a RC
        if(t_obs < t_c) # Comparando a estatistica do teste com a RC
          print("Rejeito H0: Mu_A = Mu_B")
        else 
          print("Aceito H0: Mu_A = Mu_B")
        print(paste("Intervalo de aceitação: (", round(t_c,4),";+infinito)",sep=""))
      }
      else {
        print("'H_1'deve ser 'A!=B', 'A>B' ou 'A<B'")
        return()
      }
      print(paste("Estatística do teste:" , round(t_obs,4)))
  }
}

Exemplo 13.4

Desejamos testar:

  • \(H_0:\mu_A=\mu_B\)
  • \(H_1:\mu_A>\mu_B\)

Utilizaremos então a função que teste_t programamos anteriormente:

teste_t(x_barra1=68,x_barra2=76,s1=50,s2=52,n1=12,n2=15,alpha=0.05, H_1="A<B",var_iguais=TRUE)
## [1] "Teste para igualdade de médias: Variâncias iguais e desconhecidas"
## [1] "Rejeito H0: Mu_A = Mu_B"
## [1] "Intervalo de aceitação: (-infinito ;1.7081)"
## [1] "Estatística do teste: 2.889"

Para o intervalo de confiança com 95% de confiança, a função seria:

IC_diff_medias<-function(x_barra1,x_barra2,s1,s2,n1,n2,gama=0.95){
  s_p=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2)
  t_c_ic=-qt(p = (1-gama)/2, df = n1 + n2 - 2) # calculando t com a confianca desejada
  t_obs=(x_barra2-x_barra1)/sqrt(s_p*(1/n1+1/n2)) # Calculando estatística do teste
  icl=(x_barra2-x_barra1)-t_c_ic*sqrt(s_p*(1/n1+1/n2)) # limite inferior do intervalo
  icu=(x_barra2-x_barra1)+t_c_ic*sqrt(s_p*(1/n1+1/n2))# limite superior do intervalo
  print(paste("IC(",gama,"):]",round(icl,4),";",round(icu,4),"[",sep=""))
}
IC_diff_medias(x_barra1=68,x_barra2=76,s1=50,s2=52,n1=12,n2=15,gama=0.95)
## [1] "IC(0.95):]2.2969;13.7031["

Note que este IC supõe variâncias iguais e desconhecidas.

Exemplo 13.5

Desejamos testar:

  • \(H_0:\mu_A=\mu_B\)
  • \(H_1:\mu_A\ne\mu_B\)

Mas neste caso, temos que a hipótese alternativa considera \(\mu_A\ne\mu_B\) e também não temos a suposição de igualdade de variâncias. Neste caso, então, utilizando a função teste_t:

teste_t(x_barra1=70.5,x_barra2=84.3,s1=81.6,s2=210.8,
        n1=15,n2=20,alpha=0.05, H_1="A!=B", var_iguais=FALSE)
## [1] "Teste para igualdade de médias: Variâncias desiguais e desconhecidas"
## [1] "Rejeito H0: Mu_A = Mu_B"
## [1] "Intervalo de aceitação: (-2.0369 ;2.0369)"
## [1] "Estatística do teste: 3.4522"

Exemplo 13.7

O comando para calcular a função distribuição acumulada da distribuição de Wilcoxon é pwilcox. Esta distribuição, assim como as outras, aceita os comandos rwilcox, dwilcox e qwilcox. Assim, teremos, para \(\mathbb{P}(W_s\le 87)\):

pwilcox(q=87-10*11/2,m=10,n=10)
## [1] 0.095158

Exemplo 13.8

O teste de Wilcoxon é um teste não paramétrico que utiliza os postos (ordem) dos valores amostrais para concluir sobre as hipóteses de interesse. No R, utilizaremos o comando wilcox.test:

controle<-c(1.3, 1.5, 2.1)
tratamento<-c(1.5,2.5)
pwilcox(q=7.5,m=3,n=2)
## [1] 1
wilcox.test(controle,tratamento, alternative = "less")
## Warning in wilcox.test.default(controle, tratamento, alternative = "less"):
## cannot compute exact p-value with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  controle and tratamento
## W = 1.5, p-value = 0.28
## alternative hypothesis: true location shift is less than 0

Note que o R alerta para a existência de empates nos dados e calcula então o valor-p baseado na aproximação pela distribuição Normal. Ainda assim, o valor aproximado é próximo do valor calculado de 0.3.

Exemplo 13.9

T<-c(0.6,0.63,0.83,0.85,0.91,0.95,1.01,1.03,1.03,1.16,
     1.19,1.2,1.26,1.28,1.3,1.37,1.45,1.54,1.68,2.2)
C<-c(0.52,0.77,0.79,0.79,0.81,0.81,0.89,0.98,1.01,1.18,
     1.19,1.2,1.34,1.36,1.38,1.43,1.64,1.71,2.16,2.25)
hist(T,breaks=seq(0,2.8,0.4),col="lightgray", main="", freq = FALSE,yaxt='n',ylab="")

Figura 13.4: Resistência à remoção, em kg, para o modelo C.

hist(C,breaks=seq(0,2.8,0.4),col="lightgray", main="", freq = FALSE,yaxt='n',ylab="")

Figura 13.5: Resistência à remoção, em kg, para o modelo T

Se quisermos os postos de um vetor, podemos utilizar a função rank(). Assim, se quisermos os postos para a amostra de C e T juntas, fazemos:

media<-(c(C,T))
tipo<-factor(c(rep("C",20),rep("T",20)))
posto<-rank(c(C,T))
tab13_7<-data.frame(media,tipo,posto)
tab13_7<-tab13_7[order(tab13_7$posto),]
Tabela 13.7: Postos para o Exemplo 13.9
media 0.52 0.60 0.63 0.77 0.79 0.79 0.81 0.81 0.83 0.85
tipo C T T C C C C C T T
posto 1.0 2.0 3.0 4.0 5.5 5.5 7.5 7.5 9.0 10.0
media 0.89 0.91 0.95 0.98 1.01 1.01 1.03 1.03 1.16 1.18
tipo C T T C C T T T T C
posto 11.0 12.0 13.0 14.0 15.5 15.5 17.5 17.5 19.0 20.0
media 1.19 1.19 1.20 1.20 1.26 1.28 1.30 1.34 1.36 1.37
tipo C T C T T T T C C T
posto 21.5 21.5 23.5 23.5 25.0 26.0 27.0 28.0 29.0 30.0
media 1.38 1.43 1.45 1.54 1.64 1.68 1.71 2.16 2.20 2.25
tipo C C T T C T C C T C
posto 31 32 33 34 35 36 37 38 39 40

Para o teste utilizaremos novamente a função wilcox.test. Como desta vez temos todos os dados em apenas uma coluna, a entrada para o teste é de outra forma, com a indicação da variável resposta em função do tipo de tratamento:

attach(tab13_7)
wilcox.test(media~tipo, alternative="less")
## Warning in wilcox.test.default(x = c(0.52, 0.77, 0.79, 0.79, 0.81, 0.81, :
## cannot compute exact p-value with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  media by tipo
## W = 196, p-value = 0.47
## alternative hypothesis: true location shift is less than 0

Novamente o resultado foi calculado com a aproximação pela distribuição Normal.

Exemplo 13.10

A<-c(80,72,65,78,85)
B<-c(75,70,60,72,78)
operador<-1:5
tab13_8<-data.frame(operador,A,B)

Para construir o teste t-pareado no R, utilizamos a função t.test com a opção paired=TRUE. Como queremos identificar diferenças médias grandes e positivas, i.e. \(\mu_A>\mu_B\), utilizamos alternative="greater. Então, para testar:

  • \(H_0:\mu_D=0\)
  • \(H_1:\mu_D>0\)

utilizamos:

t.test(A, B, alternative="greater", paired=TRUE, conf.level = 0.90)
## 
##  Paired t-test
## 
## data:  A and B
## t = 5.98, df = 4, p-value = 0.002
## alternative hypothesis: true difference in means is greater than 0
## 90 percent confidence interval:
##  3.7172    Inf
## sample estimates:
## mean of the differences 
##                       5

Exemplo 13.12

Queremos testar:

\(H_0: p_A=p_B\)
\(H_1: p_A\ne p_B\)

Vamos iniciar construindo uma função para realizar tal teste:

p.test<-function(p1,p2,n1,n2,alpha=0.05, H1="p1!=p2"){
  p_c=(n1*p1+n2*p2)/(n1+n2)
  z_obs=(p1-p2)/sqrt(p_c*(1-p_c)*(1/n1+1/n2))
  print("Teste de Comparação de Proporções")
  if(H1=="p1!=p2"){
    z_c=-qnorm(alpha/2)
    if(z_obs< -z_c || z_obs>z_c){
      print("Rejeito H_0!")
      print(paste("Região de Aceitação: ]",round(-z_c,4),";",round(z_c,4),"[",sep=""))
    }
  }
  else if(H1=="p1>p2"){
    z_c=qnorm(1-alpha)
    if(z_obs>z_c){
      print("Rejeito H_0!")
      print(paste("Região de Aceitação: ]-infinito;",round(z_c,4),"[",sep=""))
    }
  }
  else if(H1=="p1<p2"){
    z_c=qnorm(alpha)
    if(z_obs<z_c){
      print("Rejeito H_0!")
      print(paste("Região de Aceitação: ]",round(z_c,4),";+infinito[",sep=""))
    }
  }else {
    print("Erro: H1 deve ser 'p1!=p2', 'p1>p2' ou 'p1<p2' ")
    return()
  }
  print(paste("Z observado:",round(z_obs,4)))
}

p.test(p1=168/400, p2=180/600,n1=400,n2=600,alpha=0.05, H1="p1!=p2")
## [1] "Teste de Comparação de Proporções"
## [1] "Rejeito H_0!"
## [1] "Região de Aceitação: ]-1.96;1.96["
## [1] "Z observado: 3.9028"

Logo, rejeitamos \(H_0\) ao nível de significância 5%. O intervalo de confiança para estes dados pode ser calculado por:

IC_proporcoes<-function(p1,p2,n1,n2,gama=0.95){
  z_c=-qnorm((1-gama)/2)
  icl=(p1-p2)-z_c*sqrt((p1*(1-p1))/n1+(p2*(1-p2))/n2)
  icu=(p1-p2)+z_c*sqrt((p1*(1-p1))/n1+(p2*(1-p2))/n2)
  print(paste("IC(",gama,"):]",round(icl,4),";",round(icu,4),"[",sep=""))
}
IC_proporcoes(p1=168/400, p2=180/600,n1=400,n2=600,gama=0.95)
## [1] "IC(0.95):]0.0593;0.1807["

Para comparar proporções a partir de uma tabela de dupla entrada, podemos utilizar também a função prop.test, com sintaxe similar a que utilizamos nos demais testes do capítulo:

tab13_12A<-matrix(c(168,232,180,420),ncol=2,nrow=2,byrow = TRUE)
prop.test(x=tab13_12A, alternative = "two.sided")
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  tab13_12A
## X-squared = 14.7, df = 1, p-value = 0.00013
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.057221 0.182779
## sample estimates:
## prop 1 prop 2 
##   0.42   0.30

O teste a partir desta função, diferente da que programamos acima, utiliza a distribuição de Qui-quadrado com 1 grau de liberdade. A conclusão também e de rejeição de \(H_0\). Estes testes são equivalentes e veremos no capítulo seguinte mais detalhes sobre testes em distribuições de tabelas de dupla entrada. Os ICs calculados pela função programada e pela função prop.test são equivalentes.

Exemplo Computacional

tab13_12<-data.frame(sujeito=1:26,rbind(
  c(2.18 , 0.43),c(2.05 , 0.08),c(1.05 , 0.18),c(1.95 , 0.78),c(0.28 , 0.03),
  c(2.63 , 0.23),c(1.5  , 0.2 ),c(0.45 , 0   ),c(0.7  , 0.05),c(1.3  , 0.3 ),
  c(1.25 , 0.33),c(0.18 , 0   ),c(3.3  , 0.9 ),c(1.4  , 0.24),c(0.9  , 0.15),
  c(0.58 , 0.1 ),c(2.5  , 0.33),c(2.25 , 0.33),c(1.53 , 0.53),c(1.43 , 0.43),
  c(3.48 , 0.65),c(1.8  , 0.2 ),c(1.5  , 0.25),c(2.55 , 0.15),c(1.3  , 0.05),
  c(2.65 , 0.25)))
names(tab13_12)[2:3]<-c("antes", "depois")
tab13_12$d<-tab13_12$antes-tab13_12$depois
tab13_12$postos<-rank(abs(tab13_12$d))
Tabela 13.12: Índices de placa bacteriana.
Sujeito Antes(x_i) Depois(y_i) d_i=x_i-y_i postos de |d_i|
1 2.18 0.43 1.75 18.0
2 2.05 0.08 1.97 20.0
3 1.05 0.18 0.87 7.0
4 1.95 0.78 1.17 13.0
5 0.28 0.03 0.25 2.0
6 2.63 0.23 2.40 23.5
7 1.50 0.20 1.30 16.0
8 0.45 0.00 0.45 3.0
9 0.70 0.05 0.65 5.0
10 1.30 0.30 1.00 10.0
11 1.25 0.33 0.92 8.0
12 0.18 0.00 0.18 1.0
13 3.30 0.90 2.40 23.5
14 1.40 0.24 1.16 12.0
15 0.90 0.15 0.75 6.0
16 0.58 0.10 0.48 4.0
17 2.50 0.33 2.17 21.0
18 2.25 0.33 1.92 19.0
19 1.53 0.53 1.00 10.0
20 1.43 0.43 1.00 10.0
21 3.48 0.65 2.83 26.0
22 1.80 0.20 1.60 17.0
23 1.50 0.25 1.25 14.5
24 2.55 0.15 2.40 23.5
25 1.30 0.05 1.25 14.5
26 2.65 0.25 2.40 23.5 
xp<-list(tab13_12$antes,tab13_12$depois)
boxplot(xp,pch="-", col="lightblue", border="black", boxwex=0.3, names=c("x","y"))

Figura 13.6: Box-plot para \(x_i\)(antes) e \(y_i\)(depois). R

O teste t-pareado para a comparação da quantidade de placa bacteriana antes e depois do tratamento pode novamente ser realizado com a função t.test:

Quadro 13.1: Teste t-pareado. R.

diff_trats<-t.test(tab13_12$antes,tab13_12$depois, alternative="two.sided", paired=TRUE, conf.level = 0.95)
diff_trats
## 
##  Paired t-test
## 
## data:  tab13_12$antes and tab13_12$depois
## t = 9.29, df = 25, p-value = 1.4e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.0632 1.6691
## sample estimates:
## mean of the differences 
##                  1.3662
stripchart(tab13_12$d,method = "overplot", offset = 2, at=1,
           pch = 1, col="darkred",ylab=NA,
           xlim=c(0,3))
par(new=TRUE)
plot(x=c(0,mean(tab13_12$d)),
     y=c(0.8,0.8),xlim=c(0,3),pch=c(20,20),ylim=c(0,2),ylab="",xlab="", col="darkblue")
par(new=TRUE)
plot(x=c(diff_trats$conf.int[1],diff_trats$conf.int[2]), col="blue",
     y=c(0.8,0.8),xlim=c(0,3),pch=c("|","|"),ylim=c(0,2), type="b",ylab="",xlab="")
     text(x=c(0,mean(tab13_12$d)),y=c(0.7,0.7), labels=c("H_0","Média(d)"), col="lightblue4")

Figura 13.7: Dotplot das diferenças \(d_i\), com intervalo de confiança para \(\mu_D\); Também mostrados \(H_0:\mu_D=0\) e \(\overline{d}=1.366\).

boxplot(tab13_12$d,pch="|", col="lightblue", border="black", boxwex=0.3,horizontal = TRUE )
par(new=TRUE)
plot(x=c(0,mean(tab13_12$d)),
     y=c(0.8,0.8),xlim=c(0,3),pch=c(20,20),ylim=c(0,2),ylab="",xlab="", col="darkblue",xaxt='n', yaxt='n')
par(new=TRUE)
plot(x=c(diff_trats$conf.int[1],diff_trats$conf.int[2]), col="blue",
     y=c(0.8,0.8),xlim=c(0,3),pch=c("|","|"),ylim=c(0,2), type="b",ylab="",xlab="",
     xaxt='n', yaxt='n')
     text(x=c(0,mean(tab13_12$d)),y=c(0.7,0.7), labels=c("H_0","Média(d)"), col="lightblue4")

Figura 13.8: Box plot para as diferenças \(d_i\), com intervalo de confiança para \(\mu_D\); Também mostrados \(H_0:\mu_D=0\) e \(\overline{d}=1.366\).

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