# help(sqrt) #ajuda sobre a função raiz quadrada
# help(lm) # ajuda sobre a função linear models
# citation("ggplot2") # mostra como o pacote deve ser citadoEstatística Básica II
Prefácio
Neste curso, aprenderemos como comparar amostras quantitativas quando os dados são provenientes de distribuições normais. Utilizaremos técnicas paramétricas para a comparação de médias. Em muitos casos, ao realizar o teste de normalidade, os dados não seguem uma distribuição Normal. Assim, realizar um teste sem levar em consideração essa pressuposição pode levar o pesquisador a tomar decisões equivocadas. No entanto, há alternativas, como os testes não paramétricos, que não exigem que as amostras de dados sigam uma distribuição normal. Nesse contexto, podemos também comparar variáveis qualitativas, por exemplo, por meio de testes baseados na distribuição Qui-Quadrado. Com esses testes, também é possível estudar a relação entre variáveis qualitativas. Por fim, iniciaremos o estudo das relações entre variáveis quantitativas com a introdução à regressão linear simples.
1 Unidade I - Testes paramétricos
Nesta unidade vamos entender como comparar duas ou mais amostras a partir do testes t de Student para duas amostras independentes, com variâncias iguais e diferentes e para amostras pareadas (dependentes) a partir do programa R.
1.1 Comandos básicos do R
O
Rvem com as configurações mínimas ao ser instalado. Para realizar tarefas mais específicas pode ser necessário instalar pacotes adicionais (packages), clicando em “pacotes>>instalar pacotes”;É necessário escolher o CRAN (Comprehensive R Archive Network), que são arquivos disponibilizados na rede por diversas IES do mundo inteiro, inclusive do Brasil.
O
Ré case-sensitive, isto é, ele diferencia letras minúsculas de maiúsculas;O separador de casas decimais é o “.”;
A vírgula é usada para separar argumentos;
O prompt do
Ré o sinal de maior “>”;Após instalarmos algum pacote sempre será necessário carregá-lo com o comando: library(“nome do pacote”);
Para verificar como citar o pacote instalado digite: citation(“nome do pacote”);
Para obter ajuda dos arquivos do R use o comando help(nome da função).
Ajuda
- Exemplos:Demonstrações gráficas
- Exemplos:#demo(graphics) #pressione Enter
#demo(persp)
#demo(image)
Note
- O símbolo
#representa um comentário. Podemos usá-lo para oRnão compilar (rodar) o que estiver após este símbolo. - Um pacote amplamente utilizado para contrução de gráficos é o
ggplot.
Continhas no R (modo calculadora)
# soma
2+2-(2+6)+(2)[1] -2# multiplicação
6 + (6*78)[1] 474# Divisão
4/(3+2)[1] 0.8# potência
2**3[1] 8# potência
2^3[1] 82^(3)[1] 81.1.1 Funções do R
nome da função (argumento 1, argumento 2,... )dentro dos parênteses colocamos os argumentos da função.
os argumentos são separados por vírgula.
Exemplos:
mean(...): calcula a média aritmética.
sd(...): calcula o desvio padrão.
lm(...): estima os parâmetros de um modelo de regressão linear.
1.1.2 Criando objetos/vetores com valores numéricos
Vamos crias um vetor de notas de 10 alunos de uma turma. A função length(notas) fornece o número de observações (n) dentro do objeto.
Note
A letra c() antes dos parênteses significa concatenar (colocar junto).
notas = c(9,6,8,5,7,8,6,9,10,6)
print(notas) [1] 9 6 8 5 7 8 6 9 10 6# Criando um vetor com temperaturas diárias (em C)
temperaturas <- c(23.5, 24.8, 21.9, 25.3, 22.7)
temperaturas[1] 23.5 24.8 21.9 25.3 22.7# Criando um dataframe manualmente
# Dados sobre espécies em uma floresta
especies <- data.frame(
Especie = c("Cedro", "Ipê", "Castanheira", "Açaí"),
Altura = c(15.3, 18.2, 25.5, 12.8),
DAP = c(30, 35, 50, 20)
)
especies Especie Altura DAP
1 Cedro 15.3 30
2 Ipê 18.2 35
3 Castanheira 25.5 50
4 Açaí 12.8 20length(notas)[1] 10Para objetos com letras (variáveis qualitativas), basta colocar cada observação entre aspas.
Para listar quais objetos temos salvo utiliza-se a função ls().
Para remover um objeto: rm("nome do objeto").
letras<-c("a","b","c","d")
print(letras)[1] "a" "b" "c" "d"ls()[1] "especies" "letras" "notas" "temperaturas"1.1.3 Condicionantes
Podemos “perguntar” ao R se algo é verdade ou não a partir de operadores simples, como por exemplo:
1<5 # menor[1] TRUE2==2.0001 # igual[1] FALSEsqrt(25) != 6 # diferente[1] TRUEn1 = 8
n2 = 8.1
teste = n1>n2
teste[1] FALSEifelse(teste, "n1 maior que n2","n2 é maior que n1")[1] "n2 é maior que n1"1.1.4 Matrizes
Para inserir uma matriz no R utiliza-se a função matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL). Note que, basta colocar como argumentos o conjunto de dados e indicar o número de linhas e de colunas da matriz.
M1 = matrix(data = 1:9, nrow = 3, ncol = 3)
print(M1) [,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9Observe que a matriz foi organizada por default por colunas. Entretanto, podemos mudar isso para linhas utilizando-se o argumento byrow = TRUE na função matrix.
M2 = matrix(data = 1:9, nrow = 3, ncol = 3, byrow = TRUE)
print(M2) [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9- Índices das matrizes
Exemplo: Uma seleção de 4 firmas de ração de Minas Gerais foi obtida para avaliar a venda de rações. Cada observação bivariada forneceu a quantidade de sacos de ração vendidos e a quantidade de reais de cada venda. Os dados obtidos na forma tabular são:
\[\begin{eqnarray*} \begin{array}{l|cccc}\hline \text{Váriavel 1 (reais/vendas)} & 80 & 120 & 90 & 110 \\ \hline \text{Váriavel 2 (número de sacos de ração vendidos)} & 10 & 12 & 6 & 8\\ \hline \end{array} \end{eqnarray*}\]M_E1 = matrix(c(80,120,90,110,10,12,6,8),4,2)
print(M_E1) [,1] [,2]
[1,] 80 10
[2,] 120 12
[3,] 90 6
[4,] 110 8Pode-se acessar os elementos de uma matriz usando o nome do objeto (matriz) seguido de colchetes. Para acessar o segundo elemento da variável 1, no exemplo anterior, M_E1[2,1]. Para o terceiro elemento da variável 2, M_E1[3,2]. E, para selecionar todos os elementos da variável 1, M_E1[,1].
M_E1[2,1] # segundo elemento da variável 1[1] 120M_E1[3,2] # terceiro elemento da variável 2[1] 6M_E1[,1] # todos os elementos da variável 1[1] 80 120 90 1101.1.5 Manipulando um banco de dados (BD)
A manipulação de um banco de dados a partir o R base é feita de forma análoga às matrizes. A função data.frame tem a mesma lógica de manipulação matricial com mais algumas vantagens.
Vamos estudar alguns recursos simples e importantes dessa função a partir do famoso banco de dados Iris presente na base de dados do R.
Note
A função str mostra que o objeto df_irisé do tipo dataframe, ou seja, banco de dados, e possui 150 observações (linhas) e 5 variáveis (colunas). Tal estrutura é a mesma de uma matriz qualquer.
data("iris") # carregando o banco de dados
df_iris = iris # atribuindo um nome para este BD
str(df_iris) # verificando a estrutura deste objeto'data.frame': 150 obs. of 5 variables:
$ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
$ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
$ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
$ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...Agora vamos fazer algumas seleções com a mesma lógica feita com a matriz do exemplo anterior.
# Selecionando as linhas de 1 a 3 e todas as colunas
df_iris[1:3,] Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa# encontrar linhas onde Species=="setosa"
iris[iris$Species=="setosa",] Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
7 4.6 3.4 1.4 0.3 setosa
8 5.0 3.4 1.5 0.2 setosa
9 4.4 2.9 1.4 0.2 setosa
10 4.9 3.1 1.5 0.1 setosa
11 5.4 3.7 1.5 0.2 setosa
12 4.8 3.4 1.6 0.2 setosa
13 4.8 3.0 1.4 0.1 setosa
14 4.3 3.0 1.1 0.1 setosa
15 5.8 4.0 1.2 0.2 setosa
16 5.7 4.4 1.5 0.4 setosa
17 5.4 3.9 1.3 0.4 setosa
18 5.1 3.5 1.4 0.3 setosa
19 5.7 3.8 1.7 0.3 setosa
20 5.1 3.8 1.5 0.3 setosa
21 5.4 3.4 1.7 0.2 setosa
22 5.1 3.7 1.5 0.4 setosa
23 4.6 3.6 1.0 0.2 setosa
24 5.1 3.3 1.7 0.5 setosa
25 4.8 3.4 1.9 0.2 setosa
26 5.0 3.0 1.6 0.2 setosa
27 5.0 3.4 1.6 0.4 setosa
28 5.2 3.5 1.5 0.2 setosa
29 5.2 3.4 1.4 0.2 setosa
30 4.7 3.2 1.6 0.2 setosa
31 4.8 3.1 1.6 0.2 setosa
32 5.4 3.4 1.5 0.4 setosa
33 5.2 4.1 1.5 0.1 setosa
34 5.5 4.2 1.4 0.2 setosa
35 4.9 3.1 1.5 0.2 setosa
36 5.0 3.2 1.2 0.2 setosa
37 5.5 3.5 1.3 0.2 setosa
38 4.9 3.6 1.4 0.1 setosa
39 4.4 3.0 1.3 0.2 setosa
40 5.1 3.4 1.5 0.2 setosa
41 5.0 3.5 1.3 0.3 setosa
42 4.5 2.3 1.3 0.3 setosa
43 4.4 3.2 1.3 0.2 setosa
44 5.0 3.5 1.6 0.6 setosa
45 5.1 3.8 1.9 0.4 setosa
46 4.8 3.0 1.4 0.3 setosa
47 5.1 3.8 1.6 0.2 setosa
48 4.6 3.2 1.4 0.2 setosa
49 5.3 3.7 1.5 0.2 setosa
50 5.0 3.3 1.4 0.2 setosa# encontrar linhas onde Species=="setosa" e Sepal.Length==5
iris[iris$Species=="setosa" & iris$Sepal.Length==5,] Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5 5 3.6 1.4 0.2 setosa
8 5 3.4 1.5 0.2 setosa
26 5 3.0 1.6 0.2 setosa
27 5 3.4 1.6 0.4 setosa
36 5 3.2 1.2 0.2 setosa
41 5 3.5 1.3 0.3 setosa
44 5 3.5 1.6 0.6 setosa
50 5 3.3 1.4 0.2 setosa1.1.5.0.1 Agrupando os dados
# agrupando os dados apenas pelas variáveis Species, Petal.Width e Petal.Length
iris[, c("Species", "Petal.Width", "Petal.Length")] Species Petal.Width Petal.Length
1 setosa 0.2 1.4
2 setosa 0.2 1.4
3 setosa 0.2 1.3
4 setosa 0.2 1.5
5 setosa 0.2 1.4
6 setosa 0.4 1.7
7 setosa 0.3 1.4
8 setosa 0.2 1.5
9 setosa 0.2 1.4
10 setosa 0.1 1.5
11 setosa 0.2 1.5
12 setosa 0.2 1.6
13 setosa 0.1 1.4
14 setosa 0.1 1.1
15 setosa 0.2 1.2
16 setosa 0.4 1.5
17 setosa 0.4 1.3
18 setosa 0.3 1.4
19 setosa 0.3 1.7
20 setosa 0.3 1.5
21 setosa 0.2 1.7
22 setosa 0.4 1.5
23 setosa 0.2 1.0
24 setosa 0.5 1.7
25 setosa 0.2 1.9
26 setosa 0.2 1.6
27 setosa 0.4 1.6
28 setosa 0.2 1.5
29 setosa 0.2 1.4
30 setosa 0.2 1.6
31 setosa 0.2 1.6
32 setosa 0.4 1.5
33 setosa 0.1 1.5
34 setosa 0.2 1.4
35 setosa 0.2 1.5
36 setosa 0.2 1.2
37 setosa 0.2 1.3
38 setosa 0.1 1.4
39 setosa 0.2 1.3
40 setosa 0.2 1.5
41 setosa 0.3 1.3
42 setosa 0.3 1.3
43 setosa 0.2 1.3
44 setosa 0.6 1.6
45 setosa 0.4 1.9
46 setosa 0.3 1.4
47 setosa 0.2 1.6
48 setosa 0.2 1.4
49 setosa 0.2 1.5
50 setosa 0.2 1.4
51 versicolor 1.4 4.7
52 versicolor 1.5 4.5
53 versicolor 1.5 4.9
54 versicolor 1.3 4.0
55 versicolor 1.5 4.6
56 versicolor 1.3 4.5
57 versicolor 1.6 4.7
58 versicolor 1.0 3.3
59 versicolor 1.3 4.6
60 versicolor 1.4 3.9
61 versicolor 1.0 3.5
62 versicolor 1.5 4.2
63 versicolor 1.0 4.0
64 versicolor 1.4 4.7
65 versicolor 1.3 3.6
66 versicolor 1.4 4.4
67 versicolor 1.5 4.5
68 versicolor 1.0 4.1
69 versicolor 1.5 4.5
70 versicolor 1.1 3.9
71 versicolor 1.8 4.8
72 versicolor 1.3 4.0
73 versicolor 1.5 4.9
74 versicolor 1.2 4.7
75 versicolor 1.3 4.3
76 versicolor 1.4 4.4
77 versicolor 1.4 4.8
78 versicolor 1.7 5.0
79 versicolor 1.5 4.5
80 versicolor 1.0 3.5
81 versicolor 1.1 3.8
82 versicolor 1.0 3.7
83 versicolor 1.2 3.9
84 versicolor 1.6 5.1
85 versicolor 1.5 4.5
86 versicolor 1.6 4.5
87 versicolor 1.5 4.7
88 versicolor 1.3 4.4
89 versicolor 1.3 4.1
90 versicolor 1.3 4.0
91 versicolor 1.2 4.4
92 versicolor 1.4 4.6
93 versicolor 1.2 4.0
94 versicolor 1.0 3.3
95 versicolor 1.3 4.2
96 versicolor 1.2 4.2
97 versicolor 1.3 4.2
98 versicolor 1.3 4.3
99 versicolor 1.1 3.0
100 versicolor 1.3 4.1
101 virginica 2.5 6.0
102 virginica 1.9 5.1
103 virginica 2.1 5.9
104 virginica 1.8 5.6
105 virginica 2.2 5.8
106 virginica 2.1 6.6
107 virginica 1.7 4.5
108 virginica 1.8 6.3
109 virginica 1.8 5.8
110 virginica 2.5 6.1
111 virginica 2.0 5.1
112 virginica 1.9 5.3
113 virginica 2.1 5.5
114 virginica 2.0 5.0
115 virginica 2.4 5.1
116 virginica 2.3 5.3
117 virginica 1.8 5.5
118 virginica 2.2 6.7
119 virginica 2.3 6.9
120 virginica 1.5 5.0
121 virginica 2.3 5.7
122 virginica 2.0 4.9
123 virginica 2.0 6.7
124 virginica 1.8 4.9
125 virginica 2.1 5.7
126 virginica 1.8 6.0
127 virginica 1.8 4.8
128 virginica 1.8 4.9
129 virginica 2.1 5.6
130 virginica 1.6 5.8
131 virginica 1.9 6.1
132 virginica 2.0 6.4
133 virginica 2.2 5.6
134 virginica 1.5 5.1
135 virginica 1.4 5.6
136 virginica 2.3 6.1
137 virginica 2.4 5.6
138 virginica 1.8 5.5
139 virginica 1.8 4.8
140 virginica 2.1 5.4
141 virginica 2.4 5.6
142 virginica 2.3 5.1
143 virginica 1.9 5.1
144 virginica 2.3 5.9
145 virginica 2.5 5.7
146 virginica 2.3 5.2
147 virginica 1.9 5.0
148 virginica 2.0 5.2
149 virginica 2.3 5.4
150 virginica 1.8 5.1# Número de observações por Species
table(iris$Species)
setosa versicolor virginica
50 50 50 # mesmo resultado anterior, mas no formato de um dataframe
data.frame(table(iris$Species)) Var1 Freq
1 setosa 50
2 versicolor 50
3 virginica 501.1.5.1 Agrupando os dados por medidas descritivas
### Calculando a méda de Sepal.Length
mean(iris$Sepal.Length)[1] 5.843333# Calculando a média de Sepal.Length por Species
aggregate(Sepal.Length ~ Species, iris, mean) Species Sepal.Length
1 setosa 5.006
2 versicolor 5.936
3 virginica 6.588# calculando a média Sepal.Length por Species onde Sepal.Width >= 3
aggregate(Sepal.Length ~ Species, iris[iris$Sepal.Width>=3,],mean) Species Sepal.Length
1 setosa 5.029167
2 versicolor 6.218750
3 virginica 6.7689661.1.5.2 Medidas descritivas de uma variável
As medidas de tendência central e de dispersão são facilmente encontradas a partir de funções simples do R base.
# Média aritmética
mean(iris$Sepal.Length)[1] 5.843333# Mediana
median(iris$Sepal.Length)[1] 5.8# Variância
mean(iris$Sepal.Length)[1] 5.843333# Desvio padrão
sd(iris$Sepal.Length)[1] 0.8280661Com a função summary temos um resumo com as principais medidas descritivas de uma variável ou até mesmo de todo banco de dados.
summary(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
Median :5.800 Median :3.000 Median :4.350 Median :1.300
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
Species
setosa :50
versicolor:50
virginica :50
Um bom ótimo livro, gratuito, para mais detalhes sobre intrudução ao R e estatística básica pode ser acessado em [@introd-R]
1.2 Comparação de amostras (testes paramétricos)
1.2.1 Testes de hipóteses
Para [@mont2004] uma hipótese estatística é uma afirmativa sobre os valores dos parâmetros de uma distribuição de probabilidade, seja esta distribuição discreta ou contínua. Por exemplo, se pensamos que o diâmetro interno de uma peça é de 1,5 cm podemos expressar essa afirmativa formalmente como:
\[ H_{0}: \mu=1,5 \] \[ H_{1}: \mu\neq 1,5\]
A afirmativa \(H_{0}\) é denominada de hipótese nula e a \(H_{1}\) de hipótese alternativa. Para testar hipóteses como estas, toma-se uma amostra da população em estudo, calcula-se uma estatística de teste apropriada e, então, rejeita-se ou não a hipótese nula. O conjunto de valores que levam à rejeição de \(H_{0}\) é chamado de região crítica ou região de rejeição do teste.
A hipótese nula sempre contém uma afirmação de igualdade (\(=\), \(\leq\) ou \(\geq\)). Já a alternativa sempre contém desigualdade restrita (\(\neq\), \(<\) ou \(>\)). Há três tipos de testes: unilateral à direita, unilateral à esquerda e bilateral. O valor p, definido a seguir, é calculado com base na hipótese alternativa.
# Teste unilateral à direita
par(mfrow=c(1,2))
x<-seq(-4, 4, by =.01)
y<-dnorm(x, mean=0, sd=1, log = FALSE)
rx<-seq(1.64, 4, by =.1)
ry<-numeric(2*length(rx))
ry[1:length(rx)]<-dnorm(rx, mean=0, sd=1, log = FALSE)
rx<-c(rx, rev(rx))
plot(x, y, 'l', xlab='', ylab='', xaxt="n", yaxt="n",
main='Teste unilateral à direita')
axis(1,at=1.64,labels=expression(paste(q[alpha])),cex.axis=1.3)
polygon(rx, ry, col = "lightblue")
abline(h=0, lty=1)
text(3, 0.1, labels=expression(alpha),cex=1.3)
text(0,0.2,"R.A.",cex=1.3)
text(3,0.2,"R.R.",cex=1.3)
arrows(2.5,0.01,2.8,0.08,length = 0.1)
# Teste unilateral à esquerda
x<-seq(-4, 4, by =.01)
y<-dnorm(x, mean=0, sd=1, log = FALSE)
rx<-seq(-4, -1.64, by =.1)
ry<-numeric(2*length(rx))
ry[1:length(rx)]<-dnorm(rx, mean=0, sd=1, log = FALSE)
rx<-c(rx, rev(rx))
plot(x, y, 'l', xlab='', ylab='', xaxt="n", yaxt="n",
main='Teste unilateral à esquerda')
axis(1,at=-1.64,labels=expression(paste(-q[alpha])),cex.axis=1.3)
polygon(rx, ry, col = "lightblue")
abline(h=0, lty=1)
text(-3, 0.1, labels=expression(alpha),cex=1.3)
text(0,0.2,"R.A.",cex=1.3)
text(-3,0.2,"R.R.",cex=1.3)
arrows(-2.5,0.01,-2.8,0.08,length = 0.1)
# Teste bilateral
x<-seq(-4, 4, by =.01)
y<-dnorm(x, mean=0, sd=1, log = FALSE)
rx<-seq(-4, -1.64, by =.1)
ry<-numeric(2*length(rx))
ry[1:length(rx)]<-dnorm(rx, mean=0, sd=1, log = FALSE)
rx<-c(rx, rev(rx))
plot(x, y, 'l', xlab='', ylab='', xaxt="n", yaxt="n",
main='Teste bilateral')
axis(1,at=-1.64,labels=expression(paste(-q[alpha/2])),cex.axis=1.3)
polygon(rx, ry, col = "lightblue")
abline(h=0, lty=1)
text(-3, 0.1, labels=expression(alpha/2),cex=1.3)
text(0,0.2,"R.A.",cex=1.3)
text(-3,0.2,"R.R.",cex=1.3)
arrows(-2.5,0.01,-2.8,0.08,length = 0.1)
text(3,0.2,"R.R.",cex=1.3)
arrows(2.5,0.01,2.8,0.08,length = 0.1)
rx<-seq(1.64, 4, by =.1)
ry<-numeric(2*length(rx))
ry[1:length(rx)]<-dnorm(rx, mean=0, sd=1, log = FALSE)
rx<-c(rx, rev(rx))
polygon(rx, ry, col = "lightblue")
axis(1,at=1.64,labels=expression(paste(q[alpha/2])),cex.axis=1.3)
text(3, 0.1, labels=expression(alpha/2),cex=1.3)
Dois tipos de erros podem ser cometidos quando testamos hipóteses. Se a hipótese nula é rejeitado quando ela é verdadeira, então dizemos que ocorreu o erro tipo I. Se a hipótese nula não é rejeitada quando ela é falsa, então dizemos que ocorreu o erro tipo II. As probabilidades destes dois tipos de erro são denotadas como:
\[ \alpha=P\{\text{erro tipo I}\}=P\{\text{rejeitar} H_{0}|H_{0} \text{ é} \text{ verdadeira}\}\]
\[\beta=P\{\text{erro tipo II}\}=P\{\text{não rejeitar} H_{0}|H_{0} \text{ é falsa}\} \]
- Nível de significância (\(\alpha\))
Ao testarmos uma hipótese, a probabilidade máxima com que desejamos arriscar um erro do tipo I é chamada de nível de significância do teste. Essa probabilidade, usualmente denotada por \(\alpha\), é em geral, fixada antes da extração das amostras, de modo que os resultados obtidos não influenciam a nossa escolha. Na prática, costuma-se adotar um nível de significância de 0,05 ou 0,01, embora outros valores possam também ser usados. Se por exemplo, ao delinear um teste, se escolhe \(\alpha\)=5%, isto significa que em cerca de 5 chances em 100 rejeitaríamos a hipótese quando ela devesse ser aceita, ou seja, podemos ter 95% de confiança em que tenhamos tomado a decisão correta. Em tal situação dizemos que a hipótese foi rejeitada ao nível de significância de 5%, o que significa que podemos ter errado com uma probabilidade de 5% (Murray e Spiegel, 1978).
Probabilidade de signifiância (\(\hat{\alpha}\))
O método de construção de um teste de hipótese parte da fixação do nível de significância \(\alpha\). Pode-se argumentar que esse procedimento levar à rejeição da hipótese nula para um valor de \(\alpha\) e à não-rejeição para um valor menor. Outra maneira de proceder consiste em apresentar a probabilidade de significância ou nível descritivo ou ainda valor p (p-value) do teste, que é a probabilidade de ocorrer valores da estatística mais extremos do que o observado, sob a hipótese de \(H_{0}\) ser verdadeira.
flowchart LR
A{Valor p} --> B(> que o nível de ignificância α)
B --> C(Não rejeita a hipótese nula)
A --> D(< que o nível de ignificância α)
D --> E[Rejeita a hipótese nula]
1.2.2 Teste t de Student
O nome Student vem do pseudônimo usado pelo Estatístico inglês W. S. Gosset, que introduziu essa distribuição no início do século passado. A distribuição t de Student é importante no que se refere à inferência sobre médias populacionais. A obtenção da densidade está contida no teorema abaixo:
Teorema: Seja Z uma v.a. \(N(0,1)\) e Y uma v.a. \(\chi^2_{(1)}\), com Z e Y independentes. Então, a v.a. \(t=\frac{Z}{\sqrt{Y/\nu}}\), tem densidade dada por:
\[f_{T}(t)=\frac{\Gamma\left(\frac{\upsilon+1}{2}\right)}{\Gamma\left(\frac{\upsilon}{2}\right)}\frac{1}{(\upsilon \pi)^{1/2}}\frac{1}{(1+t^{2}/\upsilon)^{(\upsilon+1)/2}}, \quad -\infty < t <+\infty\] Fluxograma para o teste t de Student
flowchart
A[Teste de Normalidade] --> B(Dados Normais)
B --> C(Não)
B --> D(Sim)
C --> E{Teste não paramétrico}
D --> F(Variáveis dependentes) --> G{Teste t pareado}
D --> H(Variáveis independentes)
H --> I(Variâncias iguais) --> L{Teste t para variâncias iguais}
H --> J(Variâncias diferentes) --> M{Teste t para variâncias diferentes}
1.2.3 Teste T para duas Amostras Independentes (Variâncias Iguais)
Suponha que, ao testar a hipótese de igualdade de variâncias, esta não seja rejeitada, isto é, \(\sigma^2_{1}=\sigma^2_{2}\), porém essa variância comum é desconhecida. Como \(S^2_{1}=S^2_{2}\) são dois estimadores não-viesados de \(\sigma^2\), podemos combiná-los para obter um estimador comum.
\[ S^2_{p}=\frac{(n-1)S^2_{1}+(m-1)S^2_{2}}{n+m-2}=\frac{\displaystyle\sum^n_{i=1}(x_{i}-\bar{x})+\displaystyle\sum^m_{i=i}(y_{i}-\bar{y})}{n+m-2}, \]
que também é um estimador não-viesado de \(\sigma^2\). Mas ainda, cada parcela do numerador acima, quando dividida por \(\sigma^2\), terá distribuição \(\chi^2\), com (n-1) e (m-1) graus de liberdade, respectivamente. Logo, teremos que
\[ \frac{(n+m-2)S^2_{p}}{\sigma^2} \sim \chi^2_{n+m-2}, \]
as hipóteses testadas serão
\[ H_{0}:\mu_{1}=\mu_{2} \\ H_{1}:\mu_{1}\neq \mu_{2}, \]
e as estatística do teste será
\[ T=\frac{\frac{\bar{x}-\bar{y}}{\sigma\sqrt{1/n+1/m}}}{S_{p}/\sigma}=\frac{\bar{x}-\bar{y}}{S_{p}\sqrt{1/n+1/m}}, \]
terá uma distribuição t de Student, com (n+m-2) graus de liberdade,sob a hipótese \(H_{0}\), isto é, se \(\mu_{1}=\mu_{2}\).
Exemplo: Dois técnicos de controle de qualidade mediram o acabamento da superfície de uma parte de metal, obtendo os dados exibidos. Suponha que as medidas sejam normalmente distribuídas. Teste as hipóteses de que as medidas médias do acabamento da superfície feitas pelos dois técnicos são iguais. Use \(\alpha=0,05\).
| Técnico A | Técnico B |
|---|---|
| 1.45 | 1.54 |
| 1.37 | 1.41 |
| 1.21 | 1.56 |
| 1.54 | 1.37 |
| 1.48 | 1.20 |
| 1.29 | 1.31 |
| 1.34 | 1.27 |
| 1.35 |
Testaremos as seguintes hipóteses:
\[ H_{0}: \mu_{A} = \mu_{B} \;(\text{As médias das medidas dos técnicos não diferem}). \] \[ H_{1}: \mu_{A} \neq \mu_{B}\; (\text{As médias das medidas dos técnicos diferem}). \]
### dados para amostrais ##
TA = c(1.45,1.37,1.21,1.54,1.48,1.29,1.34)
TB = c(1.54,1.41,1.56,1.37,1.2,1.31,1.27,1.35)
### testando a normalidade dos dados
# H0: os dados seguem uma distribuição normal;
# H1: os dados não seguem uma distribuição normal.
shapiro.test(TA)
Shapiro-Wilk normality test
data: TA
W = 0.9816, p-value = 0.9669shapiro.test(TB)
Shapiro-Wilk normality test
data: TB
W = 0.9488, p-value = 0.6991# aceitamos a hipótese nula de que os dados são normais ### Comparação de variâcias ###
# H0: variâncias iguais;
# H1: variâncias diferentes
var.test(TA,TB)
F test to compare two variances
data: TA and TB
F = 0.84564, num df = 6, denom df = 7, p-value = 0.854
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.1652094 4.8163185
sample estimates:
ratio of variances
0.8456401 # aceitamos a hipótese nula de que as variâncias são iguais ### testando as hipóteses: ###
# H0: media_A = media_B
# H1: midia_A != media_B
# Obs.: lembrando que o símbolo "!=" significa diferente.
t.test(TA,TB,alternative = "two.sided", var.equal = TRUE)
Two Sample t-test
data: TA and TB
t = 0.10607, df = 13, p-value = 0.9171
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.1279690 0.1411833
sample estimates:
mean of x mean of y
1.382857 1.376250 # valor p > 0,05
# Portanto, não há evidências suficientes contra H0, ou seja, podemos afirmar
# com 95% de confiabilidade que não existe diferença estatisticamente
# significativa entre as medidas dos dois técnicos.1.2.4 Teste T para duas Amostras Independentes (Variâncias diferentes)
Quando a hipótese de igualdade de variâncias for rejeitada, iremos testar as seguintes hipóteses: \(H_{0}:\mu_{1}=\mu_{2}\). Onde a estatística do teste será:
\[ T=\dfrac{\bar{x}-\bar{y}}{\sqrt{S^2_{1}/n+S^2_{2}/m}} \]
Pode-se provar que, sob \(H_{0}\), a v.a. T aproxima-se de uma distribuição t de Student, com o número de graus de liberdade dado aproximadamente por
\[ \nu=\dfrac{(A+B)^2}{A^2/(n+1)+B^2/(m+1)}, \]
no qual \(A=S^2_{1}/n\) e \(B=S^2_{2}/m\). Como esse valor é geralmente fracionário, arredonde para o inteiro mais próximo para obter o número de graus de liberdade.
Exemplo: Para investigar a influência da opção profissional sobre salário inicial de recém-formados, investigaram-se dois grupos de profissionais: um de liberais em geral e outro de formados em Administração de Empresas. Com os resultados abaixo, expressos em salários mínimos, quais seriam suas conclusões?
| Profissionais Liberais |
Profissionais em Administração |
|---|---|
| 6.6 | 8.1 |
| 10.3 | 9.8 |
| 10.8 | 8.7 |
| 12.9 | 10.0 |
| 9.2 | 10.2 |
| 12.3 | 8.2 |
| 7.0 | 8.7 |
| 10.7 |
Testaremos as seguintes hipóteses:
\[ H_{0}: \mu_{Lib} = \mu_{Adm} \;(\text{Não há influência da opção prossional sobre o salário}). \] \[ H_{1}: \mu_{Lib} \neq \mu_{Adm}\; (\text{Há influência da opção prossional sobre o salário}). \]
### dados para amostrais ###
lib = c(6.6,10.3,10.8,12.9,9.2,12.3,7.0)
adm = c(8.1,9.8,8.7,10.0,10.2,8.2,8.7,10.7)
### testando a normalidade dos dados
# H0: os dados seguem uma distribuição normal;
# H1: os dados não seguem uma distribuição normal.
shapiro.test(lib)
Shapiro-Wilk normality test
data: lib
W = 0.93336, p-value = 0.5798shapiro.test(adm)
Shapiro-Wilk normality test
data: adm
W = 0.90519, p-value = 0.3214# Aceitamos a hipótese de normalidade dos dados ### Comparação de variâcias ###
# H0: variâncias iguais;
# H1: variâncias diferentes
var.test(lib,adm)
F test to compare two variances
data: lib and adm
F = 6.0223, num df = 6, denom df = 7, p-value = 0.03257
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.17655 34.29976
sample estimates:
ratio of variances
6.022287 # aceitamos a hipótese ALTERNATIVA de que as variâncias são diferentes ### testando as hipóteses: ###
# H0: media_Lib = media_Adm
# H1: midia_Lib != media_Adm
t.test(lib,adm,alternative = "two.sided", var.equal = FALSE)
Welch Two Sample t-test
data: lib and adm
t = 0.58067, df = 7.7303, p-value = 0.578
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.711745 2.854602
sample estimates:
mean of x mean of y
9.871429 9.300000 # Valor p > 0,05
# Portanto, não há evidências suficientes contra H0, ou seja, podemos afirmar
# com 95% de confiabilidade que não existe diferença estatisticamente
# significativa entre os salários de profissionais liberais e administradores.1.2.5 Teste T para duas Amostras Dependentes (Teste T Pareado)
Suponha que duas amostras \(X_{1},...,X_{n}\) e \(Y_{1},...,Y_{n}\), são pareadas,isto é, podemos considerar uma amostra de pares \((X_{1},Y_{1}),...,(X_{n},Y_{n})\). Se definirmos \(D=X-Y\), teremos a amostra \(D_{1},D_{2},...,D_{n}\), resultante das diferenças entre os valores de cada par. Observe que reduzimos a um problema com uma única população.
Supondo que a v.a. D tem distribuição normal \(N(\mu_{D},\sigma^2_{D})\), podemos deduzir que
\[ \bar{D}=\frac{1}{n}\displaystyle \sum^n_{i=1}D_{i}=\frac{1}{n}\displaystyle \sum^n_{i=1}(X_{i}-Y_{i})=\bar{X}-\bar{Y}, \] terá distribuição \(N(\mu_{D},\sigma^2_{D}/n)\).
Considere \(S^2_{D}=\frac{1}{n-1}\displaystyle \sum^n_{i=1}(D_{i}-\bar{D})^2\)
Logo \(T=\frac{\sqrt{n}(\bar{D}-\mu_{D})}{S_{D}}\) tem distribuição t de Student com (n-1) graus de liberdade.
Como \(\mu_{D}=E(D)=E(X-Y)=E(X)-E(Y)=\mu_{1}-\mu_{2}\), as hipóteses \(H_{0}:\mu_{1}=\mu_{2}; \;H_{1}:\mu_{1}\neq \mu_{2}\) são equivalentes a
\[ H_{0}:\mu_{D}=0; \]
\[ H_{1}:\mu_{D}\neq 0. \]
Exemplo: Uma pesquisa com os membros do Clube do Livro do Mês averiguou se os membros do clube gastam mais tempo vendo televisão do que lendo. Considere que uma amostra de 15 pessoas que responderam à pesquisa forneceu os seguintes dados sobre as horas semanais de ver televisão e as horas semanais de leitura. Usando um nível de significância de 0,05, pode-se concluir que os membros do Clube do Livro do Mês gastam mais horas por semana vendo televisão do que lendo?
| Horas vendo TV | Horas de leitura |
|---|---|
| 10 | 6 |
| 14 | 16 |
| 16 | 8 |
| 18 | 10 |
| 15 | 10 |
| 14 | 8 |
| 10 | 14 |
| 12 | 14 |
| 4 | 7 |
| 8 | 8 |
| 16 | 5 |
| 5 | 10 |
| 8 | 3 |
| 19 | 10 |
| 11 | 6 |
Testaremos as seguintes hipóteses:
\[H_{0}:\mu_{tv} - \mu_{liv} = 0 \leftrightarrow \mu_{tv}=\mu_{liv} \]
\[ H_{1}:\mu_{tv} - \mu_{liv} > 0 \leftrightarrow \mu_{tv} > \mu_{liv} \]
Que são equivalentes a testar:
\[ H_{0}: \mu_{D}=0; \; \text{(A direfença média entre as horas de tv e leitura é zero)} \] \[ H_{1}: \mu_{D} > 0 \; \text{(A direfença média entre as horas de tv e leitura é maior que zero)} \]
### dados amostrais ###
tv = c(10,14,16,18,15,14,10,12,4,8,16,5,8,19,11)
lei = c(6,16,8,10,10,8,14,14,7,8,5,10,3,10,6)
# testando a normalidade dos dados
shapiro.test(tv)
Shapiro-Wilk normality test
data: tv
W = 0.96625, p-value = 0.7991shapiro.test(lei)
Shapiro-Wilk normality test
data: lei
W = 0.95132, p-value = 0.5454# Aceitamos a hipótese de normalidade dos dados. ### testando as hipóteses: ###
# media_D = 0
# media_D > 0
t.test(tv,lei,paired = TRUE)
Paired t-test
data: tv and lei
t = 2.2302, df = 14, p-value = 0.04262
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
0.1148642 5.8851358
sample estimates:
mean difference
3 # Valor p < 0,05
#
# Há evidências contra H0, ou seja, a diferença entre o tempo vendo TV e lendo
# livros é estatisticamente significativa. Portanto, podemos concluir com 95% de
# confiabilidade que os membros do clube estão gastanto mais tempo vendo TV do que lendo.1.2.6 Análise de variância - ANOVA
1.2.6.1 Introdução
A Análise de Variância, ou ANOVA, é uma técnica estatística utilizada para comparar as médias de três ou mais grupos independentes. Ela avalia se existem diferenças significativas entre as médias desses grupos em relação à variação observada. É frequentemente aplicada em experimentos nos quais se deseja entender se a variabilidade observada nas medições é devida a diferenças entre os grupos ou se pode ser atribuída ao acaso. Apesar de calcular variâncias, faz comparação de médias, decompondo a variação total das unidades experimentais em duas partes: variação devida aos tratamentos; e variação devida ao acaso (ou resíduos).
Existem diferentes tipos de ANOVA, sendo os principais:
ANOVA de um fator (One-Way ANOVA): Utilizada quando há um único fator que influencia a variável dependente.
ANOVA de dois fatores (Two-Way ANOVA): Aplicada quando há dois fatores que influenciam a variável dependente, podendo ser ambos considerados de forma independente ou interagindo entre si.
ANOVA de medidas repetidas: Utilizada quando as mesmas unidades experimentais são medidas em diferentes momentos ou sob diferentes condições.
O procedimento envolve a análise das variações dentro de cada grupo e entre os grupos, comparando essas variações com base na hipótese nula de que as médias são iguais. A decisão é baseada no valor p, sendo uma rejeição da hipótese nula indicativa de diferenças significativas entre as médias. Neste texto abordamos apenas a ANOVA de um fator.
1.2.6.2 Alguns onceitos importantes
Antes de apresentar o modelo estatístico, vamos definir alguns conceitos importantes no contexto de uma ANOVA:
Tratamento: é uma condição ou conjunto de condições aplicado a uma unidade experimental no contexto de um experimento. Pode ser uma dose de um medicamento, um nível de temperatura, um método de cultivo de plantas, etc.
Variável: é uma característica ou propriedade que pode ser medida ou observada em diferentes níveis em uma unidade experimental. Variáveis podem ser classificadas como dependentes (resposta) ou independentes (fatores).
Unidade experimental: é a entidade sobre a qual a observação ou medida é realizada. Pode ser um indivíduo, um grupo de indivíduos, uma planta, uma amostra de solo, entre outros. A escolha da unidade experimental depende do contexto do experimento.
Repetição: a repetição envolve a aplicação do mesmo tratamento para múltiplas unidades experimentais. Repetições são usadas para reduzir a variabilidade experimental e melhorar a precisão das estimativas.
Parcela: refere-se a uma subdivisão da área experimental. Em experimentos agrícolas, por exemplo, uma parcela pode ser uma área específica do campo onde um tratamento é aplicado. Parcelas são frequentemente usadas em delineamentos experimentais, como o Delineamento em Blocos Casualizados (DBC).
Bloco: é uma unidade experimental que é usada para controlar a variabilidade de fatores não relacionados ao tratamento principal em estudo. Blocos são frequentemente utilizados em Delineamentos em Blocos Casualizados (DBC).
Fator: é uma variável independente que é manipulada no experimento para observar seu efeito na variável dependente. Pode haver um ou mais fatores em um experimento.
1.2.6.3 ANOVA (modelo)
Considere o modelo estatístico decomposto em duas componentes: sistemática e aleatória:
\[ Y = \mu + e. \]
Em que \(Y\) é a observação associada a uma unidade experimental; \(\mu\) é a média populacional e \(e\) é a parte aleatória. Suponha que queiramos comparar as médias de \(K\) populações, isto é:
\[ H_{0} : \mu_{1} = \mu_{2},...,\mu_{K}; \]
\[ H_{1} : \text{pelo menos uma das médias }\mu_{i} \text{ é diferente das demais}. \]
Para \(K\) amostras independentes, com \(m\) indivíduos em cada uma delas, temos
\[ Y_{ij} = \mu_{i} + e_{ij}, \; i=1,2,...,K; \; j=1,2,...,m. \]
Caso \(H_{0}\) seja verdadeira o modelo acima pode ser escrito como
\[ Y_{ij} = \mu + e_{ij}^{*}, \; i=1,2,...,K; \; j=1,2,...,m. \]
Uma maneira de quantificar a parte aleatória é a partir das somas de quadrado
\[ \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} e_{ij}^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \mu_{i} \right)^{2} \text{e} \;\; \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} (e_{ij}^{*})^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \mu \right)^{2}, \]
em que
\[ \widehat{\mu}_{i} = \dfrac{1}{m} \displaystyle \sum_{j=1}^{m} Y_{ij} = \bar{Y}_{i}, \; i = 1,...,K, \]
e
\[ \widehat{\mu} = \dfrac{1}{mK} \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} Y_{ij} = \bar{Y}. \]
Substituindo estas estimativas nas somas de quadrados, temos a soma de quadrados dentro (SQD) e soma de quadrados total (SQT) dadas respectivamente por
\[ \text{SQD} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \widehat{\mu}_{i} \right)^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \bar{Y}_{i} \right)^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} Y_{ij}^{2} - m\sum_{i=1}^{K}\bar{Y}^{2}_{i}; \]
\[ \text{SQT} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \widehat{\mu} \right)^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} \left( Y_{ij} - \bar{Y} \right)^{2} = \displaystyle \sum_{i=1}^{K}\sum_{j=1}^{m} Y_{ij}^{2} - mK\bar{Y}^{2}. \]
Temos ainda a soma dos quadrados entre (SQE) dada por
\[ \text{SQE} = \text{SQT} - \text{SQD} = \displaystyle m\sum_{i=1}^{K}\left( Y_{i} - \bar{Y} \right)^{2} = m\left( \displaystyle \sum_{i=1}^{K} \bar{Y}_{i}^{2} - K\bar{Y}^{2} \right). \]
Define-se também os quadrados médios total (QMT), dentro (QMD) e entre (QME) da seguinte maneira:
\[ \text{QMT} = \dfrac{\text{SQT}}{Km-1} \]
\[ \text{QMD} = \dfrac{\text{SQD}}{K(m-1)} \]
\[ \text{QME} = \dfrac{\text{SQE}}{m-1}. \]
Definimos então a quantidade
\[ F = \dfrac{\text{QME}}{\text{QMD}}, \]
quanto maior for o valor de \(F\), maior será \(QME\) comparado a \(QMD\) e, assim, maiores as evidências contra \(H_{0}\).
Dadas as seguintes suposições:
os resíduos \(e_{ij} = (y_{ij} - \bar{y}_{i})\) são variáveis aleatórias independentes;
os \(e_{ij} \sim N(0,\sigma^{2}_{a})\);
homocedasticidade das variâncias: todas as \(K\) populações têm variâncias iguais a \(\sigma^{2}\),
pode ser mostrado que
\[ F \sim F(K-1,K(m-1)), \]
isto é, \(F\) tem distribuição Fisher-Snedecor (ou apenas F-Snedecor) com \(K-1\) e \(K(m-1)\) graus de liberdade.
Note que, basicamente, se a distribuição dos erros for Normal, esta hipótese é equivalente a independência dos erros. Portanto, na prática verificamos a normalidade dos dados (\(y_{ij}\)), dos erros (\(e_{ij}\)) e a homoscedasticidade das variâncias (o mais importante).
| Fonte de variação |
Graus de liberdade |
Soma de quadrados |
Quadrado médio |
Estatística F | Valor p |
|---|---|---|---|---|---|
| Entre | (K - 1) | SQE | QME | QME/QMD | |
| Dentro | K(m - 1) | SQD | QMD | ||
| Total | Km - 1 | SQT | - | - | - |
1.2.6.4 Coeficiente de determinação (\(R^{2}\))
O coeficiente de determinação \(R^{2}\) é uma medida estatística que representa a proporção da variabilidade da variável dependente (Y) que é explicada pelos regressores (variáveis independentes, X) em um modelo de regressão linear. Em outras palavras, \(R^{2}\) indica a qualidade do ajuste do modelo aos dados observados. Ele é dado pela razão entre SQE (variação devida aos tratamentos) pela SQT (variação total dos valores observados).
\[ R^{2} = \dfrac{SQE}{SQT}, \]
em que \(0 \leq R^{2} \leq 1\).
Exemplo: A gerência de um depósito que armazena cargas aéreas de pequeno porte está estudando o peso das cargas que chegam ao seu terminal no interior de São Paulo. Usualmente, o terminal recebe 4 tipos de cargas: doméstica (D), administrativa (A), equipamentos industriais (E) e outros tipos (O). Deseja-se verificar se, em média, existem diferenças entre os pesos dos 4 tipos de cargas. Ao longo de um mês, cargas foram colhidas aleatoriamente e seus pesos foram anotados, fornecendo os dados (em kg):
| D | A | E | O |
|---|---|---|---|
| 24.9 | 27.9 | 38.4 | 23.8 |
| 20.4 | 28.1 | 38.6 | 25.3 |
| 24.2 | 28.4 | 41.2 | 23.5 |
| 22.3 | 25.3 | 43.9 | 27.6 |
| 20.3 | 29.3 | 40.2 | 25.5 |
| 24.0 | 28.5 | 40.2 | 23.9 |
| 23.5 | 27.9 | 37.3 | 22.6 |
- Solução
\[ \bar{Y}_{i} = (22,8; 27,9; 40,0; 24,6) \;\;\longrightarrow\text{média de cada tratamento.} \]
\[ \bar{Y} = 28,2 \;\;\longrightarrow \text{média geral.} \]
\[ \text{SQD} = (24,9-22,8)^{2} + (27,9-27,9)^{2} +...+(22,6 - 24,6)^{2} = 75,1 \]
\[ \text{SQT} = (24,9-28,2)^{2} + (27,9-28,2)^{2} +...+(22,6 - 28,2)^{2} = 1329,7 \]
\[ \text{SQE} = 1329,7 - 75,1 = 1254,6 \]
| Fonte de variação |
Graus de liberdade |
Soma de quadrados |
Quadrado médio |
Estatística F | Valor p |
|---|---|---|---|---|---|
| Entre | 3 | 1254,6 | \(\dfrac{1254,6}{3} = 418,2\) | \(\dfrac{418,2}{3,1} = 133,6\) | 0,000 |
| Dentro | 24 | 75,1 | \(\dfrac{75,1}{24}=3,1\) | ||
| Total | 27 | 1329,7 | - | - | - |
- Solução no R
### Dados amostrais ###
# variável resposta
y1 = c(24.9,20.4,24.2,22.3,20.3,24.0,23.5,
27.9,28.1,28.4,25.3,29.3,28.5,27.9,
38.4,38.6,41.2,43.9,40.2,40.2,37.3,
23.8,25.3,23.5,27.6,25.5,23.9,22.6)
# variável explicativa (categórica)
x1 = c("D","D","D","D","D","D","D",
"A","A","A","A","A","A","A",
"E","E","E","E","E","E","E",
"O","O","O","O","O","O","O")
# banco de dados
bd.anova = data.frame(y1,x1)
head(bd.anova,10) y1 x1
1 24.9 D
2 20.4 D
3 24.2 D
4 22.3 D
5 20.3 D
6 24.0 D
7 23.5 D
8 27.9 A
9 28.1 A
10 28.4 AVamos fazer uma breve análise descritiva dos dados. Note que as médias são bem diferentes entre os tratamentos. Sendo que o tipo de carga E tem valores bem mais elevados. Podemos perceber essa discrepância nos dados a pasrtir do boxplot.
# média por tratamento/grupo
aggregate(y1 ~ x1, bd.anova, mean) x1 y1
1 A 27.91429
2 D 22.80000
3 E 39.97143
4 O 24.60000boxplot(y1 ~ x1, bd.anova,
xlab = "Tipo de carga",
ylab = "Peso da carga")
Descritivamente podemos dizer que as médias entre as cargas são diferentes. Contudo, a partir da ANOVA teremos como concluir se estas diferenças são estatisticamente significativas ou não.
Os passos para a ANOVA são muito semelhantes os do teste t de Student. Vejamos.
### Verificando a normalidade dos dados ###
# H0: os dados são normais
# H1: caso contrário.
aggregate(y1 ~ x1, bd.anova, FUN = function(x)shapiro.test(x)$p.value) x1 y1
1 A 0.04503499
2 D 0.23413268
3 E 0.66417825
4 O 0.54853084# Aceitamos a hipótese de normalidade dos dadosApenas os dados da carga A apresentaram o valor-p do teste de normalidade no “limite” do nível de significância de 5%. Contudo, a ANOVA é robusta o suficiente para casos come este. Agora a homogeneidade das variâncias com função barllet.test() usada quando os dados são normais, em que as hipóteses testadas são:
\[ H_{0}: \sigma_{1}^{2} = \sigma_{2}^{2} = ... = \sigma_{n}^{2}; \text{(variâncias homogêneas)}\]
\[ H_{1}: \text{caso contrário (c.c.).}\]
Se os dados não são Normais, podemos fazer o Teste de Levene para homocedasticidade (homogeneidade) das variâncias do pacote .
### homogeneidade das variâncias ###
bartlett.test(y1 ~ x1, bd.anova)
Bartlett test of homogeneity of variances
data: y1 by x1
Bartlett's K-squared = 1.7203, df = 3, p-value = 0.6324# aceitamos a hipótese nula, ou seja, as variâncias são homogêneas.Observados os pressupostos podemos fazer a ANOVA a partir das funções aov(), lm() e anova(). Utilizaremos a função lm (linear models) para observarmos a ANOVA no formato mostrado na tabela da ANOVA.
### Testando a hipótese de igualdade das médias entre os tratamentos ###
# H0: media_D = media_A = media_E = media_O
# H1: ao menos duas das médias são diferentes
anova(lm(y1 ~ x1, bd.anova)) Analysis of Variance Table
Response: y1
Df Sum Sq Mean Sq F value Pr(>F)
x1 3 1254.6 418.19 133.64 4.136e-15 ***
Residuals 24 75.1 3.13
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Valor p < 0,05
# Podemos concluir que há fortes evidências contra H0, ou seja, as médias das
# cargas médias das cargas apresentam diferenças estatisticamente significativas.
# Essa conclusão tem 95% de confiabilidade Como há diferenças significativas entre os tratamentos, podemos seguir com a análise verificando quais tratamentos diferem entre si, dois a dois. Para tanto usa-se o teste de Tukey (Tukey’s honest significant difference) a partir da função TukeyHSD().
Podemos concluir com 95% de confiabildade que, apenas os tratamentos “O” e “D” não diferem significativamente, já que é a única difença com valor p maios que o nível de significância de 5%. Os demais tratamentos, todos apresentam diferenças estatisticamente significativas. Também podemos observar estes resultados no gráfico do teste.
### teste de Tukey
teste.tukey = TukeyHSD(aov(lm(y1 ~ x1, bd.anova)))
print(teste.tukey) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = lm(y1 ~ x1, bd.anova))
$x1
diff lwr upr p adj
D-A -5.114286 -7.7227128 -2.5058586 0.0000827
E-A 12.057143 9.4487158 14.6655699 0.0000000
O-A -3.314286 -5.9227128 -0.7058586 0.0091926
E-D 17.171429 14.5630015 19.7798557 0.0000000
O-D 1.800000 -0.8084271 4.4084271 0.2531057
O-E -15.371429 -17.9798557 -12.7630015 0.0000000plot(teste.tukey)
1.2.6.5 Exercícios
- Para avaliar dois tipos de método de ensino em matemática foram aplicados testes em duas turmas. A partir destes dados, que representam as notas dos alunos, pode-se afirmar se há diferenças entre estes métodos? Use \(\alpha\) = 0,05.
notas = c(8.2,8.5,8.6,9.5,8.7,8.6,7.4,8.6,8.1,7.6,7.3,6.8,7.9,6.9,8.5,8.9,7.1,8.0,
6.8,6.7,5.4,9.3,6.3,6.6,5.9,6.5,7.5,6.0,8.0,6.9,6.3,7.9,8.3,5.5,5.5,7.1,9.1,6.7,8.0,
6.7,9.7,9.2,8.2,5.5,7.8,7.6,7.3,8.2,8.6,6.6,6.5,5.0,6.2,9.8,9.0,8.1,6.0,7.5,8.1,8.7,
9.5,8.2,5.2,7.3,7.8,7.3,9.6,6.7,10.0,7.1,6.6,5.2,7.9,6.0,7.9,6.2,6.3,4.9,6.2,6.0)
met =rep(c("A","B"),c(40,40))
bd.notas = data.frame(notas,met)
# knitr::kable(bd.notas, col.names = c("Notas","Método"), align = "c")
knitr::kable(read.csv("C:/Users/denni/OneDrive/Documentos/Turmas 2023.1-2024/turmas2024/notas_teste.csv", dec = ","),align = "c")| Notas | Método |
|---|---|
| 8.2 | A |
| 8.5 | A |
| 8.6 | A |
| 9.5 | A |
| 8.7 | A |
| 8.6 | A |
| 7.4 | A |
| 8.6 | A |
| 8.1 | A |
| 7.6 | A |
| 7.3 | A |
| 6.8 | A |
| 7.9 | A |
| 6.9 | A |
| 8.5 | A |
| 8.9 | A |
| 7.1 | A |
| 8.0 | A |
| 6.8 | A |
| 6.7 | A |
| 5.4 | A |
| 9.3 | A |
| 6.3 | A |
| 6.6 | A |
| 5.9 | A |
| 6.5 | A |
| 7.5 | A |
| 6.0 | A |
| 8.0 | A |
| 6.9 | A |
| 6.3 | A |
| 7.9 | A |
| 8.3 | A |
| 5.5 | A |
| 5.5 | A |
| 7.1 | A |
| 9.1 | A |
| 6.7 | A |
| 8.0 | A |
| 6.7 | A |
| 9.7 | B |
| 9.2 | B |
| 8.2 | B |
| 5.5 | B |
| 7.8 | B |
| 7.6 | B |
| 7.3 | B |
| 8.2 | B |
| 8.6 | B |
| 6.6 | B |
| 6.5 | B |
| 5.0 | B |
| 6.2 | B |
| 9.8 | B |
| 9.0 | B |
| 8.1 | B |
| 6.0 | B |
| 7.5 | B |
| 8.1 | B |
| 8.7 | B |
| 9.5 | B |
| 8.2 | B |
| 5.2 | B |
| 7.3 | B |
| 7.8 | B |
| 7.3 | B |
| 9.6 | B |
| 6.7 | B |
| 10.0 | B |
| 7.1 | B |
| 6.6 | B |
| 5.2 | B |
| 7.9 | B |
| 6.0 | B |
| 7.9 | B |
| 6.2 | B |
| 6.3 | B |
| 4.9 | B |
| 6.2 | B |
| 6.0 | B |
- Medidas de um índice de placa bacteriana obtidas de 26 crianças em idade pré-escolar, antes e depois do uso de uma escova experimental de uma escova convencional. Verifique, ao nível de significânca de 1%, se há diferença os tipos de escovas. E quantos ao sexo, há diferença no índice de placa bacteriana? (Bussab e Morettin, 2017.)
db.placa = read.csv("C:/Users/denni/OneDrive/Documentos/Turmas 2023.1-2024/placa_bacte.csv", dec = ",")
knitr::kable(db.placa, align = "c")| Sexo | Antes_Hugger | Depois_Hugger | Antes_convencional | Depois_convencional |
|---|---|---|---|---|
| F | 2.18 | 0.43 | 1.20 | 0.75 |
| F | 2.05 | 0.08 | 1.43 | 0.55 |
| F | 1.05 | 0.18 | 0.68 | 0.08 |
| F | 1.95 | 0.78 | 1.45 | 0.75 |
| F | 0.28 | 0.03 | 0.50 | 0.05 |
| F | 2.63 | 0.23 | 2.75 | 1.60 |
| F | 1.50 | 0.20 | 1.25 | 0.65 |
| F | 0.45 | 0.00 | 0.40 | 0.13 |
| F | 0.70 | 0.05 | 1.18 | 0.83 |
| F | 1.30 | 0.30 | 1.43 | 0.58 |
| F | 1.25 | 0.33 | 0.45 | 0.38 |
| F | 0.18 | 0.00 | 1.60 | 0.63 |
| F | 3.30 | 0.90 | 0.25 | 0.25 |
| F | 1.40 | 0.24 | 2.98 | 1.03 |
| M | 0.90 | 0.15 | 3.35 | 1.58 |
| M | 0.58 | 0.10 | 1.50 | 0.20 |
| M | 2.50 | 0.33 | 4.08 | 1.88 |
| M | 2.25 | 0.33 | 3.15 | 2.00 |
| M | 1.53 | 0.53 | 0.90 | 0.25 |
| M | 1.43 | 0.43 | 1.78 | 0.18 |
| M | 3.48 | 0.65 | 3.50 | 0.85 |
| M | 1.80 | 0.20 | 2.50 | 1.15 |
| M | 1.50 | 0.25 | 2.18 | 0.93 |
| M | 2.55 | 0.15 | 2.68 | 1.05 |
| M | 1.30 | 0.05 | 2.73 | 0.85 |
| M | 2.65 | 0.25 | 3.43 | 0.88 |
- PlantGrowth: Este conjunto de dados está no
datasetsdorbasee contém o resultado de um experimento que investiga o efeito de diferentes tratamentos (ctrl, trt1, trt2) no crescimento de plantas. A variável de resposta é o comprimento das plantas. Você diria que há diferenças significativas nos comprimentos das plantas entre os diferentes tratamentos. Use \(alpha\) = 0,05.
db.plantas = PlantGrowth
knitr::kable(db.plantas,align = "c")| weight | group |
|---|---|
| 4.17 | ctrl |
| 5.58 | ctrl |
| 5.18 | ctrl |
| 6.11 | ctrl |
| 4.50 | ctrl |
| 4.61 | ctrl |
| 5.17 | ctrl |
| 4.53 | ctrl |
| 5.33 | ctrl |
| 5.14 | ctrl |
| 4.81 | trt1 |
| 4.17 | trt1 |
| 4.41 | trt1 |
| 3.59 | trt1 |
| 5.87 | trt1 |
| 3.83 | trt1 |
| 6.03 | trt1 |
| 4.89 | trt1 |
| 4.32 | trt1 |
| 4.69 | trt1 |
| 6.31 | trt2 |
| 5.12 | trt2 |
| 5.54 | trt2 |
| 5.50 | trt2 |
| 5.37 | trt2 |
| 5.29 | trt2 |
| 4.92 | trt2 |
| 6.15 | trt2 |
| 5.80 | trt2 |
| 5.26 | trt2 |
2 Unidade II - Teste não paramétricos
2.1 Comparação de amostras (testes não paramétricos)
2.1.1 Teste de Wilcoxon (Mann-Whitney)
A ideia básica do teste de Wilcoxon é verificar se duas amostras independentes têm a mesma distribuição de valores. Ele faz isso atribuindo uma classificação a cada valor em ambas as amostras, combinando-as em uma única lista ordenada e comparando as somas das classificações das duas amostras.
Há uma pequena diferença no contexto em que os nomes são usados:
Teste de Wilcoxon (signed-ranks): é comumente utilizado quando as duas amostras têm o mesmo tamanho e são pareadas. Neste contexto, ele é usado para comparar as diferenças entre os pares de observações em vez das observações individuais. Também é chamado de teste de Wilcoxon para amostras pareadas.
Teste de Mann-Whitney (rank sum): é comumente utilizado quando as duas amostras têm tamanhos diferentes ou não são pareadas. Neste contexto, ele é usado para comparar as distribuições de valores das duas amostras independentes. Também é chamado de teste de Wilcoxon-Mann-Whitney.
Exemplo: Avaliar se a concentração de CO2 em uma floresta aumenta após a construção de uma estrada.
- Concentração de CO2 antes da estrada (mg/m³): 400, 420, 410, 390, 430
- Concentração de CO2 antes da estrada (mg/m³): 440, 450, 430, 420, 460
\(H_{0}:\) A concentração de CO2 não aumenta significativamente após a construção de uma estrada.
\(H_{1}:\) A concentração de CO2 aumenta significativamente após a construção de uma estrada.
co2_antes <- c(400, 420, 410, 390, 430)
co2_depois <- c(440, 450, 430, 420, 460)
wilcox.test(co2_antes, co2_depois,paired = TRUE, alternative = "less",
exact = F)
Wilcoxon signed rank test with continuity correction
data: co2_antes and co2_depois
V = 0, p-value = 0.02724
alternative hypothesis: true location shift is less than 0O valor-p de 0,02724 indica que há um aumento significativa na concentração de CO2 depois da construção da estrada (nível de significância de 0,05).
Exemplo: Comparar a temperatura da água em dois rios, um com e outro sem floresta ciliar.
- Temperatura da água no rio com floresta ciliar (\(^{o}\)C): 18, 20, 21, 19, 22
- Temperatura da água no rio sem floresta ciliar (\(^{o}\)C): 24, 25, 26, 23, 27
\(H_{0}:\) Não há diferença significativa na temperatura entre os grupos.
\(H_{1}:\) Há diferença significativa na temperatura entre os grupos.
com_floresta <- c(18, 20, 21, 19, 22)
sem_floresta <- c(24, 25, 26, 23, 27)
wilcox.test(com_floresta, sem_floresta, alternative = "less")
Wilcoxon rank sum exact test
data: com_floresta and sem_floresta
W = 0, p-value = 0.003968
alternative hypothesis: true location shift is less than 0mean(com_floresta)[1] 20mean(sem_floresta)[1] 25O valor-p de 0,003968 indica que existe uma diferença significativa na temperatura da água entre os rios com e sem floresta ciliar (nível de significância de 0,05). A temperatura da água é mais alta no rio sem floresta ciliar.
Note
Até o momento, os dados (variáveis) foram inseridos individualmente no R como objetos (vetores). Quando esses dados são organizados no formato longo, ou seja, empilhados um sobre o outro, com um vetor numérico representando os dados e outro vetor qualitativo indicando a categoria correspondente a cada dado, podemos aplicar o teste de Wilcoxon seguindo este formato: wilcox.test(y ~ x, paired = TRUE, data = dados). Em que y é a variável resposta e x a categórica.
2.1.2 Teste de Kruscal-Wallis
O teste de Kruskal-Wallis é um teste estatístico que permite comparar mais de duas amostras independentes de maneira simples e sem precisar se preocupar com a normalidade dos dados. Ele utiliza a classificação dos dados para fazer a comparação. Em outras palavras, ele coloca os dados em ordem e verifica se a posição média dos valores em cada grupo é diferente.
Teste de Kruskal-Wallis vs. ANOVA:
| Característica | Teste de Kruskal-Wallis | ANOVA |
|---|---|---|
| Distribuição dos dados | Não exige normalidade | Exige normalidade |
| Número de grupos | Mais de duas | Mais de duas |
| Tipo de dados | Quantitativos | Quantitativos |
| Complexidade | Simples | Mais complexo |
| Precisão | Menos preciso para dados normalmente distribuídos | Mais preciso para dados normalmente distribuídos |
Exemplo: Comparar a largura da pétala das três espécies de flor do bando de dados iris.
\(H_{0}:\) a largura da pétala das três espécies de flor não difere de forma significativa.
\(H_{1}:\) a largura da pétala das três espécies de flor difere de forma significativa.
dados.iris = iris
aggregate(iris$Petal.Width ~ iris$Species,
FUN =function(x)shapiro.test(x)$p.value) iris$Species iris$Petal.Width
1 setosa 8.658573e-07
2 versicolor 2.727780e-02
3 virginica 8.695419e-02 kruskal.test(iris$Petal.Width ~ iris$Species, data = iris)
Kruskal-Wallis rank sum test
data: iris$Petal.Width by iris$Species
Kruskal-Wallis chi-squared = 131.19, df = 2, p-value < 2.2e-16Como o valor-p é menor que 0,05 podemos conlcuir com 95% de confiabilidade que as médias da largura da pétala das três espécies de flor aparesentam diferenças significativas.
Note
O teste de Kruskal-Wallis aceita os seguintes formatos:
kruskal.test(x1,x2,x3)kruskal.test(y,g = grupos, data = dados)kruskal.test(y ~ grupos, data = dados)
Teste post-hoc
Assim como na Anova paramétrica, a partir do teste de Tukey, existem testes post-hoc (“depois disto”) nos testes não paramétricos. Para o teste de Kruskal-Wallis o rbase já trás o teste post-hoc a partir da função pairwise.wilcox.test. Esta função realiza testes de Wilcoxon pareados entre pares de grupos em um conjunto de dados. A função fornece várias opções para ajustar os valores p para comparações múltiplas, cada uma com suas vantagens e desvantagens.
- Ajuste de Bonferroni:
- O método mais simples e conservador.
- Controla a taxa de erro familiar (FWER) em um nível estrito.
- Pode ser muito conservador e levar à perda de poder estatístico.
- Ajuste de Holm:
- Menos conservador que o Bonferroni.
- Controla a FWER de forma mais flexível.
- Pode ser mais poderoso que o Bonferroni, mas ainda pode ser conservador em alguns casos.
- Ajuste de Hochberg:
- Semelhante ao Holm, mas com maior poder estatístico.
- Controla a FWER de forma adaptativa.
- Pode ser mais poderoso que o Holm, mas pode ter um FWER ligeiramente maior.
- Ajuste de Benjamini-Hochberg:
- Controla a taxa de erro falso por descoberta (FDR).
- Permite um número maior de comparações significativas do que os métodos FWER.
- Pode ser mais poderoso que os métodos FWER, mas pode ter um FDR maior.
- Ajuste de Benjamini-Yekutieli:
- Semelhante ao Benjamini-Hochberg, mas mais poderoso em alguns casos.
- Controla o FDR de forma mais flexível.
- Pode ser mais poderoso que o Benjamini-Hochberg, mas pode ter um FDR ligeiramente maior.
Considerações:
A escolha do método de ajuste depende do contexto e das prioridades do estudo. Métodos FWER são mais conservadores e controlam o erro tipo I de forma estrita. Métodos FDR são mais poderosos, mas podem ter um maior risco de erros falsos positivos. É importante considerar o poder estatístico e o controle de erros ao escolher um método.
Recomendações:
- Para estudos com muitas comparações, os métodos FDR podem ser mais adequados.
- Para estudos com poucas comparações, os métodos FWER podem ser mais apropriados.
- É importante consultar a literatura e consultar um estatístico para escolher o método de ajuste mais adequado para cada estudo.
Como podemos observar, as três variáveis apresentam diferenças significativas na comparação 2 a 2 (valor p < 0,01). Observe que, neste caso, não houve diferença entre os métodos.
pairwise.wilcox.test(iris$Petal.Width,iris$Species,
p.adjust.method = "bonferroni")
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: iris$Petal.Width and iris$Species
setosa versicolor
versicolor < 2e-16 -
virginica < 2e-16 2.9e-16
P value adjustment method: bonferroni pairwise.wilcox.test(iris$Petal.Width,iris$Species,
p.adjust.method = "holm")
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: iris$Petal.Width and iris$Species
setosa versicolor
versicolor <2e-16 -
virginica <2e-16 <2e-16
P value adjustment method: holm pairwise.wilcox.test(iris$Petal.Width,iris$Species,
p.adjust.method = "hochberg")
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: iris$Petal.Width and iris$Species
setosa versicolor
versicolor <2e-16 -
virginica <2e-16 <2e-16
P value adjustment method: hochberg pairwise.wilcox.test(iris$Petal.Width,iris$Species,
p.adjust.method = "fdr")
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: iris$Petal.Width and iris$Species
setosa versicolor
versicolor <2e-16 -
virginica <2e-16 <2e-16
P value adjustment method: fdr 2.1.3 Teste de Friedman
O teste de Friedman é um teste estatístico não paramétrico usado para comparar várias medidas coletadas do mesmo grupo de indivíduos ao longo do tempo ou em diferentes condições. Ele utiliza a classificação dos dados para fazer a comparação. Em outras palavras, ele coloca os dados em ordem e verifica se a posição média dos valores em cada momento ou condição é diferente.
Exemplo: Em um estudo sobre emoções foram estimuladas (em ordem aleatória) em oito pacientes. As emoções analisadas foram: medo, felicidade, trizteza e calma A tabela a seguir apresenta as medidas resultantes do potencial elétrico em minivolts da pele dos pacientes.
| Medo | Felicidade | Depressão | Calma | |
|---|---|---|---|---|
| Paciente 1 | 23,1 | 22,7 | 22,5 | 22,6 |
| Paciente 2 | 57,6 | 53,2 | 53,7 | 53,1 |
| Paciente 3 | 10,5 | 9,7 | 10,8 | 8,3 |
| Paciente 4 | 23,6 | 19,6 | 21,1 | 21,6 |
| Paciente 5 | 11,9 | 13,8 | 13,7 | 13,3 |
| Paciente 6 | 54,6 | 47,1 | 39,2 | 37 |
| Paciente 7 | 21 | 13,6 | 13,7 | 14,8 |
| Paciente 8 | 20,3 | 23,6 | 16,3 | 14,8 |
\(H_{0}:\) Os dados se comportam de forma similar.
\(H_{1}:\) Pelo menos um dos grupos difere dos demais.
emocao = rep(c('medo','felicidade','depressão','calma'),8)
paciente = rep(c(1:8),each=4)
resposta=c(23.1,22.7,22.5,22.6,57.6,53.2,53.7,53.1,10.5,9.7,10.8,
8.3,23.6,19.6,21.1,21.6,11.9,13.8,13.7,13.3,54.6,47.1,
39.2,37.0,21.0,13.6,13.7,14.8,20.3,23.6,16.3,14.8)
friedman.test(resposta, groups=emocao, blocks=paciente)
Friedman rank sum test
data: resposta, emocao and paciente
Friedman chi-squared = 6.45, df = 3, p-value = 0.09166Portanto, concluímos ao nível de significância de 5% que não há evidências suficientes contra \(H_{0}\), ou seja, os grupos não apresentam diferenças estatisticamente significativas nos potenciais elétricos das emoções analisadas.
Exemplo: Numa avaliação do desempenho do veículo, seis motoristas profissionais (rotulados I, II, III, IV, V, VI) avaliaram três carros (A, B e C) em ordem aleatória. As suas notas dizem respeito apenas ao desempenho dos veículos e supostamente não são influenciadas pela marca do veículo ou por informações exógenas semelhantes. Aqui estão suas classificações na escala de 1 a 10:
| Car | I | II | III | IV | V | VI |
|---|---|---|---|---|---|---|
| A | 7 | 6 | 6 | 7 | 7 | 8 |
| B | 8 | 10 | 8 | 9 | 10 | 8 |
| C | 9 | 7 | 8 | 8 | 9 | 9 |
*[@kvam2022]
Considere que os motoristas são os blocos e o carros o tratamentos.
\(H_{0}:\) As notas não sofrem influência pela marca do veículo.
\(H_{1}:\) As notas sofrem influência por ao menos uma das marcas dos veículos.
Ao nível de significância de 5%, podemos concluir que há evidências contra \(H_{0}\), ou seja, as notas são influenciadas pela marca do veículo (valor p < 0,05).
bd.car <- matrix(c(7,8,9,6,10,7,6,8,8,7,9,8,7,10,9,8,8,9),nrow=6,
byrow=TRUE,dimnames=(list(1:6,c("A","B","C"))))
friedman.test(bd.car)
Friedman rank sum test
data: bd.car
Friedman chi-squared = 8.2727, df = 2, p-value = 0.01598carro = rep(c('A','B','C'),6)
motorista = rep(c(1:6),each=3)
resp.carro =c(7,8,9,6,10,7,6,8,8,7,9,8,7,10,9,8,8,9)
friedman.test(resp.carro,groups = carro,blocks = motorista)
Friedman rank sum test
data: resp.carro, carro and motorista
Friedman chi-squared = 8.2727, df = 2, p-value = 0.01598
Note
O teste de Friedman aceita ainda o seguinte formato:
friedman.test(resp.carro ~ carro|motorista)
Teste post-hoc
Para o teste post-hoc de Friedman, podemos usar os pacotes PMCRM e PMCMRplus. Como estes pacotes não fazem parte do rbase é necessário fazer a instalação a partir da função install.packages(nome do pacote). Feito isso, carregamos o pacote com a função library(nome do pacote).
Note que os resultados também são semelhantes. Todavia, podemos ter conclusões diferentes a depender de qual teste utilizar, isso fica a critério do pesquisador. Neste caso, resolvemos usar os resultados dos métodos “hommel” e “fdr” dadas as características descritivas os dados (média de A=6,8;média de B=8,8; média de C=8,3)
# install.packages("PMCRMplus")
library(PMCMRplus)frdAllPairsSiegelTest(bd.car,p.adjust = "bonferroni")
Pairwise comparisons using Siegel-Castellan all-pairs test for a two-way balanced complete block designdata: y A B
B 0.042 -
C 0.063 1.000
P value adjustment method: bonferronipairwise.wilcox.test(resp.carro,g = carro,p.adjust.method="bonferroni",exact=F)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: resp.carro and carro
A B
B 0.025 -
C 0.061 1.000
P value adjustment method: bonferroni pairwise.wilcox.test(resp.carro,g =carro, p.adjust.method = "hommel",exact = F)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: resp.carro and carro
A B
B 0.025 -
C 0.041 0.498
P value adjustment method: hommel pairwise.wilcox.test(resp.carro,g = carro, p.adjust.method = "fdr", exact = F)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: resp.carro and carro
A B
B 0.025 -
C 0.030 0.498
P value adjustment method: fdr 2.1.4 Testes Qui-quadrado (\(\chi^{2}\))
O teste Qui-quadrado é um teste estatístico não paramétrico utilizado para avaliar a associação entre duas variáveis categóricas. Ele compara a distribuição de frequências observadas em uma tabela de contingência com a distribuição de frequências esperadas se não houvesse associação entre as variáveis.
O teste Qui-Quadrado pode ser utilizado em diferentes contextos, como:
- Teste de Independência: Aqui, o objetivo é determinar se duas variáveis categóricas são independentes uma da outra. Por exemplo, podemos querer saber se existe uma relação entre gênero e preferência por um determinado tipo de esporte.
\(H_{0}\): As variáveis são independentes uma da outra (não há relação).
\(H_{1}\): As variáveis não são independentes uma da outra (há relação).
- Teste de Adequação (aderência): Neste caso, o teste compara a distribuição observada de uma variável categórica com uma distribuição teórica conhecida (por exemplo, uma distribuição uniforme, normal, etc.). Por exemplo, podemos querer testar se os resultados de um dado são justos, ou seja, se seguem uma distribuição uniforme.
\(H_{0}\): A distribuição observada da variável categórica é compatível com a distribuição teórica conhecida.
\(H_{1}\): A distribuição observada da variável categórica não é compatível com a distribuição teórica conhecida.
- Teste de Homogeneidade: Este teste é usado para comparar a distribuição de uma variável categórica entre dois ou mais grupos. Por exemplo, podemos querer determinar se a preferência por uma marca de refrigerante é a mesma em diferentes regiões geográficas.
\(H_{0}\): A distribuição da variável é a mesma em todos os grupos.
\(H_{1}\): A distribuição da variável categórica difere entre pelo menos dois dos grupos.
Cálculo da Estatística Qui-Quadrado:
\[ \chi^{2} = \sum_{i=1}^{n}\dfrac{(O_{i} - E_{i})^{2}}{E_{i}} \]
em que
\(O_{i}\): i-ésima frequência observada.
\(E_{i}\): i-ésima frequência esperada.
2.1.4.1 Cálculo das frequências esperadas
Teste de Independência:
No caso de um teste de independência, a frequência esperada é calculada assumindo que não há relação entre as variáveis. A frequência esperada é obtida multiplicando as marginais correspondentes para cada célula da tabela de contingência e dividindo pelo total de observações. Se temos uma tabela de contingência com linhas representando uma variável e colunas representando outra, então a frequência esperada \(E_{ij}\) para a célula (\(i,j\))é dada por:
\[ E_{ij} = \dfrac{\text{Total da linha i x Total da coluna j}}{\text{Total de observações}} \] Teste de Adequação:
No teste de adequação, a frequência esperada é calculada com base em uma distribuição de probabilidade teórica. Por exemplo, se estamos testando se os dados seguem uma distribuição uniforme, a frequência esperada para cada categoria é igualmente distribuída, ou seja, todas as categorias têm a mesma probabilidade. Portanto, a frequência esperada para cada categoria é simplesmente o total de observações dividido pelo número de categorias. Se você está testando se a distribuição de cores de flores em um jardim segue uma distribuição binomial, você precisa calcular as frequências esperadas para cada cor usando a fórmula da distribuição binomial.
Em geral, se a probabilidade para cada categoria é constante, podemos calcular as frequências esperadas da seguinte forma:
\[ E_{i} = n \times p_{i}, \]
em que \(n\) é o número de tentativas (tamanho da amostra) e \(p_{i}\) é a probabilidade assumida da i-ésima categoria.
Teste de Homogeneidade:
No teste de homogeneidade, a frequência esperada é calculada assumindo que as proporções são as mesmas em todas as categorias. Isso significa que a frequência esperada para cada categoria é a proporção dessa categoria no total de observações, multiplicada pelo total em cada grupo. Se temos uma tabela de contingência com linhas representando grupos e colunas representando categorias, então a frequência esperada \(E_{ij}\) para a célula (\(i,j\))é dada por:
\[ E_{ij} = \dfrac{\text{Total da linha i x Total da coluna j}}{\text{Total de observações}} \] Pressupostos
De uma forma geral podemos destacar os seguintes pressupostos para aplicação de um teste Qui-quadrado:
- Nível de mensuração em escala nominal.
- Tamanho da amostra (n) maior que 20 e frequências esperadas maiores que 5 no caso de tabelas 2x2.
- Se a tabela não for 2x2, por exemplo, 3x2, 3x3, 4x2 e assim por diantes, o número de células com freqüência esperada inferior a 5 deve ser menos de 20% do total de células.
Exemplo (independência): Queremos investigar se existe uma relação entre o nível de educação e a preferência por música (rock, pop, jazz). Coletamos os seguintes dados de uma amostra de 100 pessoas:
| Rock | Pop | Jazz | |
|---|---|---|---|
| Ensino Médio | 20 | 30 | 15 |
| Graduação | 10 | 25 | 20 |
| Pós-Graduação | 5 | 10 | 10 |
# Criação da tabela de contingência
ensMed = c(20,30,15)
grad = c( 10,25,20)
posGrad = c(5,10,10)
# juntando os três tipos de ensinos em um tabela por linha
edu.music = rbind(ensMed,grad,posGrad)
# atribuindo os nomes das linhas e das colunas
rownames(edu.music) <- c("Ensino Médio", "Graduação", "Pós-Graduação")
colnames(edu.music) <- c("Rock", "Pop", "Jazz")
print(edu.music) Rock Pop Jazz
Ensino Médio 20 30 15
Graduação 10 25 20
Pós-Graduação 5 10 10A partir da função chisq.test() vamos testar as seguintes hipóteses:
\(H_{0}\): Não há relação entre nível de educação e preferência por música (as variáveis são independentes).
\(H_{1}\): Há relação entre nível de educação e preferência por música (as variáveis não são independentes).
Como valor-p de 0,3077, concluímos com 95% de confiabilidade que não há evidências suficientes contra \(H_{0}\), ou seja, não há relação entre as variáveis nível de educação e preferência por música.
# Efetuando o teste Qui-quadro de independência
t.indep1 <- chisq.test(edu.music)
print(t.indep1)
Pearson's Chi-squared test
data: edu.music
X-squared = 4.8072, df = 4, p-value = 0.3077# Frequências esperadas
print(t.indep1$expected) Rock Pop Jazz
Ensino Médio 15.689655 29.13793 20.172414
Graduação 13.275862 24.65517 17.068966
Pós-Graduação 6.034483 11.20690 7.758621Exemplo (aderência): Um estudo sobre acidentes de trabalho numa indústria revelou que, em 150 acidentes, obtemos a distribuição da tabela a seguir:
| Dia | Seg. | Terça | Quarta | Quinta | Sexta | Total |
|---|---|---|---|---|---|---|
| Oi | 32 | 40 | 20 | 25 | 33 | 150 |
| Ei | 30 | 30 | 30 | 30 | 30 | 150 |
*[@bussab2017]
Vamos verificar se os acidentes ocorrem com a mesma frequência durante a semana. As hipoteses podem ser colocadas da seguinte maneira:
\(H_{0}: p_{1} = p_{2} = p_{3} = p_{4} = p_{5} = 1/5\),
\(H_{1}: p_{j} \neq 1/5,\) para pelo menos um \(j\).
Concluímos que não há evidências suficientes contra \(H_{0}\) (valor-p = 0,09405). Isto é, não há diferenças significativas entre o número de acidentes de trabalho que ocorre por dia da semana.
# dados
acidentes <- c(32,40,20,25,33)
t.adere1 <- chisq.test(acidentes)
print(t.adere1)
Chi-squared test for given probabilities
data: acidentes
X-squared = 7.9333, df = 4, p-value = 0.09405t.adere1$expected[1] 30 30 30 30 30Exemplo (homogeneidade): Derminar se a distribuição de espécies de árvores em três diferentes habitats naturais é a mesma.
| Floresta Tropical | Savana | Floresta Temperada | |
|---|---|---|---|
| Espécie A | 25 | 30 | 15 |
| Espécie B | 20 | 10 | 30 |
| Espécie C | 15 | 25 | 30 |
Hipóteses:
\(H_{0}\): A distribuição das espécies de árvores é a mesma nos três habitats naturais.
\(H_{1}\): A distribuição das espécies de árvores difere em pelo menos um dos três habitats naturais.
Há fortes evidências contra \(H_{0}\) (valor-p \(\approx 0,001\)). Portanto, podemos afirmar com 95% de confiabilidade que as espécies de árvores não se distribuem de maneira homogênea entre os tipos de floresta.
# Criação da tabela de contingência
espA = c(25, 30, 15)
espB = c(20, 10, 30)
espC = c(15, 25, 30)
# juntando os três tipos de espécies em um tabela por linha
arvores = rbind(espA,espB,espC)
# atribuindo os nomes das linhas e das colunas
colnames(arvores) <- c("Floresta Tropical", "Savana", "Floresta Temperada")
rownames(arvores) <- c("Espécie A", "Espécie B", "Espécie C")
arvores Floresta Tropical Savana Floresta Temperada
Espécie A 25 30 15
Espécie B 20 10 30
Espécie C 15 25 30# Efetuando o teste
t.homo1 = chisq.test(arvores)
print(t.homo1)
Pearson's Chi-squared test
data: arvores
X-squared = 17.717, df = 4, p-value = 0.001402t.homo1$expected Floresta Tropical Savana Floresta Temperada
Espécie A 21 22.75 26.25
Espécie B 18 19.50 22.50
Espécie C 21 22.75 26.252.1.4.2 Exercícios
- Verifique ao nível de 5% se há diferença significativa na concentração de poluentes antes e depois da implementação de uma medida de controle ambiental em um determinado local.
Concentração de poluentes antes: 26, 27, 32, 20, 26, 21, 32, 20, 32, 27, 22, 21, 22, 23, 31
Concentração de poluentes depois: 19, 21, 19, 20, 23, 20, 22, 19, 24, 23, 32, 33, 27, 27, 24
antes <- c(26, 27, 32, 20, 26, 21, 32, 20, 32, 27, 22, 21, 22, 23, 31)
depois <- c(19, 21, 19, 20, 23, 20, 22, 19, 24, 23, 32, 33, 27, 27, 24)O banco de dados
CO2dorbaseé um experimento sobre a tolerância ao frio da espécie de gramínea Echinochloa crus-galli. A absorção de CO2 por seis plantas de Quebec e seis plantas do Mississippi foi medida em vários níveis de concentração ambiente de CO2. Metade das plantas de cada tipo foram resfriadas durante a noite antes da realização do experimento. Verifique se há diferenças significativas entre as condições de tratamentos nas taxas de absorção de dióxido de carbono (uptake) ao nível de 5%.Um modelo genético especifica que animais de certa população devam estar classificados em quatro categorias, com probabilidades p1 = 0,656, p2 = 0,093, p3 = 0,093, p4 = 0,158. Dentre 197 animais, obtivemos as seguintes frequências observadas: O1 = 125, O2 = 18, O3 = 20, O4 = 34. Teste se esses dados estão de acordo com o modelo genético postulado (Use \(\alpha = 0,05\)).
*[@bussab2017]
| Categoria 1 | Categoria 2 | Categoria 3 | Categoria 4 | |
|---|---|---|---|---|
| Frequência observada | 125 | 18 | 20 | 34 |
| Probabilidades | 0,656 | 0,093 | 0,093 | 0,158 |
obs = c(125,18,20,34)
probs = c(0.656,0.093,0.093,0.158)
# chisq.test(obs,p=probs)- Uma prova básica de Estatística foi aplicada a 100 alunos de Ciências Humanas e a 100 alunos de Ciências Biológicas. As notas são classificadas segundo os graus A, B, C, D e E (em que D significa que o aluno não recebe créditos e E indica que o aluno foi reprovado).
| A | B | C | D | E | |
|---|---|---|---|---|---|
| C.Humanas | 15 | 20 | 30 | 20 | 15 |
| C.Biológicas | 8 | 23 | 18 | 34 | 17 |
*[@bussab2017]
- Qual o valor da estatística Qui-quadrado (\(\chi^{2}\))?
- Quais os valores esperados para cada casela?
- Quais as hipóteses a serem testadas? iV. Verifique se as distribuições das notas, para as diversas classes, são as mesmas para os dois grupos de alunos. Estabeleça as hipóteses e conclua apresentando o valor-p do teste.
ch = c(15,20,30,20,15)
cb = c(8,23,18,34,17)
prova = rbind(ch,cb)
teste = chisq.test(prova)- Uma companhia de seguros analisou a frequência com que 2.000 segurados (1.000 homens e 1.000 mulheres) usaram hospitais.
| Homens | Mulheres | |
|---|---|---|
| Usaram o hospital | 100 | 150 |
| Não saram o hospital | 900 | 850 |
*[@bussab2017]
- Qual o valor da estatística Qui-quadrado (\(\chi^{2}\))?
- Quais os valores esperados para cada casela?
- Quais as hipóteses a serem testadas?
- Teste a hipótese de que o uso de hospital independe do sexo do segurado.
usaram = c(100,150)
N.usaram = c(900,850)
seguro = rbind(usaram,N.usaram)
teste = chisq.test(seguro)3 Unidade III - Regressão Linear Simples
A análise de regressão linear simples é uma técnica estatística fundamental que nos permite entender e modelar a relação entre duas variáveis, chamadas de variável independente (ou explanatória) e variável dependente (ou resposta).
Exemplo: Um médico deseja verificar se a idade de 20 pacientes do sexo masculino tem alguma influência na pressão arterial diastólica (P.A.D) dos mesmos.
y1 = c(80,70,76,80,80,87,82,80,82,84,75,70,68,75,85,90,75,70,82,85)
x1 = c(35,37,45,53,67,72,40,50,60,65,54,43,39,66,67,72,40,50,60,65)
dados_reg = data.frame(y1,x1)
knitr::kable(dados_reg, col.names = c("Y (P.A.D.)","X (Idade)"), align = c("c","c")) | Y (P.A.D.) | X (Idade) |
|---|---|
| 80 | 35 |
| 70 | 37 |
| 76 | 45 |
| 80 | 53 |
| 80 | 67 |
| 87 | 72 |
| 82 | 40 |
| 80 | 50 |
| 82 | 60 |
| 84 | 65 |
| 75 | 54 |
| 70 | 43 |
| 68 | 39 |
| 75 | 66 |
| 85 | 67 |
| 90 | 72 |
| 75 | 40 |
| 70 | 50 |
| 82 | 60 |
| 85 | 65 |
Existe relação (correlação) entre as variáveis Y e X? Se afirmativo, essa relação é fraca ou forte? Positiva ou negativa? Existindo correlação entre as variáveis, pode-se estabelecer um modelo do tipo:
\[ Y = \beta_{0} + \beta_{1}X + e_{i}, \]
em que, Y é a variável dependente (ou variável resposta); X a variável independente (também chamada de explicativa ou covariável); \(\beta_{0}\) e \(\beta_{1}\) são parâmetros que devem ser estimados; e \(e_{i}\) é o erro aleatório.
A partir de métodos de estimação de mínimos quadrados (MMQ) para os parâmetros \(\beta_{0}\) e \(\beta_{1}\), que consiste em minimizar a quantidade de informação perdida pelo modelo, isto é, a soma dos quadrados dos erros (SQ):
\[SQ(\beta_{0},\beta_{1}) = \displaystyle \sum_{i=1}^{n} e_{i}^{2} = \sum_{i=1}^{n}(y_{i}-(\beta_{0}+\beta_{1}x_{i}))^{2} \]
Derivando (SQ) e igualando a zero chegamos ao seguinte modelo:
\[ \hat{y}_{i} = \hat{\beta}_{0} + \hat{\beta}_{1}x, \]
em que \(\hat{\beta}_{0}\) e \(\hat{\beta}_{1}\) são os estimadores dos parâmetros do modelo e são dados por
\[ \hat{\beta}_{0} = \bar{y} -\hat{\beta}_{1}\bar{x} \] \[ \hat{\beta}_{1} = \dfrac{\displaystyle \sum_{i = 1}^{n}x_{i}y_{i}-n\bar{x}\bar{y}}{\displaystyle \sum_{i = 1}^{n}x_{i}^{2}-n\bar{x}^{2}} \]
Note que \(\hat{\beta}_{1}\) é o coeficiente angular da regressão, ou seja, a mensuração do quanto a reta está inclinada, indicando a variação da variável Y por unidade de variação em X.
Para verificar se existe correlação entre as variáveis Y (P.A.D) e X (idade) do exemplo anterior a partir da análise de regressão linear simples, testamos as seguintes hipóteses:
\(H_{0}:\) a P.A.D dos homens não sofre alterações com a idade.
\(H_{1}:\) a P.A.D dos homens é alterada com a idade.
Ou de forma equivalente:
\(H_{0}: \beta_{1} = 0\) (não existe RLS de Y em X);
\(H_{1}: \beta_{1} \neq 0\) (existe RLS de Y em X).
3.1 Pressupostos
Os principais pressupostos da análise de regressão linear simples são:
Linearidade: A relação entre a variável dependente e as variáveis independentes é linear. Isso significa que o efeito de uma mudança em uma variável independente sobre a variável dependente é constante, independentemente dos valores das outras variáveis independentes.
Independência dos resíduos: Os resíduos (erros de previsão) não estão correlacionados entre si. Isso implica que os erros em uma observação não estão relacionados aos erros em outras observações.
Homocedasticidade: A variabilidade dos erros (ou resíduos) é constante em todas as faixas de valores das variáveis independentes. Em outras palavras, não deve haver nenhum padrão discernível nos resíduos ao longo do tempo ou em relação aos valores previstos.
Normalidade dos resíduos: Os resíduos seguem uma distribuição normal. Isso implica que os erros de previsão são distribuídos simetricamente em torno de zero e não exibem nenhum padrão discernível quando plotados em relação aos valores previstos.
Ao realizar a análise de diagnóstico de resíduos, estamos verificando se esses pressupostos estão sendo atendidos pelo modelo de regressão linear. Se os pressupostos não forem atendidos, pode ser necessário tomar medidas corretivas, como transformar variáveis ou considerar modelos alternativos.
3.2 Coeficiente de Correlação de Pearson - r
Medida estatística que descreve a correlação entre duas variáveis quantitativas. É representada por um valor entre -1 e 1, em que -1 é uma correlação negativa perfeita, 1 é uma correlação positiva perfeita e 0 indica ausência de correlação. O gráfico de dispersão nos dá uma ideia do quanto duas variáveis podem estar correlacionadas, porém, somente o coeficiente de correlação pode verificar está hipótese a partir de um teste. O coeficiente de correlação linear de Pearson é dado por:
\[ r = \dfrac{\displaystyle \sum_{i = 1}^{n}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\displaystyle \sqrt{(\sum_{i = 1}^{n}(x_{i}-\bar{x})^{2})(\sum_{i = 1}(y_{i}-\bar{y})^{2}})} \]
As estimativas dos parâmetros na RLS são dados a partir da ANOVA. Com ela, verificamos se estas estimativas são significativas bem como, a qualidade do ajuste com o coeficiente de determinação que é dados pelo quadrado do coeficiente de correlação, ou seja, \(R^{2}=r^{2}\). No , podemos testar com a função , as seguintes hipóteses:
\(H_{0}: \rho_{x,y} = 0\) (ausência de associação linear);
\(H_{1}: \rho_{x,y} \neq 0\) (presença de associação linear).
3.3 Coeficiente de determinação (\(R^{2}\))
O coeficiente de determinação \(R^{2}\) é uma medida estatística que representa a proporção da variabilidade da variável dependente (Y) que é explicada pelos regressores (variáveis independentes, X) em um modelo de regressão linear. Em outras palavras, \(R^{2}\) aos dados observados. Ele é dado pela razão entre SQE (variação devida aos tratamentos) pela SQT (variação total dos valores observados).
\[R^{2} = \dfrac{SQE}{SQT},\]
em que \(0 \leq R^{2} \leq 1\).
Exemplo: Para verificar se existe correlação entre as variáveis Y (P.A.D) e X (idade) a partir da análise de regressão linear simples, estaríamos testando as seguintes hipóteses:
\(H_{0}\): a P.A.D dos homens não sofre alterações com a idade.
\(H_{1}\): a P.A.D dos homens é alterada com a idade.
Ou de forma equivalente:
\[H_{0}: \beta_{1} = 0;\]
\[H_{1}: \beta_{1} \neq 0. \] Primeiro, vamos analisar o gráfico de dispersão. Podemos notar que, a medida que x1 aumenta, y1 também aumenta. Contudo, essa análise gráfica é subjetiva.
y1 = c(80,70,76,80,80,87,82,80,82,84,75,70,68,75,85,90,75,70,82,85)
x1 = c(35,37,45,53,67,72,40,50,60,65,54,43,39,66,67,72,40,50,60,65)
dados_reg = data.frame(y1,x1)
head(dados_reg) y1 x1
1 80 35
2 70 37
3 76 45
4 80 53
5 80 67
6 87 72plot(x1,y1)
O teste de correlação é paramétrico, logo, precisamos verificar a normalidade dos dados. A partir do teste de Shapiro-Wilk, concluímos que ambas variáveis são normais (valor-p > 0,05). Agora sim, a partir do teste de correlação, podemos afirmar com 95% de confiabilidade que existe uma correlação moderada e positiva (\(r \approx 0,69\); valor-p < 0,05), o que corrobora com a análise gráfica. Isto é, à medida que a idade aumenta, aumenta também a pressão arterial diastólica.
shapiro.test(y1)
Shapiro-Wilk normality test
data: y1
W = 0.95303, p-value = 0.4154shapiro.test(x1)
Shapiro-Wilk normality test
data: x1
W = 0.92002, p-value = 0.09918cor.test(x1,y1)
Pearson's product-moment correlation
data: x1 and y1
t = 4.0286, df = 18, p-value = 0.000788
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3538909 0.8669348
sample estimates:
cor
0.6885776 Ajuste do modelo:
A partir da função lm() obtemos as estimativas dos coeficientes do modelo RLS, em que, \(\beta_{0} = 60,25\) e \(\beta_{1} = 0,34\) são significativas (valor-p < 0,05). O coeficiente de determinação (\(R^{2}=0,4741\)) indica que 47,41% da variabilidade da pressão arterial diastólica é explicada pela idade dos pacientes, o que nos leva a concluir que a qualidade do ajuste é boa. Note que se calcularmos a raiz quadrada do \(R^{2}\), temos o valor do coeficiente de correlação \(r \approx 0,69\).
Modelo ajustado:
\[y (P.A.D) = 60,25 + 0,34\times(Idade)\]
A cada unidade de X, temos um aumento de 0,34 em Y. Ou seja, a cada unidade da idade, temos um aumento de 0,34 na pressão arterial diastólica destes pacientes.
ajuste1 = lm(y1 ~ x1, dados_reg) # ajuste do modelo
summary(ajuste1) # resumo das estimativas do modelo
Call:
lm(formula = y1 ~ x1, data = dados_reg)
Residuals:
Min 1Q Median 3Q Max
-7.922 -3.399 1.139 2.118 8.009
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 60.25010 4.71942 12.766 1.85e-10 ***
x1 0.34352 0.08527 4.029 0.000788 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.628 on 18 degrees of freedom
Multiple R-squared: 0.4741, Adjusted R-squared: 0.4449
F-statistic: 16.23 on 1 and 18 DF, p-value: 0.000788plot(x1,y1,
ylab = "(P.A.D.)",xlab = "Idade")
abline(ajuste1, col = "blue")
3.4 Análise de resíduos
Nesta seção verificamos os pressupostos do modelo RLS.
- Normalidade dos resíduos: como o valor-p do teste de Shapiro-Wilk é maior que 5%, concluímos que os resíduos são normalmente distribuídos.
shapiro.test(ajuste1$residuals)
Shapiro-Wilk normality test
data: ajuste1$residuals
W = 0.94252, p-value = 0.2675- Homocedasticidade: como o valor-p do teste de Breush-Pagan é maior que 5%, concluímos que os resíduos são homocedasticos.
library(lmtest) # carregando o pacote para o teste de Breush-PaganCarregando pacotes exigidos: zoo
Anexando pacote: 'zoo'Os seguintes objetos são mascarados por 'package:base':
as.Date, as.Date.numericbptest(ajuste1)
studentized Breusch-Pagan test
data: ajuste1
BP = 1.8255, df = 1, p-value = 0.1767Resíduos em função dos valores estimados (gráfico 1) : se os resíduos se distribuem de maneira razoavelmente aleatória e com mesma amplitude em torno do zero dizemos que os pressupostos de independência e a homocedasticidade foram aceitos.
Normalidade dos resíduos (gráfico 2): se pontos estão próximos à reta teórica, o pressuposto de normalidade dos resíduos é aceito.
(Gráfico 3) O mesmo que o primeiro, mas, com os resíduos padronizados.
Distâncias de Cook (gráfico 4): medida de influência que pode indicar a presença de outliers quando possui valor maior do que 1
par(mfrow = c(2,2))
plot(ajuste1, which = (1:4), pch = 20)
3.4.1 Exercícios
- trees: Este conjunto de dados fornece medições de diâmetro, altura e volume de madeira em 31 cerejeiras pretas derrubadas. Ajuste um modelo de RLS e responda:
- As variáveis são correlacionadas?
- Escreva o modelo ajustado da altura em funcão do volume.
- Faça a análise dos testes de normalidade e de homocedasticidade dos resíduos.
- Interprete os gráficos do modelo ajustado.
- Os dados abaixo correspondem a altura de uma planta com o tempo de exposição à luz solar.
| Tempo | Altura |
|---|---|
| 1 | 1.439524 |
| 2 | 3.769823 |
| 3 | 7.558708 |
| 4 | 8.070508 |
| 5 | 10.129288 |
| 6 | 13.715065 |
| 7 | 14.460916 |
| 8 | 14.734939 |
| 9 | 17.313147 |
| 10 | 19.554338 |
Faça a análise de cada saída do r para este experimento.
### normalidade dos dados
shapiro.test(Tempo)
Shapiro-Wilk normality test
data: Tempo
W = 0.97016, p-value = 0.8924shapiro.test(Altura)
Shapiro-Wilk normality test
data: Altura
W = 0.96356, p-value = 0.8256### correlação
# gráfico de dispersão
plot(planta$Tempo,planta$Altura,
ylab = "Altura",xlab = "Tempo")
# teste de correlação
cor.test(Tempo,Altura)
Pearson's product-moment correlation
data: Tempo and Altura
t = 17.835, df = 8, p-value = 1e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.9468098 0.9971814
sample estimates:
cor
0.9876575 ### Ajustar o modelo de regressão linear simples
modelo <- lm(Altura ~ Tempo, data = planta)
# Visualizar o resumo do modelo
summary(modelo)
Call:
lm(formula = Altura ~ Tempo, data = planta)
Residuals:
Min 1Q Median 3Q Max
-1.1348 -0.5624 -0.1393 0.3854 1.6814
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5255 0.6673 0.787 0.454
Tempo 1.9180 0.1075 17.835 1e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9768 on 8 degrees of freedom
Multiple R-squared: 0.9755, Adjusted R-squared: 0.9724
F-statistic: 318.1 on 1 and 8 DF, p-value: 1e-07### Plotar o gráfico da regressão
plot(planta$Tempo,planta$Altura,
ylab = "Altura",xlab = "Tempo")
abline(modelo, col = "blue")
### Testes para os pressupostos
# normalidade dos resíduos
shapiro.test(modelo$residuals)
Shapiro-Wilk normality test
data: modelo$residuals
W = 0.92954, p-value = 0.4434### homocedasticidade dos resíduos
library(lmtest) # carregando o pacote para o teste de Breush-Pagan
bptest(Altura ~ Tempo, data = planta)
studentized Breusch-Pagan test
data: Altura ~ Tempo
BP = 0.24532, df = 1, p-value = 0.6204### verificação gráfica dos pressupostos
par(mfrow = c(2,2))
plot(modelo, which = (1:4), pch = 20)
- Seja a reta de regressão \(y = 4 + 1,3\times X\), com r = 0,73 (valor-p<0,01). Sendo a variável resposta o tempo de recuperação de pacientes com lombalgia e a variável independente a idade
- Qual a interpretação para o coeficiente \(\beta_{1}\)?
- Há correlação entre o tempo de recuperação e a idade? Justifique.
- Qual o percentual da variabilidade do tempo de recuperação é explicada pela idade dos pacientes?
- Estamos interessados em estudar a relação entre a dose de um medicamento administrado e o tempo necessário para que o efeito desejado seja alcançado. Abaixo temos o resumo de um modelo de regressão linear simples ajustado aos dados. Com base nestes resultados, responda as seguintes questões:
- Escreva a reta de regressão do modelo ajustado.
- Qual o valor do coeficiente de correlação?
- Interprete o coeficiente de determinação.
- Você aceita ou rejeita a hipótese de regressão?
- Escreva um pequeno texto com suas conclusões sobre este experimento. Lembre (releia) das conclusões feitas durante as aulas e que estão na apostila. Sempre que fizer uma conclusão, apresente no texto, entre parênteses ou não, estatísticas que provem o que você está comentando, ou seja, coloque os coeficientes, estatísticas, valor-p dos testes.
Call:
lm(formula = Tempo ~ Dose, data = dados_farmacia)
Residuals:
Min 1Q Median 3Q Max
-5.6738 -2.8121 -0.6962 1.9268 8.4071
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.6273 3.3364 3.785 0.005352 **
Dose 0.7180 0.1075 6.677 0.000156 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.884 on 8 degrees of freedom
Multiple R-squared: 0.8478, Adjusted R-squared: 0.8288
F-statistic: 44.58 on 1 and 8 DF, p-value: 0.0001564`geom_smooth()` using formula = 'y ~ x'