Primeiro Exemplo

Sign Test for Two-sample Paired Data

O Teste de Sinal de duas amostras avalia o número de observações em um grupo que são maiores que as observações pareadas no outro grupo, sem levar em conta a magnitude da diferença. O teste é similar em propósito ao teste de Wilcoxon de duas amostras, mas examina especificamente o valor mediano das diferenças (se os valores forem numéricos) e não é afetado pela distribuição dos dados.

A função SIGN.test no pacote BSDA requer que os dados sejam separados em duas variáveis, cada uma ordenada para que a primeira observação de cada uma esteja emparelhada, e assim por diante. Informações sobre opções para a função podem ser exibidas com? SIGN.test. A função SignTest no pacote DescTools é semelhante.

Hipóteses:No caso da hipótese nula (para dados numéricos), a mediana das diferenças pareadas na população da qual a amostra foi retirada é igual a zero. Quanto a hipótese alternativa (bilateral): (para dados numéricos), a mediana das diferenças pareadas na população da qual a amostra foi sorteada não é igual a zero.Quanto a interpretação dos resultados eles podem ser significativos e podem ser relatados do tipo: “Houve uma diferença significativa nos valores entre o grupo A e o grupo B.”

Neste exemplo vamos utilizar os dados abaixo.

if(!require(psych)){install.packages("psych")}
## Loading required package: psych
if(!require(BSDA)){install.packages("BSDA")}
## Loading required package: BSDA
## Loading required package: lattice
## 
## Attaching package: 'BSDA'
## The following object is masked from 'package:datasets':
## 
##     Orange
if(!require(DescTools)){install.packages("DescTools")}
## Loading required package: DescTools
## 
## Attaching package: 'DescTools'
## The following objects are masked from 'package:psych':
## 
##     AUC, ICC, SD
Input =("
 Speaker  Time  Student  Likert
 Pooh      1     a        1       
 Pooh      1     b        4       
 Pooh      1     c        3        
 Pooh      1     d        3       
 Pooh      1     e        3       
 Pooh      1     f        3       
 Pooh      1     g        4       
 Pooh      1     h        3        
 Pooh      1     i        3       
 Pooh      1     j        3       
 Pooh      2     a        4
 Pooh      2     b        5
 Pooh      2     c        4 
 Pooh      2     d        5
 Pooh      2     e        4
 Pooh      2     f        5
 Pooh      2     g        3
 Pooh      2     h        4 
 Pooh      2     i        3
 Pooh      2     j        4      
")

Data = read.table(textConnection(Input),header=TRUE)

library(psych)
library(BSDA)
library(DescTools)

headTail(Data)
##     Speaker Time Student Likert
## 1      Pooh    1       a      1
## 2      Pooh    1       b      4
## 3      Pooh    1       c      3
## 4      Pooh    1       d      3
## ...    <NA>  ...    <NA>    ...
## 17     Pooh    2       g      3
## 18     Pooh    2       h      4
## 19     Pooh    2       i      3
## 20     Pooh    2       j      4
str(Data)
## 'data.frame':    20 obs. of  4 variables:
##  $ Speaker: Factor w/ 1 level "Pooh": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Time   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Student: Factor w/ 10 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Likert : int  1 4 3 3 3 3 4 3 3 3 ...
summary(Data)
##  Speaker        Time        Student      Likert    
##  Pooh:20   Min.   :1.0   a      :2   Min.   :1.00  
##            1st Qu.:1.0   b      :2   1st Qu.:3.00  
##            Median :1.5   c      :2   Median :3.50  
##            Mean   :1.5   d      :2   Mean   :3.55  
##            3rd Qu.:2.0   e      :2   3rd Qu.:4.00  
##            Max.   :2.0   f      :2   Max.   :5.00  
##                          (Other):8

Two-sample sign test with BSDA package 

Time.1 = Data$Likert [Data$Time == 1]
Time.2 = Data$Likert [Data$Time == 2]

library(BSDA)

SIGN.test(x = Time.1,
          y = Time.2,
          alternative = "two.sided",
          conf.level = 0.95)
## 
##  Dependent-samples Sign-Test
## 
## data:  Time.1 and Time.2
## S = 1, p-value = 0.03906
## alternative hypothesis: true median difference is not equal to 0
## 95 percent confidence interval:
##  -2.0000000 -0.3244444
## sample estimates:
## median of x-y 
##            -1 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt  U.E.pt
## Lower Achieved CI     0.8906     -2 -1.0000
## Interpolated CI       0.9500     -2 -0.3244
## Upper Achieved CI     0.9785     -2  0.0000

Two-sample sign test with DescTools package 

Time.1 = Data$Likert [Data$Time == 1]
Time.2 = Data$Likert [Data$Time == 2]

library(DescTools)

SignTest(x = Time.1,
         y = Time.2)
## 
##  Dependent-samples Sign-Test
## 
## data:  Time.1 and Time.2
## S = 1, number of differences = 9, p-value = 0.03906
## alternative hypothesis: true median difference is not equal to 0
## 97.9 percent confidence interval:
##  -2  0
## sample estimates:
## median of the differences 
##                        -1

 

Segundo Exemplo

Sign tests for paired data in R

O teste de sinal é um teste não paramétrico que pode ser usado para testar se há mais valores negativos ou positivos em uma amostra de dados. A hipótese nula é que existe um número igual de valores negativos e positivos.

Por exemplo, se você tiver dados pareados (por exemplo, ‘antes’ e ‘depois’) para os indivíduos da amostra, você pode usar um teste de sinais para testar se as diferenças entre os valores ‘antes’ e ‘depois’ tendem a ser positivas ou negativo (a hipótese nula é que existe um número igual de diferenças negativas e positivas). 

x1 <- c(488, 478, 480, 426, 440, 410, 458, 460)
x2 <- c(484, 478, 492, 444, 436, 398, 464, 476)
difference <- x1 - x2

library("BSDA")
SIGN.test(difference, md=0, alternative = "less")
## 
##  One-sample Sign-Test
## 
## data:  difference
## s = 3, p-value = 0.5
## alternative hypothesis: true median is less than 0
## 95 percent confidence interval:
##  -Inf    4
## sample estimates:
## median of x 
##          -3 
## 
## Achieved and Interpolated Confidence Intervals: 
## 
##                   Conf.Level L.E.pt U.E.pt
## Lower Achieved CI     0.8555   -Inf      4
## Interpolated CI       0.9500   -Inf      4
## Upper Achieved CI     0.9648   -Inf      4
s = 3

Aqui três das diferenças são positivas (4, 4, 12), então a estatística de teste é 3. Existem 7 diferenças diferentes de zero. Assumindo que os sinais + e - tenham a mesma probabilidade, a distribuição do número de sinais + de 7 é B (7, 0,5). Podemos calcular o valor P para um teste unilateral como a probabilidade de observar 3 ou menos sinais positivos em 7: 

pbinom(3, size=7, p=0.5)
## [1] 0.5

Isso concorda com o resultado de \(SIGN.test\) ().Para este teste unilateral, a hipótese nula é que o número de sinais + e - é igual, e a hipótese alternativa é que há menos + sinais que - sinais.

 

Terceiro Exemplo

Paired Samples Wilcoxon Test in R

O teste Wilcoxon de amostras pareadas (também conhecido como teste de postos sinalizados de Wilcoxon) é uma alternativa não paramétrica ao teste t pareado usado para comparar dados pareados. É usado quando seus dados não são normalmente distribuídos.

Aqui, usaremos um conjunto de dados de exemplo, que contém o peso de 10 camundongos antes e depois do tratamento. Queremos saber se existe alguma diferença significativa nos pesos medianos antes e depois do tratamento? 

# Data in two numeric vectors

# Weight of the mice before treatment
before <-c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)

