Na aula do dia 25/08/2023, vimos como carregar um pacote, ler um banco de dados e excluir conlunas.
OBS: Um asterisco antes e depois de um trecho faz com que ele fique em itálico. Já dois asteriscos, antes e depois de um trecho, faz com que ele fique em negrito.
# Lendo o pacote:
library(titanic)Warning: package 'titanic' was built under R version 4.2.3# atribuindo titanic_train para o objeto dados
dados <- titanic_train
# excluindo colunas
dados$PassengerId <- NULL
dados$Ticket <- NULL
dados$Cabin <- NULL
dados$Name <- NULLNa aula do dia 01/09/2023, aprendemos a ler um banco de dados externo.
library(readxl)Warning: package 'readxl' was built under R version 4.2.3dados_titanic <- read_excel("C:/Users/18304197731/Downloads/dados_titanic.xlsx")
dados_titanic$Nome<-NULL# Vendo a estrutura dos dados
str(dados_titanic)tibble [891 × 8] (S3: tbl_df/tbl/data.frame)
$ Sobreviveu : num [1:891] 0 1 1 1 0 0 0 0 1 1 ...
$ Classe : num [1:891] 3 1 3 1 3 3 1 3 3 2 ...
$ Sexo : chr [1:891] "male" "female" "female" "female" ...
$ Idade : num [1:891] 22 38 26 35 35 NA 54 2 27 14 ...
$ N_irmaos_conjuges: num [1:891] 1 1 0 1 0 0 0 3 0 1 ...
$ N_pais_filhos : num [1:891] 0 0 0 0 0 0 0 1 2 0 ...
$ Tarifa : num [1:891] 7.25 71.28 7.92 53.1 8.05 ...
$ Porto : chr [1:891] "S" "C" "S" "S" ...# Transformando para fator
dados_titanic$Sobreviveu<- as.factor(dados_titanic$Sobreviveu)
dados_titanic$Classe<-as.factor(dados_titanic$Classe)
dados_titanic$Porto<-as.factor(dados_titanic$Porto)
dados_titanic$Sexo<-as.factor(dados_titanic$Sexo)levels(dados_titanic$Sobreviveu)[1] "0" "1"levels(dados_titanic$Sobreviveu) <- c("Não", "Sim")
levels(dados_titanic$Classe)[1] "1" "2" "3"levels(dados_titanic$Classe)<- c("Primeira", "Segunda", "Terceira")
levels(dados_titanic$Sexo)[1] "female" "male" levels(dados_titanic$Sexo)<- c("Feminino", "Masculino")colnames(dados_titanic)[1] "Sobreviveu" "Classe" "Sexo"
[4] "Idade" "N_irmaos_conjuges" "N_pais_filhos"
[7] "Tarifa" "Porto" colnames(dados_titanic)[8] <- "Porto_de_Embarque"# substiruir {r} por {r, output=F} faz com que a saída do código não apareça no relatório, pois é muito grande.
dados_titanic$Faixa_Etaria<-cut(dados_titanic$Idade,c(0,18,65,200))
dados_titanic$Faixa_Etaria
dados_titanic$Faixa_Etaria
levels(dados_titanic$Faixa_Etaria)<-c("Até 18 anos", "Maior que 18 anos e até 65 anos", "Maior que 65 anos")Vamos trabalhar com a variável classe econômica. Vamos construir uma tabela de distribuição de frequências. Conclui-se que a maioria dos passageiros pertencia à terceira classe (55,1%) e a minoria à segunda classe (20,7%). Ver figura 1.
library(summarytools) #chamando o pacoteWarning: package 'summarytools' was built under R version 4.2.3freq(dados_titanic$Classe) #como queremos fazer uma tabela de distribuição de frequência, usamos a função "freq"Frequencies
dados_titanic$Classe
Type: Factor
Freq % Valid % Valid Cum. % Total % Total Cum.
-------------- ------ --------- -------------- --------- --------------
Primeira 216 24.24 24.24 24.24 24.24
Segunda 184 20.65 44.89 20.65 44.89
Terceira 491 55.11 100.00 55.11 100.00
<NA> 0 0.00 100.00
Total 891 100.00 100.00 100.00 100.00# % Valid= frequência relativa (porcentagem)
# % Valis cum.= frequencia relativa acumulada (soma das porcenagens)
# diferença entre % Valid e % Total: a % valid considera apenas as pessoas que tem informação; a % total considera os passageiros que não possuem informação (os N/A). Ou seja, só haverá diferença quando houver presença de N/A.# importante deixar o cursor dentro desse chuck e clicar em addins, clicar em ggplot2 builder, se pedir pra instalar algo instale
library(ggplot2)Warning: package 'ggplot2' was built under R version 4.2.3ggplot(dados_titanic) +
aes(x = Classe) +
geom_bar(fill = "#ED73E3") +
labs(x = "Classe Econômica",
y = "Frequência", title = "Figura 1", caption = "Fonte: Autora") +
theme_minimal() +
theme(plot.title = element_text(face = "bold.italic", hjust = 0.5))
# CTRL + SHIFT + C= Torna toda a área selecionada como comentário
#Para mudar a ordem das categorias de classe econômica:
# dados_titanic$Classe<-factor(dados_titanic$Classe, levels=c("Terceira", "Segunda", "Primeira"))
# factor(dados_titanic)
#eu não preciso fazer todo o caminho através do Addins para mudar algo do grafico. eu posso alterar apenas ali dentro do codigo que gerou o grafico.
# para mudar a cor, vou no console e escrevo colors(), aparecerá uma lista de cores.
# para saber como é cada cor (pois n da pra saber pelo nome), copie todas as opçoes de cores que apareceu e cole em um arquivo novo, irá aparecer todas os nomes das cores coloridos A idade média dos passageiros foi de 29,70 anos (dp=14,53 anos), com mediana igual a 28 anos. A variável pode ser classificada como assimetrica à direita (positiva) e heterogênea (CV>0,30). Por ser uma variável univariada quantitativa, podemos plotar em histograma, em gráfico de densidade e em boxplot (o gráfico em linha não é recomendado pois a idade nao e algo que varia com o tempo nesse caso).
library(summarytools)#ferramenta pra acessar a "biblioteca" de dados
descr(dados_titanic$Idade)#chamando o pacote $ chamando a variávelDescriptive Statistics
dados_titanic$Idade
N: 891
Idade
----------------- --------
Mean 29.70
Std.Dev 14.53
Min 0.42
Q1 20.00
Median 28.00
Q3 38.00
Max 80.00
MAD 13.34
IQR 17.88
CV 0.49
Skewness 0.39
SE.Skewness 0.09
Kurtosis 0.16
N.Valid 714.00
Pct.Valid 80.13library(dplyr)Warning: package 'dplyr' was built under R version 4.2.3
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(ggplot2)
dados_titanic %>%
filter(!is.na(Idade)) %>%
ggplot() +
aes(x = Idade) +
geom_histogram(bins = 30L, fill = "#ED73E3") +
labs(x = "Idade (em anos)",
y = "Frequência Absoluta", title = "Distribuição da Variável Idade") +
theme_gray()
library(dplyr)
library(ggplot2)
dados_titanic %>%
filter(!is.na(Idade)) %>%
ggplot() +
aes(x = Idade) +
geom_density(fill = "#ED73E3") +
labs(x = "Idade (em anos)",
y = "Densidade", title = "Figura 3") +
theme_gray()
library(dplyr)
library(ggplot2)
dados_titanic %>%
filter(!is.na(Idade)) %>%
ggplot() +
aes(y = Idade) +
geom_boxplot(fill = "#ED73E3", width = 0.35) +
labs(x = "",
y = "Idade (em anos)", title = "Figura 4") +
theme_gray() + xlim(c(-1,1))
Será que existe relação entre a classe econômica e a sobrevivência ao desastre do Titanic? - Variável explicativa: classe econômica - Variável resposta: sobrevivência A tabela abaixo mostra que, dos passageiros da primeira classe, 37% deles não sobreviveram, sendo este percentual maior na segunda (52,7%) e terceira (75,8%) classes.
library(summarytools)
ctable(dados_titanic$Classe, dados_titanic$Sobreviveu)Cross-Tabulation, Row Proportions
Classe * Sobreviveu
Data Frame: dados_titanic
---------- ------------ ------------- ------------- --------------
Sobreviveu Não Sim Total
Classe
Primeira 80 (37.0%) 136 (63.0%) 216 (100.0%)
Segunda 97 (52.7%) 87 (47.3%) 184 (100.0%)
Terceira 372 (75.8%) 119 (24.2%) 491 (100.0%)
Total 549 (61.6%) 342 (38.4%) 891 (100.0%)
---------- ------------ ------------- ------------- --------------library(ggplot2)
ggplot(dados_titanic) +
aes(x = Classe, fill = Sobreviveu) +
geom_bar(position = "dodge") +
scale_fill_brewer(palette = "Set1",
direction = 1) +
labs(x = "Classe econômica", y = "Frequência absoluta", title = "Figura 5: Relação entre Classe econômica e desfecho do passageiro") +
theme_minimal()
library(ggplot2)
ggplot(dados_titanic) +
aes(x = Classe, fill = Sobreviveu) +
geom_bar(position = "fill") +
scale_fill_brewer(palette = "Set1",
direction = 1) +
labs(x = "Classe econômica", y = "Frequência relativa", title = "Figura 6: Relação entre classe econômica e desfecho do passageiro") +
theme_minimal()
Existe correlação entre o número de filhos? Pais a bordo e a idade do passageiro? Ao analisar o coeficiente de correlação de Spearman, observamos uma correlção fraca e negativa (rho = -0.25).
library(ggplot2)
ggplot(dados_titanic) +
aes(x = Idade, y = N_pais_filhos) +
geom_point(shape = "circle", size = 1.5,
colour = "#0C4C8A") +
geom_smooth(span = 1L) +
labs(x = "Idade (anos)", y = "Número de Pais - Filhos",
title = "Figura 7: Correlação entre o número de Pais - Filhos e idade dos passageiros") +
theme_gray()`geom_smooth()` using method = 'loess' and formula = 'y ~ x'Warning: Removed 177 rows containing non-finite values (`stat_smooth()`).Warning: Removed 177 rows containing missing values (`geom_point()`).
cor(dados_titanic$Idade , dados_titanic$N_pais_filhos , use = "complete.obs" , method = "spearman")[1] -0.2542121Existe relação entre a idade do passageiro e sua sobrevivência? A partir da estatística, podemos observar que ambos os grupos apresentaram medianas iguais a 28 anos (IQR = 18 anos para grupo “Não” e 17 anos para grupo “Sim”). Quanto à assimetria, ambas as distribuiçãos foram assimétricas positivas. A idade possui distribuição heterogênea em ambos os defechos. Cv > 0.30
library(summarytools)
with(dados_titanic, stby(Idade, Sobreviveu, descr))Descriptive Statistics
Idade by Sobreviveu
Data Frame: dados_titanic
N: 549
Não Sim
----------------- -------- --------
Mean 30.63 28.34
Std.Dev 14.17 14.95
Min 1.00 0.42
Q1 21.00 19.00
Median 28.00 28.00
Q3 39.00 36.00
Max 74.00 80.00
MAD 11.86 13.34
IQR 18.00 17.00
CV 0.46 0.53
Skewness 0.58 0.18
SE.Skewness 0.12 0.14
Kurtosis 0.25 -0.10
N.Valid 424.00 290.00
Pct.Valid 77.23 84.80library(dplyr)
library(ggplot2)
dados_titanic %>%
filter(!is.na(Idade)) %>%
ggplot() +
aes(x = "", y = Idade, fill = Sobreviveu) +
geom_boxplot() +
scale_fill_brewer(palette = "Set1",
direction = 1) +
labs(x = " ", y = "Idade (anos)", title = "Figura 8: Distribuição da idade segundo o desfecho dos passageiros ") +
theme_gray()
# Aula 4
# Tabela Descritiva com dados filtrados
A tabela 1 abaixo apresenta as características dos passageiros a bordo do Titanic. Podemos observar que 62% dos passageiros não sobreviveram ao acidente. Quanto a classe ecônomica, a maioria pertencia a terceira classe, representando 55% dos passageiros. 65% eram do sexo masculino e a idade mediana doi de 28 anos (IQR = 20- 38 anos). Quanto ao número de pais ou filhos a bordo, 76% dos passageiros não tinham esse tipo de acompanhante.
#Criando a tabela com os dados filtrados
# Filtrando as variáveis que usamos nas nossas perguntas-foco
library(gtsummary)Warning: package 'gtsummary' was built under R version 4.2.3dados2<- dados_titanic[, c("Sobreviveu",
"Classe",
"Sexo",
"Idade",
"N_pais_filhos", "Tarifa")]
tbl_summary(dados2, label = list (Sobreviveu ~ "sobreviventes",
Classe ~ "Classe econômica",
N_pais_filhos ~ "Número de pais/filhos"),
missing_text = "Sem informação") %>%
modify_header(label ~ "**variável**") %>%
modify_caption("Tabela 1. Características dos passageiros do Titanic.") %>% bold_labels() %>% italicize_levels()| variável | N = 8911 |
|---|---|
| sobreviventes | |
| Não | 549 (62%) |
| Sim | 342 (38%) |
| Classe econômica | |
| Primeira | 216 (24%) |
| Segunda | 184 (21%) |
| Terceira | 491 (55%) |
| Sexo | |
| Feminino | 314 (35%) |
| Masculino | 577 (65%) |
| Idade | 28 (20, 38) |
| Sem informação | 177 |
| Número de pais/filhos | |
| 0 | 678 (76%) |
| 1 | 118 (13%) |
| 2 | 80 (9.0%) |
| 3 | 5 (0.6%) |
| 4 | 4 (0.4%) |
| 5 | 5 (0.6%) |
| 6 | 1 (0.1%) |
| Tarifa | 14 (8, 31) |
| 1 n (%); Median (IQR) | |
#bold_labels é para deixar em negrito as categorias e italicize_levels é para deixar os resultados/opções das categorias em itálico# Aula 5
## Testes de hipóteses para dois grupos independentes
Vamos agora avaliar se existe diferença estatisticamente significativa entre as idades dos passageiros que sobreviveram e morreram no acidente do Titanic.
vimos dois tipos de testes de hipóteses:
-Teste t (paramétrico)
Exige normalidade em ambos os grupos . Testar com o teste de Shapiro-Wilk: se pelo menos um dos grupos não apresentar normalidade, usar teste de Mann-Whitney
Exige homocedasticidade . Testar com teste de Leavene . Se não tivermos homocedasticidade, usar teste t com welch
-Teste de Mann-Whitney (não paramétrico) Usar esse quando pelo menos uma variável não é normal (e vemos se é normal através do teste de shapiro) (se todas as variaveis são normais, usar o test t)
# avaliar a normalidade apenas da variavel quantitativa (a qualitativa é sobreviveu)
# dividir em dois grupos a variável quantitativa para realizar o teste (no caso, é a idade)
# o que nos interessa de ftao é o valor p (a probabilidade)
# o valor máximo de uma probabilidade é 1; logo, o valor p=7.816e^-08 (elevado a -8)
# se p<0.001, escrever apenas isso pois não precisa escrever o valor dado pelo teste
# nosso alfa é de 0.05; logo, é maior que o valor de p. Por conta disso, deve-se rejeitar a hipotese nula e ficar com a hipotese alternativa
shapiro.test(dados2$Idade[dados2$Sobreviveu=="Não"])
Shapiro-Wilk normality test
data: dados2$Idade[dados2$Sobreviveu == "Não"]
W = 0.96894, p-value = 7.816e-08Uma vez que o valor-p foi menor que 0,0001 para a idade no grupo de que não sobreviveram, temos que a variável não apresenta distribuição normal para este grupo, logo, teremos que utilizar o teste de Mann-Whitney.
Dessa forma, as hipóteses em estudo são:
H0: Não há diferença de idade entre os que morreram e os que sobreviveram H1: Há diferença de idade entre os que morreram e os que sobreviveram
#colocar a variável quantitativa em função da qualitativa
wilcox.test(Idade ~ Sobreviveu,
data = dados2,
paired= F)
Wilcoxon rank sum test with continuity correction
data: Idade by Sobreviveu
W = 65278, p-value = 0.1605
alternative hypothesis: true location shift is not equal to 0Uma vez que o teste de Mann-whitney apresentou valor = 0.16, nossos dados sugerem que não há diferença de idade entre os que não sobreviveram e os que morreram.
para avaliar se o segundo grupo tem distribuição normal só mudar o “Não” para “Sim” no teste de shapiro-wilk (normalidade acima de 0,05)
Sendo os dois grupos normais, precisamos avaliar a homocedasticidade.
library(car)Warning: package 'car' was built under R version 4.2.3Carregando pacotes exigidos: carDataWarning: package 'carData' was built under R version 4.2.3
Attaching package: 'car'The following object is masked from 'package:dplyr':
recodeleveneTest(Idade ~ Sobreviveu, data = dados2)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 1.1954 0.2746
712 Como o valores- p do teste de Levene foi igual a 0,27 temos que os grupos são homocedásticos e, portanto, faremos o teste t (nesta simulação).
H0: não há diferença entre as médias das idades dos que não morreram e sobreviveram H1; há diferença entre as médias das idades dos que não morreram e sobreviveram
t.test(Idade~Sobreviveu,
data= dados2,
var.equal = T)
Two Sample t-test
data: Idade by Sobreviveu
t = 2.0667, df = 712, p-value = 0.03912
alternative hypothesis: true difference in means between group Não and group Sim is not equal to 0
95 percent confidence interval:
0.114181 4.450798
sample estimates:
mean in group Não mean in group Sim
30.62618 28.34369 #vemos que as conclusoes dos dois testes feitos são incoerentes entre si. Logo, cada caso exige um tipo de teste diferente.
#Se o intervalo de confiança nao contem o zero, hã diferença. Se contém o zero, não há diferença entre as idades.Aplicando o teste t encontramos um valor-p igual a 0.039, indicando haver uma diferença estatisticamente significativa entre as médias.
Supondo que temos heterocedasticidade:
#Mesmo código do teste t, porém alterando T por F
t.test(Idade~Sobreviveu,
data= dados2,
var.equal = F)
Welch Two Sample t-test
data: Idade by Sobreviveu
t = 2.046, df = 598.84, p-value = 0.04119
alternative hypothesis: true difference in means between group Não and group Sim is not equal to 0
95 percent confidence interval:
0.09158472 4.47339446
sample estimates:
mean in group Não mean in group Sim
30.62618 28.34369 #Os dois testes levam á mesma conclusão de rejeitar a hipótese nulaAplicando o teste t encontramos um valor-p igual a 0.04, indicando haver uma diferença estatisticamente significativa entre as médias, sendo o grupo dos que não sobreviveram o grupo com maior média (igual a 30.6 anos)
Temos dois tipos de testes para avaliar grupos pareados (é usado grupo pareado pois os dados vieram das mesmas pessoas): - teste t pareado (paramétrico) . A diferença de valores precisa ser normal: avaliar comk o teste de Shapiro-wilk - teste de Wilcoxon ( não paramétrico) . Usar caso a diferença não seja normal
# Foi baixado o pacote BSDA e o pacote blood foi chamado
library(BSDA)Warning: package 'BSDA' was built under R version 4.2.3Carregando pacotes exigidos: lattice
Attaching package: 'BSDA'The following objects are masked from 'package:carData':
Vocab, WoolThe following object is masked from 'package:datasets':
Orangedata(Blood)Blood$dif <- Blood$machine - Blood$expertshapiro.test(Blood$dif)
Shapiro-Wilk normality test
data: Blood$dif
W = 0.92609, p-value = 0.2383# objeitov do teste é avaliar a normalidade
# resultado é 0.23
# conclusão é que é maior que 0,05
# hipotese nula = há normalidade ( não rejeita a hipotese nula logo há distribuição nomrmal)Uma vez que o teste de shapiro-wilk apresentou valor-p igual a 0,24, podemos inferir que os dados seguem distribuição normal, logo, aplicaremos o teste t pareado. Hipóteses:
H0: não há diferença entre as médias das pressões mensuradas por máquina e homem H1: há diferença entre as médias das pressões mensuradas por máquina e homem
t.test(Blood$machine,
Blood$expert,
paired= T)
Paired t-test
data: Blood$machine and Blood$expert
t = 0.68162, df = 14, p-value = 0.5066
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-2.146615 4.146615
sample estimates:
mean difference
1 Ao analisar a saída acima, qual foi o teste utilizado a partir das informações? Justifique. Foi o teste t variado, por conta do intervalo de confiança (algo q so aparece p testes parametricos), pela presença da média (se uso média, é um teste paramétrico).
A partir dessa saída, a pessoa rejeitou ou não a hipotese nula? Justifique. O intervalo de confiança contém o valor zero (-2 a 4), se contem o zero, não há diferença, logo, a pessoa não rejeitou a hipótese nula (pois é a hipotese nula que afirma que não há diferença).
É importante atrelar normalidade com média sempre!!
Sempre quando for falar o resultado do teste tem que ter o valor-p atrelado a resposta A partir do teste t pareado, podemos inferir que não há diferença estatisticamente significativa entre as médias das pressões mensuradas pela máquina e pelo homem (valor-p = 0,51)
H0: Não há diferença entre as pressões mensuradas por máquina e homem H1: Há diferença entre as pressões mensuradas por máquina e homem
wilcox.test(Blood$machine,
Blood$expert,
paired= T)
Wilcoxon signed rank test with continuity correction
data: Blood$machine and Blood$expert
V = 64, p-value = 0.489
alternative hypothesis: true location shift is not equal to 0A partir do teste de wilcoxon, podemos inferur qye não há diferença estatisticamente significativa entre as pressões mensuradas pela máquina e pelo homem ( valor-p = 0,49)
Existe diferença de idade em relação à classe economica? (quali X quanti: a palavra chave é diferença) H0: há normalidade HA: não há normalidade (normalidade é so pra variavel quantitativa) se p>o,05 não rejeita
Para analisar esta pergunta, temos 3 possibilidade: ANOVA, ANOVA com Welch e teste de Kruskal-Walls.
Para começar, precisamos avaliar a base normalidade dos 3 grupos:
shapiro.test(dados2$Idade[dados2$Classe
== "Primeira"])
Shapiro-Wilk normality test
data: dados2$Idade[dados2$Classe == "Primeira"]
W = 0.99169, p-value = 0.3643shapiro.test(dados2$Idade[dados2$Classe
== "Segunda"])
Shapiro-Wilk normality test
data: dados2$Idade[dados2$Classe == "Segunda"]
W = 0.97695, p-value = 0.005648Como o valor-p do teste de Shapiro-wil foi menor que 0,05 para a idade dos passageiros da segunda classe, teremos que aplicar o teste de Krusakai-wallis
H0: Não há diferença entre as idades segundo a classe econômica H1: Há diferença entre as idades segundo a classe econômica
kruskal.test(Idade ~ Classe, data = dados2)
Kruskal-Wallis rank sum test
data: Idade by Classe
Kruskal-Wallis chi-squared = 95.995, df = 2, p-value < 2.2e-16Uma vez que o teste de Kruskal-wallis (não possui pressupostos, ele é usado quando a normalidade é violada, porém não posso calcular a média pois ele n tem normalidade então pode existir outliers) apresentou valor-p < 0,001. Sabemos que pelo menos um dos grupos difere. Portanto, precisamos realizar o teste post-hoc de Dunn para verificar quais grupos diferem entre si.
library(DescTools)Warning: package 'DescTools' was built under R version 4.2.3
Attaching package: 'DescTools'The following object is masked from 'package:car':
RecodeDunnTest( Idade ~ Classe, data = dados2)
Dunn's test of multiple comparisons using rank sums : holm
mean.rank.diff pval
Segunda-Primeira -104.50956 3.2e-06 ***
Terceira-Primeira -182.19861 < 2e-16 ***
Terceira-Segunda -77.68906 4.8e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Após aplicar o teste de Dunn (saber quem difere de quem), foi possível observar que os três grupos diferem, sendo a primeira classe o grupo com maiores valores de idade, seguido pela segunda classe.
Precisamos avaliar se os grupos paresentam homocedasticidade.
library(car)
leveneTest(Idade ~ Classe, data = dados2)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 5.6202 0.003787 **
711
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p>0,05: há homocedasticidade p<0,05: não há homocedasticidade
O teste de Levene indicou um valor-p menor que o nível de significância, o teste correto a ser aplicado seria a ANOVA com Welch
H0: Não há diferença entre as médias das idades segundo a Classe econômica H1 Há diferença entre as médias das idades segundo a Classe econômica
oneway.test(Idade ~ Classe, data = dados2)
One-way analysis of means (not assuming equal variances)
data: Idade and Classe
F = 53.355, num df = 2.00, denom df = 359.44, p-value < 2.2e-16Visto que a ANOVA com Welch (requer normalidade) indicou que pelo menos 1 dos grupos difere (O valor de p é menor que 0,05; logo a hipótese nula foi rejeitada), precisamos identificar quais diferem entre eles. Para isso, usamos o teste de Tukey.
Obs: Caso tivéssemos homocedasticidade, faríamos a ANOVA, esse teste calcula a média entre as categorias de classes econômicas. Existe diferença nos perfis de idade em cada classe econômica? A NOVA requer normalidade e homocedasticidade.
anova<- aov(Idade ~ Classe, data = dados2)
summary(anova) Df Sum Sq Mean Sq F value Pr(>F)
Classe 2 20930 10465 57.44 <2e-16 ***
Residuals 711 129527 182
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
177 observations deleted due to missingnessRealizando o teste de Tukey, temos:
TukeyHSD(anova) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Idade ~ Classe, data = dados2)
$Classe
diff lwr upr p adj
Segunda-Primeira -8.355811 -11.704133 -5.007489 0.0000000
Terceira-Primeira -13.092821 -15.962198 -10.223445 0.0000000
Terceira-Segunda -4.737010 -7.676279 -1.797742 0.0004884Diff e upr é o intervalo de confiança pra cada diferença das médias. Em média a diferença entre as classes Segunda e Primeira é de 8 anos (tomando como exemplo a primeira linha).
A partir do teste de Tukey (aplicado quando há diferença entre os grupos para saber qual grupo difere) foi possível observar uma diferença significativa entre as médias de todas as comparações, indicando que a primeira classe apresenta média maior que a segunda, seguida da terceira classe.
Existe associação entre a classe econômica e o desfecho do passageiro? HO: não existe associação entre a classe econômica e o desfecho do passageiro H1: existe associação entre a classe econômica e o desfecho do passageiro
chisq.test(dados2$Classe, dados2$Sobreviveu)
Pearson's Chi-squared test
data: dados2$Classe and dados2$Sobreviveu
X-squared = 102.89, df = 2, p-value < 2.2e-16O valor de p é menor que 0.001. Logo, a hipótese nula é rejeitada e, consequentemente, existe associação entre a classe e o desfecho.
#para achar as porcentagens para saber os valores é necessario fazer uma tabela de contingência
library(summarytools)
ctable(dados2$Classe, dados2$Sobreviveu)Cross-Tabulation, Row Proportions
Classe * Sobreviveu
Data Frame: dados2
---------- ------------ ------------- ------------- --------------
Sobreviveu Não Sim Total
Classe
Primeira 80 (37.0%) 136 (63.0%) 216 (100.0%)
Segunda 97 (52.7%) 87 (47.3%) 184 (100.0%)
Terceira 372 (75.8%) 119 (24.2%) 491 (100.0%)
Total 549 (61.6%) 342 (38.4%) 891 (100.0%)
---------- ------------ ------------- ------------- --------------Visto que o teste qui-quadrado apresentou valor-p < 0,001, os dados sugerem que há uma assosiação estatisticamente significativa entre a classe econômica e o desfecho dos passageiros, indicando que passageiros da primeira classe sobreviveram mais, enquanto passageiros da terceira classe morreram mais
##CASO TIVESSE APARECIDO UMA MENSAGEM DE ALERTA/AVISO…
usariamos o teste exato de Fisher:
fisher.test(dados2$Classe, dados2$Sobreviveu)
Fisher's Exact Test for Count Data
data: dados2$Classe and dados2$Sobreviveu
p-value < 2.2e-16
alternative hypothesis: two.sidedExiste correlação entre a idade dos passageiros com o número de pais ou filhos a bordo?
H0: não há correlação entre a idade dos passageiros com o número de pais ou filhos a bordo H1: há correlação entre a idade dos passageiros com o número de pais ou filhos a bordo
Vamos começar analisando a normalidade da idade (depois, se necessário, analiso a normalidade do num de pais ou filhos)
shapiro.test(dados2$Idade)
Shapiro-Wilk normality test
data: dados2$Idade
W = 0.98146, p-value = 7.337e-08Para relembrar: A correlaçao, a associação e a diferença sao os testes finais. Antes, analisamos a homocedasticidade e a normalidade, para que, a partir do resultado destes, seja definido qual teste final prosseguir.
Visto que a variável idade não apresenta normalidade (p<0,001, hipotese nula foi rejeitada), o teste correto a ser realizado a ser realizado é o teste de correlação de Spearman.
cor.test(dados2$Idade, dados2$N_pais_filhos, method = "spearman")Warning in cor.test.default(dados2$Idade, dados2$N_pais_filhos, method =
"spearman"): Impossível calcular o valor exato de p com empates
Spearman's rank correlation rho
data: dados2$Idade and dados2$N_pais_filhos
S = 76087537, p-value = 5.409e-12
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2542121 O valor que rho (-0.25) representa é o coeficiente de correlação. 1- Ver de há uma correlação significativa 2- Ver se essa correlação é positiva ou negativa
Ao realizar o teste de Correlação de Spearman foi obtido um valor p<0,001, indicando uma correlação estatisticamente significativa entre as variáveis. O coeficiene de correlação rho = -0.25 indica que esta correlação é negativa, porém fraca. Dessa forma, quanto maior a idade do passageiro, menor o número de pais/filhos a bordo.
cor.test(dados2$Idade, dados2$N_pais_filhos, method = "pearson")
Pearson's product-moment correlation
data: dados2$Idade and dados2$N_pais_filhos
t = -5.1391, df = 712, p-value = 3.57e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2588990 -0.1173757
sample estimates:
cor
-0.1891193 O valor “cor” ou chamado de valor “r” (no trabalho, chamar de valor r) é referente ao valor do coeficiente de correlação. Perceba que ele está menor que o anterior pois esse teste não é o ideal para essa situação pois uma das variáveis não apresenta normalidade, logo, o valor não se adequa perfeitamente à variação real dos valores
##AULA 7 ##
vamos começar nossa análise contruindo um modelo de regressão linear. Para isso, nossa variável resposta precisa ser quantitativa, de preferência contínua. Vale lembrar que nossos resíduos devem ter distribuição normal em torno da média zero, ou seja, a maioria dos resíduos estarão proximos de zero (residuo = zero é aquele em que o ponto está exatamente em cima da reta). Obs.1: Valor significativo é aquele que o valor p é menor que 0,05%. Variáveis não significativas devem ser retiradas. Obs.2: Quando possuimos n categorias, serão geradas n-1 variáveis. Obs.3: A análise feita aqui sempre se dá através de comparações das variáveis em relação a uma variável-referência (o R escolhe a primeira variável que aparece nos levels. Caso eu queira escolher qual será a referência, preciso mudar a primeira nos levels).
rl<-lm(Tarifa~ Sexo+ Classe+Idade+N_pais_filhos, data = dados2)
summary(rl)
Call:
lm(formula = Tarifa ~ Sexo + Classe + Idade + N_pais_filhos,
data = dados2)
Residuals:
Min 1Q Median 3Q Max
-89.63 -9.45 -0.46 4.97 430.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 99.1396 5.6557 17.529 < 2e-16 ***
SexoMasculino -5.5850 3.3220 -1.681 0.09316 .
ClasseSegunda -69.0954 4.3887 -15.744 < 2e-16 ***
ClasseTerceira -78.7481 4.0287 -19.547 < 2e-16 ***
Idade -0.3369 0.1155 -2.916 0.00366 **
N_pais_filhos 11.5859 1.8631 6.219 8.57e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40.49 on 708 degrees of freedom
(177 observations deleted due to missingness)
Multiple R-squared: 0.4185, Adjusted R-squared: 0.4144
F-statistic: 101.9 on 5 and 708 DF, p-value: < 2.2e-16Ao observar os resultados do modelo, verificamos que apenas a variável sexo não foi significativa (valor-p > 0,05) Passageiros da segunda classe pagam, em média, 69 “moedas” a menos que os passageiros da primeira classe, enquanto os da terceira classe pagam, em média 79 moedas a menos comparados a essa mesma classe. Em relação a idade observamos que a cada um anos a mais reduz-se em média, a tarifa em 0.33 moedas. Quanto ao número de pais/filhos o aumento de uma unidade nesta variável aumenta a tarifa em 11.6 moedas. Em relação ao coeficiente de determinação (R2), este foi abaixo, explicando apenas 44% da variabilidade da variável resposta.
rl2<-lm(Tarifa~ + Classe+Idade+N_pais_filhos, data = dados2)
summary(rl2)
Call:
lm(formula = Tarifa ~ +Classe + Idade + N_pais_filhos, data = dados2)
Residuals:
Min 1Q Median 3Q Max
-88.19 -9.26 -0.96 4.71 428.23
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.6890 5.4717 17.671 < 2e-16 ***
ClasseSegunda -69.4477 4.3893 -15.822 < 2e-16 ***
ClasseTerceira -80.0308 3.9609 -20.205 < 2e-16 ***
Idade -0.3598 0.1149 -3.133 0.0018 **
N_pais_filhos 12.3098 1.8150 6.782 2.5e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40.55 on 709 degrees of freedom
(177 observations deleted due to missingness)
Multiple R-squared: 0.4162, Adjusted R-squared: 0.4129
F-statistic: 126.4 on 4 and 709 DF, p-value: < 2.2e-16Ao observar os resultados do modelo, sem a variável sexo, os passageiros da segunda classe pagam, em média, 69 “moedas” a menos que os passageiros da primeira classe, enquanto os da terceira classe pagam, em média 79 moedas a menos comparados a essa mesma classe. Em relação a idade observamos que a cada um anos a mais reduz-se em média, a tarifa em 0.33 moedas. Quanto ao número de pais/filhos o aumento de uma unidade nesta variável aumenta a tarifa em 11.6 moedas. Em relação ao coeficiente de determinação (R2), este foi abaixo, explicando apenas 44% da variabilidade da variável resposta.
Vamos avalias os pressupostos acima:
plot(rl2)
Através do gráfico 1, pode-se observar que na tarifa de 0 a 40, os valores estão próximos de 0, o que significa que as tarifas estão bem ajustadas nesses valores, ou seja, bem preditos pelo modelo. Não há normalidade pois há muitos valores fugindo de valores próximos de zero (vê-se nos valores de 60 a 120). Contudo, a linearidade talvez tenha sido respeitada pois não há padrões curvilíneos.
Já no gráfico 2, avaliamos a normalidade. O ideal é que todos os pontinhos estivem alinhados na reta tracejada. Os pontos fora da reta tracejada indicam uma cauda no gráfico de normalidade.
O gráfico 3 é indicado pra avaliar a homocedasticidade, vê-se um padrão de crescimento relacionado a um funil, indicando a violação da homocedasticidade (variancias iguais). Caso os pontos estivessem na linha, a variancia seria igual e a homocedasticidade seria ok.
O gráfico 4 demonstra o ponto de alavancagem.
Ao analisar os gráficos do modelos, observamos que os pressupostos de normalidade e homocedasticidade foram violados. Dessa forma, o modelo de regressão linear não é o mais indicado para modelar a variabilidade da tarifa com base na classe, idade e número de pais/filhos a bordo.
rlog <- glm(Sobreviveu ~ ., data= dados2 , family= binomial(link = "logit"))
summary(rlog)
Call:
glm(formula = Sobreviveu ~ ., family = binomial(link = "logit"),
data = dados2)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6966 -0.6750 -0.4061 0.6323 2.4481
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.809152 0.473171 8.050 8.26e-16 ***
ClasseSegunda -1.244302 0.314836 -3.952 7.74e-05 ***
ClasseTerceira -2.480356 0.331877 -7.474 7.80e-14 ***
SexoMasculino -2.583233 0.214780 -12.027 < 2e-16 ***
Idade -0.037821 0.007834 -4.828 1.38e-06 ***
N_pais_filhos -0.164902 0.119644 -1.378 0.168
Tarifa 0.001356 0.002404 0.564 0.573
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 964.52 on 713 degrees of freedom
Residual deviance: 645.25 on 707 degrees of freedom
(177 observations deleted due to missingness)
AIC: 659.25
Number of Fisher Scoring iterations: 5Visto que número de pais e filhos e tarifa não foram significativas, vamos fazer um segundo modelo sem estas variáveis.
rlog2 <- glm(Sobreviveu ~ . -N_pais_filhos -Tarifa, data= dados2 , family= binomial(link = "logit"))
summary(rlog2)
Call:
glm(formula = Sobreviveu ~ . - N_pais_filhos - Tarifa, family = binomial(link = "logit"),
data = dados2)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.7303 -0.6780 -0.3953 0.6485 2.4657
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.777013 0.401123 9.416 < 2e-16 ***
ClasseSegunda -1.309799 0.278066 -4.710 2.47e-06 ***
ClasseTerceira -2.580625 0.281442 -9.169 < 2e-16 ***
SexoMasculino -2.522781 0.207391 -12.164 < 2e-16 ***
Idade -0.036985 0.007656 -4.831 1.36e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 964.52 on 713 degrees of freedom
Residual deviance: 647.28 on 709 degrees of freedom
(177 observations deleted due to missingness)
AIC: 657.28
Number of Fisher Scoring iterations: 5library(gtsummary)
tbl_regression(rlog2, exponentiate=T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| Classe | |||
| Primeira | — | — | |
| Segunda | 0.27 | 0.16, 0.46 | <0.001 |
| Terceira | 0.08 | 0.04, 0.13 | <0.001 |
| Sexo | |||
| Feminino | — | — | |
| Masculino | 0.08 | 0.05, 0.12 | <0.001 |
| Idade | 0.96 | 0.95, 0.98 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
Por não termos uma varíavel explicativa quantitativa, aqui as análises não serão por comparação, mas sim por probabilidade (no caso, probabilidade de sobreviver). Na tabela, a variável com - é a nossa referência. Toda vez que tivermos um valor menor que 1, por exemplo, na segunda classe, se faz o seguinte: 0,27=0,73x100=73% (ou seja, o valor dado de 0,27 é igual a o que falta para chegar em 100, vezes 100 para chegar na porcentagem). Ou posso falar 0,27 vezes menor, mas o entendimento fica melhor se passar para porcentagem.
Contudo, quando ultrapassa 1, os valores podem ser ditos sem passar para porcentagem (o entendimento fica mais fácil do que quando é menor que 1). Mas também pode apresentar o valor em porcentagem (ex: 1,4: teve um aumento de 1,4 vezes ou teve um aumento de 40%).
A partir do modelo final, observamos que passageiros da segunda classe apresentam 73% menos chances de sobreviver do que os da primeira classe, já os da terceira classe apresentam 92% menos chances de sobreviver. Em relação ao sexo os homens tambem apresentam 92% menos chances de sobreviver quando comparados ás mulheres, o aumento de um ano na idade reduz as chacnes de sobrevivência em 4%
#install.packages("effects")
library(effects)Warning: package 'effects' was built under R version 4.2.3Use the command
lattice::trellis.par.set(effectsTheme())
to customize lattice options for effects plots.
See ?effectsTheme for details.plot(allEffects(rlog2))