Teste T Student nas Notas de Ciências Naturais do enem no Colegio Y
Luis Henrique
2024-06-06
Teste estatístico para averiguar diferenças nas médias de grupos por nota de Ciências Naturais
Um professor queria entender em qual categoria de sexos deveria tentar chamar mais a atenção para suas aulas
Então ele decidiu tirar as médias dos grupos por sua matéria e verificou que os homens mantinham mais pontos, em média, por sua amostra, que as mulheres. Mas o quão confiável é essa diferença visto que vieram de uma amostra?
library(dplyr)
library(car)
library(RVAideMemoire)
library(ggplot2)
library(ggfun)
library(ggthemes)
library(ggpubr)
library(webr)
library(gridExtra)
setwd("C:/Users/Luis Henrique/Desktop/Docs/Estudos R")
enem2019_tratado <- read.csv('enem2019_tratado.csv', sep = ",")
# Objetivo:
# Analisar a diferença entre as notas de Homens e mulheres do colegio estudado 35151506 em CN
# Criando o dataframe de interesse
colegioy <- enem2019_tratado %>% filter(CO_ESCOLA == "35151506")
#DIFERENÇA ENCONTRADA DAS AMOSTRAS ENTRE OS SEXOS
group_by(colegioy, TP_SEXO) %>% get_summary_stats(NOTA_CN)## # A tibble: 2 × 14
## TP_SEXO variable n min max median q1 q3 iqr mad mean sd
## <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 F NOTA_CN 15 480. 614. 566. 557. 574. 17.4 12.8 564. 32.2
## 2 M NOTA_CN 15 552. 671. 601 578. 628 50.3 42.0 604. 37.1
## # ℹ 2 more variables: se <dbl>, ci <dbl># HOMOGENEIDADE DAS VARIÂNCIAS (HOMOCEDASTICIDADE)
# Variabilidade dos erros constante.
# Ho = variâncias homogêneas : p > 0.05
# Ha = variâncias não homogêneas : p <= 0.05
leveneTest(NOTA_CN ~ TP_SEXO, colegioy, center=mean) #PRECISA COLOCAR CENTER - MEAN PQ O DEFAULT É MEDIAN## Levene's Test for Homogeneity of Variance (center = mean)
## Df F value Pr(>F)
## group 1 1.3411 0.2566
## 28#NORMALIDADE
#Nível de significância (alfa) : 0,05
#Ho = distribuição normal : p > 0.05
#Ha = distribuição != normal : p <= 0.05
byf.shapiro(NOTA_CN ~ TP_SEXO, colegioy) #teste da nota_cn quebrado por tp_sexo##
## Shapiro-Wilk normality tests
##
## data: NOTA_CN by TP_SEXO
##
## W p-value
## F 0.9032 0.1065
## M 0.9603 0.6974teste_shap <- shapiro.test(colegioy$NOTA_CN)
#GRÁFICO DA NORMALIDADE DOS DADOS
ggplot(colegioy, aes(sample=NOTA_CN))+
stat_qq(shape=21, size=2.5,col='darkblue',fill='cornflowerblue')+
stat_qq_line(lwd=2.0, col='coral1')+
labs(title = 'Distribuição dos valores - Ciências Naturais',
x='Quantis Teóricos (Distribuição Normal',
y='Quantis da Amostra')+
theme_gdocs()+
annotate('text', x= 1.5, y= 620, size=3 ,label=paste("Shapiro Wilk: ", round(teste_shap$statistic,2)
, "P_Value: "))
##PODEMOS CONCLUIR QUE OS DADOS POSSUEM DISTRIBUIÇÕES NORMAIS E SÃO HOMOGÊNEAS ##########
##PORTANTO ESTÁ VÁLIDA A UTILIZAÇÃO DO TESTE T DE DUAS AMOSTRAS SEM A CORREÇÃO DE WELCH #Os dados são normais e com distribuições normais.
Levando isso em conta, como desejamos verificar a significância em diferenças das médias, vamos abordar um teste de média de grupos independentes paramétrico o famoso Teste T de Student
############# APLICANDO O TESTE T
# Grupo Experimental (Nível 1) : Homens
# Grupo de Controle (Nível 2) : Mulheres
# TESTE t PARA AMOSTRAS INDEPENDENTES
# Ho = A DIFERENÇA DE MÉDIA É IGUAL A 0 OU MENOR QUE 0 : p > 0.05
# Ha = HÁ DIFERENÇAA É MAIOR QUE 0 : p <= 0.05
#PRIMEIRO VAMOS COLOCAR UM FACTOR PARA PODERMOS DEFINIR OS NÍVEIS A SEREM AVALIADOS NO TESTE
colegioy$TP_SEXO <- as.factor(colegioy$TP_SEXO)
levels(colegioy$TP_SEXO)## [1] "F" "M"colegioy$TP_SEXO <- factor(colegioy$TP_SEXO, levels = c("M", "F")) #definindo níveis corretos
#teste com valor p menor que 0.05 e ponto dentro da região critica (vermelha) :
# significa que podemos rejeitar a hipótese nula nesse caso e aceitar Ha
teste1 <- t.test(NOTA_CN ~ TP_SEXO, colegioy, var.equal=TRUE,alternative="greater") #variâncias homogêneas (var.equal=TRUE).
teste1##
## Two Sample t-test
##
## data: NOTA_CN by TP_SEXO
## t = 3.1809, df = 28, p-value = 0.001787
## alternative hypothesis: true difference in means between group M and group F is greater than 0
## 95 percent confidence interval:
## 18.78201 Inf
## sample estimates:
## mean in group M mean in group F
## 604.3067 563.9333plot(teste1) + theme_gdocs() + annotate('text',
x=teste1$statistic,
y=0.02, size = 4,
label=paste("p_value: ",
round(teste1$p.value,2)))
Após aplicado o Teste T foi compreendido que o t_statistic não foi encontrado ao acaso, ou em outras palavaras, a diferença entre as médias não foi encontrada ao acaso, portanto podemos rejeitar a hipótese nula
Porém, devemos tomar cuidado com o Erro Tipo II, um erro mais subjetivo. Para isso vamos montar as curvas de Poder e OC
Erro Tipo II é entendido como a probabilidade da diferença esperada menos a diferença das amostras dividido pelo desvio das amostras (uma das formas de ser calculado) Ele indica a probabilidade de se aceitar a hipótese nula quando a hipótese alternativa está correta
###### AGORA PRECISAMOS CALCULAR O PODER DO TESTE T
#Para isso precisamos calcular a probabilida de B e poder 1 - B
group1 <- colegioy %>% filter(TP_SEXO == "M")
sigma1 <- sd(group1$NOTA_CN)
mi1 <- mean(group1$NOTA_CN)
group2 <- colegioy %>% filter(TP_SEXO == "F")
sigma2 <- sd(group2$NOTA_CN)
mi2 <- mean(group2$NOTA_CN)
delta <- mi1 - mi2
alpha <- 0.05
n1 <- nrow(group1)
n2 <- nrow(group2)
df <- n1 + n2 -2
t_alpha <- qt(1 - alpha, df)
SE <- sqrt((sigma1^2 / n1) + (sigma2^2 / n2))
t_beta <- delta / SE
##Probabilidade de erro tipo II
beta <- pt(t_alpha - t_beta, df)
beta## [1] 0.07504665#Poder estatistico
power <- 1 - beta
power## [1] 0.9249533#Quantidade de valores de diferenças - 0 até 70, de 0.5 a 0.5
delta_range <- seq(0, 70, by = 0.5)
# Inicializar vetores para armazenar os valores de poder e OC
power_values <- numeric(length(delta_range))
OC_values <- numeric(length(delta_range))
# Calcular o poder para cada valor de delta
for (i in seq_along(delta_range)) {
delta <- delta_range[i]
SE <- sqrt((sigma1^2 / n1) + (sigma2^2 / n2))
t_beta <- delta / SE
beta <- pt(t_alpha - t_beta, df)
power <- 1 - beta
power_values[i] <- power
OC_values[i] <- beta
}
# Criar um dataframe com os resultados
df_power <- data.frame(delta = delta_range, power = power_values, OC = OC_values)
# Plotar as curvas de poder e OC
power <- ggplot(df_power, aes(x = delta)) +
geom_line(aes(y = power), lwd=1.5, color='tomato') +
labs(title = "Curva de Poder",
subtitle = "Poder do teste dado as diferenças dos grupos",
x = "Delta (Diferença de Médias)",
y = "Probabilidade")
CO <- ggplot(df_power, aes(x = delta)) +
geom_line(aes(y = OC), col='coral1',lwd=1.5) +
labs(title = "Curva de OC",
subtitle = 'Curva de Probabilidade de erro tipo II dado as diferenças dos grupos ',
x = "Delta (Diferença de Médias)",
y = "Probabilidade") +
annotate('text', x = 10, y = 0.25, size = 4,label= paste("Prob.β:", round(beta,2)))
ggarrange(power + theme_gdocs(), CO + theme_gdocs(), ncol = 1)
Agora que definimos nosso poder para as diversas diferenças, vamos ao passo final e plotar dois boxplot para verificar isso mais de perto
############# BOXPLOT PARA ANALISAR OS GRUPOS
delta_medias <- mean(group1$NOTA_CN) - mean(group2$NOTA_CN)
delta_medias## [1] 40.37333p <- ggplot(data = colegioy, aes(x=TP_SEXO, y=NOTA_CN))+
geom_boxplot(aes(fill=TP_SEXO))+
scale_fill_manual(values = c("tomato", "lightblue"), name = 'SEXO')+
stat_summary(fun = mean, geom = "point", shape=21, size=3, fill="white")+
labs(title = "Notas de Ciências Naturais entre M e F",
subtitle="Após conferência e análise do teste t, vamos inferir sobre a diferença média dos grupos\nA diferença média é em média mais de 40 pontos a mais para o grupo Homens",
x="Sexo",
y="Nota de Ciências Naturais")+
theme_gdocs() +
theme(
legend.position = c(0.95, 1),
legend.justification = c("right", "top"),
legend.background = element_blank()
)
posicao_media_x <- mean(as.numeric(colegioy$TP_SEXO))
p <- p + stat_compare_means(method = "t.test", comparisons = list(c("M", "F")))
p <- p + annotate('text', x = posicao_media_x, y= 660, label = paste("Diferença média: ", round(delta_medias,2)))
p