# Weight of the mice after treatment
after <-c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)

# Create a data frame
my_data <- data.frame( 
                group = rep(c("before", "after"), each = 10),
                weight = c(before,  after)
                )
print(my_data)
##     group weight
## 1  before  200.1
## 2  before  190.9
## 3  before  192.7
## 4  before  213.0
## 5  before  241.4
## 6  before  196.9
## 7  before  172.2
## 8  before  185.5
## 9  before  205.2
## 10 before  193.7
## 11  after  392.9
## 12  after  393.2
## 13  after  345.1
## 14  after  393.0
## 15  after  434.0
## 16  after  427.9
## 17  after  422.0
## 18  after  383.9
## 19  after  392.3
## 20  after  352.2

As estatísticas de resumo de cálculo (mediana e intervalo interquartil (IQR)) por grupos usando o pacote dplyr podem ser usadas.

Visualize seus dados usando gráficos do tipo box plots. Aqui, usaremos o pacote R ggpubr para facilitar a visualização de dados baseada em ggplot2. Instale a versão mais recente do ggpubr do GitHub da seguinte forma (recomendado): 

library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
group_by(my_data, group) %>%
  summarise(
    count = n(),
    median = median(weight, na.rm = TRUE),
    IQR = IQR(weight, na.rm = TRUE)
  )
## # A tibble: 2 x 4
##   group  count median   IQR
##   <fct>  <int>  <dbl> <dbl>
## 1 after     10   393.  28.8
## 2 before    10   195.  12.6

Ou então, instale o CRAN: 

# Plot weight by group and color by group
library("ggpubr")
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## Loading required package: magrittr
#library("ggplot2")
library("magrittr")
ggboxplot(my_data, x = "group", y = "weight", 
          color = "group", palette = c("#00AFBB", "#E7B800"),
          order = c("before", "after"),
          ylab = "Weight", xlab = "Groups")

Box plots mostra o aumento, mas perde a informação emparelhada. Você pode usar a função plot.paired () [no pacote pairedData] para traçar dados emparelhados (gráfico “antes - depois”). 

library("MASS")
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# Subset weight data before treatment
before <- subset(my_data,  group == "before", weight,
                 drop = TRUE)
# subset weight data after treatment
after <- subset(my_data,  group == "after", weight,
                 drop = TRUE)
# Plot paired data
library(PairedData)
## Loading required package: gld
## Loading required package: mvtnorm
## 
## Attaching package: 'PairedData'
## The following object is masked from 'package:base':
## 
##     summary
pd <- paired(before, after)
plot(pd, type = "profile") + theme_bw()

Pergunta: Existe alguma alteração significativa nos pesos dos ratos antes do tratamento?

1) Compute paired Wilcoxon test - Method 1:Os dados são salvos em dois vetores numéricos diferentes. 

res <- wilcox.test(before, after, paired = TRUE)
res
## 
##  Wilcoxon signed rank test
## 
## data:  before and after
## V = 0, p-value = 0.001953
## alternative hypothesis: true location shift is not equal to 0

2) Compute paired Wilcoxon-test - Method 2: Os dados são salvos em um data frame. 

# Compute t-test
res <- wilcox.test(weight ~ group, data = my_data, paired = TRUE)
res
## 
##  Wilcoxon signed rank test
## 
## data:  weight by group
## V = 55, p-value = 0.001953
## alternative hypothesis: true location shift is not equal to 0
# print only the p-value
res$p.value
## [1] 0.001953125

Como se pode ver, os dois métodos fornecem os mesmos resultados. O p-valor do teste é 0,001953, que é menor que o nível de significância alfa = 0,05. Podemos concluir que o peso mediano dos camundongos antes do tratamento é significativamente diferente do peso mediano após o tratamento com um valor de p = 0,001953.

Note que: se você quiser testar se a mediana do peso antes do tratamento é menor que a mediana do peso após o tratamento, digite o seguinte: 

wilcox.test(weight ~ group, data = my_data, paired = TRUE, alternative = "less")
## 
##  Wilcoxon signed rank test
## 
## data:  weight by group
## V = 55, p-value = 1
## alternative hypothesis: true location shift is less than 0

Ou então, se você quiser testar se a mediana do peso antes do tratamento é maior que a mediana do peso após o tratamento, digite o seguinte: 

wilcox.test(weight ~ group, data = my_data, paired = TRUE, alternative = "greater")
## 
##  Wilcoxon signed rank test
## 
## data:  weight by group
## V = 55, p-value = 0.0009766
## alternative hypothesis: true location shift is greater than 0

 

Quarto Exemplo

Wilcoxon signed rank test in R

Teste de hipótese estatística não paramétrica, para comparação das médias entre duas amostras pareadas. O prefeito de uma cidade quer ver se os níveis de poluição são reduzidos, fechando as ruas para o tráfego de carros. Isso é medido pela taxa de poluição a cada 60 minutos (8h 22h: total de 15 medições) em um dia quando o tráfego é aberto, e em um dia de fechamento do tráfego. Aqui os valores da poluição do ar: 

# com tráfico: 214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234
# sem tráfico: 159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112

É claro que os dois grupos estão emparelhados, porque existe um vínculo entre as leituras, consistindo no fato de estarmos considerando a mesma cidade (com suas peculiaridades climáticas, ventilação, etc.) embora em dois dias diferentes. Não sendo capaz de assumir uma distribuição gaussiana para os valores registrados, devemos prosseguir com um teste não-paramétrico, o teste de classificação sinalizada de Wilcoxon. 

a <- c(214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234)
b <- c(159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112)

wilcox.test(a,b, paired=TRUE)
## 
##  Wilcoxon signed rank test
## 
## data:  a and b
## V = 80, p-value = 0.2769
## alternative hypothesis: true location shift is not equal to 0

Como o valor-p é maior que 0,05, concluímos que as médias permaneceram essencialmente inalteradas (aceitamos a hipótese nula H0), então bloquear o tráfego para um único dia não levou a nenhuma melhoria em termos de poluição da cidade. O valor V = 80 corresponde à soma das classificações atribuídas às diferenças com sinal positivo.

Podemos calcular manualmente a soma das classificações atribuídas às diferenças com sinal positivo e a soma das classificações atribuídas às diferenças com sinal negativo, para comparar esse intervalo com o intervalo tabulado nas tabelas de Wilcoxon para amostras pareadas e confirmar nossa estatística decisão. Veja como calcular as duas somas. 

a <- c(214, 159, 169, 202, 103, 119, 200, 109, 132, 142, 194, 104, 219, 119, 234)
b <- c(159, 135, 141, 101, 102, 168, 62, 167, 174, 159, 66, 118, 181, 171, 112)

 diff <- c(a - b) #calculating the vector containing the differences
 
 diff <- diff[ diff!=0 ] #delete all differences equal to zero
 diff.rank <- rank(abs(diff)) #check the ranks of the differences, taken in absolute
 diff.rank.sign <- diff.rank * sign(diff) #check the sign to the ranks, recalling the signs of the values of the differences
 ranks.pos <- sum(diff.rank.sign[diff.rank.sign > 0]) 
 #calculating the sum of ranks assigned to the differences as a positive, ie greater than zero
 ranks.neg <- -sum(diff.rank.sign[diff.rank.sign < 0]) #calculating the sum of ranks assigned to the differences as a negative, ie less than zero
 
 ranks.pos #it is the value V of the wilcoxon signed rank test
## [1] 80
 ranks.neg
## [1] 40

O intervalo computado é (40, 80). O intervalo tabelado em tabelas de amostras pareadas de Wilcoxon, com 15 diferenças, é (25, 95). Como o intervalo calculado está contido na tabela, aceitamos a hipótese nula H0 de igualdade das médias. Como previsto pelo valor p, o fechamento de estradas para o tráfego não trouxe qualquer melhoria em termos de taxa de poluição.

 

Quinto Exemplo

Non Parametric tests

Um método estatístico é chamado de não-paramétrico se não fizer nenhuma suposição sobre a distribuição da população ou o tamanho da amostra. Isto está em contraste com a maioria dos métodos paramétricos em estatísticas elementares que assumem que os dados são quantitativos, a população tem uma distribuição normal e o tamanho da amostra é suficientemente grande.Em geral, as conclusões tiradas de métodos não paramétricos não são tão poderosas quanto as paramétricas. No entanto, como os métodos não paramétricos fazem menos suposições, eles são mais flexíveis, mais robustos e aplicáveis a dados não quantitativos.

 

Exemplo com o Sign Test

Uma empresa de refrigerantes inventou uma nova bebida e gostaria de descobrir se será tão popular quanto a bebida favorita. Para este propósito, seu departamento de pesquisa organiza 18 participantes para testes de sabor. Cada participante experimenta as duas bebidas em ordem aleatória antes de dar sua opinião.Acontece que 5 dos participantes gostam mais da nova bebida e o resto prefere a antiga. No nível de significância de 0,05, podemos rejeitar a noção de que as duas bebidas são igualmente populares?

Obs.:A hipótese nula é que as bebidas são igualmente populares. Aqui nós aplicamos a função binom.test. Como o valor p é 0,096525 e é maior que o nível de significância de 0,05, não rejeitamos a hipótese nula 

binom.test(5, 18) 
## 
##  Exact binomial test
## 
## data:  5 and 18
## number of successes = 5, number of trials = 18, p-value = 0.09625
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.09694921 0.53480197
## sample estimates:
## probability of success 
##              0.2777778

No nível de significância de 0,05, não rejeitamos a noção de que as duas bebidas são igualmente populares.

 

Exemplo Mann-Whitney-Wilcoxon Test

Aqui, usaremos um conjunto de dados de exemplo, que contém as vendas antes e depois do tratamento (desconto).Queremos saber se existe alguma diferença significativa nas medianas de vendas antes e depois do tratamento? 

# Weight of the mice before treatment
before <-c(200.1, 190.9, 192.7, 213, 241.4, 196.9, 172.2, 185.5, 205.2, 193.7)
# Weight of the mice after treatment
after <-c(392.9, 393.2, 345.1, 393, 434, 427.9, 422, 383.9, 392.3, 352.2)
# Create a data frame
my_data <- data.frame( 
                group = rep(c("before", "after"), each = 10),
                sales = c(before,  after)
                )

group_by(my_data, group) %>%
  summarise(
    count = n(),
    median = median(sales, na.rm = TRUE),
    IQR = IQR(sales, na.rm = TRUE)
  )
## # A tibble: 2 x 4
##   group  count median   IQR
##   <fct>  <int>  <dbl> <dbl>
## 1 after     10   393.  28.8
## 2 before    10   195.  12.6

Pergunta: Há alguma mudança significativa nas vendas antes do tratamento? 

res <- wilcox.test(sales ~ group, data = my_data, paired = TRUE)
res
## 
##  Wilcoxon signed rank test
## 
## data:  sales by group
## V = 55, p-value = 0.001953
## alternative hypothesis: true location shift is not equal to 0

O p-valor do teste é 0,001953, que é menor que o nível de significância alfa = 0,05. Podemos concluir que a mediana das vendas antes do tratamento é significativamente diferente da mediana de vendas após o tratamento com um valor de p = 0,001953.

 

Kruskal-Wallis Test

O teste de Kruskal-Wallis por classificação é uma alternativa não paramétrica ao teste one-way ANOVA, que amplia o teste de duas amostras de Wilcoxon na situação em que há mais de dois grupos. É recomendado quando as suposições do teste ANOVA unidirecional não são atendidas.

Novo Exemplo

Aqui, usaremos o conjunto de dados R incorporado denominado PlantGrowth. Contém o peso das plantas obtidas sob um controle e duas condições de tratamento diferentes. 

my_data <- PlantGrowth
head(my_data)
##   weight group
## 1   4.17  ctrl
## 2   5.58  ctrl
## 3   5.18  ctrl
## 4   6.11  ctrl
## 5   4.50  ctrl
## 6   4.61  ctrl

Resumindo as informações: 

group_by(my_data, group) %>%
  summarise(
    count = n(),
    mean = mean(weight, na.rm = TRUE),
    sd = sd(weight, na.rm = TRUE),
    median = median(weight, na.rm = TRUE),
    IQR = IQR(weight, na.rm = TRUE)
  )
## # A tibble: 3 x 6
##   group count  mean    sd median   IQR
##   <fct> <int> <dbl> <dbl>  <dbl> <dbl>
## 1 ctrl     10  5.03 0.583   5.15 0.743
## 2 trt1     10  4.66 0.794   4.55 0.662
## 3 trt2     10  5.53 0.443   5.44 0.467

Queremos saber se existe alguma diferença significativa entre os pesos médios das plantas nas três condições experimentais. 

kruskal.test(weight ~ group, data = my_data)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  weight by group
## Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842

Como o valor de p é menor que o nível de significância de 0,05, podemos concluir que existem diferenças significativas entre os grupos de tratamento. A partir do resultado do teste de Kruskal-Wallis, sabemos que há uma diferença significativa entre os grupos, mas não sabemos quais pares de grupos são diferentes.É possível usar a função \(pairwise.wilcox.test ()\) para calcular comparações entre pares entre níveis de grupo com correções para vários testes. 

# pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group, p.adjust.method = "BH")

#ctrl  trt1 
#trt1 0.199 -    
#trt2 0.095 0.027

O R avisa que não é possível computar o valor de p exato com o de desempate.A comparação pareada mostra que apenas trt1 e trt2 são significativamente diferentes (p <0,05).

 

Sexto Exemplo

Wilcoxon Signed-Rank Test

Duas amostras de dados são combinadas se vierem de observações repetidas do mesmo assunto. Usando o Teste de Classificação Assinada de Wilcoxon, podemos decidir se as distribuições de população de dados correspondentes são idênticas, sem pressupor que elas sigam a distribuição normal.

Exemplo

No conjunto de dados interno denominado immer, o rendimento de cevada nos anos 1931 e 1932 do mesmo campo é registrado. Os dados de rendimento são apresentados nas colunas do quadro de dados Y1 e Y2. 

head(immer) 
##   Loc Var    Y1    Y2
## 1  UF   M  81.0  80.7
## 2  UF   S 105.4  82.3
## 3  UF   V 119.7  80.4
## 4  UF   T 109.7  87.2
## 5  UF   P  98.3  84.2
## 6   W   M 146.6 100.4

A hipótese nula é que os rendimentos de cevada dos dois anos de amostragem são populações idênticas. Para testar a hipótese, aplicamos a função wilcox.test para comparar as amostras correspondentes. Para o teste emparelhado, definimos o argumento “pareado” como TRUE. Como o valor p é 0,005318 e é menor que o nível de significância de 0,05, rejeitamos a hipótese nula. 

#wilcox.test(immer$Y1, immer$Y2, paired=TRUE) 

#Wilcoxon signed rank test with continuity correction

#data:  immer$Y1 and immer$Y2
#V = 368.5, p-value = 0.005318
#alternative hypothesis: true location shift is not equal to 0

#Warning message: 
#In wilcox.test.default(immer$Y1, immer$Y2, paired = TRUE) : 
  #cannot compute exact p-value with ties

O R avisa que não é possível computar o valor de p exato com o de desempate.No nível de significância de 0,05, concluímos que os rendimentos de cevada de 1931 e 1932 do conjunto de dados são populações não-idênticas.

 

Sétimo Exemplo

Kruskal-Wallis Test

Uma coleção de amostras de dados é independente se elas vêm de populações não relacionadas e as amostras não afetam umas às outras. Usando o Teste de Kruskal-Wallis, podemos decidir se as distribuições populacionais são idênticas, sem presumir que elas sigam a distribuição normal.

Exemplo

No conjunto de dados interno denominado airquality, as medições diárias da qualidade do ar em Nova York, de maio a setembro de 1973, são registradas. A densidade do ozono é apresentada na coluna do quadro de dados Ozone. 

 head(airquality) 
##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6

Problema:Sem assumir que os dados tenham distribuição normal, teste em um nível de significância de 0,05 se a densidade mensal de ozônio em Nova York tiver distribuições de dados idênticas de maio a setembro de 1973.

Solução:A hipótese nula é que a densidade mensal de ozônio é de populações idênticas. Para testar a hipótese, aplicamos a função kruskal.test para comparar os dados mensais independentes. O valor p é quase zero (6.901e-06). Por isso rejeitamos a hipótese nula. 

kruskal.test(Ozone ~ Month, data = airquality) 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ozone by Month
## Kruskal-Wallis chi-squared = 29.267, df = 4, p-value = 6.901e-06

Resposta:No nível de significância de 0,05, concluímos que a densidade de ozônio mensal em Nova York de maio a setembro de 1973 é de populações não-identitárias.

 

Oitavo Exemplo

Exemplo utilizando testes Paramétricos

Testes estatísticos paramétricos estão entre os mais comuns que você encontrará. Eles incluem teste t, análise de variância e regressão linear. Eles são usados quando a variável dependente é uma variável de dados de intervalo / razão. Isso pode incluir variáveis medidas na ciência, como comprimento do peixe, altura da criança, peso de produção ou concentração de poluentes na água.Uma vantagem de usar testes estatísticos paramétricos é que seu público provavelmente estará familiarizado com as técnicas e interpretação dos resultados. Esses testes são também mais flexíveis e mais poderosos que seus análogos não paramétricos. Sua principal desvantagem é que todos os testes paramétricos pressupõem algo sobre a distribuição dos dados subjacentes. Se essas suposições forem violadas, as estatísticas de teste resultantes não serão válidas e os testes não serão tão poderosos quanto para os casos em que as suposições forem cumpridas.

Dados de contagem podem não ser apropriados para testes paramétricos comuns. Um erro freqüente é usar modelos paramétricos comuns e testes com dados de contagem para a variável dependente. Em vez disso, os dados de contagem podem ser analisados usando testes de dados nominais ou usando métodos de regressão apropriados para dados de contagem. Estes incluem regressão de Poisson, regressão binomial negativa e regressão de Poisson inflada a zero. Veja o capítulo Regressão para dados de contagem.

Quando usar testes paramétricos para dados de contagem? Às vezes, é permitido usar testes paramétricos comuns para dados de contagem ou outros dados discretos. Elas podem ser usadas nos casos em que as contagens são usadas como um tipo de medição de alguma propriedade dos indivíduos, desde que 1) a distribuição de dados ou resíduos da análise satisfaçam aproximadamente as premissas do teste; e 2) há poucas ou nenhumas contagens próximas ou próximas de zero, ou perto de um máximo, se existir uma. Exemplos permitidos podem incluir pontuações de testes, idade ou número de etapas tomadas durante o dia. Tecnicamente, cada uma dessas medidas é limitada por zero e são medidas discretas em vez de contínuas. No entanto, se outras condições forem atendidas, é razoável tratá-las como se fossem variáveis de medição contínuas.Esse tipo de dados de contagem às vezes precisará ser transformado para atender às premissas da análise paramétrica. Raiz quadrada e transformações logarítmicas são comuns. No entanto, se houver muitas contagens iguais ou próximas a zero, é improvável que a transformação ajude. Geralmente, não vale a pena tentar forçar a contagem de dados para atender aos pressupostos da análise paramétrica com transformações, uma vez que existem métodos mais apropriados disponíveis.

Suposições em estatísticas paramétricas: Todas as análises paramétricas têm hipóteses sobre os dados subjacentes, e essas suposições devem ser confirmadas ou assumidas com um bom motivo ao usar esses testes. Se essas suposições forem violadas, as estatísticas e conclusões resultantes não serão válidas e os testes poderão não ter energia em relação aos testes alternativos. 

Input = ("
Student  Sex     Teacher  Steps  Rating
a        female  Catbus    8000   7
b        female  Catbus    9000  10
c        female  Catbus   10000   9
d        female  Catbus    7000   5
e        female  Catbus    6000   4
f        female  Catbus    8000   8
g        male    Catbus    7000   6
h        male    Catbus    5000   5
i        male    Catbus    9000  10
j        male    Catbus    7000   8
k        female  Satsuki   8000   7
l        female  Satsuki   9000   8
m        female  Satsuki   9000   8
n        female  Satsuki   8000   9
o        male    Satsuki   6000   5
p        male    Satsuki   8000   9
q        male    Satsuki   7000   6
r        female  Totoro   10000  10
s        female  Totoro    9000  10
t        female  Totoro    8000   8
u        female  Totoro    8000   7
v        female  Totoro    6000   7
w        male    Totoro    6000   8
x        male    Totoro    8000  10
y        male    Totoro    7000   7
z        male    Totoro    7000   7
")

Data = read.table(textConnection(Input),header=TRUE)

Data
##    Student    Sex Teacher Steps Rating
## 1        a female  Catbus  8000      7
## 2        b female  Catbus  9000     10
## 3        c female  Catbus 10000      9
## 4        d female  Catbus  7000      5
## 5        e female  Catbus  6000      4
## 6        f female  Catbus  8000      8
## 7        g   male  Catbus  7000      6
## 8        h   male  Catbus  5000      5
## 9        i   male  Catbus  9000     10
## 10       j   male  Catbus  7000      8
## 11       k female Satsuki  8000      7
## 12       l female Satsuki  9000      8
## 13       m female Satsuki  9000      8
## 14       n female Satsuki  8000      9
## 15       o   male Satsuki  6000      5
## 16       p   male Satsuki  8000      9
## 17       q   male Satsuki  7000      6
## 18       r female  Totoro 10000     10
## 19       s female  Totoro  9000     10
## 20       t female  Totoro  8000      8
## 21       u female  Totoro  8000      7
## 22       v female  Totoro  6000      7
## 23       w   male  Totoro  6000      8
## 24       x   male  Totoro  8000     10
## 25       y   male  Totoro  7000      7
## 26       z   male  Totoro  7000      7
library(psych)
library(DescTools)
headTail(Data)
##     Student    Sex Teacher Steps Rating
## 1         a female  Catbus  8000      7
## 2         b female  Catbus  9000     10
## 3         c female  Catbus 10000      9
## 4         d female  Catbus  7000      5
## ...    <NA>   <NA>    <NA>   ...    ...
## 23        w   male  Totoro  6000      8
## 24        x   male  Totoro  8000     10
## 25        y   male  Totoro  7000      7
## 26        z   male  Totoro  7000      7
str(Data)
## 'data.frame':    26 obs. of  5 variables:
##  $ Student: Factor w/ 26 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Sex    : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Teacher: Factor w/ 3 levels "Catbus","Satsuki",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Steps  : int  8000 9000 10000 7000 6000 8000 7000 5000 9000 7000 ...
##  $ Rating : int  7 10 9 5 4 8 6 5 10 8 ...
summary(Data)
##     Student       Sex        Teacher       Steps           Rating      
##  a      : 1   female:15   Catbus :10   Min.   : 5000   Min.   : 4.000  
##  b      : 1   male  :11   Satsuki: 7   1st Qu.: 7000   1st Qu.: 7.000  
##  c      : 1               Totoro : 9   Median : 8000   Median : 8.000  
##  d      : 1                            Mean   : 7692   Mean   : 7.615  
##  e      : 1                            3rd Qu.: 8750   3rd Qu.: 9.000  
##  f      : 1                            Max.   :10000   Max.   :10.000  
##  (Other):20
rm(Input)

Amostragem aleatória

Todos os testes estatísticos assumem que os dados capturados na amostra são escolhidos aleatoriamente da população como um todo. O viés de seleção obviamente afetará a validade do resultado da análise.

Os testes também pressupõem que as observações são independentes umas das outras, exceto quando a análise leva em conta a não-independência. Um caso comum de observações não independentes é em experimentos de medidas repetidas, em que o mesmo sujeito é observado ao longo do tempo. Se você estivesse medindo as pontuações dos testes dos alunos ao longo do tempo, você pode esperar que os alunos com uma pontuação alta no teste em uma data tenham uma pontuação alta nas próximas datas. Neste caso, a observação em uma data não seria independente das observações em outras datas.A independência da observação é frequentemente assumida a partir de um bom desenho experimental. Além disso, dados ou resíduos podem ser plotados, por exemplo, para ver se as observações de uma data estão correlacionadas com as de outra data.

Distribuição normal de dados ou resíduos

Testes paramétricos assumem que os dados provêm de uma população de distribuição conhecida. Em particular, os testes discutidos nesta seção assumem que a distribuição dos dados é condicionalmente normal na distribuição. Ou seja, os dados são normalmente distribuídos uma vez que os efeitos das variáveis no modelo são levados em consideração. Em termos práticos, isso significa que os resíduos da análise devem ser distribuídos normalmente. Isso geralmente será avaliado com um histograma de resíduos, um gráfico de densidade, como mostrado abaixo, ou com um quantil-quantil plot.Um número selecionado de testes exigirá que os dados sejam distribuídos normalmente. Isso será limitado ao teste t de uma amostra, teste t de duas amostras e teste t pareado. Para outros testes, a distribuição dos resíduos será investigada.

Residuais de uma análise também são comumente chamados de erros. Eles são a diferença entre as observações e o valor previsto pelo modelo. Por exemplo, se a média calculada de uma amostra for 10 e uma observação for 12, o residual para essa observação será 2. Se outra observação for 7, o residual para essa observação será –3.Tenha cuidado para não ficar confuso sobre essa suposição. Você pode ver a discussão sobre como os dados devem ser normalmente distribuídos para testes paramétricos. Isso geralmente é errado. O teste t assume que as observações para cada grupo são normalmente distribuídas, mas se houver uma diferença nos grupos, podemos esperar uma distribuição bimodal, não uma distribuição normal simples, para os dados combinados. É por isso que, na maioria dos casos, analisamos a distribuição dos resíduos, não os dados brutos.

Opcional: Considerando as distribuições de dados para grupos e seus resíduos: O código a seguir serve apenas para ilustrar esse princípio e não é normalmente usado na análise. Primeiro, as distribuições de Passos por Sexo no conjunto de dados Catbus são examinadas, e então a distribuição dos resíduos - médias dos grupos subtraídos de cada observação - é examinada. 

ggplot(Data,
       aes(Steps, fill = Sex)) +
       geom_density(position="dodge",
                    alpha = 0.6)
## Warning: Width not defined. Set with `position_dodge(width = ?)`

O código a seguir cria uma variável Mean for Steps para cada Sex no conjunto de dados Catbus e subtrai Steps de cada média, chamando o resultado Residual. Isto é apenas para fins ilustrativos; Você normalmente não precisará manipular os dados dessa maneira. 

M1 = mean(Data$Steps[Data$Sex=="female"])
M2 = mean(Data$Steps[Data$Sex=="male"])

Data$Mean[Data$Sex=="female"] = M1
Data$Mean[Data$Sex=="male"]   = M2

Data$Residual = Data$Steps - Data$Mean

### Density plot of residuals for mean for each sex for the Catbus data set

ggplot(Data, 
       aes(Residual)) +
         geom_density(fill = "gray")

Homogeneidade de variância: Análises paramétricas também assumirão uma homogeneidade de variância entre os grupos. Ou seja, para o teste t de Student que compara dois grupos, cada grupo deve ter a mesma variação.Uma boa abordagem para avaliar essa suposição é plotar os valores residuais em relação aos valores previstos. Os resíduos devem ter aproximadamente o mesmo spread em todos os valores previstos.Homogeneidade de variância também é chamada de homocedasticidade. O oposto é heterocedasticidade. 

### Residuals from mean for each sex vs. mean for Catbus data set 

plot(jitter(Residual) ~ Mean, 
     data = Data)

Como nota técnica, por padrão R conduz uma variante do teste t chamado teste t de Welch. Este teste não assume homogeneidade de variância e, portanto, pode ser usado para comparar dois grupos com variâncias desiguais.

Aditividade dos efeitos do tratamento

Modelos para análise de variância de duas vias e análises similares são construídos como modelos lineares nos quais a variável dependente é predita como uma combinação linear das variáveis independentes.Uma violação dessa suposição é às vezes indicada quando um gráfico de valores residuais versus valores preditos exibe um padrão curvo.

Outliers

Outliers são observações cujo valor está muito além do esperado. Eles podem atrapalhar as análises paramétricas, pois afetam a distribuição dos dados e influenciam fortemente a média. Há uma variedade de testes formais para detectar outliers, mas eles não serão discutidos aqui. A melhor abordagem é aquela que considera os resíduos após uma análise. Boas ferramentas são o gráfico “Residuais vs. alavancagem” e outros gráficos na seção “Outros diagramas de diagnóstico” abaixo. É minha opinião que outliers não devem ser removidos dos dados a menos que haja um bom motivo, geralmente quando um valor é impossível ou um erro de medição é suspeito.

Testes paramétricos são de certa forma robustos

Alguns testes paramétricos são de certa forma robustos a violações de certas suposições. Por exemplo, o teste t é razoavelmente robusto para violações da normalidade para distribuições simétricas, mas não para amostras com variâncias desiguais (a menos que o teste t de Welch seja usado). Uma análise unidirecional da variância é também razoavelmente robusta a violações na normalidade.

O resultado é que as premissas do modelo devem sempre ser verificadas, mas você pode tolerar pequenas violações na distribuição de resíduos ou homocedasticidade. Grandes violações invalidarão o teste. É importante ser honesto com suas avaliações ao verificar as premissas do modelo. É melhor transformar dados, alterar seu modelo, usar um método robusto ou usar um teste não-paramétrico do que não confiar em sua análise.

Avaliando os pressupostos do modelo

Para este exemplo, revisitaremos os dados do Catbus. Em seguida, definiremos um modelo linear em que Steps é a variável dependente e Sex and Teacher são as variáveis independentes. 

Input = ("
Student  Sex     Teacher  Steps  Rating
a        female  Catbus    8000   7
b        female  Catbus    9000  10
c        female  Catbus   10000   9
d        female  Catbus    7000   5
e        female  Catbus    6000   4
f        female  Catbus    8000   8
g        male    Catbus    7000   6
h        male    Catbus    5000   5
i        male    Catbus    9000  10
j        male    Catbus    7000   8
k        female  Satsuki   8000   7
l        female  Satsuki   9000   8
m        female  Satsuki   9000   8
n        female  Satsuki   8000   9
o        male    Satsuki   6000   5
p        male    Satsuki   8000   9
q        male    Satsuki   7000   6
r        female  Totoro   10000  10
s        female  Totoro    9000  10
t        female  Totoro    8000   8
u        female  Totoro    8000   7
v        female  Totoro    6000   7
w        male    Totoro    6000   8
x        male    Totoro    8000  10
y        male    Totoro    7000   7
z        male    Totoro    7000   7
")

Data = read.table(textConnection(Input),header=TRUE)


###  Check the data frame

library(psych)

headTail(Data)
##     Student    Sex Teacher Steps Rating
## 1         a female  Catbus  8000      7
## 2         b female  Catbus  9000     10
## 3         c female  Catbus 10000      9
## 4         d female  Catbus  7000      5
## ...    <NA>   <NA>    <NA>   ...    ...
## 23        w   male  Totoro  6000      8
## 24        x   male  Totoro  8000     10
## 25        y   male  Totoro  7000      7
## 26        z   male  Totoro  7000      7
str(Data)
## 'data.frame':    26 obs. of  5 variables:
##  $ Student: Factor w/ 26 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Sex    : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 2 2 2 2 ...
##  $ Teacher: Factor w/ 3 levels "Catbus","Satsuki",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Steps  : int  8000 9000 10000 7000 6000 8000 7000 5000 9000 7000 ...
##  $ Rating : int  7 10 9 5 4 8 6 5 10 8 ...
summary(Data)
##     Student       Sex        Teacher       Steps           Rating      
##  a      : 1   female:15   Catbus :10   Min.   : 5000   Min.   : 4.000  
##  b      : 1   male  :11   Satsuki: 7   1st Qu.: 7000   1st Qu.: 7.000  
##  c      : 1               Totoro : 9   Median : 8000   Median : 8.000  
##  d      : 1                            Mean   : 7692   Mean   : 7.615  
##  e      : 1                            3rd Qu.: 8750   3rd Qu.: 9.000  
##  f      : 1                            Max.   :10000   Max.   :10.000  
##  (Other):20
### Remove unnecessary objects

rm(Input)


#Define a linear model

model = lm(Steps ~ Sex + Teacher,
           data = Data)

Usando testes formais para avaliar a normalidade dos resíduos

Existem testes formais para avaliar a normalidade dos resíduos. Testes comuns incluem Shapiro-Wilk, Anderson-Darling, Kolmogorov-Smirnov e D’Agostino-Pearson. Estes são apresentados na seção “Análises opcionais: testes formais para a normalidade”.Em geral, não recomendo usar esses testes porque os resultados deles dependem do tamanho da amostra. Quando o tamanho da amostra é grande, os testes podem indicar uma partida estatisticamente significativa da normalidade, mesmo que essa partida seja pequena. E quando os tamanhos das amostras são pequenos, eles não detectam desvios da normalidade.O artigo do blog Fells Stats, na seção “Referências”, contém exemplos bastante convincentes desses problemas.

Inclinação e curtose

Não há diretrizes definitivas sobre qual faixa de inclinação ou curtose são aceitáveis para considerar os resíduos normalmente distribuídos.Em geral, eu não recomendaria depender de cálculos de inclinação e curtose, mas sim usar histogramas e outros gráficos.Se eu fosse forçado a dar conselhos para cálculos de assimetria, eu diria, seja cauteloso se o valor absoluto for> 0,5, e considere que não é normalmente distribuído se o valor absoluto for> 1.0. Alguns autores usam o 2.0 como um ponto de corte para a normalidade, e outros usam um limite maior para a curtose. 

library(psych)

x = residuals(model)

describe(x,
         type=2)
##    vars  n mean      sd median trimmed     mad     min     max   range
## X1    1 26    0 1132.26 -39.58   11.07 1114.18 -2202.4 2123.25 4325.65
##     skew kurtosis     se
## X1 -0.13    -0.11 222.05

Usando inspeção visual para avaliar a normalidade dos resíduos*: Normalmente, o melhor método para ver se os resíduos do modelo atendem às premissas de distribuição normal e homocedasticidade é plotá-los e inspecionar visualmente os gráficos.

Histograma com curva normal*

Um histograma dos resíduos deve ser aproximadamente normal, sem inclinação excessiva ou curtose. Adicionando uma curva normal com a mesma média e desvio padrão que os dados ajudam a avaliar o histograma. 

x = residuals(model)


library(rcompanion)

#normal.histogram(x)

plotNormalHistogram(x, prob = FALSE, col = "gray", main = "",
  linecol = "blue", lwd = 2, length = 1000)

Gráfico de densidade do kernel com curva normal

Um gráfico de densidade do kernel é semelhante a um histograma, mas é suavizado em uma curva. Às vezes, um gráfico de densidade fornece uma representação melhor da distribuição de dados, porque a aparência do histograma depende de quantas caixas são usadas. A função \(plotNormalDensity\) irá produzir este gráfico. As opções incluem aquelas para a função de plotagem, bem como ajustar, bw e kernel que são passados para a função de densidade. col1, col2 e col3 alteram as cores do gráfico e lwd altera a espessura da linha. 

x = residuals(model)

plotNormalDensity(x,
                  adjust = 1)  ### Decrease this number 

                                 ###  to make line less smooth

Plot of residuals vs. fitted values

Padrões no gráfico de resíduos versus valores ajustados podem indicar uma falta de homocedasticidade ou que os erros não são independentes dos valores ajustados.

plot(fitted(model),
     residuals(model))

Outros diagramas de diagnóstico

A utilização da função de plotagem para os tipos mais comuns de modelos mostrará plotagens de diagnóstico ou outros gráficos úteis para o modelo. Para modelos lineares, a função plot produz gráficos de 1) valores residuais vs. valores ajustados; 2) quantil normal - quantil plot; 3) raiz quadrada de resíduos padronizados versus valores ajustados; e 4) resíduos padronizados versus alavancagem.

 

Embora esses sejam diagramas de diagnóstico úteis, sua interpretação não será discutida aqui.

plot(model)

 

Análises opcionais: Testes formais para normalidade

Código para a realização de testes formais de normalidade estão incluídos aqui. Eu não recomendo usá-los, mas eles podem ser instrutivos em treinar o olho para interpretar histogramas e gráficos Q-Q, ou para responder a um revisor menos experiente.Em cada caso, a hipótese nula é que a distribuição de dados não é diferente da normal. Ou seja, um valor de p significativo (p <0,05) sugere que os dados não são normalmente distribuídos. 

#Define a linear model

model = lm(Steps ~ Sex + Teacher,
           data = Data)


#Shapiro–Wilk normality test
x = residuals(model)
shapiro.test(x)
## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.96256, p-value = 0.4444
#Anderson-Darling normality test
 if(!require(nortest)){install.packages("nortest")}
## Loading required package: nortest
library(nortest)
x = residuals(model)
ad.test(x)
## 
##  Anderson-Darling normality test
## 
## data:  x
## A = 0.39351, p-value = 0.3506
#One-sample Kolmogorov-Smirnov test
#x = residuals(model)
#ks.test(x,"pnorm", mean = mean(x), sd   = sd(x))

# O R emite a seguinte mensagem:"ties should not be present for the Kolmogorov-Smirnov test"


#D'Agostino Normality Test
if(!require(fBasics)){install.packages("fBasics")}
## Loading required package: fBasics
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:psych':
## 
##     outlier
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:psych':
## 
##     tr
library(fBasics)
x = residuals(model)
dagoTest(x)
## 
## Title:
##  D'Agostino Normality Test
## 
## Test Results:
##   STATISTIC:
##     Chi2 | Omnibus: 0.1071
##     Z3  | Skewness: -0.3132
##     Z4  | Kurtosis: 0.0948
##   P VALUE:
##     Omnibus  Test: 0.9479 
##     Skewness Test: 0.7541 
##     Kurtosis Test: 0.9245 
## 
## Description:
##  Tue Jun 25 18:26:45 2019 by user: 1545-6_IRON

Como mencionado anteriormente, uma limitação ao uso desses testes é que, conforme o número de observações é aumentado, a capacidade do teste de retornar um valor de p significativo aumenta, mesmo para pequenos desvios de uma distribuição normal.

 

Testes de homogeneidade de Variância

A homogeneidade da variância é uma suposição subjacente tanto aos testes t quanto aos testes F (ANOVAs) nos quais as variâncias populacionais (ou seja, a distribuição ou “spread” das pontuações em torno da média) de duas ou mais amostras são consideradas iguais . Em correlações e regressões, o termo “homogeneidade de variância em matrizes”, também chamado de “homocedasticidade”, refere-se à suposição de que, dentro da população, a variância de Y para cada valor de X é constante. Esta entrada concentra-se na homogeneidade da variância no que se refere aos testes t e ANOVAs.

 

Brown–Forsythe or robust Levene’s test

A modificação de Brown-Forsythe do teste de Levene torna-o mais robusto para partidas na normalidade dos dados. 

if(!require(car)){install.packages("car")} 
## Loading required package: car
## Loading required package: carData
## 
## Attaching package: 'carData'
## The following objects are masked from 'package:BSDA':
## 
##     Vocab, Wool
## 
## Attaching package: 'car'
## The following object is masked from 'package:fBasics':
## 
##     densityPlot
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:DescTools':
## 
##     Recode
## The following object is masked from 'package:psych':
## 
##     logit
library(car)
x = residuals(model)
leveneTest(x ~ Sex * Teacher, data=Data)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  0.4015  0.842
##       20
if(!require(lawstat)){install.packages("lawstat")} 
## Loading required package: lawstat
## 
## Attaching package: 'lawstat'
## The following object is masked from 'package:car':
## 
##     levene.test
library(lawstat)
x = residuals(model)
levene.test(x, interaction(Data$Sex, Data$Teacher))
## 
##  Modified robust Brown-Forsythe Levene-type test based on the
##  absolute deviations from the median
## 
## data:  x
## Test Statistic = 0.4015, p-value = 0.842

 

Bartlett’s test for homogeneity of variance

O teste de Bartlett é conhecido por ser sensível à não normalidade nas amostras. Ou seja, amostras não normais podem resultar em um teste significativo devido à não normalidade. 

x = residuals(model)

bartlett.test(x ~ interaction(Sex, Teacher), 
              data=Data)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  x by interaction(Sex, Teacher)
## Bartlett's K-squared = 3.7499, df = 5, p-value = 0.586

 

Fligner-Killeen test

O teste de Fligner-Killeen é outro teste de homogeneidade de variâncias que é robusto para desvios na normalidade dos dados. 

x = residuals(model)

fligner.test(x ~ interaction(Sex, Teacher), data=Data)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  x by interaction(Sex, Teacher)
## Fligner-Killeen:med chi-squared = 2.685, df = 5, p-value = 0.7484