Análise de Tabela de Contingência em R
MSP4112: Estatística Aplicada à Medicina
09 novembro 2025 00:15h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(asht, warn.conflicts=FALSE))
suppressMessages(library(boot, warn.conflicts=FALSE))
suppressMessages(library(coin, warn.conflicts=FALSE))
suppressMessages(library(confintr, warn.conflicts=FALSE))
suppressMessages(library(corrplot, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(diffdepprop, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(eiras2x2, warn.conflicts=FALSE))
suppressMessages(library(epiR, warn.conflicts=FALSE))
suppressMessages(library(exact2x2, warn.conflicts=FALSE))
suppressMessages(library(fishmethods, warn.conflicts=FALSE))
suppressMessages(library(FSA, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(haven, warn.conflicts=FALSE))
suppressMessages(library(irr, warn.conflicts=FALSE))
suppressMessages(library(irrCAC, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(misty, warn.conflicts=FALSE))
suppressMessages(library(plyr, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(rcompanion, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(RVAideMemoire, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
suppressMessages(library(tidyverse, warn.conflicts=FALSE))
suppressMessages(library(ufs, warn.conflicts=FALSE))
suppressMessages(library(vcd, warn.conflicts=FALSE))
suppressMessages(library(vcdExtra, warn.conflicts=FALSE))# pacote eiras2x2: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/5QTMBW
source("eiras.bartitle.R")
source("eiras.chisqgof.R")
source("eiras.cramerV.R")
source("eiras.matrix2dataframe.R")
source("eiras.show.MCSTAR.R")
source("friendlycolor.R")Adicionar
jmv::logLinear
Material
- HTML de Rmarkdown em
RPubs
Objetivos
- Identificar e realizar estatística descritiva de variáveis nominais e ordinais.
- Identificar distribuição amostral \(\chi^2\) (qui-quadrado).
- Formular hipótese estatística na análise de dados qualitativos.
- Identificar, aplicar e interpretar teste qui-quadrado de aderência.
- Definir e determinar frequências observada e esperada.
- Identificar e analisar tabela de contingência \(l\times c\) (número de linhas versus de colunas).
- Analisar tabela de contingência com variáveis ordinal e intervalar.
- Analisar tabela de contingência tridimensional \(l\times c \times k\), com a inclusão de variável de estratificação.
- Realizar os testes estatísticos em R.
Roteiro de funções R para análise de tabela de contingência
Este roteiro facilita a localização das principais funções R usadas neste texto, de acordo com o tipo de delineamento do estudo.
Delineamento entre participantes
Teste de independência
- Tabela unidimensional
- Uma variável nominal
- Teste de aderência:
chisq.test(simulate.p.value=TRUE, B=1e6,...)
- Teste de aderência:
- Uma variável nominal
- Tabela bidimensional
- Duas variáveis nominais dicotômicas
- Teste qui-quadrado de Pearson:
chisq.test(simulate.p.value=TRUE, B=1e6,...),confintr::ci_cramersv,effectsize::effectsize - Teste exato de Fisher de OR bruto:
exact2x2::fisher.exact,coin::chisq_test(..., distribution="exact"),confintr::ci_cramersv - Teste kappa de Cohen:
epiR::epi.kappa
- Teste qui-quadrado de Pearson:
- Duas variáveis nominais com pelo menos uma variável politômica
- Teste qui-quadrado de Pearson:
chisq.test(simulate.p.value=TRUE, B=1e6,...),coin::chisq_test,confintr::ci_cramersv,effectsize::effectsize - Teste qui-quadrado de Cochran-Armitage generalizado:
coin::chisq_test - Teste kappa de Cohen não ponderado:
psych::cohen.kappa - Uma variável nominal dicotômica e uma ordinal
- Teste qui-quadrado de Cochran-Armitage (trend):
stats::prop.trend.test,DescTools::CochranArmitageTest
- Teste qui-quadrado de Cochran-Armitage (trend):
- Teste qui-quadrado de Pearson:
- Uma variável nominal politômica e uma ordinal
- Teste qui-quadrado de Cochran-Armitage generalizado:
coin::chisq_test,coin::cmh_test,confintr::ci_cramersv
- Teste qui-quadrado de Cochran-Armitage generalizado:
- Duas variáveis ordinais
- Teste qui-quadrado de Mantel-Haenszel:
coin::lbl_test,coin::cmh_test,DescTools::MHChisqTest,vcdExtra::GKgamma
- Teste qui-quadrado de Mantel-Haenszel:
- Duas variáveis nominais dicotômicas
- Tabela tridimensional
- Duas variáveis nominais dicotômicas e uma variável de estratificação
nominal politômica
- Teste de Mantel-Haenszel de OR comum:
mantelhaen.test(exact=TRUE,...),epiR::2by2,coin::cmh_test
- Teste de Mantel-Haenszel de OR comum:
- Duas variáveis nominais politômicas e uma variável de estratificação
nominal politômica
- Teste de Mantel-Haenszel de OR comum generalizado:
coin::cmh_test
- Teste de Mantel-Haenszel de OR comum generalizado:
- Duas variáveis ordinais e uma variável de estratificação nominal
politômica
- Teste qui-quadrado ordinal estratificado (linear-by-linear
association model):
coin::lbl_test
- Teste qui-quadrado ordinal estratificado (linear-by-linear
association model):
- Duas variáveis nominais dicotômicas e uma variável de estratificação
nominal politômica
Teste de concordância
- Dois ou mais avaliadores e variável de desfecho nominal politômica
- Teste kappa de Fleiss:
irr::kappam.fleiss
- Teste kappa de Fleiss:
- Dois ou mais avaliadores e variável de desfecho intervalar
- Teste do coeficiente de correlação intraclasse (ICC):
irr::icc
- Teste do coeficiente de correlação intraclasse (ICC):
- Dois ou mais avaliadores e variável de desfecho nominal ou ordinal
politômica
- Testes AC1 (duas nominais) e AC2 (duas ordinais) de Gwet:
rrCAC::gwet.ac1.table,irrCAC::gwet.ac1.raw,eiras2x2::AC1
- Testes AC1 (duas nominais) e AC2 (duas ordinais) de Gwet:
- Dois ou mais avaliadores e duas variáveis intervalares politômicas
- Teste qui-quadrado de simetria de McNemar (pesos lineares):
FSA::ageBias,fishmethods::compare2 - Teste de simetria de Evans-Hoenig (projetado diagonalmente) (pesos
quadráticos):
FSA::ageBias,fishmethods::compare2 - Teste de simetria de Bowker (unpooled) (pesos cúbicos):
FSA::ageBias
- Teste qui-quadrado de simetria de McNemar (pesos lineares):
Delineamento intraparticipantes
- Tabela bidimensional
- Duas variáveis nominais dicotômicas
- Teste qui-quadrado de mudança de McNemar (\(H_0\): não mudança ou simetria):
coin::mh_test,exact2x2::mcnemar.exact - Intervalo de confiança para diferença de proporção em dados
pareados:
diffdepprop::diffpci
- Teste qui-quadrado de mudança de McNemar (\(H_0\): não mudança ou simetria):
- Duas variáveis nominais dicotômicas
Introdução
Observamos eventos que podem estar, de alguma forma, interligados. Por exemplo:
Pessoas chegam todos os dias ao pronto-socorro de um certo hospital e temos a impressão que o movimento no fim de semana não é igual ao dos outros dias ou que a segunda-feira tem menos procura:
- Hitzek, J et al (2022) Influence of Weekday and Seasonal Trends on Urgency and In-hospital Mortality of Emergency Department Patients. Frontiers in public health 10: 711235. https://doi.org/10.3389/fpubh.2022.711235
Alguns acreditam que o número de partos coincide com a mudança da fase da Lua (Stabolidou et al., 2008, Wake et al., 2010, Gudziunaite & Moshammer, 2022).
Pharma Hoje (2017) Vinte alimentos nutritivos que ajudam a manter o bom funcionamento do organismo:
A linhaça […] aumenta o HDL — conhecido também como o bom colesterol, fundamental para a formação de diversas estruturas em nossas células — e ainda reduz os triglicérides e o LDL (colesterol ruim).
O azeite de oliva extravirgem traz diversos benefícios. Ele atua inibindo os radicais livres e é um verdadeiro aliado quando o assunto é combater doenças. Também combate a oxidação do LDL, sendo um dos alimentos nutritivos que ajudam o bom funcionamento do organismo. Não é à toa que algumas nações, como os italianos e os franceses, são reconhecidamente saudáveis em todo o mundo. O hábito de utilizar o azeite de oliva em seu dia a dia, em vez de outros óleos menos benéficos, tem muito a ver com a disposição e a qualidade de vida desses povos.
Se gosta desse tipo de alimento [vegetais crucíferos], temos uma boa notícia: ao consumi-los, você está melhorando sua saúde! Alimentos como a couve-flor, a couve de Bruxelas, o agrião, o rabanete, o repolho e a mostarda são bons exemplos de vegetais crucíferos, ricos em vitaminas e minerais. Eles também ajudam na prevenção de várias doenças, inclusive do câncer — na medida em que são ricos em glicosinolatos.
Como saber se tais ocorrências não são apenas coincidências?
Contingência
Em sentido geral, contingência pode significar qualquer relação de dependência entre eventos ambientais ou entre eventos comportamentais e ambientais (Catania, 1993; Skinner, 1953, 1969; Todorov, 1985). Embora possa ser encontrado nos dicionários com diferentes significados, esse termo é empregado, na análise do comportamento, como termo técnico para enfatizar como a probabilidade de um evento pode ser afetada ou causada por outros eventos. (Catania, 1999, p. 368)
Teste qui-quadrado
Os testes baseados na estatística qui-quadrado de Pearson podem ser divididos em dois tipos:
Teste qui-quadrado de aderência de uma variável nominal com duas ou mais categorias: teste de proporções;
Teste qui-quadrado de independência entre duas variáveis nominais com duas ou mais categorias.
A decisão estatística destes testes utiliza a distribuição de mesmo nome (veja o Apêndice 1 para mais detalhes).
Teste qui-quadrado de aderência: uma variável nominal
O teste qui-quadrado de aderência é também conhecido por teste multinomial.
Em inglês é chamado também de Goodness of Fit (GoF).
Este teste permite descobrir se um conjunto de frequências observadas difere de um outro conjunto de frequências hipotetizadas.
Suposição:
Conforme Siegel & Castellan Jr. (1988, p. 111), independência das observações, i.e., cada observação da unidade observacional deve ser alocada em apenas uma categoria da variável nominal, é uma suposição necessária para validade do teste estatístico.
O teste está implementado em uma função criada por nós e chamada
chisqgof, disponível no arquivo eiras.chisqgof.R.
|
O uso de funções pode parecer coisa de programador, mas traz vantagens. Uma vez que uma função está organizada, os códigos que as utilizam são muito mais curtos, legíveis e fáceis de adaptar para casos particulares. Também tem a vantagem de que, uma vez bem desenvolvidas, podemos deixar de nos preocupar com os detalhes internos de uma função. Neste exemplo, a função foi definida emeiras.chisqgof.R.
|
Esta função recebe um dataframe com três colunas contendo, respectivamente, os nomes dados ao fator, as frequências observadas e as frequências hipotetizadas; por default, faz o teste robusto por bootstrapping com \(10^5\) reamostragens. O principal desta função é o teste em si. Encontre a linha que traz:
# Exemplo biológico: H-W com simulação
# Dados observados (AA, Aa, aa)
frq_obs <- c(50, 30, 20)
n <- sum(frq_obs)
# Frequências alélicas e H0 Hardy–Weinberg
p <- (2*frq_obs[1] + frq_obs[2]) / (2*n)
q <- 1 - p
probH0 <- c(p^2, 2*p*q, q^2)
set.seed(123)
gof <- chisq.test(x = frq_obs,
p = probH0,
simulate.p.value = TRUE,
B = 100000)
print(gof)
Chi-squared test for given probabilities with simulated p-value (based
on 1e+05 replicates)
data: frq_obs
X-squared = 11.605, df = NA, p-value = 0.0031Esta é a função nativa do R que faz o teste de aderência, recebendo
os valores em x e as probabilidades em probs
como dois vetores de igual tamanho (preparados, a partir do
dataframe, pelas linhas de código anteriores ao teste). Você
pode experimentar usar a função sozinha, mas como sua saída não é muito
instrutiva, aproveitamos para reorganizar os resultados no restante da
função e acrescentamos o cálculo do tamanho de efeito utilizando o
V de Cramér.
Exemplo 1: Preferência por marca de chocolate
Setenta e nove pessoas foram solicitadas a manifestar suas preferências com respeito às marcas de chocolate Dream Brown, Best Cocoa e Wonka. Queremos verificar se alguma das marcas é a preferida em detrimento de outras. No experimento, as marcas foram ocultadas dos participantes (estudo cego) e para os pesquisadores (duplo cego).
Caso os chocolates fossem igualmente preferidos pelas pessoas, deveríamos esperar aproximadamente a mesma proporção de pessoas com preferência para cada uma das marcas, i.e., \(1/3\) para cada marca de chocolate (\(H_0\)).
Na prática, dificilmente haverá exatamente o mesmo número de participantes em cada categoria, mas os números deveriam ter uma distribuição com razoável proximidade a 33.33% das pessoas com preferência para cada tipo de chocolate. O teste avalia tal proximidade, permitindo uma decisão estatística.
O teste está implementado em TesteQuiQuadradoAderencia_Chocolate.R,
resultando:
# TesteQuiQuadradoAderencia_Chocolate.R
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.chisq.gof.R")
Dados <- data.frame(
marca = c("Dream Brown", "Best Cocoa", "Wonka"),
prefreq = c(12, 23, 44),
probH0 = c(1/3, 1/3, 1/3),
stringsAsFactors = FALSE
)
seed <- 123
set.seed(seed)
alpha <- 0.05
B <- 1e3
res_gof <- chisq.gof(Dados,
alpha = alpha,
simulate.p.value = TRUE,
B = B)
cv <- CramerV.gof(Dados$prefreq,
Dados$probH0,
conf.level = 1-alpha,
B = B,
seed = seed)
cat(sprintf("\nCramer V: V=%.4f \nIC95%% analítico=[%.4f, %.4f]",
cv$V, cv$ci_analytic[1], cv$ci_analytic[2]))
if (all(!is.na(cv$ci_bca)))
cat(sprintf("\nIC95%% BCa=[%.4f, %.4f]\n", cv$ci_bca[1], cv$ci_bca[2])) else cat("\n")Goodness of Fit, robusto (p simulado)
Dados
label prefreq probH0
Dream Brown 12 0.3333333
Best Cocoa 23 0.3333333
Wonka 44 0.3333333
Análise de significância
X² = 20.0759 df = 2 p = 0.000999
p simulado baseado em B = 1000 reamostragens
Heurística: rejeita-se H0 se X²/df > 2; aqui X²/df = 10.038
Resíduos padronizados
Dream Brown Best Cocoa Wonka
-3.421 -0.796 4.216
Análise de significância prática
V de Cramér = 0.3565
IC95% analítico (não central): [0.1883, 0.5049]
IC95% bootstrap BCa: [0.2007, 0.5079]
Análise post-hoc (Bonferroni)
Pairwise comparisons using exact binomial tests
data: frq_obs
Dream Brown Best Cocoa
Best Cocoa 0.269 -
Wonka 6.3e-05 0.042
P value adjustment method: bonferroni
Inferência populacional, entre pares:
- Best Cocoa e Dream Brown: pode não haver diferença de proporção.
- Wonka e Dream Brown: diferença de proporção.
- Wonka e Best Cocoa: diferença de proporção.
Cramer V: V=0.3565
IC95% analítico=[0.1883, 0.5049]
IC95% BCa=[0.2007, 0.5079]Com \(p < 0.001\), rejeita-se \(H_0\). Portanto, temos evidência de que a preferência pelas marcas de chocolate não é homogênea.
Em outras palavras, o teste indica que as diferenças entre a frequência observada e a hipotetizada de \(1/3\) de preferência em cada marca é significante, improvável, o que é dado pelo valor p.
Observando a tabela, vemos que o desequilíbrio ocorreu pela maior preferência pelo chocolate da marca Wonka.
Exemplo 2: codominância × epistasia recessiva
As probabilidades relacionadas com \(H_0\) não precisam ser sempre todas iguais. Por exemplo, em genética, é comum precisarmos avaliar como é a herança de um genótipo a partir dos fenótipos observados.
Duas linhagens puras de plantas foram cruzadas: uma com pétalas amarelas e outra com pétalas vermelhas. A primeira geração, \(F_1\), tem pétalas cor de laranja. Então \(F_1\) é cruzada entre si, gerando \(F_2\) com 320 plantas, das quais 182 têm pétalas cor de laranja, 61 amarelas e 77 vermelhas.
Há possibilidades para explicar este resultado, e queremos saber qual é o caso para o fenótipo das pétalas das flores desta espécie:
1. codominância
Ocorre quando dois alelos ocupam um determinado locus gênico e combinam seus resultados. Suponha os seguintes alelos:
- Y - permite a produção de pigmento amarelo.
- R - permite a produção de pigmento vermelho.
As duas linhagens iniciais são puras (RR e YY) e, portanto, \(F_1\) é inteiramente composta por plantas RY (laranja). Ao cruzar as plantas \(F_1\) entre si, para \(F_2\) os seguintes genótipos e frequências esperadas são:
- YR: laranja (\(1/2\))
- YY: amarela (\(1/4\))
- RR: vermelha (\(1/4\))
Portanto seguindo as proporções 2:1:1, respectivamente 160 plantas com flores cor de laranja, 80 amarelas e 80 vermelhas.
Ocorre quando genes interagem para determinar a cor das pétalas. Os alelos Y estão em um locus gênico e os alelos R em outro, com segregação independente durante a meiose.
Vamos supor os seguintes genes e alelos:
- Gene Y: Responsável pela presença ou ausência de um pigmento
amarelo.
- Y (dominante) permite a produção de pigmento amarelo.
- y (recessivo) impede a produção de pigmento amarelo.
- Gene R: Modifica a cor do pigmento.
- R (dominante) transforma o pigmento amarelo em laranja.
- r (recessivo) deixa a cor original (amarelo, se houver pigmento, ou vermelho que é a cor de base).
Neste caso, para obter uma geração F1 homogênea (todas com pétalas laranja), as linhagens parentais precisam ser homozigotas para os genes envolvidos:
- Linhagem amarela (P1): O fenótipo amarelo indica que esta linhagem tem o alelo Y dominante, mas é homozigota recessiva para o gene R, o que resulta em flores amarelas. O genótipo será YYrr.
- Linhagem vermelha (P2): O fenótipo vermelho indica que esta linhagem é homozigota recessiva para o gene Y (não produz pigmento amarelo), independentemente do gene R. O genótipo será yyRR ou yyRr (neste exemplo, assumiremos yyRR para simplificação).
Assim, \(F_1\) é composta inteiramente por plantas YyRr (laranja). Ao serem cruzadas entre si, para \(F_2\) os seguintes genótipos e frequências esperadas são:
- Y-R-: laranja (\(9 \over 16\)) - o alelo Y permite a produção de pigmento amarelo, e o alelo R modifica para laranja.
- Y-rr: amarela (\(3 \over 16\)) - o alelo Y permite a produção de pigmento amarelo, mas rr não altera essa cor, permanecendo amarelo.
- yy–: vermelha (\(4 \over 16\)) - o genótipo yy não codifica pigmento amarelo, e a flor fica vermelha independentemente do outro gene.
Seguindo as proporções 9:3:4, respectivamente, 180 plantas com flores cor de laranja, 60 amarelas e 80 vermelhas.
A qual das duas hipóteses a herança da cor das flores desta espécie adere?
CorPetalas Frequencia Codominancia Epistasis
laranja 182 0.50 0.5625
vermelha 77 0.25 0.2500
amarela 61 0.25 0.1875As duas primeiras são as cores das flores e os números observados de
plantas, respectivamente similares às marcas e número de pessoas. A
terceira coluna tem as probabilidades \(0.50,
0.25, 0.25\), que são a \(H_0\)
da codominância. A quarta coluna tem \(\dfrac{9}{16}, \dfrac{4}{16},
\dfrac{3}{16}\), a \(H_0\) da
epistasia recessiva. Por isso, chamamos chisq.gof duas
vezes, passando as colunas c(1,2,3) e c(1,2,4)
para testarmos as duas possibilidades de herança genética.
O código TesteQuiQuadradoAderencia_Flores.R
produz a seguinte saída:
# Comparação de modelos de transmissão genética (Flores)
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))
# Carrega funções
source("eiras.chisq.gof.R") # deve conter chisq.gof, CramerV.gof e auxiliares
# Lê dados
if (file.exists("Flores.rds")) {
Dados <- readRDS("Flores.rds")
} else {
Dados <- data.frame(readxl::read_excel("flores.xlsx"))
Dados$CorPetalas <- factor(Dados$CorPetalas)
}
# Espera-se estrutura com:
# Coluna 1 = fator (cores)
# Coluna 2 = frequências observadas
# Coluna 3 = probH0 codominância (0.50, 0.25, 0.25)
# Coluna 4 = probH0 epistasia recessiva (9/16, 4/16, 3/16)
set.seed(123)
alpha <- 0.05
B <- 1e5
# Monta data.frames específicos (1,2,3) e (1,2,4)
Dados_cod <- data.frame(Dados[[1]], Dados[[2]], Dados[[3]])
names(Dados_cod) <- c(names(Dados)[1], names(Dados)[2], "probH0")
Dados_epi <- data.frame(Dados[[1]], Dados[[2]], Dados[[4]])
names(Dados_epi) <- c(names(Dados)[1], names(Dados)[2], "probH0")
# Roda GOF para cada hipótese
cat("\n=== Codominância (H0: 0.50, 0.25, 0.25) ===\n")
res_cod <- chisq.gof(Dados_cod, alpha = alpha, simulate.p.value = TRUE, B = B)
cv_cod <- CramerV.gof(Dados_cod[[2]], Dados_cod$probH0, conf.level = 1 - alpha, B = B, seed = 123)
cat("\n=== Epistasia recessiva (H0: 9/16, 4/16, 3/16) ===\n")
res_epi <- chisq.gof(Dados_epi, alpha = alpha, simulate.p.value = TRUE, B = B)
cv_epi <- CramerV.gof(Dados_epi[[2]], Dados_epi$probH0, conf.level = 1 - alpha, B = B, seed = 123)
# Resumo comparativo
sumario <- data.frame(
Modelo = c("Codominância", "Epistasia recessiva"),
X2 = c(res_cod$X2, res_epi$X2),
df = c(res_cod$df, res_epi$df),
p = c(res_cod$p_value, res_epi$p_value),
V = c(res_cod$V, res_epi$V),
V_L_analit = c(res_cod$ciV_analytic[1], res_epi$ciV_analytic[1]),
V_U_analit = c(res_cod$ciV_analytic[2], res_epi$ciV_analytic[2]),
V_L_BCa = c(cv_cod$ci_bca[1], cv_epi$ci_bca[1]),
V_U_BCa = c(cv_cod$ci_bca[2], cv_epi$ci_bca[2])
)
cat("\n=== Resumo comparativo (V e IC95%) ===\n")
num.cols <- vapply(sumario, is.numeric, TRUE)
sumario[num.cols] <- lapply(sumario[num.cols], function(x) round(x, 4))
print(sumario, row.names = FALSE)
=== Codominância (H0: 0.50, 0.25, 0.25) ===
Goodness of Fit, robusto (p simulado)
Dados
label prefreq probH0
laranja 182 0.50
vermelha 77 0.25
amarela 61 0.25
Análise de significância
X² = 7.6500 df = 2 p = 0.021580
p simulado baseado em B = 100000 reamostragens
Heurística: rejeita-se H0 se X²/df > 2; aqui X²/df = 3.825
Resíduos padronizados
laranja vermelha amarela
2.460 -0.387 -2.453
Análise de significância prática
V de Cramér = 0.1093
IC95% analítico (não central): [0.0108, 0.1812]
IC95% bootstrap BCa: [0.0323, 0.1719]
Análise post-hoc (Bonferroni)
Pairwise comparisons using exact binomial tests
data: frq_obs
laranja vermelha
vermelha 0.6 -
amarela 1.1e-14 1.7e-10
P value adjustment method: bonferroni
Inferência populacional, entre pares:
- vermelha e laranja: pode não haver diferença de proporção.
- amarela e laranja: diferença de proporção.
- amarela e vermelha: diferença de proporção.
=== Epistasia recessiva (H0: 9/16, 4/16, 3/16) ===Goodness of Fit, robusto (p simulado)
Dados
label prefreq probH0
laranja 182 0.5625
vermelha 77 0.2500
amarela 61 0.1875
Análise de significância
X² = 0.1514 df = 2 p = 0.931971
p simulado baseado em B = 100000 reamostragens
Heurística: rejeita-se H0 se X²/df > 2; aqui X²/df = 0.076
Resíduos padronizados
laranja vermelha amarela
0.225 -0.387 0.143
Análise de significância prática
V de Cramér = 0.0154
IC95% analítico (não central): [0.0000, 0.0589]
IC95% bootstrap BCa: [0.0000, 0.0230]
Análise post-hoc (Bonferroni)
Não aplicável (omnibus p = 0.931971)
=== Resumo comparativo (V e IC95%) ===
Modelo X2 df p V V_L_analit V_U_analit V_L_BCa
Codominância 7.6500 2 0.0216 0.1093 0.0108 0.1812 0.0323
Epistasia recessiva 0.1514 2 0.9320 0.0154 0.0000 0.0589 0.0000
V_U_BCa
0.1719
0.0230Neste exemplo, buscamos não rejeitar \(H_0\). Estes resultados sugerem que a herança das cores das pétalas desta espécie é explicada por epistasia recessiva. Esta é a conclusão porque, para o habitual \(\alpha=0.05\), \(H_0\) é rejeitada com as proporções esperadas pela codominância, com \(p \approx 0.02\), mas não é rejeitada com as proporções esperadas para epistasia recessiva, com \(p \approx 0.93\).
| Em teste de aderência, o número de graus de liberdade corresponde ao número de categorias menos 1. Portanto, para o exemplo da preferência pelos chocolates com 3 marcas ou para as flores com três cores possíveis, temos 2 graus de liberdade. Um teste para o número de pacientes atendidos de acordo com o dia da semana, de segunda a domingo, tem 6 graus de liberdade. |
Duas variáveis nominais dicotômicas
- tabela 2×2
- independência
- dependência positiva
- dependência negativa
Associação
O teste qui-quadrado de Pearson testa associação entre duas variáveis nominais; não testa causação.
- causa e efeito positivo?
- causa e efeito negativo?
- causa e efeito positivo?
- causa e efeito?
“Quando alguém lhe traz uma opinião coincidente com sua crença, você a toma como um fato. Quando alguém lhe traz um fato discordante de sua crença, você o toma como uma opinião.” Prof. Eduardo Massad
Exemplo 1: Gato e Alergia
- gato causa alergia …
- estresse
- mecanismo imune
- desinfetantes
- alergia e estresse
- alergia é a causa do gato
Exemplo 2: Vinho e Doença Coronariana
Conforme Pharma Hoje (2017) Vinte alimentos nutritivos que ajudam a manter o bom funcionamento do organismo:
“O suco de uva tinto integral - ou até mesmo o consumo moderado de vinho - pode fazer muito bem. Como reduz a pressão arterial e é rico em flavonoides, ele traz uma série de benefícios para quem o ingere. O resultado disso é a diminuição do risco de algumas doenças, como a aterosclerose.”
- ambiente inseguro
- moradia precária
- má nutrição
- boa moradia
- segurança
- boa alimentação
- outros cuidados
- é o vinho?
Doença coronariana: quais são as evidências?
Conforme Pandit K (2017) (grifos nossos):
“A educação do paciente é o processo pelo qual os profissionais de saúde transmitem informações aos pacientes e aos cuidadores que irão alterar seus comportamentos de saúde ou melhorar seu estado de saúde. Uma revisão recente da Cochrane nos diz que, atualmente, há poucas evidências de que a educação do paciente, como parte de um programa de reabilitação cardíaca, reduza os eventos relacionados ao coração e melhore a qualidade de vida relacionada à saúde. Na verdade, não há informações suficientes no momento para compreender completamente os benefícios ou malefícios da educação do paciente para pessoas com doenças cardíacas. No entanto, os resultados sugerem provisoriamente que as pessoas com doenças cardíacas devem receber reabilitação abrangente que inclui educação.
Um segundo componente da reabilitação é o exercício. Expliquei ao meu paciente a importância dos exercícios para seus resultados de curto e longo prazo. Outra revisão recente da Cochrane mostrou que a reabilitação cardíaca baseada em exercícios reduz o risco de morte por causa cardiovascular, diminui a duração da internação hospitalar e melhora a qualidade de vida relacionada à saúde quando comparada com aqueles que não praticam exercícios.
Ataques cardíacos e cirurgia cardíaca podem ser assustadores e traumáticos, podendo causar aumento do sofrimento psicológico. Uma terceira revisão recente da Cochrane avaliou a eficácia das intervenções psicológicas para doença coronariana. Embora as evidências tenham mostrado que as intervenções psicológicas, como parte da reabilitação cardíaca, não reduzem a mortalidade total de pessoas com doença coronariana, há algumas evidências de que podem aliviar os sintomas de depressão, ansiedade e estresse relatados pelo paciente [3]. Para prevenir e controlar esses sintomas no futuro, recomendei que meu paciente fizesse terapia psicológica regular. Mesmo que não tenha nenhum papel na redução do risco de mortalidade por doença coronariana, pode melhorar outros sintomas psicológicos importantes.”
Dieta de estilo mediterrâneo para a prevenção de doenças cardiovasculares
Conforme Rees, K et al. (2019) (grifos nossos):
“Esta revisão avaliou os efeitos de fornecer aconselhamento dietético para seguir uma dieta de estilo mediterrâneo ou fornecimento de alimentos relevantes para a dieta (ou ambos) para adultos saudáveis, pessoas com risco aumentado de doença cardiovascular e aquelas com doença cardiovascular, a fim de prevenir o ocorrência ou recorrência de doença cardiovascular e reduzir os fatores de risco a ela associados.
As definições de um padrão alimentar mediterrâneo variam e incluímos apenas ensaios clínicos randomizados de intervenções que relataram ambos os seguintes componentes principais: uma alta proporção de gordura monoinsaturada / saturada (uso de azeite de oliva como ingrediente principal para cozinhar e / ou consumo de outro alimentos tradicionais ricos em gorduras monoinsaturadas, como nozes, e uma alta ingestão de alimentos vegetais, incluindo frutas, vegetais e legumes. Componentes adicionais incluídos: baixo a moderado consumo de vinho tinto; alto consumo de grãos inteiros e cereais; baixo consumo de carnes e derivados e aumento do consumo de peixes; consumo moderado de leite e produtos lácteos.
[…]
Apesar do número relativamente grande de estudos incluídos nesta revisão, ainda há alguma incerteza em relação aos efeitos de uma dieta de estilo mediterrâneo nos desfechos clínicos e nos fatores de risco de doença cardiovascular para prevenção primária e secundária. A qualidade da evidência para os benefícios modestos sobre os fatores de risco de DCV na prevenção primária é baixa ou moderada, com um pequeno número de estudos relatando danos mínimos. Existem poucas evidências de prevenção secundária. Os estudos em andamento podem fornecer mais certeza no futuro.”
Exemplo 3: Tabagismo e Câncer de Pulmão
The Cigarette Century
Como se poderia atribuir o câncer de pulmão ao tabagismo?
“A impressão pública de incerteza científica e médica que se seguiu foi um elemento crucial para manter o mercado dos fumantes atuais, bem como para recrutar novos. A literatura da indústria, por exemplo, frequentemente apontava o fato de que não fumantes também podiam desenvolver câncer de pulmão. Portanto, argumentavam eles, como se poderia atribuir o câncer de pulmão ao tabagismo?”
O conhecimento médico é provisório.
“O conhecimento médico sempre foi provisório e contingente. Assim como medicamentos considerados “eficazes” não funcionam em todos os casos, uma causa de doença também não resulta sempre na doença. Como Richard Doll explicaria mais tarde, os epidemiologistas haviam identificado uma causa do câncer de pulmão (e de outras doenças), não a causa.”
Fisher era contra a associação entre o tabagismo e o câncer de pulmão.
“Outro crítico veemente das conclusões sobre o câncer de pulmão foi Sir Ronald Fisher, o principal biométrico e geneticista da Grã-Bretanha durante a primeira metade do século XX e um homem profundamente comprometido em aplicar a análise estatística à genética e à experimentação agrícola. Seu livro de 1925, Statistical Methods for Research Workers, rapidamente se tornou um clássico, levando-o a cargos acadêmicos no University College London e na Universidade de Cambridge. As críticas de Fisher eram semelhantes às de Berkson. A impossibilidade ética de realizar um experimento randomizado levou-o a questionar os resultados dos estudos epidemiológicos. Como crente em noções genéticas de causalidade do câncer, Fisher especulou que havia algum fator constitucional que levava certos indivíduos tanto a se tornarem fumantes quanto a desenvolverem câncer de pulmão, mesmo que fumar e o câncer de pulmão não estivessem causalmente relacionados. Doll e Hill refutaram repetidamente essa teoria, retornando à questão crucial de como explicar o aumento dos casos de câncer de pulmão durante o século XX se a doença fosse simplesmente “constitucional”.”
“Uma das inovações de [A. Bradford] Hill foi o primeiro ensaio clínico randomizado e duplo-cego, projetado para reduzir o viés do pesquisador na avaliação dos desfechos clínicos. […] Esse método, que se baseava na experimentação agrícola de Fisher em genética, tornou-se uma nova ferramenta essencial para a avaliação de tratamentos médicos.”
“Por volta da década de 1930, a relação entre o tabagismo e o câncer era um tema de debate não resolvido. Foi a indústria de seguros de vida, que, assim como as corporações de tabaco, cresceu rapidamente na primeira metade do século XX, que assumiu a liderança na compreensão dos efeitos do tabagismo sobre a saúde.
No início da década de 1950, porém, estava claramente evidente que as provas que implicavam o cigarro como risco à saúde pertenciam agora a uma outra ordem. Em primeiro lugar, a ligação entre o tabagismo e a doença era categórica, fora do campo do julgamento clínico individual.
A indústria do tabaco respondeu com uma nova e sem precedentes estratégia de relações públicas. Seu objetivo era produzir e sustentar o ceticismo e a controvérsia científica a fim de interromper o consenso emergente sobre os danos do tabagismo. Essa estratégia exigia intrusões no processo e nos procedimentos científicos. […] Enquanto houvesse aparência de dúvida, enquanto a indústria pudesse afirmar “não comprovado”, […] o futuro do cigarro dependeria agora da produção bem-sucedida de uma controvérsia científica.”
Teste qui-quadrado de Pearson de independência
O teste qui-quadrado de Pearson permite o teste da existência de associação ou efeito de interação entre duas variáveis qualitativas nominais com duas ou mais categorias.
Suposição:
Conforme Siegel & Castellan Jr. (1988, p. 111), o teste qui-quadrado de Pearson de independência necessita que cada observação da unidade observacional deve ser alocada em apenas uma categoria de cada uma das duas variáveis nominais.
Exemplo: capacete e trauma craniano
Considere um estudo que investiga a efetividade dos capacetes de segurança de bicicleta na prevenção de traumas cranianos. Os dados consistem de uma amostra aleatória de 793 indivíduos envolvidos em acidentes ciclísticos durante um certo período. Deseja-se testar, com \(\alpha=0.05\), se o uso do capacete tem funcionado como um fator de proteção.
Formulamos as hipóteses nula e alternativa:
\[ \begin{cases} H_0: \text{não há associação entre uso de capacete e trauma craniano}\\ H_1: \text{há associação entre uso de capacete e trauma craniano} \end{cases}\\ \alpha = 5\% \]
Os dados são:
tabela <- as.table(matrix(c( 17, 138,
130, 508),
nrow = 2,
byrow = TRUE))
colnames(tabela) <- c("Trauma+","Trauma-")
rownames(tabela) <- c("Capacete+","Capacete-")
print(tabela) Trauma+ Trauma-
Capacete+ 17 138
Capacete- 130 508Operacionalização do teste qui-quadrado de Pearson
Sob a hipótese nula (\(H_0\)), i.e., ausência de ssociação entre o uso de capacete e o trauma craniano nos acidentes, os valores esperados sob a hipótese nula (i.e., supondo-se que não existe associação entre o uso de capacete e o trauma craniano) são:
Trauma+ Trauma-
Capacete+ 28.7 126.3
Capacete- 118.3 519.7O teste é baseado na comparação entre os valores observados e esperados.
O resultado do teste, já executado com chisq.test e
armazenado em chi2, pode ser visto com:
Pearson's Chi-squared test
data: tabela
X-squared = 7.3099, df = 1, p-value = 0.006858alpha <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=tabela,
probs=c(alpha/2, 1-alpha/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.09601
Confidence interval:
2.5% 97.5%
0.000 0.149 print(effectsize::effectsize(chi2,
type="cramers_v",
ci=1-alpha,
alternative="two.sided"),
digits=4)Cramer's V (adj.) | 95% CI
------------------------------------
0.0892 | [0.0000, 0.1619]É mais fácil entender a interpretação deste teste, representando-o
graficamente (capacetetrauma_testeX2.R):
(para melhor visibilidade, computamos \(0 \le X^2 \le 9\) e truncamos o eixo das ordenadas em 1.2)
O interesse no teste qui-quadrado é saber a partir de qual valor de \(X^2\) a discrepância entre os valores observados e esperados sob \(H_0\) é menos provável que o valor de \(\alpha\) escolhido. Isto é encontrar o valor que separa a área de \(5\%\) sob a curva de distribuição em sua cauda direita. Este valor é chamado de qui-quadrado crítico, \(\chi^2_c\).
|
O raciocínio é bicaudal, buscando associação positiva ou negativa, mas a operacionalização do teste só utiliza a cauda superior. Isto acontece porque \(X^2=0\) correponde à igualdade numérica entre os valores observados e esperados sob \(H_0\). Como os valores de \(X^2\) aumentam com o aumento da discrepância entre os valores observados e esperados, utilizamos \(\alpha=0.05\) somente na cauda direita. |
Podemos encontrar \(\chi^2_c\) pelo complemento (\(1-\alpha=0.95\)) com
alfa <- 0.05
chi2.critico <- qchisq(p=1 - alfa, df=1)
cat("qui-quadrado crítico (df=1, alfa=5%): ", chi2.critico, "\n", sep="")qui-quadrado crítico (df=1, alfa=5%): 3.841459pois este é o valor de \(\chi^2\) que deixa \(95\%\) da área sob a curva da distribuição \(\chi^2\) à sua esquerda e \(\alpha=5\%\) à sua direita (área hachurada em azul). Equivalentemente, valores de \(X^2\) à esquerda de \(\chi^2_c\) estão na faixa de não rejeição de \(H_0\) e à direita está a faixa de rejeição de \(H_0\).
O teste obteve a estatística de teste com o valor \(X^2\)=7.3098819, com um grau de liberdade
(df=1), correspondendo ao valor p=0.0068576.
Comparando estes valores com \(X^2\) ou
\(\alpha\), respectivamente, podemos
decidir pela rejeição ou não-rejeição de \(H_0\).
|
Os graus de liberdade em uma tabela com \(L\) linhas e \(C\) colunas são dados por \[df = (L-1) \times (C-1)\] Portanto, uma tabela \(2\text{ x }2\) tem somente um grau de liberdade, e uma tabela \(4\text{ x }3\) tem seis graus de liberdade. |
Decisão pela estatística de teste
Para 1 grau de liberdade e \(\alpha=0.05\), computamos \(\chi^2_c(95\%)=\) 3.8414588. Neste exemplo rejeitamos \(H_0\) porque calculamos \(X^2=\) 7.3098819, valor maior que o valor crítico.
Decisão pelo valor p
À direita do valor calculado \(X^2=\) 7.3098819 define-se uma área, que
corresponde ao valor p (área hachurada em laranja), acessível
em chi2$p.value=0.0068576. Como \(p < \alpha\), concluimos pela rejeição
de \(H_0\).
|
As duas formas de decisão, pela comparação de \(\chi_c^2\) com \(X^2\) ou pela comparação entre \(\alpha\) e \(p\), são igualmente válidas e equivalentes. Exibir o valor p é a forma mais moderna e recomendada atualmente. |
Teste qui-quadrado de Pearson robusto
Há várias condições para a aplicação do teste qui-quadrado
(comentadas adiante) que podem, muitas vezes, restringir seu uso. Uma
alternativa é computar o teste em sua versão robusta. A função
chisq.test pode ser utilizada com:
alpha <- 0.05
seed <- 123
B <- 1e4
set.seed(123)
chi2.robusto <- chisq.test(tabela,
simulate.p.value=TRUE,
B=B)
print(chi2.robusto)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tabela
X-squared = 7.3099, df = NA, p-value = 0.007499set.seed(123)
print(confintr::ci_cramersv(tabela,
probs = c(alpha/2, 1-alpha/2),
type = "bootstrap",
R = B,
seed = seed),
digits=4)
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.09601
Confidence interval:
2.5% 97.5%
0.000 0.149 Com os parâmetros simulate.p.value=TRUE e
B=1e5, o teste foi feito com bootstrapping,
executando \(1 \times 10^5 = 100000\)
reamostragens.
|
Como a amostra é grande, os valores de \(X^2\) e \(p\) não mudam muito entre os dois métodos.
Há um defeito na implementação da versão robusta, que não calcula os
graus de liberdade (sabemos que |
A decisão, neste caso, não mudou. Tanto pelo critério do \(\chi^2_c(1-\alpha)\) quanto pelo \(\alpha\) rejeita-se \(H_0\).
Risco relativo
Há dois tipos principais de risco relativo: razão de riscos (RR, risk ratio) e razão de chances (OR, odds ratio).
Vamos considerar a seguinte tabela:
Com desfecho Sem desfecho
Com exposição a b a+b
Sem exposição c d c+d
a+c b+d totalEsta tabela de contingência genericamente relaciona exposição (e.g., hábito de fumar) com um efeito ou desfecho (e.g., câncer de pulmão). Traz as contagens dos:
- expostos que apresentaram o desfecho (\(a\)),
- expostos que não apresentaram o desfecho (\(b\)),
- não expostos que apresentaram o desfecho (\(c\)) e
- não expostos que não apresentaram o desfecho (\(d\)).
As marginais totalizam
- total de expostos (\(a+b\)),
- total de não expostos (\(c+d\)),
- total dos que apresentaram o desfecho (\(a+c\)) e
- total dos que não apresentaram o desfecho (\(b+d\)).
Razão de riscos (RR, risk ratio)
Risco é probabilidade. A razão de riscos é dada por:
\[ \Large{\text{RR} = {\dfrac{\dfrac{a}{a+b}}{{\dfrac{c}{c+d}}}}} \]
Repare que \(\dfrac{a}{a+b}\) é a probabilidade de um indivíduo exposto ter o desfecho (e.g., tabagista ter câncer); \(\dfrac{c}{c+d}\) é a probabilidade de um indivíduo não exposto ter o desfecho (e.g., não-tabagista ter câncer).
\(\text{RR}\), portanto, é um odds: quantas vezes mais é mais provável o desfecho (ocorrência de câncer) em quem é exposto (tabagista) em comparação com quem não é exposto (não tabagista). O foco do RR, portanto, está na exposição.
|
Odds e probabilidade Odds (traduzido habitualmente por chance) e probabilidade são formas equivalentes para expressar possibilidades e relacionadas por: \[ \textit{Odds} = {{\text{Probabilidade}} \over{1-\text{Probabilidade}}} \] e \[ \text{Probabilidade} = {{\textit{Odds}} \over{1+\textit{Odds}}} \] Por exemplo, uma probabilidade de \(80\%\) corresponde a \[ \textit{Odds} = {\dfrac{0.8}{1-0.8}} = \dfrac{0.8}{0.2} = 4 \] indicando que dizer “\(80\%\) de probabilidade” é equivalente a dizer “4 vezes mais chance de ocorrer do que não ocorrer” um determinado evento. Resersamente, Odds de 2 é: \[ \text{Probabilidade} = \dfrac{2}{1+2} = \dfrac{2}{3} \approx 66.67\% \] e \(\textit{Odds}=1\) resulta em: \[ \text{Probabilidade} = \dfrac{1}{1+1} = \dfrac{1}{2} = 50\% \] Associa-se o máximo de incerteza com probabilidade de \(50\%\); por este motivo, \(Odds=1\) será um valor necessário para decisões estatísticas, adiante. |
Razão de chances (OR, odds ratio)
Odds ratio, como diz o nome, é uma razão entre dois odds. Ainda considerando a tabela:
Com desfecho Sem desfecho
Com exposição a b a+b
Sem exposição c d c+d
a+c b+d total\(\text{OR}\) é dado por:
\[ \Large{\text{OR} = \dfrac{ \dfrac{ \dfrac{a}{a+c} }{ \dfrac{c}{a+c} } }{ \dfrac{ \dfrac{b}{b+d} }{ \dfrac{d}{b+d} } } } \]
O numerador de \(\text{OR}\)
- \(\dfrac{\dfrac{a}{a+c}}{\dfrac{c}{a+c}}\) é odds entre os que tiveram o desfecho em exposto e não exposto; por exemplo, quantas vezes é mais provável ter câncer para tabagista em relação ao não tabagista.
e o denominador de \(OR\)
- \(\dfrac{\dfrac{b}{b+d}}{\dfrac{d}{b+d}}\) é odds entre os que não tiveram o desfecho em exposto e não exposto; por exemplo, quantas vezes é mais provável não ter câncer para tabagista em relação a não tabagista.
Para os dois grupos, consideram-se os odds da exposição (quantas vezes é mais provável ser exposto) e verifica quantas vezes maior é o odds de ser exposto dos que tiveram o desfecho em relação aos que não tiveram o desfecho. O foco do OR, portanto, está nos odds (nas chances) do desfecho. Sua interpretação, especialmente para quem não está acostumado a pensar em odds, é convoluta.
No entanto, é interessante (com alguma álgebra) verificar que \(\text{OR}\) reduz-se a:
\[ \Large{\text{OR} = \dfrac{a \cdot d}{b \cdot c}} \]
fácil de lembrar e assim conhecido como “produto cruzado” (com o defeito de ocultar o motivo pelo qual é um odds ratio). Como produto cruzado, percebe-se que é uma medida de quanto pesa a:
- diagonal principal, do que é concordante: \(\text{desfecho entre expostos}~~\text{E}~~\text{não desfecho entre não expostos}\)
em relação à
- diagonal secundária, do que é discordante: \(\text{desfecho entre não expostos}~~\text{E}~~\text{não desfecho entre expostos}\).
Esta é uma forma mais fácil de interpretar \(\text{OR}\).
|
Os valores de \(\text{RR}\) e \(\text{OR}\) são próximos quando menos que \(5\%\) dos não expostos têm desfecho. |
No R nativo, existe a função fisher.test para o cálculo
de \(\text{OR}\). No entanto, esta
função tem defeito no cálculo do intervalo de confiança. Usaremos a
função exact2x2::fisher.exact. No exemplo de capacete e
trauma, temos:
library(exact2x2)
tabela <- as.table(matrix(c(17, 138,
130, 508),
nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma+","Trauma-")
rownames(tabela) <- c("Capacete+","Capacete-")
print(tabela)
ft <- exact2x2::fisher.exact(tabela)
print(ft, digits=3)Carregando pacotes exigidos: exactciCarregando pacotes exigidos: ssanvCarregando pacotes exigidos: testthat Trauma+ Trauma-
Capacete+ 17 138
Capacete- 130 508
Two-sided Fisher's Exact Test (usual method using minimum likelihood)
data: tabela
p-value = 0.006
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.279 0.833
sample estimates:
odds ratio
0.482
|
A função nativa Existe uma alternativa, implementada em
Para ver qual é o problema, usando outra tabela com menos valores, compare: Neste caso, a função nativa produz um intervalo de confiança que inclui o valor unitário em contradição com o valor p (aliás, igual nas duas funções). Como |
Há também implementações de OR no pacote DescTools. No
entanto, a implementação usando método mle não é
recomendável, pois tem o mesmo problema da função nativa
fisher.test.
library(exact2x2)
# Crude RR+ & OR (Wald)
# OR por Wald: not because it is the best, but because it is the most commonly used
alpha <- 0.05
tabela <- as.table(matrix(c(2, 14,
12, 15), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma+","Trauma-")
rownames(tabela) <- c("Capacete+","Capacete-")
print(tabela) Trauma+ Trauma-
Capacete+ 2 14
Capacete- 12 15odds ratio lwr.ci upr.ci
0.1786 0.0338 0.9436 odds ratio lwr.ci upr.ci
0.1945 0.0244 0.9018 odds ratio lwr.ci upr.ci
0.1854 0.0172 1.0645 qui2 <- chisq.test(tabela, correct=FALSE)
print(effectsize::effectsize(qui2,
type="oddsratio",
ci=1-alpha,
alternative="two.sided"),
digits=4)Odds ratio | 95% CI
-----------------------------
0.1786 | [0.0338, 0.9436]
Fisher's Exact Test for Count Data
data: tabela
p-value = 0.04471
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.01722368 1.06448721
sample estimates:
odds ratio
0.1854297
Two-sided Fisher's Exact Test (usual method using minimum likelihood)
data: tabela
p-value = 0.04471
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.0255 0.9252
sample estimates:
odds ratio
0.1854297 Decisão pelo intervalo de confiança
A estimativa pontual de OR está em
ft$estimate=0.4817645, mas para a decisão é obrigatório
verificar seu intervalo de confiança de 95%, com limites inferior e
superior respectivamente guardados em ft$conf.int[1]=0.2791
e ft$conf.int[2]=0.8334, pois este é o intervalo onde
confiamos estar o valor populacional (que nunca saberemos).
Como a ausência de efeito é dada por \(\text{OR}=1\) (correspondendo a \(H_0\)) e este valor unitário está fora do intervalo, rejeita-se \(H_0\). Além disto, como o intervalo de confiança de 95% está à esquerda do \(1\), o uso de capacete é protetor, associado com redução dos traumas cranianos.
|
Caso a tabela fosse construída com as colunas em posições trocadas (compare com a tabela anterior), obter-se-ia o seguinte: OR 2.0757031 é agora maior que o valor unitário, o qual está fora do intervalo (1.2 e 3.5833), rejeitando-se a hipótese nula e a associação é positiva: a conclusão é a mesma, o uso de capacete está associado com o aumento de “Trauma-”, portanto é protetor. Aliás, observe que são os mesmos valores obtidos anteriormente e guardados no objeto
|
|
Lembre-se de nunca decidir a respeito de OR sem considerar o intervalo de confiança. Por isso, evite fazer o cálculo manualmente: use R! |
Risco relativo em epiR::epi.2by2
O pacote epiR produz estimativas pontuais e por
intervalo de confiança de OR e RR para alguns delineamentos.
Incidence OR (
cohort.count) e Prevalence OR (POR) (cross.sectional)Incidence RR (
cohort.count) e Prevalence RR (PRR ou PR) (cross.sectional)
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
tabela <- as.table(matrix(c(17, 138,
130, 508),
nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma+","Trauma-")
rownames(tabela) <- c("Capacete+","Capacete-")
print(tabela)
print(epiR::epi.2by2(dat=tabela,
method="cohort.count",
units=100,
digits=3,
outcome="as.columns")) Trauma+ Trauma-
Capacete+ 17 138
Capacete- 130 508
Outcome+ Outcome- Total Inc risk *
Exposure+ 17 138 155 10.968 (6.520 to 16.979)
Exposure- 130 508 638 20.376 (17.315 to 23.714)
Total 147 646 793 18.537 (15.891 to 21.420)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 0.538 (0.335, 0.865)
Inc odds ratio 0.481 (0.281, 0.826)
Attrib risk in the exposed * -9.408 (-15.237, -3.580)
Attrib fraction in the exposed (%) -85.783 (-199.767, -17.562)
Attrib risk in the population * -1.839 (-5.972, 2.294)
Attrib fraction in the population (%) -9.920 (-10.709, -8.967)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units Para chamar a função epiR::epi.2by2 precisamos obedecer
seu padrão, para que adote os valores de exposição e desfecho
corretamente:
- a primeira linha e primeira coluna são tratados como eventos
positivos. Tabelas que não estão nesta ordem podem ser alteradas com o
uso de
DescTools::Rev(tabela,margin=…) - o parâmetro
outcome=“as.columns”indica que, nesta tabela, as colunas são as categorias do desfecho - o delinamento dado por
method=“cohort.count”informa que a tabela tem contagens e, por isso, calcula \(\text{RR}\) e \(\text{OR}\) units=100pede para que a saída seja dada em porcentagem (para as estimativas que forem relevantemente apresentadas em porcentagem)
As duas primeiras linhas da saída desta função trazem o \(\text{RR}\) (risk ratio) e o \(\text{OR}\) (odds ratio).
O raciocínio pode ficar confuso com a terminologia “Capacete+” e “Capacete-”; vamos simplificar pensando na falta do uso do capacete como um fator de exposição, aumentando o risco de trauma craniano em acidentes.
(note o uso de DescTools::Rev(tabela,margin=1) para
alterar a ordem das linhas, ajeitando a tabela para
epiR::epi.2by2)
Este \(\text{RR}\) tradicional é o \(\text{RR}_+\), visto acima. \(\text{RR}_+\) refere-se ao desfecho positivo (a ocorrência do desfecho entre os expostos em comparação com a ocorrência do desfecho entre os não expostos). Podemos conceituar o \(\text{RR}_-\) para o desfecho negativo (a não ocorrência do desfecho entre os expostos em comparação com a não ocorrência do desfecho entre os não expostos). Algebricamente, como
\[ \text{OR} = \dfrac{\text{RR}_+}{\text{RR}_-} \]
O risco relativo negativo é dado por:
\[ \text{RR}_- = \dfrac{\text{RR}_+}{\text{OR}} \]
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
tabela <- as.table(matrix(c(17, 138,
130, 508),
nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma+","Trauma-")
rownames(tabela) <- c("Capacete+","Capacete-")
print(tabela <- DescTools::Rev(tabela, margin=1))
print(epiR::epi.2by2(dat=tabela,
method="cohort.count",
units=100,
digits=3,
outcome="as.columns"))
a <- tabela[1,1]; b <- tabela[1,2]; c <- tabela[2,1]; d <- tabela[2,2]
RRp <- (a/(a+b))/(c/(c+d))
RRm <- (b/(a+b))/(d/(c+d))
OR <- (a*d)/(b*c)
cat(sprintf("RR+ = %.6f\nRR- = %.6f\nOR = %.6f\nRR+/RR- = %.6f\nRR+/OR = %.6f\n",
RRp, RRm, OR, RRp/RRm, RRp/OR))
all.equal(OR, RRp/RRm) # TRUE
all.equal(RRm, RRp/OR) # TRUE
# RR+ e OR na tabela original
resP <- epiR::epi.2by2(dat=tabela,
method="cohort.count",
units=100,
digits=3,
outcome="as.columns")
print(resP)
RRp <- as.numeric(resP$massoc.summary$est[1])
OR <- as.numeric(resP$massoc.summary$est[2])
# RR− invertendo as linhas (Capacete- como “exposto”)
tabelaR <- DescTools::Rev(tabela, margin=2)
resM <- epiR::epi.2by2(dat=tabelaR,
method="cohort.count",
units=100,
digits=3,
outcome="as.columns")
print(resM)
RRm <- as.numeric(resM$massoc.summary$est[1])
# Checagem das identidades
print(sprintf("RR+ = %.6f RR- = %.6f OR = %.6f", RRp, RRm, OR))
cat(sprintf("OR ≈ RR+/RR- : %.6f\n", RRp / RRm))
cat(sprintf("RR- ≈ RR+/OR : %.6f\n", RRp / OR)) Trauma+ Trauma-
Capacete- 130 508
Capacete+ 17 138
Outcome+ Outcome- Total Inc risk *
Exposure+ 130 508 638 20.376 (17.315 to 23.714)
Exposure- 17 138 155 10.968 (6.520 to 16.979)
Total 147 646 793 18.537 (15.891 to 21.420)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 1.858 (1.156, 2.985)
Inc odds ratio 2.077 (1.211, 3.562)
Attrib risk in the exposed * 9.408 (3.580, 15.237)
Attrib fraction in the exposed (%) 46.174 (14.938, 66.641)
Attrib risk in the population * 7.569 (1.956, 13.183)
Attrib fraction in the population (%) 40.834 (20.733, 58.967)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units
RR+ = 1.857828
RR- = 0.894326
OR = 2.077351
RR+/RR- = 2.077351
RR+/OR = 0.894326
Outcome+ Outcome- Total Inc risk *
Exposure+ 130 508 638 20.376 (17.315 to 23.714)
Exposure- 17 138 155 10.968 (6.520 to 16.979)
Total 147 646 793 18.537 (15.891 to 21.420)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 1.858 (1.156, 2.985)
Inc odds ratio 2.077 (1.211, 3.562)
Attrib risk in the exposed * 9.408 (3.580, 15.237)
Attrib fraction in the exposed (%) 46.174 (14.938, 66.641)
Attrib risk in the population * 7.569 (1.956, 13.183)
Attrib fraction in the population (%) 40.834 (20.733, 58.967)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units
Outcome+ Outcome- Total Inc risk *
Exposure+ 508 130 638 79.624 (76.286 to 82.685)
Exposure- 138 17 155 89.032 (83.021 to 93.480)
Total 646 147 793 81.463 (78.580 to 84.109)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 0.894 (0.836, 0.957)
Inc odds ratio 0.481 (0.281, 0.826)
Attrib risk in the exposed * -9.408 (-15.237, -3.580)
Attrib fraction in the exposed (%) -11.816 (-19.010, -3.553)
Attrib risk in the population * -7.569 (-13.183, -1.956)
Attrib fraction in the population (%) -9.292 (-11.140, -5.652)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units
[1] "RR+ = 1.857828 RR- = 0.894326 OR = 2.077351"
OR ≈ RR+/RR- : 2.077351
RR- ≈ RR+/OR : 0.894326Riscos relativos por regressão
Há críticas sobre estas formas tradicionais para o cálculo dos riscos relativos:
Apesar das críticas deste artigo, as estimativas de \(\text{OR}\), \(\text{RR}_+\) e \(\text{RR}_-\) com seus intervalos de confiança podem ser obtidos por meio de regressões.
Modelo linear generalizado (GzLM): glm e
logbin::logbin
- Função de Distribuição (FD) da VD
- Família de distribuições exponenciais
- E.g.: normal ou gaussiana, binomial (sucesso/fracasso), Poisson (eventos raros: acidente, incidência de doença), gama (valores positivos com assimetria positiva: duração de vida), gaussiana inversa (valores positivos com assimetria positiva: duração de casamento) etc.
- Função de ligação (FL)
- Função não-linear que estabelece uma relação linear entre a média condicionada da VD e a combinação linear das VIs
A combinação entre a distribuição da VD e a função de ligação resulta em diferentes modelos:
- Linear geral (GLM) (ANOVA & Regressão)
- FD: normal
- FL: identidade
- Regressão logística binária
- FD: binomial
- FL: logit
- Probit
- FD: binomial
- FL: probit
- Loglinear
- FD: Poisson
- FL: log
- Regressão ordinal
- FD: multinomial
- FL: logit
OR_RR_Regressao_Logistica_Poisson_Logbinomial.R
Trauma+ Trauma-
Capacete+ 17 138
Capacete- 130 508 Trauma+ Trauma-
Capacete- 130 508
Capacete+ 17 138 Outcome+ Outcome- Total Prev risk *
Exposure+ 130 508 638 20.38 (17.32 to 23.71)
Exposure- 17 138 155 10.97 (6.52 to 16.98)
Total 147 646 793 18.54 (15.89 to 21.42)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Prev risk ratio 1.86 (1.16, 2.98)
Prev odds ratio 2.08 (1.21, 3.56)
Attrib prev in the exposed * 9.41 (3.58, 15.24)
Attrib fraction in the exposed (%) 46.17 (14.94, 66.64)
Attrib prev in the population * 7.57 (1.96, 13.18)
Attrib fraction in the population (%) 40.83 (20.73, 58.97)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units Trauma- Trauma+
Capacete- 508 130
Capacete+ 138 17 Outcome+ Outcome- Total Prev risk *
Exposure+ 508 130 638 79.62 (76.29 to 82.68)
Exposure- 138 17 155 89.03 (83.02 to 93.48)
Total 646 147 793 81.46 (78.58 to 84.11)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Prev risk ratio 0.89 (0.84, 0.96)
Prev odds ratio 0.48 (0.28, 0.83)
Attrib prev in the exposed * -9.41 (-15.24, -3.58)
Attrib fraction in the exposed (%) -11.82 (-19.01, -3.55)
Attrib prev in the population * -7.57 (-13.18, -1.96)
Attrib fraction in the population (%) -9.29 (-11.14, -5.65)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 7.310 Pr>chi2 = 0.007
Fisher exact test that OR = 1: Pr>chi2 = 0.006
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units
Converte a tabela em data frame: Capacete Trauma
1 1 1
2 1 1
3 1 1
4 1 1
5 1 1
6 1 1 Capacete Trauma
788 0 0
789 0 0
790 0 0
791 0 0
792 0 0
793 0 0
Codificando (conferindo):
1=Capacete-
0=Capacete+
1=Trauma+
0=Trauma- Trauma
Capacete 1 0
1 130 508
0 17 138
Logistic Regression: OR
Call:
glm(formula = Trauma ~ Capacete, family = binomial(link = "logit"),
data = dt_captrm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.0940 0.2570 -8.147 3.73e-16 ***
Capacete 0.7311 0.2752 2.657 0.00789 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 760.40 on 792 degrees of freedom
Residual deviance: 752.32 on 791 degrees of freedom
AIC: 756.32
Number of Fisher Scoring iterations: 4
OR = 2.08, IC95 = [1.24, 3.68]
Poisson Regression: RR+
Call:
glm(formula = Trauma ~ Capacete, family = poisson(link = "log"),
data = dt_captrm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.2102 0.2425 -9.113 <2e-16 ***
Capacete 0.6194 0.2579 2.402 0.0163 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 495.50 on 792 degrees of freedom
Residual deviance: 488.76 on 791 degrees of freedom
AIC: 786.76
Number of Fisher Scoring iterations: 6Waiting for profiling to be done...
RR+ = 1.86, IC95(RR+) = [1.15, 3.19]
Crude RR-: Poisson regression
Call:
glm(formula = Trauma_rev ~ Capacete, family = poisson, data = dt_captrm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.11617 0.08513 -1.365 0.172
Capacete -0.11169 0.09599 -1.163 0.245
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 264.89 on 792 degrees of freedom
Residual deviance: 263.57 on 791 degrees of freedom
AIC: 1559.6
Number of Fisher Scoring iterations: 4Waiting for profiling to be done...
Regressao Poisson
Crude RR- = 0.89IC95(RR-) = [0.7433, 1.0833]
log-binomial Regression: RR+logbin parameterisation 1/1
Deviance = 752.3221 Iterations - 6
Call:
logbin::logbin(formula = Trauma ~ Capacete, data = dt_captrm,
method = "em", accelerate = "squarem", trace = 1)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.6751 -0.6751 -0.6751 -0.4820 2.1025
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.2102 0.2288 -9.658 <2e-16 ***
Capacete 0.6194 0.2419 2.561 0.0104 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Null deviance: 760.40 on 792 degrees of freedom
Residual deviance: 752.32 on 791 degrees of freedom
AIC: 756.32
AIC_c: 756.34
Number of iterations: 6
Crude RR+ = 1.86IC95(RR+) = [1.16, 2.98]
Crude RR-: log-binomial regressionlogbin parameterisation 1/1
Deviance = 752.3221 Iterations - 8
Call:
logbin::logbin(formula = Trauma_rev ~ Capacete, data = dt_captrm,
method = "em", accelerate = "squarem", trace = 1)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1025 0.4820 0.6751 0.6751 0.6751
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.11617 0.02819 -4.121 3.78e-05 ***
Capacete -0.11169 0.03458 -3.230 0.00124 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Null deviance: 760.40 on 792 degrees of freedom
Residual deviance: 752.32 on 791 degrees of freedom
AIC: 756.32
AIC_c: 756.34
Number of iterations: 8
Regressao log-binomial
Crude RR- = 0.89IC95(RR-) = [0.8357, 0.957]
OR (estimativa, IC95%): medida metodo est LI LS
OR Logística (crude) 2.077 1.242 3.677
OR epiR (tabela, crude) 2.077 1.211 3.562
RR+ (estimativa, IC95%): medida metodo est LI LS
RR+ Poisson (crude) 1.858 1.154 3.193
RR+ Log-binomial (crude) 1.858 1.156 2.985
RR+ epiR (tabela, crude) 1.858 1.156 2.985
RR- (estimativa, IC95%): medida metodo est LI LS
RR- epiR (tabela, crude) 0.894 0.836 0.957
RR- Poisson (crude, desfecho invertido) 0.894 0.743 1.083
RR- Log-binomial (crude, desfecho invertido) 0.894 0.836 0.957Saída em arquivo-texto com sink
Em tabelacontingencia2x2_capacetetrauma.R,
que você pode adaptar para as suas necessidades, processamos o mesmo
exemplo, utilizando a função sink que desvia o resultado da
tela para um arquivo texto. Além da estatística \(X^2\), computamos outros valores
interessantes. Execute este código R:
# Criando tabela de contingencia 2x2 e realizando
# testes qui-quadrado de Pearson e de razao de chances (OR)
alfa <- 0.05
source("eiras.show.MCSTAR.R")
# para guardar o resultado em um arquivo txt
# substitua por NA para emitir o resultado na tela
file_text <- "tabelacontingencia2x2_capacetetrauma.txt"
# file_text <- NA
# Usar apenas espaco em branco para formatar tabela de contigencia
fator_exposicao <- "Uso do capacete"
Tabela <- ("
Capacete Trauma Nao_trauma
Sim 17 138
Nao 130 508
")
# Nao edite daqui em diante
TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE, row.names=1))
if (!is.na(file_text))
{
sink(file_text)
}
cat("\n")
cat("------------------------------\n")
cat("Testes qui-quadrado de Pearson\n")
cat("------------------------------\n")
cat("\n")
cat("Fator de exposicao: ",fator_exposicao,"\n")
print(TC)
# GARCIA-PEREZ, MA & NUNEZ-ANTON, V (2003)
# Cellwise residual analysis in two-way contingency tables.
# Educational and Psychological Measurement 63(5): 825-839
cat("\nX^2 critico de 95% =",qchisq(p=0.95,df=1),"\n")
cat("\n----------\nTeste qui-quadrado assintotico sem correcao de Yates\n")
test <- chisq.test(TC,correct=FALSE)
print (test)
cat("\n----------\nTeste qui-quadrado assintotico com correcao de Yates\n")
test <- chisq.test(TC,correct=TRUE)
print (test)
cat("\n----------\nTeste qui-quadrado exato (robusto)\n")
res <- chisq.test(TC,correct=FALSE,simulate.p.value = TRUE, B=1e5)
print(res)
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
show.MCSTAR(MCSTAR=MCSTAR,alpha=alfa,df=df)
X2 <- res$statistic # estatistica de teste qui-quadrado
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
cat("V de Cramer =",V,"\n")
ncp <- X2 # parametro de nao-centralidade
poder <- 1-pchisq(qchisq(0.95,df),df,ncp)
cat("Poder observado =",poder,"\n")
cat("\nTeste de razao de chances (OR) robusto\n")
resft <- fisher.test(TC) # Teste de OR robusto
print(resft)
OR <- resft$estimate
cat("Razao de chances (OR) =",OR,"\n")
preval <- 0.05
RR <- DescTools::ORToRelRisk(OR,preval)
cat("Razao de riscos (RR) a partir de OR\n(",
preval, "de prevalencia do desfecho na populacao nao-exposta)\n",
"RR =",RR,"\n")
# Personality Project - Revelle – 2014
# Covariance, Regression, and Correlation in R - Chapter 4
d <- abs(log(OR)/1.81)
cat("abs(d de Cohen) =",d,"\n")
r <- d/sqrt(d^2+4)
cat("abs(correl) =",r,"\n\n")
# EDWARDES, MD & BALTZAN, M (2000) The generalization of the odds ratio,
# risk ratio and risk difference to rxk tables.
# Statistics in Medicine 19:1901-14.
gama <- DescTools::GoodmanKruskalGamma(TC)
ORg <- (1+gama)/(1-gama)
cat("Gama de Goodman-Kruskal =", gama,"\n")
cat("OR generalizado = (1+",gama,")/(1-",gama,") =", ORg,"\n")
lnORg <- log(ORg)
cat("ln(OR generalizado) =",lnORg,"\n")
RRg <- DescTools::ORToRelRisk(ORg,preval)
cat("Razao de riscos generalizada (RRg) a partir de ORg\n",
"(ZHANG and YU,1998)\n",
"(",preval, "de prevalencia do desfecho na populacao nao-exposta)\n",
"RRg =",RRg,"\n")
dg <- abs(log(ORg)/1.81)
cat("abs(d de Cohen) =",dg,"\n")
rg <- dg/sqrt(dg^2+4)
cat("abs(correl) =",rg,"\n")
# fecha o arquivo, se abriu antes
if (!is.na(file_text))
{
sink()
cat ("\nResultados armazenados em ",file_text,"\n", sep="")
}e observe sua saída em tabelacontingencia2x2_capacetetrauma.txt.
Teste qui-quadrado interativo por bootstrapping
Em Qui2.R também
podemos processar o mesmo exemplo ou qualquer outro. Veja o código para
aprender a recolher respostas do usuário. Esta versão aplica
bootstrapping de um jeito mais manual e produz gráficos
associados aos conceitos que discutimos.
Execute com:
# Simulacao: Qui-quadrado
# rodando este script np RStudio, observe em Plots
# suppress warnings
options(warn=-1)
cat ("**** TABELA DE CONTINGÊNCIA ****\n")
cat ("\nConsidere:\n")
output <- data.frame(
c1 = c("", "Linha +", "Linha -", "" ),
c2 = c("Coluna+", "a", "c", "Total col. +"),
c3 = c("Coluna-", "b", "d", "Total col. -"),
c4 = c("", "Total lin. +", "Total lin. -", "Total geral")
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n\nInforme:")
nomelin <- readline(prompt="Nome nas linhas (ate 10 letras): ")
nomecol <- readline(prompt="Nome nas colunas (ate 10 letras): ")
cat ("Definido:\n")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), "a", "c", paste("Total de",nomelin,"+")),
c3 = c(paste(nomecol,"-"), "b", "d", paste("Total de",nomelin,"-")),
c4 = c("", paste("Total de",nomecol,"+"), paste("Total de",nomecol,"-"), "Total geral")
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n")
cat ("Forneca os valores de a, b, c, d:\n")
a <- readline(prompt="a:")
a <- as.numeric(a)
if (is.na(a)) {a <- 0}
b <- readline(prompt="b:")
b <- as.numeric(b)
if (is.na(b)) {b <- 0}
c <- readline(prompt="c:")
c <- as.numeric(c)
if (is.na(c)) {c <- 0}
d <- readline(prompt="d:")
d <- as.numeric(d)
if (is.na(d)) {d <- 0}
# caso tenha fornecido em proporcoes
total_abcd = a+b+c+d
if (total_abcd == 1)
{
cat ("\n")
cat ("Qual o tamanho de sua amostra ( = Total geral)?\n")
cat ("(entre um numero inteiro; default n=10000\n")
total_abcd <- readline(prompt="n: ")
total_abcd <- as.integer(total_abcd)
if (is.na(total_abcd)) {total_abcd <- 10000}
a <- a*total_abcd
b <- b*total_abcd
c <- c*total_abcd
d <- d*total_abcd
}
cat ("\n")
cat ("Defina alfa = prob. erro do tipo I).\n")
cat ("(número entre 0 e 1, default = 0.05).\n")
alfa <-readline(prompt="alfa: ")
alfa <- as.numeric(alfa)
if (is.na(alfa)) {alfa <- 0.05; cat(alfa);}
cat ("\n")
cat ("Quantas iteracoes para simular?\n")
cat ("(entre um numero inteiro; default n=3e3)\n")
num_iteracoes <- readline(prompt="iteracoes: ")
num_iteracoes <- as.integer(num_iteracoes)
if (is.na(num_iteracoes)) {num_iteracoes <- 3e3; cat(num_iteracoes);}
cat ("\n")
cat ("Exibir tabelas? (exibir lentifica a simulacao)\n")
exibetabela <- readline(prompt="0=nao, 1=sim; default eh 0: ")
exibetabela <- as.integer(exibetabela)
if (is.na(exibetabela)) {exibetabela <- 0; cat(exibetabela);}
# H1
a_obs <- a/total_abcd
b_obs <- b/total_abcd
c_obs <- c/total_abcd
d_obs <- d/total_abcd
# H0
a_esp <- (a_obs+b_obs)*(a_obs+c_obs)
b_esp <- (a_obs+b_obs)*(b_obs+d_obs)
c_esp <- (a_obs+c_obs)*(c_obs+d_obs)
d_esp <- (c_obs+d_obs)*(b_obs+d_obs)
cat ("\n\nIniciando a simulacao\n")
cat ("\n\nDefinidos:\n")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(a,3), round(c,3), round((a+c),3)),
c3 = c(paste(nomecol,"-"), round(b,3), round(d,3), round((b+d),3)),
c4 = c("", round(a+b,3), round(c+d,3), round(total_abcd,3))
)
names(output) <- NULL
print(output, row.names = FALSE)
# Exibe H0 e H1
cat ("\n")
cat ("Precisamos estabelecer se",nomecol,"depende de",nomelin,".\n")
cat("\n")
cat ("\nH0 (em proporcoes, valores esperados):\n")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(a_esp,3), round(c_esp,3), round((a_esp+c_esp),3)),
c3 = c(paste(nomecol,"-"), round(b_esp,3), round(d_esp,3), round((b_esp+d_esp),3)),
c4 = c("", round(a_esp+b_esp,3), round(c_esp+d_esp,3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat("\n")
cat ("\nH1 (em proporcoes, valores observados):\n")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(a_obs,3), round(c_obs,3), round((a_obs+c_obs),3)),
c3 = c(paste(nomecol,"-"), round(b_obs,3), round(d_obs,3), round((b_obs+d_obs),3)),
c4 = c("", round(a_obs+b_obs,3), round(c_obs+d_obs,3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat("\n")
cat("\namostra=",total_abcd)
cat ("\nalfa =",alfa)
cat ("\niteracoes =",num_iteracoes)
# com que frequencia exibe o grafico (20 graficos)
exibegrf <- round(num_iteracoes/10);
exibepto <- round(num_iteracoes/100);
# valores de referencia
h0_ref <- c(a_esp,b_esp,c_esp,d_esp)
h0_probs <- c()
acm <- 0
for (idx in 1:4)
{
acm <- acm+h0_ref[idx]
h0_probs <- c(h0_probs, acm)
}
h1_ref <- c(a_obs,b_obs,c_obs,d_obs)
h1_probs <- c()
acm <- 0
for (idx in 1:4)
{
acm <- acm+h1_ref[idx]
h1_probs <- c(h1_probs, acm)
}
# iterando
h0_qui <- c()
h0_p <- c()
h0_OR <- c()
h1_qui <- c()
h1_p <- c()
h1_OR <- c()
h1_VCramer <- c()
for (i in 1:num_iteracoes)
{
# distribuindo os individuos sob H0
aloca_h0 <- c(0,0,0,0)
for (a in 1:total_abcd)
{
rnd=runif(1,min=0,max=1)
for (idx in 1:4)
{
if (rnd <= h0_probs[idx])
{
aloca_h0[idx] <- aloca_h0[idx]+1
break
}
}
}
# converte em % ... aloca_h0 na ordem a,b,c,d
for (idx in 1:4)
{
aloca_h0[idx] <- aloca_h0[idx]/total_abcd
}
h0_sml_obs <- c(0,0,0,0)
h0_sml_obs[1] <- aloca_h0[1]
h0_sml_obs[2] <- aloca_h0[2]
h0_sml_obs[3] <- aloca_h0[3]
h0_sml_obs[4] <- aloca_h0[4]
h0_sml_esp <- c(0,0,0,0)
h0_sml_esp[1] <- (aloca_h0[1]+aloca_h0[2])*(aloca_h0[1]+aloca_h0[3])
h0_sml_esp[2] <- (aloca_h0[1]+aloca_h0[2])*(aloca_h0[2]+aloca_h0[4])
h0_sml_esp[3] <- (aloca_h0[1]+aloca_h0[3])*(aloca_h0[3]+aloca_h0[4])
h0_sml_esp[4] <- (aloca_h0[3]+aloca_h0[4])*(aloca_h0[2]+aloca_h0[4])
phi_h0 <- 0
for (idx in 1:4)
{
phi_h0 <- phi_h0 + ((h0_sml_obs[idx]-h0_sml_esp[idx])^2)/h0_sml_esp[idx]
}
qui2_h0 <- phi_h0 * total_abcd
h0_qui <- c(h0_qui,qui2_h0)
h0_p <- c(h0_p,1-pchisq(qui2_h0,1))
OR_h0 <- (aloca_h0[1]*aloca_h0[4])/(aloca_h0[2]*aloca_h0[3])
if (is.finite(OR_h0))
{
h0_OR <- c(h0_OR,OR_h0)
}
# distribuindo os individuos sob H1
aloca_h1 <- c(0,0,0,0)
for (a in 1:total_abcd)
{
rnd=runif(1,min=0,max=1)
for (idx in 1:4)
{
if (rnd <= h1_probs[idx])
{
aloca_h1[idx] <- aloca_h1[idx]+1
break
}
}
}
# converte em % ... aloca_h1 na ordem a,b,c,d
for (idx in 1:4)
{
aloca_h1[idx] <- aloca_h1[idx]/total_abcd
}
h1_sml_obs <- c(0,0,0,0)
h1_sml_obs[1] <- aloca_h1[1]
h1_sml_obs[2] <- aloca_h1[2]
h1_sml_obs[3] <- aloca_h1[3]
h1_sml_obs[4] <- aloca_h1[4]
h1_sml_esp <- c(0,0,0,0)
h1_sml_esp[1] <- (aloca_h1[1]+aloca_h1[2])*(aloca_h1[1]+aloca_h1[3])
h1_sml_esp[2] <- (aloca_h1[1]+aloca_h1[2])*(aloca_h1[2]+aloca_h1[4])
h1_sml_esp[3] <- (aloca_h1[1]+aloca_h1[3])*(aloca_h1[3]+aloca_h1[4])
h1_sml_esp[4] <- (aloca_h1[3]+aloca_h1[4])*(aloca_h1[2]+aloca_h1[4])
phi_h1 <- 0
for (idx in 1:4)
{
phi_h1 <- phi_h1 + ((h1_sml_obs[idx]-h1_sml_esp[idx])^2)/h1_sml_esp[idx]
}
qui2_h1 <- phi_h1 * total_abcd
h1_qui <- c(h1_qui,qui2_h1)
h1_p <- c(h1_p,1-pchisq(qui2_h1,1))
OR_h1 <- (aloca_h1[1]*aloca_h1[4])/(aloca_h1[2]*aloca_h1[3])
if (is.finite(OR_h1))
{
h1_OR <- c(h1_OR,OR_h1)
}
# V de Cramer
nL <- 2 #nrow(TC) # ou dim(TC)[1]
nC <- 2 #ncol(TC) # ou dim(TC)[2]
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V_h1 <- sqrt(phi_h1/Dim) # V ou C de Cramer
h1_VCramer <- c(h1_VCramer,V_h1)
if (exibetabela == 1)
{
cat("\n")
cat ("\n----------------------------------------------\n")
cat ("\nValores simulados sob H0:\n")
cat ("\n- Individuos alocados:")
cat ("\n--- observados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(aloca_h0[1]*total_abcd,0), round(aloca_h0[3]*total_abcd,0), round((aloca_h0[1]+aloca_h0[3])*total_abcd,0)),
c3 = c(paste(nomecol,"-"), round(aloca_h0[2]*total_abcd,0), round(aloca_h0[4]*total_abcd,0), round((aloca_h0[2]+aloca_h0[4])*total_abcd,0)),
c4 = c("", round((aloca_h0[1]+aloca_h0[2])*total_abcd,0), round((aloca_h0[3]+aloca_h0[4])*total_abcd,0), total_abcd )
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n--- esperados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h0_sml_esp[1]*total_abcd,3), round(h0_sml_esp[3]*total_abcd,3), round((h0_sml_esp[1]+h0_sml_esp[3])*total_abcd,3)),
c3 = c(paste(nomecol,"-"), round(h0_sml_esp[2]*total_abcd,3), round(h0_sml_esp[4]*total_abcd,3), round((h0_sml_esp[2]+h0_sml_esp[4])*total_abcd,3)),
c4 = c("", round((h0_sml_esp[1]+h0_sml_esp[2])*total_abcd,3), round((h0_sml_esp[3]+h0_sml_esp[4])*total_abcd,3), total_abcd)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n- Em probabilidades:")
cat ("\n--- observados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h0_sml_obs[1],3), round(h0_sml_obs[3],3), round((h0_sml_obs[1]+h0_sml_obs[3]),3)),
c3 = c(paste(nomecol,"-"), round(h0_sml_obs[2],3), round(h0_sml_obs[4],3), round((h0_sml_obs[2]+h0_sml_obs[4]),3)),
c4 = c("", round(h0_sml_obs[1]+h0_sml_obs[2],3), round(h0_sml_obs[3]+h0_sml_obs[4],3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n--- esperados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h0_sml_esp[1],3), round(h0_sml_esp[3],3), round((h0_sml_esp[1]+h0_sml_esp[3]),3)),
c3 = c(paste(nomecol,"-"), round(h0_sml_esp[2],3), round(h0_sml_esp[4],3), round((h0_sml_esp[2]+h0_sml_esp[4]),3)),
c4 = c("", round(h0_sml_esp[1]+h0_sml_esp[2],3), round(h0_sml_esp[3]+h0_sml_esp[4],3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\nPhi^2 de Cohen:",phi_h0)
cat ("\nQui-quadrado:",qui2_h0,", p:",round(1-pchisq(qui2_h0,1),5))
cat ("\nOdds Ratio:",OR_h0,"\n")
cat("\n")
cat ("\nValores simulados sob H1:\n")
cat ("\n- Individuos alocados:")
cat ("\n--- observados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(aloca_h1[1]*total_abcd,0), round(aloca_h1[3]*total_abcd,0), round((aloca_h1[1]+aloca_h1[3])*total_abcd,0)),
c3 = c(paste(nomecol,"-"), round(aloca_h1[2]*total_abcd,0), round(aloca_h1[4]*total_abcd,0), round((aloca_h1[2]+aloca_h1[4])*total_abcd,0)),
c4 = c("", round((aloca_h1[1]+aloca_h1[2])*total_abcd,0), round((aloca_h1[3]+aloca_h1[4])*total_abcd,0), total_abcd )
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n--- esperados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h1_sml_esp[1]*total_abcd,3), round(h1_sml_esp[3]*total_abcd,3), round((h1_sml_esp[1]+h1_sml_esp[3])*total_abcd,3)),
c3 = c(paste(nomecol,"-"), round(h1_sml_esp[2]*total_abcd,3), round(h1_sml_esp[4]*total_abcd,3), round((h1_sml_esp[2]+h1_sml_esp[4])*total_abcd,3)),
c4 = c("", round((h1_sml_esp[1]+h1_sml_esp[2])*total_abcd,3), round((h1_sml_esp[3]+h1_sml_esp[4])*total_abcd,3), total_abcd)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n- Em probabilidades:")
cat ("\n--- observados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h1_sml_obs[1],3), round(h1_sml_obs[3],3), round((h1_sml_obs[1]+h1_sml_obs[3]),3)),
c3 = c(paste(nomecol,"-"), round(h1_sml_obs[2],3), round(h1_sml_obs[4],3), round((h1_sml_obs[2]+h1_sml_obs[4]),3)),
c4 = c("", round(h1_sml_obs[1]+h1_sml_obs[2],3), round(h1_sml_obs[3]+h1_sml_obs[4],3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\n--- esperados:")
output <- data.frame(
c1 = c("", paste(nomelin,"+"), paste(nomelin,"-"), "" ),
c2 = c(paste(nomecol,"+"), round(h1_sml_esp[1],3), round(h1_sml_esp[3],3), round((h1_sml_esp[1]+h1_sml_esp[3]),3)),
c3 = c(paste(nomecol,"-"), round(h1_sml_esp[2],3), round(h1_sml_esp[4],3), round((h1_sml_esp[2]+h1_sml_esp[4]),3)),
c4 = c("", round(h1_sml_esp[1]+h1_sml_esp[2],3), round(h1_sml_esp[3]+h1_sml_esp[4],3), 1)
)
names(output) <- NULL
print(output, row.names = FALSE)
cat ("\nPhi^2 de Cohen:",phi_h1)
cat ("\nQui-quadrado:",qui2_h1,", p:",round(1-pchisq(qui2_h1,1),5))
cat ("\nOdds Ratio:",OR_h1,"\n")
cat("\n")
cat ("\nV de Cramer:",V_h1,"\n")
}
if (i>50 & (i %% exibepto) == 0)
{
cat(".")
}
if ( length(h0_qui)>=2 & ((i %% exibegrf) == 0 || (i>=10 & i <= 50)) )
{
if(i>50)
{
cat("\n")
}
# densidades
h0_qui <- sort(h0_qui)
h1_qui <- sort(h1_qui)
dens_h0 <- density(h0_qui, na.rm=TRUE)
dens_h1 <- density(h1_qui, na.rm=TRUE)
# xlim, ylim
df <- 1
# quantil 95% (empirico)
q95 <- quantile(h0_qui,probs=1-alfa, na.rm=TRUE)
# media da dist. qui^2
mediah0 <- df
dph0 <- sqrt(2*df)
mediah1 <- mean(h1_qui)
dph1 <- sqrt(2*df+mediah1)
# sugestao para os eixos, evitando cauda muito longa
truncar_x <- 0.05
maxx <- mediah1+3*dph1
maxy <- max(dens_h0$y, dens_h1$y)
cor_H0 = "#1965B0"
cor_alfa_transparencia = paste(cor_H0,"88",sep="")
cor_H1 = "#a30b1b"
cor_beta_transparencia = paste(cor_H1,"55",sep="")
cor_poder_transparencia = "#F7F056aa"
titulo = paste(nomelin," x ",nomecol," (n=",total_abcd,", ",i," iteracoes)", sep="")
# H0
plot (dens_h0$x[dens_h0$x>=0],dens_h0$y[dens_h0$x>=0],
main=titulo,
xlab="Valor da Estatistica Qui-Quadrado",
ylab="Densidade",
col=cor_H0, lwd=3, type="l",
xlim=c(truncar_x,maxx), ylim=c(0,maxy))
# H1
lines (dens_h1, col=cor_H1, lwd=2, lty=2)
# beta
d_beta_x <- c(dens_h1$x[dens_h1$x<=q95])
d_beta_x <- c(min(d_beta_x),d_beta_x,max(d_beta_x))
d_beta_y <- c(dens_h1$y[dens_h1$x<=q95])
d_beta_y <- c(0,d_beta_y,0)
polygon(d_beta_x,d_beta_y,col=cor_beta_transparencia, border=cor_beta_transparencia)
# poder
d_poder_x <- c(dens_h1$x[dens_h1$x>=q95])
d_poder_x <- c(min(d_poder_x),d_poder_x,max(d_poder_x))
d_poder_y <- c(dens_h1$y[dens_h1$x>=q95])
d_poder_y <- c(0,d_poder_y,0)
polygon(d_poder_x,d_poder_y,col=cor_poder_transparencia, border=cor_poder_transparencia)
# alfa
q95 <- quantile(h0_qui,probs=1-alfa, na.rm=TRUE)
d_alfa_x <- c(dens_h0$x[dens_h0$x>=q95]);
d_alfa_x <- c(min(d_alfa_x),d_alfa_x,max(d_alfa_x))
d_alfa_y <- c(dens_h0$y[dens_h0$x>=q95])
d_alfa_y <- c(0,d_alfa_y,0)
polygon(d_alfa_x,d_alfa_y,col=cor_alfa_transparencia, border=cor_alfa_transparencia)
# cutoff
lines(c(q95,q95),c(0,maxy), lty=3)
txtH0 <- paste("H0")
txtH1 <- paste("H1")
legend ("topright",c(txtH0,txtH1,"alfa","beta","poder"),
lwd=c(3,2,15,15,15),
lty=c(1,2,1,1,1),
col=c(cor_H0,cor_H1,cor_alfa_transparencia,cor_beta_transparencia,cor_poder_transparencia),
box.lwd=0,
bg="transparent")
# pausa para ver o grafico
if (i>=10 & i <= 100)
{
Sys.sleep(0.2)
}
}
} # iteracoes
# distribuicao de OR
h1_OR <- sort(h1_OR)
h1_ICOR <- quantile(h1_OR,probs=c(alfa/2,1-(alfa/2)))
d_OR <- density(h1_OR)
medianaOR <- median(h1_OR)
cor_OR = "#ac4d12"
cor_mediana = "#000000"
cor_IC = "#0d5092"
maxx = max(d_OR$x,1)
maxy <- max(d_OR$y)
plot (d_OR, main="Distribuicao do Odds-ratio",
xlab = "OR", ylab = "Densidade",
xlim = c(0,maxx), ylim = c(0,maxy),
col=cor_OR, lwd=3, lty=1
)
legend ("topright",
c("OR simul.",
"mediana",
paste("IC",round((1-alfa)*100,1),"%",sep=""),
"referencia"),
pch=c(NA,NA,NA,19),
lwd=c(3,1,2,NA),
lty=c(1,2,1,NA),
col=c(cor_OR,cor_mediana,cor_IC,"#000000"),
box.lwd=0,
bg="transparent")
# faixas de OR
# protecao
# 0.1) = grande
# [0.1 a 0.3) = intermediaria
# [0.3 a 0.7) = pequena
# [0.7 a 1) = muito pequena
# neutra
# [1] =
# risco
# (1 a 1.5] = muito pequeno
# (1.5 a 3.5] = pequeno
# (3.5 a 9] = intermediario
# (9 = grande
maxx <- max(d_OR$x)
maxy <- max(d_OR$y)
if (maxx<9) {maxx = 9}
faixa_lim = c(
0,
0.1,
0.3,
0.7,
1,
1.5,
3.5,
9,
maxx
)
faixa_txt = c(
"ef-.gde", # verde musgo
"ef-.int", # verde folha
"ef-.peq", # verde claro 1
"ef-.min", # verde claro 3
"ef+.min", # tijolo
"ef+.peq", # terra
"ef+.int", # marrom claro
"ef+.gde" # marrom médio
)
# green
# [12] = "#507052"; # verde musgo
# [14] = "#4EB265"; # verde folha
# [15] = "#90C987"; # verde claro 1
# [17] = "#CAE0AB"; # verde claro 3
# red
# [29] = "#f43328"; # tijolo
# [28] = "#a3261b"; # terra
# [26] = "#A5170E"; # marrom claro
# [25] = "#72190E"; # marrom médio
faixa_cor = c(
"#50705288", # verde musgo
"#4EB26588", # verde folha
"#90C98788", # verde claro 1
"#CAE0AB88", # verde claro 3
"#f4332888", # tijolo
"#a3261b88", # terra
"#A5170E88", # marrom claro
"#72190E88" # marrom médio
)
for (f in 1:(length(faixa_lim)-1))
{
polygon(c(faixa_lim[f],faixa_lim[f],faixa_lim[f+1],faixa_lim[f+1]),c(0,maxy/10,maxy/10,0),col=faixa_cor[f], border=faixa_cor[f])
}
# legenda das faixas
leg_txt <- c()
leg_col <- c()
for (f in 1:length(faixa_txt))
{
if (faixa_lim[f] <= maxx)
{
leg_txt <- c(leg_txt,faixa_txt[f])
leg_col <- c(leg_col,faixa_cor[f])
}
}
maxx <- max(d_OR$x)
if (medianaOR <= 3.5*maxx/10 | medianaOR >= 6.5*maxx/10)
{
# fora do centro
posleg <- "top"
} else # mediana no centro
{
if (medianaOR > 5*maxx/10) {posleg <- "topleft"}
if (medianaOR < 5*maxx/10) {posleg <- "right"}
}
legend (posleg,
leg_txt,
lwd=10,
lty=1,
col=leg_col,
box.lwd=0,
bg="transparent")
# IC, mediana e referencia
lines(c(h1_ICOR[1],h1_ICOR[1]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
lines(c(h1_ICOR[2],h1_ICOR[2]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
lines(c(h1_ICOR[1],h1_ICOR[2]),c(maxy/10/2,maxy/10/2),lty=1,lwd=2,col=cor_IC)
lines(c(medianaOR,medianaOR),c(0,maxy),lty=2,lwd=1,col=cor_mediana)
points(c(1),c(maxy/10/2), pch=19,cex=1,col="#000000")
# # distribuicao de log(OR)
# h1_OR <- sort(h1_OR)
# h1_ICOR <- quantile(h1_OR,probs=c(alfa/2,1-(alfa/2)))
# d_OR <- density(h1_OR)
# medianaOR <- median(h1_OR)
#
# cor_OR = "#ac4d12"
# cor_mediana = "#000000"
# cor_IC = "#0d5092"
#
# minx = min(d_OR$x[d_OR$x!=0],1e-4)
# maxx = max(d_OR$x,1)
# maxy <- max(d_OR$y)
# plot (d_OR, main="Distribuicao do ln(Odds-ratio)",
# xlab = "ln(OR)", ylab = "Densidade",
# xlim=c(minx,maxx), ylim = c(0,maxy), log = "x",
# col=cor_OR, lwd=3, lty=1
# )
# legend ("topleft",
# c("ln(OR) simul.",
# "mediana",
# paste("IC",round((1-alfa)*100,1),"%",sep="")
# ),
# pch=c(NA,NA,NA),
# lwd=c(3,1,2),
# lty=c(1,2,1),
# col=c(cor_OR,cor_mediana,cor_IC),
# box.lwd=0,
# bg="transparent")
# # faixas de OR
# # protecao
# # 0.1) = grande
# # [0.1 a 0.3) = intermediaria
# # [0.3 a 0.7) = pequena
# # [0.7 a 1) = muito pequena
# # neutra
# # [1] =
# # risco
# # (1 a 1.5] = muito pequeno
# # (1.5 a 3.5] = pequeno
# # (3.5 a 9] = intermediario
# # (9 = grande
# maxx <- max(d_OR$x)
# minx <- max(d_OR$x)
# if (minx >= 0.1) {minx = 0.001}
# maxy <- max(d_OR$y)
# if (maxx<9) {maxx = 9}
# faixa_lim = c(
# minx,
# 0.1,
# 0.3,
# 0.7,
# 1,
# 1.5,
# 3.5,
# 9,
# maxx
# )
# faixa_txt = c(
# "ef-.gde", # verde musgo
# "ef-.int", # verde folha
# "ef-.peq", # verde claro 1
# "ef-.min", # verde claro 3
# "ef+.min", # tijolo
# "ef+.peq", # terra
# "ef+.int", # marrom claro
# "ef+.gde" # marrom médio
# )
# # green
# # [12] = "#507052"; # verde musgo
# # [14] = "#4EB265"; # verde folha
# # [15] = "#90C987"; # verde claro 1
# # [17] = "#CAE0AB"; # verde claro 3
# # red
# # [29] = "#f43328"; # tijolo
# # [28] = "#a3261b"; # terra
# # [26] = "#A5170E"; # marrom claro
# # [25] = "#72190E"; # marrom médio
# faixa_cor = c(
# "#50705288", # verde musgo
# "#4EB26588", # verde folha
# "#90C98788", # verde claro 1
# "#CAE0AB88", # verde claro 3
# "#f4332888", # tijolo
# "#a3261b88", # terra
# "#A5170E88", # marrom claro
# "#72190E88" # marrom médio
# )
# for (f in 1:(length(faixa_lim)-1))
# {
# polygon(c(faixa_lim[f],faixa_lim[f],faixa_lim[f+1],faixa_lim[f+1]),c(0,maxy/10,maxy/10,0),col=faixa_cor[f], border=faixa_cor[f])
# }
# # legenda das faixas
# leg_txt <- c()
# leg_col <- c()
# for (f in 1:length(faixa_txt))
# {
# if (faixa_lim[f] <= maxx)
# {
# leg_txt <- c(leg_txt,faixa_txt[f])
# leg_col <- c(leg_col,faixa_cor[f])
# }
# }
# legend ("left",
# leg_txt,
# lwd=10,
# lty=1,
# col=leg_col,
# box.lwd=0,
# bg="transparent")
# # IC, mediana e referencia
# lines(c(h1_ICOR[1],h1_ICOR[1]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
# lines(c(h1_ICOR[2],h1_ICOR[2]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
# lines(c(h1_ICOR[1],h1_ICOR[2]),c(maxy/10/2,maxy/10/2),lty=1,lwd=2,col=cor_IC)
# lines(c(medianaOR,medianaOR),c(0,maxy),lty=2,lwd=1,col=cor_mediana)
# # points(c(0),c(maxy/10/2), pch=19,cex=1,col="#000000")
# points(minx,maxy/10/2,pch=19,col="black",bg="black")
# distribuicao de V de Cramer
h1_VCramer <- sort(h1_VCramer)
h1_ICV <- quantile(h1_VCramer,probs=c(alfa/2,1-(alfa/2)))
d_V <- density(h1_VCramer)
medianaV <- median(h1_VCramer)
cor_V = "#994F88"
cor_mediana = "#000000"
cor_IC = "#0d5092"
plot (d_V, main="Distribuicao do V de Cramer",
xlab = "V de Cramer", ylab = "Densidade",
col=cor_V, lwd=3, lty=1
)
maxy <- max(d_V$y)
legend ("topright",
c("V simul.",
"mediana",
paste("IC",round((1-alfa)*100,1),"%",sep="")
),
pch=c(NA,NA,NA),
lwd=c(3,1,2),
lty=c(1,2,1),
col=c(cor_V,cor_mediana,cor_IC),
box.lwd=0,
bg="transparent")
# faixas de Cramer
# protecao
# 0.1) = minima
# [0.1 a 0.3) = pequena
# [0.3 a 0.5) = intermediaria
# [0.5 = grande
maxx <- max(d_V$x)
maxy <- max(d_V$y)
if (maxx<0.5) {maxx = 0.6}
faixa_lim = c(
0, # verde musgo
0.1, # verde folha
0.3, # verde claro 1
0.5, # verde claro 3
maxx
)
faixa_txt = c(
"corr.min", # verde musgo
"corr.peq", # verde folha
"corr.int", # verde claro 1
"corr.gde" # verde claro 3
)
# green
# [12] = "#507052"; # verde musgo
# [14] = "#4EB265"; # verde folha
# [15] = "#90C987"; # verde claro 1
# [17] = "#CAE0AB"; # verde claro 3
faixa_cor = c(
"#CAE0AB88", # verde claro 3
"#90C98788", # verde claro 1
"#4EB26588", # verde folha
"#50705288" # verde musgo
)
for (f in 1:(length(faixa_lim)-1))
{
polygon(c(faixa_lim[f],faixa_lim[f],faixa_lim[f+1],faixa_lim[f+1]),c(0,maxy/10,maxy/10,0),col=faixa_cor[f], border=faixa_cor[f])
}
# legenda das faixas
maxx <- max(d_V$x)
leg_txt <- c()
leg_col <- c()
for (f in 1:length(faixa_txt))
{
if (faixa_lim[f] <= maxx)
{
leg_txt <- c(leg_txt,faixa_txt[f])
leg_col <- c(leg_col,faixa_cor[f])
}
}
maxx <- max(d_V$x)
if (medianaV <= 3.5*maxx/10 | medianaV >= 6.5*maxx/10)
{
# fora do centro
posleg <- "top"
} else # mediana no centro
{
if (medianaV > 5*maxx/10) {posleg <- "topleft"}
if (medianaV < 5*maxx/10) {posleg <- "right"}
}
legend (posleg,
leg_txt,
lwd=10,
lty=1,
col=leg_col,
box.lwd=0,
bg="transparent")
# IC e mediana
lines(c(h1_ICV[1],h1_ICV[1]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
lines(c(h1_ICV[2],h1_ICV[2]),c(0,maxy/10),lty=1,lwd=2,col=cor_IC)
lines(c(h1_ICV[1],h1_ICV[2]),c(maxy/10/2,maxy/10/2),lty=1,lwd=2,col=cor_IC)
maxy <- max(d_V$y)
lines(c(medianaV,medianaV),c(0,maxy),lty=2,lwd=1,col=cor_mediana)
# # distribuicao de p
h0_p <- sort(h0_p)
h0_ICp <- quantile(h0_p,probs=c(alfa/2,1-(alfa/2)))
d_p_h0 <- density(h0_p)
medianap_h0 <- median(h0_p)
d_p_h1 <- density(h1_p)
pexato_h1 <- median(h1_p, na.rm = TRUE)
#
# cor_h1 = "#EE8026"
# cor_h0 = "#507052"
# cor_mediana = "#000000"
# cor_IC = "#0d5092"
#
# maxx = max(d_p_h0$x,d_p_h1$x)
# maxy <- max(d_p_h0$y,d_p_h1$y)
# plot (d_p_h1, main="Distribuicao do valor-p",
# xlab = "p", ylab = "Densidade",
# xlim = c(0,maxx), ylim = c(0,maxy),
# col=cor_h1, lwd=3, lty=1
# )
# lines (d_p_h0, col=cor_h0, lwd=3, lty=2)
# # IC e mediana
# lines(c(pexato_h1,pexato_h1),c(0,maxy),lty=2,lwd=1,col=cor_mediana)
# legend ("topright",
# c("H0","H1",
# "valor-p exato"
# ),
# pch=c(NA,NA,NA),
# lwd=c(3,3,1),
# lty=c(2,1,2),
# col=c(cor_h0,cor_h1,cor_mediana),
# box.lwd=0,
# bg="transparent")
# Grafico isolando H0, alfa e p
# densidades p sob H0
h0_qui <- sort(h0_qui)
dens_h0 <- density(h0_qui, na.rm=TRUE)
# xlim, ylim
df <- 1
# media da dist. qui^2
mediah0 <- df
dph0 <- sqrt(2*df)
# sugestao para os eixos, evitando cauda muito longa
truncar_x <- 0.05
q95 <- quantile(h0_qui,probs=(1-pexato_h1), na.rm=TRUE)
maxx <- mediah0+3*dph0
if (maxx < as.vector(q95)) {maxx = as.vector(q95)}
maxx <- maxx+1
maxy <- max(dens_h0$y, dens_h1$y)
# quantil 95% (empirico)
q95 <- quantile(h0_qui,probs=1-alfa, na.rm=TRUE)
cor_H0 = "#1965B0"
cor_alfa_transparencia = paste(cor_H0,"88",sep="")
cor_p = "#a30b1b"
titulo = paste(nomelin," x ",nomecol," (n=",total_abcd,", ",i," iteracoes)", sep="")
# H0
plot (dens_h0$x[dens_h0$x>=0],dens_h0$y[dens_h0$x>=0],
main=titulo,
xlab="Valor da Estatistica Qui-Quadrado",
ylab="Densidade",
col=cor_H0, lwd=3, type="l",
xlim=c(truncar_x,maxx), ylim=c(0,maxy))
# alfa
q95 <- quantile(h0_qui,probs=1-alfa, na.rm=TRUE)
d_alfa_x <- c(dens_h0$x[dens_h0$x>=q95]);
d_alfa_x <- c(min(d_alfa_x),d_alfa_x,max(d_alfa_x))
d_alfa_y <- c(dens_h0$y[dens_h0$x>=q95])
d_alfa_y <- c(0,d_alfa_y,0)
polygon(d_alfa_x,d_alfa_y,col=cor_alfa_transparencia, border=cor_alfa_transparencia)
# p
q_h1 <- quantile(h0_qui,probs=(1-pexato_h1), na.rm=TRUE)
d_p_x <- c(dens_h0$x[dens_h0$x>=q_h1]);
d_p_x <- c(min(d_p_x),d_p_x,max(d_p_x))
d_p_y <- c(dens_h0$y[dens_h0$x>=q_h1])
d_p_y <- c(0,d_p_y,0)
polygon(d_p_x,d_p_y,col=cor_p, border=cor_p)
lines(c(q_h1,q_h1),c(0,maxy/10), lwd=3, col=cor_p, lty=4)
# cutoff alfa
lines(c(q95,q95),c(0,maxy), lty=3)
txtH0 <- paste("H0")
legend ("topright",c(txtH0,"alfa","p"),
lwd=c(3,15,15),
lty=c(1,1,1),
col=c(cor_H0,cor_alfa_transparencia,cor_p),
box.lwd=0,
bg="transparent")
# calculos
q95 <- quantile(h0_qui,probs=1-alfa, na.rm=TRUE)
cat("\n\nTerminado:\n")
cat("alfa: ",round(alfa*100,3),"%\n", sep="")
cat("\nvalor-p exato:", pexato_h1,"\n")
cat("\nOR:\n")
cat("exato:", medianaOR,"\n")
cat("IC",round((1-alfa)*100,1),"%: [",round(h1_ICOR[1],3),", ",round(h1_ICOR[2],3),"]\n", sep="")
cat("\nV de Cramer:\n")
cat("exato:", medianaV,"\n")
cat("IC",round((1-alfa)*100,1),"%: [",round(h1_ICV[1],3),", ",round(h1_ICV[2],3),"]\n", sep="")
# enable warnings
options(warn=0)fornecendo os dados à medida que são pedidos, como por exemplo:
**** TABELA DE CONTINGÊNCIA ****
Considere:
Coluna+ Coluna-
Linha+ a b Total lin.+
Linha- c d Total lin.-
Total col. + Total col. - Total geral
Informe:
Nome nas linhas (ate 10 letras): Capacete
Nome nas colunas (ate 10 letras): Trauma
Definido:
Trauma+ Trauma-
Capacete+ a b Total de Trauma+
Capacete- c d Total de Trauma-
Total de Capacete + Total de Capacete - Total geral
Forneca os valores de a, b, c, d:
a: 17
b: 138
c: 130
d: 508
Defina alfa = prob. erro do tipo I).
(número entre 0 e 1, default = 0.05).
alfa: 0.05
Quantas iteracoes para simular?
(entre um numero inteiro; default n=1e4)
iteracoes: 1e4
Exibir tabelas? (exibir lentifica a simulacao)
0=nao, 1=sim; default eh 0: 0
Observe a simulação em andamento na área de Plots. Quando terminar, procure os gráficos finalizados. Os principais são:
- as distribuições \(\chi^2\) para a hipótese nula e para o valor observado de \(H_1\):
| O valor de \(\beta\) e o poder associados, aqui, são os valores a posteriori ou retrospectivos, portanto inúteis para qualquer decisão. |
|
Flashback: Recorde a tomada de decisão com a distribuição binomial:
A estatística e o formato das distribuições mudam (em função, também, do tipo de variável), mas as áreas de \(\alpha\), \(\beta\) e o raciocínio são os mesmos. Lembre-se, porém, que \(\beta\) a posteriori não tem valor para a decisão estatística. |
- a distribuição, com intervalo de confiança, do odds ratio simulado (compare com o obtido acima), indicando a rejeição de \(H_0\) (o valor 1, de referência, está fora do intervalo) e que seu efeito está entre mínimo a intermediário:
|
Odds ratio é uma medida da intensidade de associação entre duas variáveis nominais dicotômicas de exposição e desfecho. De acordo com Sharpe (2015), OR é, também, uma medida de tamanho de efeito [sic] (não depende do tamanho da amostra, é adimensional e adirecional, mas não tem limite superior (o limite inferior é zero)). |
- a distribuição do tamanho de efeito (V de Cramer), mostrando que o efeito, apesar de significante, está entre mínimo e pequeno:
|
O programa já indica o tamanho de efeito de acordo com
V de Cramer é a correlação de Pearson absoluta implícita entre duas variáveis nominais:
|
- novamente a distribuição \(\chi^2\), indicando as áreas de \(\alpha\) e do valor p, concluindo-se pela rejeição de \(H_0\):
Teste qui-quadrado de Pearson em tabela \(L\times C\)
O teste qui-quadrado opera em tabelas maiores, com \(L\) linhas e \(C\) colunas.
Por exemplo, três quimioterápicos foram comparados quanto ao efeito colateral (náusea). Considere \(\alpha=0.05\). Há associação entre náusea e drogas?
Os dados obtidos foram:
| Droga | Náusea | Sem náusea |
|---|---|---|
| A | 3 | 5 |
| B | 7 | 2 |
| C | 6 | 3 |
Em Qui2_LxC.R
implementamos este exemplo (versão robusta):
alfa <- 0.05
B <- 1e5
tbl <- ("
Droga ComNausea SemNausea
A 3 5
B 7 2
C 6 3
")
tcrxc <- as.matrix(read.table(textConnection(tbl),header=TRUE,row.names=1))
print(tcrxc)
res <- chisq.test(tcrxc,simulate.p.value=TRUE,B=B)
print(res)
x2 <- res$statistic # estatistica de teste qui-quadrado
n <- sum(res$observed)
r <- nrow(tcrxc)
c <- ncol(tcrxc)
df <- (r-1)*(c-1) # graus de liberdade de tabela de contingencia rxc
qui_critico <- qchisq(1-alfa,df)
phi2 <- x2/n # phi ou w de Cohen
dim <- min(r,c)-1 # dimensao de tcrxc: dim = max(phi^2)
V <- sqrt(phi2/dim)
cat("\nV de Cramer = ", round(V,4), "\n", sep="")
print(effectsize::interpret_cramers_v(V))
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(tcrxc),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
# distribuicao qui quadrado sob H0
qui <- rchisq(1e6, df)
dqui <- density(qui)
# distribuicao qui quadrado sob H1
quih1 <- rchisq(1e6, df, x2)
dquih1 <- density(quih1)
# sugestao para os eixos, evitando cauda muito longa
# media da dist. qui^2
mediah0 <- df
dph0 <- sqrt(2*df)
mediah1 <- df+x2
dph1 <- sqrt(2*(df+2*x2))
maxx <- mediah1+3*dph1
maxy <- max(c(dqui$y,dquih1$y))
maxyXcal <- max(dquih1$y)
plot(dqui$x[dqui$x>0],dqui$y[dqui$x>0],xlim=c(0,maxx),
ylim=c(0,maxy), xlab="qui-quadrado", ylab="densidade", type="l")
lines(c(qui_critico,qui_critico),c(0,maxy), lty=2)
lines(c(x2,x2),c(0,maxyXcal), lwd=2, lty=2)
legend ("topright",
c("H0",
paste("X^2 critico=",round(qui_critico,4),sep=""),
paste("X^2=",round(x2,4),sep="")),
pch=c(NA,NA,NA),
lwd=c(1,1,2),
lty=c(1,2,2),
col=c("#000000","#000000","#000000"),
box.lwd=0,
bg="transparent") ComNausea SemNausea
A 3 5
B 7 2
C 6 3
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: tcrxc
X-squared = 3.0559, df = NA, p-value = 0.275
V de Cramer = 0.3428
X-squared
"large"
(Rules: funder2019)
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.3428
Confidence interval:
2.5% 97.5%
0.0000 0.6095 
Observe que se rejeita \(H_0\) se (os dois critérios são sempre equivalentes se o método analítico é usado):
- \(X^2 > \chi^2_c(1-\alpha)\)
- \(p > \alpha\)
Estude o código R para ver como foi implementado este teste.
Análise post hoc
Esta análise permite, quando se rejeita \(H_0\), localizar quais células da tabela de contingência mais contribuíram para esta rejeição.
Por exemplo, para descobrir se existe uma relação entre tabagismo e etilismo a partir de um estudo realizado com 100 estudantes universitários, encontrou-se:
| Categoria | Tabagista | Não tabagista |
|---|---|---|
| Etilista | 50 | 15 |
| Não etilista | 20 | 25 |
| A manipulação de dados em R é poderosa, mas nem sempre fácil. Os dados precisam ser ajustados de acordo com a forma que cada função os recebe: dataframes, tabelas e matrizes, por exemplo. Veja o Apêndice 3 para mais detalhes. |
Neste exemplo, o objetivo é executar o teste qui-quadrado mas, no
caso de rejeitar-se \(H_0\), executar
uma outra alternativa para análise post hoc, utilizando
MCSTAR proposto por García-Pérez & Núñez-Antón
(2003).
O código de TabelaContingencia2x2_TabagismoEtilismo.R
utiliza os dados da planilha tabagismo_e_etilismo_2x2.xlsx
e que produz a seguinte saída:
# TabelaContingencia2x2_TabagismoEtilismo.R
library(exact2x2)
source("eiras.show.MCSTAR.R")
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- data.frame(readxl::read_excel("tabagismo_e_etilismo_2x2.xlsx"))
TC <- data.matrix(TCtmp)
rownames(TC) <- as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
gplots::balloonplot(t(as.table(TC)), label = FALSE,
xlab = "", ylab = "",
show.margins = TRUE)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
B <- 1e5
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=B))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
invisible(show.MCSTAR(MCSTAR=MCSTAR,alpha=alfa,df=df))
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
cat("\nV de Cramer = ", round(V,4), "\n", sep="")
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(TC),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
cat("\nTeste de razao de chances (OR) robusto:\n")
resft <- exact2x2::fisher.exact(TC,conf.level = 1-alfa) # Teste de OR robusto
print(resft)New names:
Registered S3 method overwritten by 'gplots': method from reorder.factor
DescTools
• `` -> `...1`
------------------------------------------
Dados:
------------------------------------------
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: TC
X-squared = 12.121, df = NA, p-value = 0.00067
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 12.12149
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standardizados corrigidos por momento (MCSTAR):
|MCSTAR critico| (alphaBonferroni=5%/1) = 1.959964
Tabagista NaoTabagista
Etilista #6.96 -6.96
NaoEtilista -6.96 #6.96
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.332
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.332
Confidence interval:
2.5% 97.5%
0.1388 0.5003
Teste de razao de chances (OR) robusto:
Two-sided Fisher's Exact Test (usual method using minimum likelihood)
data: TC
p-value = 0.0006368
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.7922 9.7333
sample estimates:
odds ratio
4.107342 A novidade está na sessão Resíduos ajustados standardizados corrigidos por momento (MCSTAR). Analise pelos valores positivos acima do valor assinalado como [MCSTAR crítico]: são escores \(z\) relativos a quanto cada célula contribuiu para o qui-quadrado estimado. Neste exemplo, a diagonal principal é a responsável pela rejeição de \(H_0\), indicando associação concordante entre etilismo e tabagismo.
|
A única “sujeira” deste código é:
Isto desaparece se colocarmos algo, que será desprezado, na célula A1
da planilha. Experimente rodar o código com a planilha Caso tenha interesse no código
que são as responsáveis por inibir os possíveis warnings (só os iniba depois que seu código estiver correto; warnings são importantes para evitar erros enquanto desenvolve um script). |
Comentários sobre os métodos tradicionais
Antes das implementações computacionais e a facilidade do uso do R, a partir de uma tabela 2×2 era necessário construir a tabela dos valores esperados sob \(H_0\).
- Valores observados:
| Trauma+ | Trauma- | Total | |
|---|---|---|---|
| Capacete+ | 17 | 138 | 155 |
| Capacete- | 130 | 508 | 638 |
| Total | 147 | 646 | 793 |
- Valores esperados sob:
| Trauma+ | Trauma- | |
|---|---|---|
| Capacete+ | 28.7 | 126.3 |
| Capacete- | 118.3 | 519.7 |
Cada valor esperado é obtido pelo produto das marginais, dividido pelo total geral:
- \((147 \times 155) / 793 \approx 28.7\)
- \((147 \times 638) / 793 \approx 118.3\)
- \((646 \times 155) / 793 \approx 126.3\)
- \((646 \times 638) / 793 \approx 519.7\)
Condições para aplicar o teste qui-quadrado de Pearson
Avaliando os dados observados e os valores esperados sob \(H_0\), era necessário verificar se o teste qui-quadrado podia ser aplicado. Eram as seguintes condições:
Siegel & Castellan Jr., 1988, p. 123:
- O teste exato de Fisher para testar independência é adequado apenas para tabela de contingência 2×2 se \(n < 20\).
- Se \(n\) entre 20 e 40, o teste qui-quadrado de Pearson pode ser usado se todas as frequências esperadas sob \(H_0\) são maiores que 5.
- Se \(n > 40\), o teste qui-quadrado de Pearson com correção de Yates pode ser usado.
Cochran, 1954:
- Um valor esperado mínimo de 1 em alguma célula é permitido, desde que não mais que 20% das células tenham valor abaixo de 5.
Então, o valor de \(X^2\) era dado por
\[ X^2 = \sum{ ({\text{observado}-\text{esperado}})^2 \over {\text{esperado}} } \]
neste exemplo calculado por:
\[ \begin{align} X^2 &= { {{(17-28.7)^2} \over {28.7}} + {{(138-126.3)^2} \over {126.3}} +\\ {{(130-118.3)^2} \over {118.3}} + {{(508-519.7)^2} \over {519.7}} }\\ X^2&\approx 7.31 \end{align} \]
ou, em tabelas 2×2, era recomendado usar a correção de continuidade de Yates, dada por:
\[ X^2 = \sum{ (|{\text{observado}-\text{esperado}}|-0.5)^2 \over {\text{esperado}} } \]
Neste exemplo:
\[ \begin{align} X^2 &= \dfrac{(|17-28.7|-0.5)^2}{28.7} + \dfrac{(|138-126.3|-0.5)^2}{126.3} + \\ &\quad\dfrac{(|130-118.3|-0.5)^2}{118.3} + \dfrac{(|508-519.7|-0.5)^2}{519.7}\\ X^2&\approx 6.7 \end{align} \]
Então, era necessário consultar uma tabela com os valores críticos de \(\chi^2_c\) previamente calculados para a tomada de decisão (sem obter, portanto, o valor p). A tabela permitia escolher o valor de \(\alpha\) e, sabendo-se os graus de liberdade (neste caso \(\nu = 1\)), localizar \(\chi^2_c\) (neste exemplo igual a 3.841):
Cálculo manual do teste exato de Fisher
Caso \(\chi^2\) não atendesse as exigências para ser aplicado, a alternativa (se fosse tabela 2×2) era o Teste Exato de Fisher. Este envolve calcular o valor p por fatoriais de tabelas progressivamente mais extremas. Por exemplo:
Neste caso o valor exato de \(p\) era calculado e comparado diretamente com \(\alpha\) para a tomada de decisão.
O que diz a literatura mais recente sobre estes métodos tradicionais e não robustos?
Conforme Ludbrook (2011, p. 925), em seu artigo intitulado “Ainda há lugar para o teste qui-quadrado de Pearson e o teste exato de Fisher na pesquisa cirúrgica?”
“Conclusões:
Não posso concluir outra coisa senão que o teste qui-quadrado de Pearson e o teste exato de Fisher quase nunca devem ser usados para analisar os resultados de estudos cirúrgicos. Os fundamentos para essa conclusão são os seguintes:
Esses testes têm pré-requisitos que raramente, se é que algum dia, podem ser cumpridos em estudos da vida real (ver Tabela 3).
Ambos os procedimentos testam hipóteses muito não específicas (ver Tabela 3). É muito melhor usar procedimentos que são projetados para detectar desigualdade de proporções, uma razão de proporções que difere da unidade ou uma tendência nas proporções.
É melhor usar testes de permutação exatos em vez de testes aproximados baseados nas distribuições qui-quadrado ou normal.”
Métodos avançados: delineamento entre participantes
No teste qui-quadrado de Pearson tradicional ambas as variáveis são tratadas como nominais e todos os indivíduos são distribuídos em células, de acordo com a contagem. Há situações em que informações podem ser estratificadas, ou as categorias podem ser ordinais. Os métodos descritos aqui servem para aproveitar estas informações.
Como nos casos anteriores, o delineamento pode ser:
- entre participantes ou
- intraparticipantes
Teste qui-quadrado de Cochran-Mantel-Haenszel
Appleton et al. (1996) apresentam um estudo no qual mulheres foram avaliadas quanto à sobrevida e o tabagismo.
| Morta | Viva | |
|---|---|---|
| Tabagista | 139 | 443 |
| Não tabagista | 230 | 502 |
Podemos aplicar o teste qui-quadrado de Pearson ou o teste exato de
Fisher nesta tabela 2×2, utilizando o código R TabelaContingencia2x2_TabagismoMorte.R
e os dados em fumo_e_faixaetaria_total.xlsx,
obtendo
# TabelaContingencia2x2_TabagismoMorte.R
# paradoxo de Simpson
library(exact2x2)
source("eiras.show.MCSTAR.R")
TCtmp <- data.frame(readxl::read_excel("fumo_e_faixaetaria_total.xlsx"))
TC <- data.matrix(TCtmp)
rownames(TC) <- as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
gplots::balloonplot(t(as.table(TC)), label = FALSE,
xlab = "", ylab = "",
show.margins = TRUE)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(ftable(TC))
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
B <- 1e5
cat("\nTeste qui-quadrado de Pearson exato:\n")
print(res <- chisq.test(TC,simulate.p.value=TRUE,B=B))
cat("X^2 critico de 95% = ",qchisq(p=1-alfa,df=1), "\n", sep="")
N <- sum(res$observed)
nL <- nrow(TC) # ou dim(TC)[1]
nC <- ncol(TC) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
show.MCSTAR(MCSTAR=MCSTAR,alpha=alfa,df=df)
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica: V de Cramer")
cat("\n------------------------------------------\n")
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
V <- sqrt(phi2/Dim) # V ou C de Cramer
cat("\nV de Cramer = ", round(V,4), "\n", sep="")
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(TC),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
cat("\nTeste de razao de chances (OR):\n")
resft <- exact2x2::fisher.exact(TC,conf.level=1-alfa) # Teste de OR
print(resft)
------------------------------------------
Dados:
------------------------------------------
Morta Viva
Tabagista 139 443
Não tabagista 230 502
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: TC
X-squared = 9.1209, df = NA, p-value = 0.00289
X^2 critico de 95% = 3.841459
Graus de liberdade (nao fornecidos por bootstrapping): 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 9.120903
------------------------------------------
Analise post hoc:
------------------------------------------
Residuos ajustados standardizados corrigidos por momento (MCSTAR):
|MCSTAR critico| (alphaBonferroni=5%/1) = 1.959964
Morta Viva
Tabagista -6.04 #6.04
Não tabagista #6.04 -6.04
------------------------------------------
Analise de significancia pratica: V de Cramer
------------------------------------------
V de Cramer = 0.0833
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.08331
Confidence interval:
2.5% 97.5%
0.02822 0.13470
Teste de razao de chances (OR):
Two-sided Fisher's Exact Test (usual method using minimum likelihood)
data: TC
p-value = 0.002989
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.5338 0.8786
sample estimates:
odds ratio
0.6850392 A conclusão é pela rejeição de \(H_0\). Observando os MCSTAR e o \(\text{OR}\) fornecidos, verificamos que o fumo tem efeito protetor: tabagismo, portanto está associado com vida e não tabagismo com morte. Algo parece estranho?
Neste exemplo temos a contagem de mortes entre fumantes e não fumantes estratificadas por faixa etária e utilizaremos o teste conhecido, em inglês, como Cochran-Mantel-Haenszel Chi-Squared Test for Count Data. Esta versão de qui-quadrado faz o teste global e, depois, estratifica pela terceira variável (faixa etária, neste exemplo).
| Tabagista | Desfecho | 18-24 | 25-34 | 35-44 | 45-54 | 55-64 | 65-74 | 75- |
|---|---|---|---|---|---|---|---|---|
| Sim | Morta | 2 | 3 | 14 | 27 | 51 | 29 | 13 |
| Viva | 53 | 121 | 95 | 103 | 64 | 7 | 0 | |
| Não | Morta | 1 | 5 | 7 | 12 | 40 | 101 | 64 |
| Viva | 61 | 152 | 114 | 66 | 81 | 28 | 0 |
Para aproveitar a estratificação por idade, os dados da tabela
precisam ser reorganizados para que construamos uma tabela
tridimensional requerida pela funções relacionadas a este teste. Em
nossa solução, verifique como os dados foram acomodados na planilha fumo_e_faixaetaria.xlsx.
O código está implementado em Teste qui-quadrado de Mantel-Haenszel.R,
que lê esta planilha, rearranja os dados em uma tabela tridimensional
adequada e aplica o teste de Mantel-Haenszel obtendo:
# Teste qui-quadrado de Mantel-Haenszel.R
# Tabela de contingencia 2x2 segmentada
# Teste robusto da razao de chances (OR) de Mantel-Haenszel
# Relacao entre tabagismo e sobrevivência em 20 anos (1974-1994)
# em 1.134 mulheres adultas do Reino Unido
# Delineamento: coorte
# Fonte: APPLETON, D. R. et al. (1996) Ignoring a covariate:
# An example of Simpson's paradox. The American Statistician,
# 50(4): 340-1.
TCtmp <- readxl::read_excel("fumo_e_faixaetaria.xlsx")
TCtmp <- tidyr::gather(TCtmp, var, Contagem, -Idade)
TCtmp <- tidyr::separate(TCtmp, var, c('Exposicao','Desfecho'))
TCS <- TCtmp %>% xtabs(Contagem ~ factor(Exposicao, levels=c("Tabagista","NaoTabagista")) + factor(Desfecho, levels=c("Morta", "Viva")) + Idade, .)
print(TCS)
# Mantel-Haenszel test of the null that two nominal variables
# are conditionally independent in each stratum,
# assuming that there is no three-way interaction.
# Woolf test for homogeneity of ORs across strata on 2 x 2 x k tables:
# If significant, M-H test is not appropriate.
cat("\nTabela global horizontalizada:\n")
prmatrix(ftable(TCS)) # Display a flattened table (tabela horizontalizada)
cat("\nHomogeneidade dos OR:\n")
print(vcd::woolf_test(TCS))
# Teste para tabela 2x2xK
cat("\nQui-quadrado de Mantel-Haenszel:\n")
print(mh <- mantelhaen.test(TCS, exact=TRUE)) # Teste exato de OR de Mantel-Haenszel
# Teste para tabela LxCxK
cat("\nTeste de Cochran-Mantel-Haenszel (coin::cmh_test):\n")
print(coin::cmh_test(TCS, distribution = coin::approximate(nresample = 1e5)))
cat("\nAnalise estratificada por idade:\n")
print(rcompanion::groupwiseCMH(TCS,
group = 3,
fisher = TRUE,
gtest = FALSE,
chisq = FALSE,
method = "bonferroni",
correct = "none",
digits = 3))
cat("\n")
print(ors <- vcd::oddsratio(TCS, log=FALSE)) # Show OR for each 2x2
cat("\n")
lnors <- vcd::oddsratio(TCS, log=TRUE) # Show ln(OR) for each 2x2
print(lnors)
print(summary(lnors), digits=3)
print(confint(lnors), digits=3)
plot(lnors, main="Intervalos de confiança por faixa etária",
xlab = "Faixa Etaria", ylab = "ln(OR)"), , Idade = 18-24
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 2 53
NaoTabagista 1 61
, , Idade = 25-34
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 3 121
NaoTabagista 5 152
, , Idade = 35-44
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 14 95
NaoTabagista 7 114
, , Idade = 45-54
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 27 103
NaoTabagista 12 66
, , Idade = 55-64
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 51 64
NaoTabagista 40 81
, , Idade = 65-74
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 29 7
NaoTabagista 101 28
, , Idade = 75+
factor(Desfecho, levels = c("Morta", "Viva"))
factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) Morta Viva
Tabagista 13 0
NaoTabagista 64 0
Tabela global horizontalizada:
18-24 25-34 35-44 45-54 55-64 65-74 75+
Tabagista_Morta 2 3 14 27 51 29 13
Tabagista_Viva 53 121 95 103 64 7 0
NaoTabagista_Morta 1 5 7 12 40 101 64
NaoTabagista_Viva 61 152 114 66 81 28 0
Homogeneidade dos OR:
Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
data: TCS
X-squared = 3.2061, df = 6, p-value = 0.7826
Qui-quadrado de Mantel-Haenszel:
Exact conditional test of independence in 2 x 2 x k tables
data: TCS
S = 139, p-value = 0.01591
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
1.068889 2.203415
sample estimates:
common odds ratio
1.530256
Teste de Cochran-Mantel-Haenszel (coin::cmh_test):
Approximative Generalized Cochran-Mantel-Haenszel Test
data: factor.Desfecho..levels...c..Morta....Viva... by
factor.Exposicao..levels...c..Tabagista....NaoTabagista... (Tabagista, NaoTabagista)
stratified by Idade
chi-squared = 5.8443, p-value = 0.01585
Analise estratificada por idade:
Group Test p.value adj.p
1 18-24 Fisher 0.6000 1.000
2 25-34 Fisher 1.0000 1.000
3 35-44 Fisher 0.0706 0.494
4 45-54 Fisher 0.3650 1.000
5 55-64 Fisher 0.0830 0.581
6 65-74 Fisher 1.0000 1.000
7 75+ Fisher 1.0000 1.000
odds ratios for factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) and factor(Desfecho, levels = c("Morta", "Viva")) by Idade
18-24 25-34 35-44 45-54 55-64 65-74 75+
1.9158879 0.7987280 2.3179756 1.4135266 1.6067566 1.1044335 0.2093023
log odds ratios for factor(Exposicao, levels = c("Tabagista", "NaoTabagista")) and factor(Desfecho, levels = c("Morta", "Viva")) by Idade
18-24 25-34 35-44 45-54 55-64 65-74
0.65018114 -0.22473479 0.84069420 0.34608770 0.47421763 0.09933253
75+
-1.56397554
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
18-24 0.6502 1.0496 0.62 0.536
25-34 -0.2247 0.6945 -0.32 0.746
35-44 0.8407 0.4706 1.79 0.074 .
45-54 0.3461 0.3756 0.92 0.357
55-64 0.4742 0.2681 1.77 0.077 .
65-74 0.0993 0.4606 0.22 0.829
75+ -1.5640 2.0223 -0.77 0.439
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2.5 % 97.5 %
18-24 -1.4070 2.71
25-34 -1.5859 1.14
35-44 -0.0817 1.76
45-54 -0.3900 1.08
55-64 -0.0513 1.00
65-74 -0.8035 1.00
75+ -5.5276 2.40
Procure a sessão “Qui-quadrado de Mantel-Haenszel” nesta saída. Para
\(\alpha=0.05\), a associação entre
fumo e morte é significante (\(p \approx
0.016\)) mas o \(\text{OR}\)
generalizado (mh$estimate=1.5302557) deste teste deixa o
valor unitário à esquerda do intervalo
[mh$conf.int[1:2]=1.068889, 2.2034151], indicando que
Tabagismo está associado com Morte, diferente do que foi obtido com o
teste de Qui-quadrado tradicional em uma tabela 2×2. Esta inversão de
conclusão é conhecida como paradoxo de Simpson. Verifique,
também, o trecho da saída que mostra odds ratios for
factor, mostrando que duas das idades têm odds ratio
menores que 1 (correspondendo, na saída gráfica, ao logarítmo de \(\text{OR}\) abaixo de 0). O grupo de 25 a
34 anos é numeroso e o de maiores que 75 anos têm zeros; possivelmente
foram estes dois grupos que distorceram o teste qui-quadrado, invertendo
a conclusão quando a covariável idade não foi considerada.
Teste qui-quadrado mínimo de Serra
Segundo este autor, a correção de Yates não deve ser usada:
Note a relação entre esta forma de calcular o qui-quadrado e o odds-ratio:
Teste gama de Goodman-Kruskal para variáveis ordinais
O teste qui-quadrado de Pearson tradicional trata os dois fatores com variáveis nominais. Há muitos casos, porém, em que estes fatores são ordinais. Embora não seja incorreto utilizar o qui-quadrado de Pearson (toda variável ordinal é também nominal), alguma informação está sendo perdida.
Existe uma medida, gama de Goodman-Kruskal, que considera ambas as variáveis como categóricas ordinais, potencialmente melhorando o resultado do teste. É uma estatistica que mede a força de associação (também de tamanho de efeito) em tabela de contingência, com as hipóteses:
\[ \begin{cases} H_0: \gamma = 0\\ H_1: \gamma \ne 0 \end{cases}\\ \alpha=5\% \]
Além disto, existe o \(\text{OR}\) generalizado proposto por Edwards & Baltzan (2000).
|
Especificamente para tabelas 2×2 existe o Q de Yule, que corresponde ao gama de Goodman-Kruskal para este caso particular, dado por \[ Q = \dfrac{ad-bc}{ad+bc} \] Uma função está implementada em Assim como o gama de Goodman-Kruskal, varia de \(-1\) a \(+1\), mas como se aplica a tabelas 2×2, relaciona-se com o odds-ratio: \[ Q = \dfrac{\text{OR}-1}{\text{OR}+1} \] Consequentemente, \(\text{OR}\) é conceitualmente ligado ao gama de Goodman-Kruskal e, portanto, mede a associação entre variáveis ordinais dicotômicas, ou equivalentemente, entre duas variáveis nominais dicotômicas. A generalização nos permite tratar \(\text{OR}\) como medida de associação para variáveis ordinais politômicas dispostas em tabelas de contigência com qualquer número de linhas e colunas, nem mesmo necessitando ser uma matriz quadrada: \[ \text{OR} = {{1+\gamma} \over {1-\gamma}} \] |
Para o \(\text{OR}\) generalizado as hipóteses são:
\[ \begin{cases} H_0: &\text{OR} = 1\\ H_1: &\text{OR} \ne 1 \end{cases} \]
Um exemplo de aplicação é fornecido pelo próprio autor (Goodman, 1979), desenvolvendo a análise com variáveis ordinais.
A tabela a seguir mostra a classificação cruzada de uma amostra de
3498 homens britânicos segundo a categoria de status ocupacional de cada
indivíduo e a categoria de status ocupacional de seu pai, utilizando
oito categorias de status. Os dados desta tabela são, originalmente, de
Duncan (1979) e estão disponíveis em
datasets::occupationalStatus.
As categorias de status são as seguintes:
I. Profissional e alta administração
Gerencial e executiva
Inspetoria, supervisão e outras ocupações não manuais (nível superior)
Inspetoria, supervisão e outras ocupações não manuais (nível inferior)
V(a). Ocupações rotineiras não manuais
V(b). Trabalho manual qualificado
Trabalho manual semiqualificado
Trabalho manual não qualificado
Portanto, a variável de status ocupacional é, claramente, ordinal.
O repositório do R tem mais que pacotes. Há vários bancos de dados
armazenados nele, para uso em exemplos. Aqui usaremos os dados deste
estudo, disponíveis no pacote datasets.
|
Implementado em TCrXc_StatusOcupacional.R, a saída
é:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source ("eiras.cramerV.R")
source ("eiras.show.MCSTAR.R")
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
# Duncan, O.D. (1979), How Destination Depends on Origin in the
# Occupational Mobility Table. American Journal of Sociology 84, 793-803.
tcrxc <- datasets::occupationalStatus
# para recolocar as categorias originais
categories <- c("I","II","III","IV","Va","Vb","VI","VII")
rownames(tcrxc) <- categories
colnames(tcrxc) <- categories
print(tcrxc, row.names=FALSE)
gplots::balloonplot(t(as.table(tcrxc)), label = FALSE,
xlab = "StatusFilho", ylab = "StatusPai",
show.margins = TRUE)
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato\n")
cat("(assumindo as variáveis como nominais):\n")
res <- chisq.test(tcrxc,simulate.p.value=TRUE,B=1e5)
aux <- chisq.test(tcrxc)
res$parameter <- aux$parameter
print(res)
N <- sum(res$observed)
nL <- nrow(tcrxc) # ou dim(TC)[1]
nC <- ncol(tcrxc) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
cat("X^2 critico de ",round((1-alfa)*100,1),"% = ",qchisq(p=1-alfa,df=df), "\n", sep="")
X2 <- res$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade (nao fornecidos por bootstrapping): ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\nTestes Gama de Goodman-Kruskal e OR generalizado\n")
cat("(assumindo as variáveis como ordinais):\n")
# Goodman, L. A. (1979) Simple Models for the Analysis of Association
# in Cross-Classifications having Ordered Categories.
# J. Am. Stat. Assoc., 74 (367), 537–552.
# gama <- DescTools::GoodmanKruskalGamma(tcrxc)
gama <- vcdExtra::GKgamma(tcrxc,level=1-alfa)
cat("\nGama de Goodman-Kruskal = ", round(gama$gamma,2),
"\n\tIC",round((1-alfa)*100,1),"%(gama G-K) = [",round(gama$CI[1],2),", ",round(gama$CI[2],2),"]\n", sep="")
# EDWARDES, MD & BALTZAN, M (2000) The generalization of the odds ratio,
# rik ratio and risk difference to rxk tables.
# Statistics in Medicine 19:1901-14.
ORg <- (1+gama$gamma)/(1-gama$gamma)
ORg.LI <- (1+gama$CI[1])/(1-gama$CI[1])
ORg.LS <- (1+gama$CI[2])/(1-gama$CI[2])
cat("OR bruto generalizado = (1 + ",round(gama$gamma,2),")/(1 - ",round(gama$gamma,2),") = ",round(ORg,2),sep="")
cat("\n\tIC",round((1-alfa)*100,1),"%(OR bruto generalizado) = [",round(ORg.LI,2), ", ", round(ORg.LS,2),"]\n",sep="")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("(assumindo as variáveis como nominais):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
show.MCSTAR(MCSTAR,alfa,df)
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
cat("(assumindo as variáveis como nominais):\n")
V <- sqrt(phi2/Dim) # V ou C de Cramer
cat("\nV de Cramer = ", round(V,4), "\n", sep="")
print(effectsize::interpret_cramers_v(V))
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(tcrxc),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
------------------------------------------
Dados:
------------------------------------------
destination
origin I II III IV Va Vb VI VII
I 50 19 26 8 7 11 6 2
II 16 40 34 18 11 20 8 3
III 12 35 65 66 35 88 23 21
IV 11 20 58 110 40 183 64 32
Va 2 8 12 23 25 46 28 12
Vb 12 28 102 162 90 554 230 177
VI 0 6 19 40 21 158 143 71
VII 0 3 14 32 15 126 91 106
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Pearson exato
(assumindo as variáveis como nominais):
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: tcrxc
X-squared = 1416, df = 49, p-value = 1e-05
X^2 critico de 95% = 66.33865
Graus de liberdade (nao fornecidos por bootstrapping): 49
Heuristica de significancia (rej. H0 se X^2/gl > 2): 28.89877
Testes Gama de Goodman-Kruskal e OR generalizado
(assumindo as variáveis como ordinais):
Gama de Goodman-Kruskal = 0.42
IC95%(gama G-K) = [0.39, 0.45]
OR bruto generalizado = (1 + 0.42)/(1 - 0.42) = 2.45
IC95%(OR bruto generalizado) = [2.29, 2.64]
------------------------------------------
Analise post hoc:
------------------------------------------
(assumindo as variáveis como nominais):
Residuos ajustados standardizados corrigidos por momento (MCSTAR):
|MCSTAR critico| (alphaBonferroni=5%/49) = 3.284839
I II III IV Va Vb VI VII
I #28.02 #6.47 #4.85 -2.71 -0.8 -7.09 -4.34 -4.28
II #6.54 #15.19 #6.48 -0.48 0.2 -6.22 -4.43 -4.44
III 0.71 #6.01 #7.2 #3.98 2.78 -3.97 -6.13 -4.13
IV -1.37 -0.93 1.7 #6.77 0.83 0.85 -3.45 -5.13
Va -1.44 0.41 -0.87 0.7 #5.19 -1.36 0.39 -1.98
Vb -6.55 -6.4 -3.51 -1.86 -0.7 #7.93 0.03 1.55
VI -4.57 -4.08 -4.74 -3.41 -2.46 0.33 #9.98 2.72
VII -4.15 -4.32 -4.74 -3.43 -2.9 -0.68 #4.17 #11.15
------------------------------------------
Analise de significancia pratica:
------------------------------------------
(assumindo as variáveis como nominais):
V de Cramer = 0.2405
X-squared
"medium"
(Rules: funder2019)
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.2405
Confidence interval:
2.5% 97.5%
0.2164 0.2575 É praticamente o mesmo código utilizado para o exemplo de tabagismo e etilismo, com algumas diferenças:
- para efeito de comparação, exibimos o resultado do qui-quadrado de Pearson que trata as variáveis como nominais.
- em uma tabela com 8 linhas e 8 colunas não existe definição para o odds ratio. No lugar deste, aparecem o gama de Goodman-Kruskal e o \(\text{OR}\) generalizado, aqui aproveitando a informação de ordem das variáveis. As decisões são tomadas com base nos respectivos intervalos de confiança; neste caso, as \(H_0\) são rejeitadas e, portanto, há associação entre a categoria profissionais de pais e filhos.
- a análise post hoc e de tamanho de efeito tratam as variáveis como
nominais. A função implementada em
eiras.show.MCSTARexibe a matriz resultante de MCSTAR com menor número de casas decimais e símbolos adicionais (“#”) para que os valores acima de \(z\) crítico sejam mais facilmente localizados. Alternativamente, utilizamoscorrplot::corrplotpara exibir o gráfico correspondente.
Teste qui-quadrado de Mantel-Haenszel para duas variáveis ordinais
A seguir é apresentado um exemplo com o uso de funções R para testar a dissociação (independência) entre duas variáveis ordinais.
O qui-quadrado de Mantel-Haenszel para variáveis ordinais utiliza a
implementação disponível em DescTools::MHChisqTest. Existe
uma alternativa, o Approximative Linear-by-Linear Association
Test feito com bootstrapping e implementado em
coin::chisq_test.
Também comparamos o resultado do teste qui-quadrado de Pearson por bootstrapping, mas que considera as duas variáveis como nominais.
Ambos estão implementados em Teste_2Ordinais.R. A saída é:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source ("eiras.show.MCSTAR.R")
Job <- matrix(c(1,2,1,0,
3,3,6,1,
10,10,14,9,
6,7,12,11),
4, 4,
dimnames = list(income = c("<15k", "15-25k",
"25-40k", ">40k"),
satisfaction = c("VeryD", "LittleD",
"ModerateS", "VeryS")))
print(ftable(Job))
gplots::balloonplot(t(as.table(Job)), label = FALSE,
xlab = "satisfaction", ylab = "income",
show.margins = TRUE)
nL <- nrow(Job) # ou dim(TC)[1]
nC <- ncol(Job) # ou dim(TC)[2]
df <- (nL-1)*(nC-1)
alfa <- 0.05
cat("\nTeste qui-quadrado de Pearson exato\n")
cat("(assumindo as variáveis como nominais):\n")
res <- chisq.test(Job,
simulate.p.value=TRUE,
B=1e6)
aux <- chisq.test(Job)
res$parameter <- aux$parameter
X2 <- res$statistic
N <- sum(res$observed)
print(res)
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
cat("(assumindo as variáveis como nominais):\n")
STAR <- res$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
show.MCSTAR(MCSTAR,alfa,df)
phi2 <- X2/N # phi ou w de Cohen
Dim <- min(nL,nC)-1 # dimensoes da TC: dim = max(phi^2)
# The Mantel-Haenszel chi-square statistic tests the alternative hypothesis
# that there is a linear association between the row variable and the
# column variable.
# Both variables must lie on an ordinal scale.
# 7.5 = (15-0)/2
# 20 = (25-15)/2
# 32.5 = (40-25)/2
# 60 = (160-40)/2
#
# > 20-7.5
# [1] 12.5
# > 32.5-20
# [1] 12.5
# > 60-32.5
# [1] 27.5
# > (160-40)/2
# [1] 60
mhtest <- DescTools::MHChisqTest(Job,
srow=c(7.5,20,32.5,60))
print(mhtest)
# X-squared = 3.8075, df = 1, p-value = 0.05102
mhtest <- DescTools::MHChisqTest(Job)
print(mhtest)
# X-squared = 2.983, df = 1, p-value = 0.08414
# A diferenca está na suposicao que o valor máximo de income é 160.
# Approximative Linear-by-Linear Association Test
cointest <- coin::chisq_test(as.table(Job),
distribution =
coin::approximate(nresample = 1e6),
scores =
list(income = 1:ncol(Job),
satisfaction = 1:nrow(Job)))
print(cointest)
cat("\n------------------------------------------\n")
cat("Analise de significancia pratica:")
cat("\n------------------------------------------\n")
cat("Duas variáveis ordinais:\n")
print(vcdExtra::GKgamma(as.table(Job)))
cat("\nDuas variáveis nominais:\n")
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(Job),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4) satisfaction VeryD LittleD ModerateS VeryS
income
<15k 1 3 10 6
15-25k 2 3 10 7
25-40k 1 6 14 12
>40k 0 1 9 11
Teste qui-quadrado de Pearson exato
(assumindo as variáveis como nominais):
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: Job
X-squared = 5.9655, df = 9, p-value = 0.7702
------------------------------------------
Analise post hoc:
------------------------------------------
(assumindo as variáveis como nominais):
Residuos ajustados standardizados corrigidos por momento (MCSTAR):
|MCSTAR critico| (alphaBonferroni=5%/9) = 2.772921
VeryD LittleD ModerateS VeryS
<15k 0.28 0.29 0.7 -1.04
15-25k 1.76 0.02 0.09 -0.84
25-40k -0.54 1.28 -0.45 -0.22
>40k -1.44 -1.77 -0.27 2.12
Mantel-Haenszel Chi-Square
data: Job
X-squared = 3.8075, df = 1, p-value = 0.05102
Mantel-Haenszel Chi-Square
data: Job
X-squared = 2.983, df = 1, p-value = 0.08414
Approximative Linear-by-Linear Association Test
data: satisfaction (ordered) by income (<15k < 15-25k < 25-40k < >40k)
Z = 1.7362, p-value = 0.09343
alternative hypothesis: two.sided
------------------------------------------
Analise de significancia pratica:
------------------------------------------
Duas variáveis ordinais:
gamma : 0.221
std. error : 0.117
CI : -0.009 0.451
Duas variáveis nominais:
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.1439
Confidence interval:
2.5% 97.5%
0.0000 0.1768 Observe que a consequência do aproveitamento da informação de ordinalidade das variáveis reflete-se no valor p dos dois primeiros testes em comparação com o obtido por meio do qui-quadrado tradicional.
Teste qui-quadrado para tendência de Cochran-Armitage: uma variável ordinal (dose) e nominal dicotômica (resposta)
Para situações em que uma variável é ordinal (dose) e a outra é nominal (resposta), por exemplo quando há um modelo de dose-resposta com um desfecho dicotômico associado a uma exposição que tem mais que dois níveis, é preferível utilizar o teste conhecido como Cochran-Armitage test for trend.
Exemplo: espessura da dobra cutânea e menarca
Os pesquisadores Beckles et al. (1985) estudaram a associação entre a variável de a faixa etária da menarca (\(<\) 12 anos e \(\ge\) 12 anos) e o fator a espessura da dobra cutânea do tríceps (pequena, intermediária e grande).
|
|
|
Caso ambas as variáveis sejam consideradas nominais, o teste adequado é o qui-quadrado de Pearson. Caso a variável considerada como exposição seja ordinal e o desfecho seja dicotômico (tabela \(k\times 2\)), o teste qui-quadrado de Cochran-Armitage pode ser usado.
Para comparação, ambos os testes são executados em Teste qui-quadrado de Cochran-Armitage.R
que produz:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
TC <- matrix(c(15,156,
29,197,
36,186),
nrow = 3, byrow = TRUE,
dimnames = list(
c("Pequena","Intermediária","Grande"),
c("<12",">=12")
))
gplots::balloonplot(t(as.table(TC)),
xlab = "Menarca",
ylab = "Espessura da dobra cutânea",
label = TRUE, show.margins = TRUE)
# ------------------------------------------------------------
# Qui-quadrado de Pearson (com simulação, se quiser)
# ------------------------------------------------------------
alfa <- 0.05
res_chi <- suppressWarnings(chisq.test(TC, simulate.p.value = TRUE, B = 1e5))
N <- sum(TC)
nL <- nrow(TC); nC <- ncol(TC)
df_chi <- (nL - 1) * (nC - 1)
X2 <- unname(res_chi$statistic)
cat("\nTeste qui-quadrado de Pearson (simulado):\n")
print(res_chi)
cat(sprintf("gl = %d\n", df_chi))
cat(sprintf("X^2 crítico (%.0f%%) = %.4f\n", (1-alfa)*100, qchisq(p=1-alfa, df=df_chi)))
cat(sprintf("Heurística X^2/gl = %.3f\n", X2/df_chi))
# ------------------------------------------------------------
# Cochran–Armitage (requer 2×k ou k×2)
# Escolha orientação automaticamente.
# ------------------------------------------------------------
if(nL == 2 || nC == 2) {
# DescTools detecta orientação; pesos padrão = 1:k (tendência linear)
res_ca <- DescTools::CochranArmitageTest(TC)
x2ca <- unname(res_ca$statistic)^2 # estatística ao quadrado ~ χ^2_1
p_ca <- res_ca$p.value
cat("\nTeste de tendência de Cochran–Armitage:\n")
print(res_ca)
cat(sprintf("\nX^2(trend) = %.4f, gl = 1, p-valor = %.5f\n", x2ca, p_ca))
} else {
cat("\n[Aviso] Cochran–Armitage não aplicado: a tabela não é 2×k nem k×2.\n")
}
# ------------------------------------------------------------
# Tamanho de efeito (V de Cramér) + IC bootstrap
# ------------------------------------------------------------
phi2 <- X2 / N
Dim <- min(nL, nC) - 1
V <- sqrt(phi2 / Dim)
cat("\nTamanho de efeito:\n")
cat(sprintf("V de Cramér = %.4f\n", V))
print(effectsize::interpret_cramers_v(V))
set.seed(123)
ciV <- confintr::ci_cramersv(
x = as.table(TC),
probs = c(alfa/2, 1 - alfa/2),
type = "bootstrap",
R = 1e4
)
print(ciV, digits = 4)
Teste qui-quadrado de Pearson (simulado):
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: TC
X-squared = 4.7594, df = NA, p-value = 0.09238
gl = 2
X^2 crítico (95%) = 5.9915
Heurística X^2/gl = 2.380
Teste de tendência de Cochran–Armitage:
Cochran-Armitage test for trend
data: TC
Z = 2.1783, dim = 3, p-value = 0.02938
alternative hypothesis: two.sided
X^2(trend) = 4.7449, gl = 1, p-valor = 0.02938
Tamanho de efeito:
V de Cramér = 0.0877
[1] "very small"
(Rules: funder2019)
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.08769
Confidence interval:
2.5% 97.5%
0.0000 0.1493 Compare os valores \(p\) de ambos os testes. Cochran-Armitage, ao considerar a espessura da prega cutânea como “dose”, ordinal, indica significância estatística onde o teste convencional não detectava.
Teste de Cochran-Armitage generalizado: uma variável ordinal politômica (dose) e uma variável nominal politômica (resposta)
Para desfecho com três ou mais categorias (variável de desfecho
nominal) a função no pacote coin::chisq.test pode ser
parametrizada. Apresenta-se a opção de teste post hoc com
rcompanion::pairwiseOrdinalIndependence.
Um exemplo é investigar se o salário e a satisfação dos empregados
com o trabalho estão associados (neste exemplo optamos por tratar apenas
o nível salarial como ordinal, mas deixamos o grau de satisfação como
nominal - caso queira tratar as duas como ordinais, veja Mantel-Haenszel
para variáveis ordinais, acima). O exemplo está implementado em Teste qui-quadrado de Cochran-Armitage generalizado.R,
que produz:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
# The chisq_test function in the coin package can be used to conduct a test of association for a contingency table with one ordered nominal variable and one non-ordered nominal variable.
# The Cochran–Armitage test is a special case of this where the non-ordered variable has only two categories.
# The scores option is used to indicate which variable should be treated as ordered, and the spacing of the levels of this variable.
Job <- matrix(c(1,2,1,0,
3,3,6,1,
10,10,14,9,
6,7,12,11),
4, 4,
dimnames = list(income = c("< 15k", "15-25k",
"25-40k", "> 40k"),
satisfaction = c("VeryD", "LittleD",
"ModerateS", "VeryS")))
Job <- as.table(Job)
ftable(Job)
gplots::balloonplot(t(as.table(Job)), label = FALSE,
xlab = "satisfaction", ylab = "income",
show.margins = TRUE)
cat("\nDuas variáveis ordinais:\n")
print(out <- vcdExtra::GKgamma(as.table(Job)))
cat("\nDuas variáveis nominais:\n")
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(Job),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
# teste tradicional
cat("\nQui-quadrado de Pearson\n")
chisqtest <- chisq.test(Job, simulate.p.value = TRUE, B=1e6)
chisqtest$parameter <- (nrow(Job)-1)*(ncol(Job)-1)
print(chisqtest)
cat("\nQui-quadrado de Cochran-Armitage generalizado\n")
print(coin::chisq_test(Job,
scores = list("income"= 1:4)))
cat("\nTeste post hoc\n")
print(ph <- rcompanion::pairwiseOrdinalIndependence(Job,
compare="column"))
sig <- NA
try(sig <- rcompanion::cldList(p.value ~ Comparison,
data = ph,
threshold = 0.05),
silent = TRUE)
if (sum(is.na(sig))==0)
{
print(sig)
}
Duas variáveis ordinais:
gamma : 0.221
std. error : 0.117
CI : -0.009 0.451
Duas variáveis nominais:
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.1439
Confidence interval:
2.5% 97.5%
0.0000 0.1768
Qui-quadrado de Pearson
Pearson's Chi-squared test with simulated p-value (based on 1e+06
replicates)
data: Job
X-squared = 5.9655, = 9, p-value = 0.7697
Qui-quadrado de Cochran-Armitage generalizado
Asymptotic Generalized Pearson Chi-Squared Test
data: satisfaction by income (< 15k < 15-25k < 25-40k < > 40k)
chi-squared = 3.1362, df = 3, p-value = 0.3711
Teste post hoc
Comparison p.value p.adjust
1 VeryD : LittleD 0.451 0.541
2 VeryD : ModerateS 0.351 0.526
3 VeryD : VeryS 0.161 0.526
4 LittleD : ModerateS 0.698 0.698
5 LittleD : VeryS 0.242 0.526
6 ModerateS : VeryS 0.271 0.526
Group Letter MonoLetter
1 VeryD a a
2 LittleD a a
3 ModerateS a a
4 VeryS a aTeste qui-quadrado de Agresti para duas variáveis ordinais politômicas e estratificado por variável nominal politômica
Agresti (2002) apresenta uma generalização do teste de Mantel-Haenszel estratificado, aqui aproveitando a informação de ordem das duas variáveis principais, analisadas em estratos (variável nominal).
São usadas as funções coin::lbl_test e
coin::cmh_test.
Usaremos o exemplo de associação entre renda e satisfação com o trabalho, estratificando por sexo.
O exemplo está implementado em Teste qui-quadrado de Agresti.R, que
produz:
# stratified ordinal chi-square test
## Agresti (2002), p. 287 and 297
## Job Satisfaction example
Satisfaction <-
as.table(array(c(1, 2, 0, 0, 3, 3, 1, 2,
11, 17, 8, 4, 2, 3, 5, 2,
1, 0, 0, 0, 1, 3, 0, 1,
2, 5, 7, 9, 1, 1, 3, 6),
dim = c(4, 4, 2),
dimnames =
list(Income =
c("<5000", "5000-15000",
"15000-25000", ">25000"),
"Job Satisfaction" =
c("V_D", "L_S", "M_S", "V_S"),
Gender = c("Female", "Male"))))
print(ftable(Satisfaction))
names(dimnames(Satisfaction)) <- c("Income","JobSatisfaction","Gender")
vcd::mosaic(~ Income + JobSatisfaction + Gender,
data = Satisfaction, shade = TRUE, legend = TRUE)
print(coin::lbl_test(Satisfaction,
distribution =
coin::approximate(nresample=1e6))) Gender Female Male
Income Job Satisfaction
<5000 V_D 1 1
L_S 3 1
M_S 11 2
V_S 2 1
5000-15000 V_D 2 0
L_S 3 3
M_S 17 5
V_S 3 1
15000-25000 V_D 0 0
L_S 1 0
M_S 8 7
V_S 5 3
>25000 V_D 0 0
L_S 2 1
M_S 4 9
V_S 2 6
Approximative Linear-by-Linear Association Test
data: JobSatisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01106
alternative hypothesis: two.sided\[\Diamond\]
Efeito da codificação do parâmetro scores das funções
coin::lbl_test e coin::cmh_test
No script abaixo é mostrado o efeito da codificação do parâmetro
scores das funções coin::lbl_test e
coin::cmh_test. Note que os valores relativos atribuídos
aos níveis dos fatores ordinais afetam o valor p, pois o que
importa são as distâncias relativas entre os níveis. Os valores
absolutos atribuídos aos níveis dos fatores ordinais não afetam o valor
p. Portanto, se o fator ordinal é dicotômico, há apenas uma
distância entre estes níveis, não importando assim os valores absolutos
distintos atribuídos aos dois níveis do fator ordinal. Desta forma, o
valor p é o mesmo se os dois fatores são considerados como
nominais ou ordinais (com qualquer codificação em
scores).
Gender Female Male
Income Job.Satisfaction
<5000 Very Dissatisfied 1 1
A Little Satisfied 3 1
Moderately Satisfied 11 2
Very Satisfied 2 1
5000-15000 Very Dissatisfied 2 0
A Little Satisfied 3 3
Moderately Satisfied 17 5
Very Satisfied 3 1
15000-25000 Very Dissatisfied 0 0
A Little Satisfied 1 0
Moderately Satisfied 8 7
Very Satisfied 5 3
>25000 Very Dissatisfied 0 0
A Little Satisfied 2 1
Moderately Satisfied 4 9
Very Satisfied 2 6 'table' int [1:4, 1:4, 1:2] 1 2 0 0 3 3 1 2 11 17 ...
- attr(*, "dimnames")=List of 3
..$ Income : chr [1:4] "<5000" "5000-15000" "15000-25000" ">25000"
..$ Job.Satisfaction: chr [1:4] "Very Dissatisfied" "A Little Satisfied" "Moderately Satisfied" "Very Satisfied"
..$ Gender : chr [1:2] "Female" "Male"# Exemplo do Help da funcao coin::lbl_test
# Asymptotic linear-by-linear association test (Agresti, p. 297)
# Note: coin::lbl_test: 'Job.Satisfaction' and 'Income' sempre são ordinais
cat("\ncoin::lbl_test\n")
coin::lbl_test
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.4812, p-value = 0.01309
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.4406, p-value = 0.01466
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sidedprint(lt <- coin::lbl_test(TC,
scores=list("Job.Satisfaction"=c(0.1,0.2,0.3,0.4),
"Income"=c(1,2,3,4))))
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided# Note: coin::cmh_test: 'Job.Satisfaction' and 'Income' podem ser ordinais
cat("\ncoin::cmh_test\n")
coin::cmh_test
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.4812, p-value = 0.01309
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.4406, p-value = 0.01466
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided
Asymptotic Generalized Cochran-Mantel-Haenszel Test
data: Job.Satisfaction (ordered) by
Income (<5000, 5000-15000, 15000-25000, >25000)
stratified by Gender
chi-squared = 9.2259, df = 3, p-value = 0.02643
Asymptotic Generalized Cochran-Mantel-Haenszel Test
data: Job.Satisfaction by
Income (<5000, 5000-15000, 15000-25000, >25000)
stratified by Gender
chi-squared = 10.2, df = 9, p-value = 0.3345print(lt <- coin::cmh_test(TC,
scores=list("Job.Satisfaction"=c(0.1,0.2,0.3,0.4),
"Income"=c(1,2,3,4))))
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sided
Asymptotic Linear-by-Linear Association Test
data: Job.Satisfaction (ordered) by
Income (<5000 < 5000-15000 < 15000-25000 < >25000)
stratified by Gender
Z = 2.5736, p-value = 0.01006
alternative hypothesis: two.sidedAnálise loglinear: teste qui-quadrado e modelo paramétrico
Conforme Field (2012, p. 833):
“OK, tudo isso é muito bom, mas o título desta seção realmente insinuava que eu mostraria como o teste qui-quadrado pode ser conceituado como um modelo linear. Bem, basicamente, o teste qui-quadrado verifica se duas variáveis são independentes; portanto, não tem interesse no efeito combinado das duas variáveis, apenas no seu efeito único. Assim, podemos conceituar o qui-quadrado de maneira muito semelhante ao modelo saturado, exceto que não incluímos o termo de interação.”
Streiner & Lin, 1998
Por mais de quatro décadas, a General Social Survey (GSS){{target=“_blank”}} estudou a crescente complexidade da sociedade americana. É a única pesquisa de probabilidade total com entrevista pessoal projetado para monitorar mudanças em ambas as características sociais e atitudes que estão sendo conduzidas atualmente nos Estados Unidos.
Como exemplo usaremos alguns dados obtidos em 2000. O arquivo está em
formato *.sav (do SPSS). O SPSS guarda a codificação e o
rótulo correspondente em cada entrada (demo_GSS2000_01_data.R):
# For more than four decades, the General Social Survey (GSS)
# has studied the growing complexity of American society.
# It is the only full-probability, personal-interview survey
# designed to monitor changes in both social characteristics
# and attitudes currently being conducted in the United States.
# https://gss.norc.org/
dt_gss <- haven::read_sav("GSS2000.sav")
print(head(dt_gss))
cat("\n\nCodificacao:")
for (c.aux in 2:ncol(dt_gss))
{
cat("\n")
print(unique(dt_gss[c.aux]))
}# A tibble: 6 × 5
id renda feliz casado sexo
<dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl>
1 1 2 [ClasseB] 1 [Feliz] 2 [Descasado] 2 [Mulher]
2 2 3 [ClasseA] 1 [Feliz] 1 [Casado] 1 [Homem]
3 3 2 [ClasseB] NA 2 [Descasado] 2 [Mulher]
4 4 3 [ClasseA] NA 2 [Descasado] 2 [Mulher]
5 5 2 [ClasseB] NA 2 [Descasado] 1 [Homem]
6 6 3 [ClasseA] 1 [Feliz] 2 [Descasado] 1 [Homem]
Codificacao:
# A tibble: 4 × 1
renda
<dbl+lbl>
1 2 [ClasseB]
2 3 [ClasseA]
3 1 [ClasseC]
4 NA
# A tibble: 3 × 1
feliz
<dbl+lbl>
1 1 [Feliz]
2 NA
3 2 [Infeliz]
# A tibble: 3 × 1
casado
<dbl+lbl>
1 2 [Descasado]
2 1 [Casado]
3 3 [Nuncasado]
# A tibble: 2 × 1
sexo
<dbl+lbl>
1 2 [Mulher]
2 1 [Homem] Neste exemplo podemos começar verificando se existe associação entre sexo (categorias: Homem e Mulher) e felicidade (categorias: Feliz e Infeliz).
Para comparação, executamos o teste qui-quadrado para a tabela 2×2
(demo_GSS2000_00_x2.R),
source("eiras.show.MCSTAR.R")
dt_gss <- haven::read_sav("GSS2000.sav")
print(head(dt_gss[,c("sexo","feliz")]))
cat("\nTeste de associacao:\n")
alpha <- 0.05
tb <- xtabs(~sexo+feliz,data=dt_gss)
out_x2 <- chisq.test(tb,correct=FALSE)
print(out_x2)
cat("\nobservado:\n")
print(out_x2$observed)
cat("\nesperado sob H0:\n")
print(out_x2$expected)
cat("\ncodificacao: Sexo (1=Homem, 2=Mulher); Feliz (1=Sim, 2=Nao)\n")
cat("\n------------------------------------------\n")
cat("Analise post hoc:")
cat("\n------------------------------------------\n")
n <- sum(out_x2$observed)
nL <- nrow(tb)
nC <- ncol(tb)
df <- (nL-1)*(nC-1)
X2 <- out_x2$statistic # estatistica de teste qui-quadrado
cat("Graus de liberdade: ", df, "\n", sep="")
cat("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
STAR <- out_x2$stdres # STandardized Adjusted Residual (STAR)
MCSTAR <- STAR/(sqrt((1-1/nL)*(1-1/nC))) # Moment-correct (MCSTAR)
show.MCSTAR(MCSTAR=MCSTAR,alpha=alfa,df=df)# A tibble: 6 × 2
sexo feliz
<dbl+lbl> <dbl+lbl>
1 2 [Mulher] 1 [Feliz]
2 1 [Homem] 1 [Feliz]
3 2 [Mulher] NA
4 2 [Mulher] NA
5 1 [Homem] NA
6 1 [Homem] 1 [Feliz]
Teste de associacao:
Pearson's Chi-squared test
data: tb
X-squared = 10.34, df = 1, p-value = 0.001302
observado:
feliz
sexo 1 2
1 588 61
2 611 109
esperado sob H0:
feliz
sexo 1 2
1 568.4083 80.59167
2 630.5917 89.40833
codificacao: Sexo (1=Homem, 2=Mulher); Feliz (1=Sim, 2=Nao)
------------------------------------------
Analise post hoc:
------------------------------------------
Graus de liberdade: 1
Heuristica de significancia (rej. H0 se X^2/gl > 2): 10.3397
Residuos ajustados standardizados corrigidos por momento (MCSTAR):
|MCSTAR critico| (alphaBonferroni=5%/1) = 1.959964
1 2
1 #6.43 -6.43
2 -6.43 #6.43
O gráfico de mosaico representa este tipo de situação (demo_GSS2000_06_mosaic.R):
# mosaic plot
dt_gss$sexoTxt <- plyr::mapvalues(as.character(dt_gss$sexo),
from=c("1","2"),
to=c("Homem","Mulher"))
dt_gss$felizTxt <- plyr::mapvalues(as.character(dt_gss$feliz),
from=c("1","2"),
to=c("Feliz","Infeliz"))
tbtxt <- xtabs(~sexoTxt+felizTxt,data=dt_gss)
mosaicplot(tbtxt,shade=TRUE,main="GSS: sexo e felicidade")
cat("\nobservado:\n")
print(out_x2$observed)
cat("\nesperado sob H0:\n")
print(out_x2$expected)
cat("\nresiduos:\n")
print(out_x2$residuals)
cat("\ncodificacao: Sexo (1=Homem, 2=Mulher); Feliz (1=Sim, 2=Nao)\n")
observado:
feliz
sexo 1 2
1 588 61
2 611 109
esperado sob H0:
feliz
sexo 1 2
1 568.4083 80.59167
2 630.5917 89.40833
residuos:
feliz
sexo 1 2
1 0.8217530 -2.1823602
2 -0.7801847 2.0719656
codificacao: Sexo (1=Homem, 2=Mulher); Feliz (1=Sim, 2=Nao)Este gráfico é obtido com mosaicplot. Para ser
interpretado:
- o tamanho de cada retângulo representa o número de ocorrências.
- a cor e o limite de cada retângulo nos informa sobre os resíduos
padronizados.
- um retângulo com borda sólida representa resíduo positivo; borda tracejada para negativo.
- um retângulo azul representa resíduo maior que 2 e vermelho para menor que -2 (heurística de significância estatística).
Na forma de uma regressão linear múltipla podemos incluir o efeito de interação (\(\text{sexo:feliz}\)) tornando o modelo saturado:
\[ \ln(O_{ij}) = { \beta_0 + \beta_1 \;\text{sexo}_i + \beta_2 \;\text{feliz}_j + \beta_3 \;\text{sexo:feliz}_{ij} } \]
em que \(O_{ij}\) é a frequência absoluta observada em uma casela da tabela de contingência.
Note que usaremos o logaritmo das frequências observadas (a análise adotada aqui é loglinear).
|
Usamos no código R, abaixo, a função
Esta forma de escrever é economica; em sua forma mais explícita é o mesmo que
em que Esta interação representa a dependência entre as duas variaveis nominais. Há, ainda, outra possibilidade (com resultado idêntico):
|
O modelo fica não saturado quando não consideramos o efeito de interação:
\[\ln\left(O_{ij}\right) = { \beta_0 + \beta_1\; \text{sexo}_i + \beta_2 \;\text{feliz}_j + \varepsilon_{ij} }\]
Veremos adiante como o modelo saturado se comporta, pois ele será o modelo referência para os testes estatísticos do efeito de interação entre sexo e felicidade que investigamos neste exemplo. A ideia é aninhar o modelo não saturado no modelo saturado, de forma que, a partir da diferença entre os dois, restará o que é relativo ao efeito de interação entre sexo e felicidade sob investigação (i.e., a hipótese nula é a ausência do efeito de interação destas duas variáveis), o que corresponde à hipótese nula do teste qui-quadrado de Pearson.
Observe que a felicidade foi codificada como:
- 1 … Feliz
- 2 … Infeliz
e sexo como:
- 1 … Homem
- 2 … Mulher
O efeito de interação é o produto destas variáveis nominais. Criamos
uma nova coluna com (demo_GSS2000_02_interacao.R):
dt_gss$interacao <- dt_gss$sexo * dt_gss$feliz
print(head(dt_gss[,c("sexo","feliz","interacao")]))# A tibble: 6 × 3
sexo feliz interacao
<dbl+lbl> <dbl+lbl> <dbl>
1 2 [Mulher] 1 [Feliz] 2
2 1 [Homem] 1 [Feliz] 1
3 2 [Mulher] NA NA
4 2 [Mulher] NA NA
5 1 [Homem] NA NA
6 1 [Homem] 1 [Feliz] 1Para a execução da regressão linear múltipla que será feita adiante,
precisamos criar uma tabela de contingência e, com base nela, repetir a
contagem das combinações possíveis de sexo e
feliz como a variável dependente, e, então, calcular seu
valor em logaritmo, o que se obtém com (demo_GSS2000_03_vd.R):
cat("\nTabela de contingencia\n")
tb <- xtabs(~sexo+feliz,data=dt_gss)
print(tb)
# contagem dos valores observados (VD)
dt_gss$Observed <- NA
rnames <- rownames(tb) # Training
cnames <- colnames(tb) # Dance
for (r.aux in 1:length(rnames))
{
for (c.aux in 1:length(cnames))
{
dt_gss$Observed[as.character(dt_gss$sexo)==rnames[r.aux]&
as.character(dt_gss$feliz)==cnames[c.aux]] <- tb[r.aux,c.aux]
}
}
dt_gss$LnObserved <- log(dt_gss$Observed)
print(head(dt_gss[,c("sexo","feliz","interacao","Observed","LnObserved")]))
Tabela de contingencia
feliz
sexo 1 2
1 588 61
2 611 109
# A tibble: 6 × 5
sexo feliz interacao Observed LnObserved
<dbl+lbl> <dbl+lbl> <dbl> <int> <dbl>
1 2 [Mulher] 1 [Feliz] 2 611 6.42
2 1 [Homem] 1 [Feliz] 1 588 6.38
3 2 [Mulher] NA NA NA NA
4 2 [Mulher] NA NA NA NA
5 1 [Homem] NA NA NA NA
6 1 [Homem] 1 [Feliz] 1 588 6.38Com este procedimento, o valor 588 e seu logaritmo 6.38 são repetidos 588 vezes, em todas as colunas nas quais sexo tem o valor 1 e feliz tem o valor 1, o valor 61 e seu logaritmo 4.11 são repetidos 61 vezes, em todas as colunas nas quais sexo tem o valor 1 e feliz tem o valor 2, o valor 611 e seu logaritmo 6.42 são repetidos 611 vezes, em todas as colunas nas quais sexo tem o valor 2 e feliz tem o valor 1 e o valor 109 e seu logaritmo 4.69 são repetidos 109 vezes, em todas as colunas nas quais sexo tem o valor 2 e feliz tem o valor 2.
O modelo saturado é computado com a função lm desta
maneira (demo_GSS2000_04_sat.R):
cat("\nRegressao linear com interacao (modelo saturado):\n")
out_s <- lm(LnObserved~sexo*feliz,data=dt_gss)
print(summary(out_s))
b0 <- as.numeric(out_s$coefficients[1])
b1 <- as.numeric(out_s$coefficients[2])
b2 <- as.numeric(out_s$coefficients[3])
b3 <- as.numeric(out_s$coefficients[4])
dt_gss$SaturatedHat <- b0 + b1*dt_gss$sexo + b2*dt_gss$feliz + b3*dt_gss$interacao
dt_gss$SaturatedRes <- dt_gss$LnObserved-dt_gss$SaturatedHat
print(head(dt_gss[,c("sexo","feliz","LnObserved","SaturatedHat","SaturatedRes")]))
cat("\nResíduos encontrados na coluna SaturatedRes:\n")
cat("\t",unique(round(dt_gss$SaturatedRes,5)))
Regressao linear com interacao (modelo saturado):
Call:
lm(formula = LnObserved ~ sexo * feliz, data = dt_gss)
Residuals:
Min 1Q Median 3Q Max
-1.609e-13 -4.300e-17 0.000e+00 7.700e-17 1.859e-13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.146e+00 2.268e-15 4.032e+15 <2e-16 ***
sexo -5.037e-01 1.359e-15 -3.708e+14 <2e-16 ***
feliz -2.808e+00 1.986e-15 -1.414e+15 <2e-16 ***
sexo:feliz 5.421e-01 1.171e-15 4.630e+14 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.886e-15 on 1365 degrees of freedom
(1396 observations deleted due to missingness)
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.927e+30 on 3 and 1365 DF, p-value: < 2.2e-16
# A tibble: 6 × 5
sexo feliz LnObserved SaturatedHat SaturatedRes
<dbl+lbl> <dbl+lbl> <dbl> <dbl> <dbl>
1 2 [Mulher] 1 [Feliz] 6.42 6.42 8.88e-15
2 1 [Homem] 1 [Feliz] 6.38 6.38 -2.66e-15
3 2 [Mulher] NA NA NA NA
4 2 [Mulher] NA NA NA NA
5 1 [Homem] NA NA NA NA
6 1 [Homem] 1 [Feliz] 6.38 6.38 -2.66e-15
Resíduos encontrados na coluna SaturatedRes:
0 NAInicia-se com o modelo saturado, que é um modelo completo, consegue-se o ajuste perfeito.
Note que a regressão linear múltipla forneceu os coeficientes:
(Intercept) sexo feliz sexo:feliz
9.146314 -0.503734 -2.807957 0.542104 e, portanto, a coluna SaturatedHat foi computada com a
equação
\(\ln(\hat{O_{ij}^{s}}) =\) 9.14631 -0.50373 \(\text{sexo}_i\) -2.80796 \(\text{feliz}_j +\) 0.5421 \(\text{interacao}_{ij}\)
mostrando que o modelo saturado faz o ajuste perfeito aos dados,
tanto que os resíduos calculados na coluna SaturatedRes
(dada pela diferença entre LnObserved e
SaturatedHat) são praticamente nulos.
Observe, ainda, na saída da regressão, que \(R^2\) é igual a \(1\) porque \(\ln\left(O_{ij}\right) =\ln\left(\hat{O}_{ij}^{s}\right)\).
|
A forma como calculamos a coluna Existe uma forma mais automatizada, usando a função
Uma solução possível para usar Outra alternativa é preparar o data frame deixando apenas as
colunas de interesse (neste exemplo, descartando as colunas
O que fizemos até aqui é equivalentemente obtido com: |
O modelo não saturado é computado com (demo_GSS2000_05_naosat.R):
cat("\nRegressao linear com interacao (modelo nao saturado):\n")
out_ns <- lm(LnObserved~sexo+feliz,data=dt_gss)
print(summary(out_ns))
b0 <- as.numeric(out_ns$coefficients[1])
b1 <- as.numeric(out_ns$coefficients[2])
b2 <- as.numeric(out_ns$coefficients[3])
dt_gss$NotSaturatedHat <- b0 + b1*dt_gss$sexo + b2*dt_gss$feliz
dt_gss$NotSaturatedRes <- dt_gss$LnObserved-dt_gss$NotSaturatedHat
print(head(dt_gss[,c("sexo","feliz","LnObserved","NotSaturatedHat","NotSaturatedRes")]))
cat("\nResíduos encontrados na coluna NotSaturatedRes:\n")
cat("\t",unique(round(dt_gss$NotSaturatedRes,5)))
Regressao linear com interacao (modelo nao saturado):
Call:
lm(formula = LnObserved ~ sexo + feliz, data = dt_gss)
Residuals:
Min 1Q Median 3Q Max
-0.3075 -0.0307 0.0319 0.0319 0.1721
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.170376 0.010504 777.84 <2e-16 ***
sexo 0.100961 0.004687 21.54 <2e-16 ***
feliz -1.926505 0.007097 -271.45 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.08627 on 1366 degrees of freedom
Multiple R-squared: 0.9818, Adjusted R-squared: 0.9818
F-statistic: 3.684e+04 on 2 and 1366 DF, p-value: < 2.2e-16
# A tibble: 6 × 5
sexo feliz LnObserved NotSaturatedHat NotSaturatedRes
<dbl+lbl> <dbl+lbl> <dbl> <dbl> <dbl>
1 2 [Mulher] 1 [Feliz] 6.42 6.45 -0.0307
2 1 [Homem] 1 [Feliz] 6.38 6.34 0.0319
3 1 [Homem] 1 [Feliz] 6.38 6.34 0.0319
4 1 [Homem] 2 [Infeliz] 4.11 4.42 -0.307
5 2 [Mulher] 1 [Feliz] 6.42 6.45 -0.0307
6 2 [Mulher] 2 [Infeliz] 4.69 4.52 0.172
Resíduos encontrados na coluna NotSaturatedRes:
-0.03069 0.0319 -0.30745 0.17206Observe no código a sintaxe do modelo, que agora utiliza
lm(LnObserved~sexo + feliz,data=dt_gss)
fazendo com que a regressão múltipla tenha uma linha a menos e o
ajustamento do modelo não seja mais perfeito. Como há quatro valores
diferentes na coluna LnObserved, quatro valores diferentes
aparecem nos resíduos.
Para conferência, vamos verificar uma ocorrência de cada combinação existente de sexo e felicidade:
for (s in sort(unique(dt_gss$sexo)))
{
for (f in sort(unique(dt_gss$feliz)))
{
dt_tmp <- na.omit(dt_gss[dt_gss$sexo==s & dt_gss$feliz==f,
c("sexo","feliz",
"LnObserved",
"SaturatedHat","NotSaturatedHat")])
prmatrix(dt_tmp[1,],rowlab = "")
}
} sexo feliz LnObserved SaturatedHat NotSaturatedHat
1 1 6.376727 6.376727 6.344831
sexo feliz LnObserved SaturatedHat NotSaturatedHat
1 2 4.110874 4.110874 4.418327
sexo feliz LnObserved SaturatedHat NotSaturatedHat
2 1 6.415097 6.415097 6.445792
sexo feliz LnObserved SaturatedHat NotSaturatedHat
2 2 4.691348 4.691348 4.519287A coluna NotSaturatedHat foi computada a partir dos
coeficientes obtidos pela função lm com o modelo não
saturado com a equação:
\(\ln\left(\hat{O}_{ij}^{ns}\right) =\) 8.17038\(\;+\;\) 0.10096 \(\text{sexo}_i\;\) -1.9265 \(\text{feliz}_j\)
O pacote MASS tem a função MASS::loglm que resolve e
calcula o teste estatístico para estas regressões loglineares
trabalhosas que fizemos acima. O modelo saturado é escrito como (demo_GSS2000_07_mass.R):
tb <- xtabs(~sexo+feliz,data=dt_gss)
cat("\n\tmodelo saturado:\n")
out_mass_s <- MASS::loglm(~sexo+feliz+sexo:feliz,data=tb,fit=TRUE)
print(out_mass_s)
modelo saturado:
Call:
MASS::loglm(formula = ~sexo + feliz + sexo:feliz, data = tb,
fit = TRUE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0 0 1
Pearson 0 0 1Observe a linha “Pearson”, a qual (correspondendo ao ajuste perfeito
do modelo saturado que obtivemos com lm) verifica se todas
as inclinações (\(\beta_i\), exceto
\(i=0\)) são nulas.
O modelo não saturado é feito com (demo_GSS2000_08_mass.R):
cat("\n\tmodelo nao saturado:\n")
out_mass_ns <- MASS::loglm(~sexo+feliz,data=tb,fit=TRUE)
print(out_mass_ns)
modelo nao saturado:
Call:
MASS::loglm(formula = ~sexo + feliz, data = tb, fit = TRUE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 10.49635 1 0.001196107
Pearson 10.33970 1 0.001301989Caso subtraíssemos este resultado do modelo saturado (análise hieráquica loglinear), o qual mostrou \(X^2=0\) e \(df=0\), a diferença são os próprios valores do modelo não saturado neste caso simples. Compare com o resultado que vimos, acima, com o teste qui-quadrado de Pearson:
Pearson's Chi-squared test
data: tb
X-squared = 10.34, df = 1, p-value = 0.001302Mais formalmente, a comparação entre os dois modelos (subtrair o
modelo saturado do modelo não saturado) utiliza a função nativa
anova que, executada desta forma (demo_GSS2000_09_mass.R):
cat("\nTeste do modelo nao saturado aninhado no modelo saturado:\n")
print(anova(out_mass_s,out_mass_ns))
Teste do modelo nao saturado aninhado no modelo saturado:
LR tests for hierarchical log-linear models
Model 1:
~sexo + feliz
Model 2:
~sexo + feliz + sexo:feliz
Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
Model 1 10.49635 1
Model 2 0.00000 0 10.49635 1 0.0012
Saturated 0.00000 0 0.00000 0 1.0000mostra a diferença entre os modelos saturado e não saturado (exceto que usa os valores de usando LR em vez de Pearson, próximos neste caso; LR é mais robusto para amostras pequenas).
Portanto, o teste qui-quadrado de Pearson pode ser expresso na forma de um teste paramétrico.
A vantagem desta abordagem paramétrica é que podemos adicionar mais variáveis nominais dicotômicas e politômicas.
Descrição do modelo
- Análise de associações entre 3 ou mais variáveis nominais: extensão do teste qui-quadrado de Pearson
- Análise hierárquica de dependências (interações)
- Modelo linear geral (GLM):
- VD: logaritmo da frequência observada da casela = \(\ln(O)\)
- VI: 3 ou mais variáveis nominais
- UE: observações independentes
- Efeitos: principais e interações
- Suposições: as mesmas do teste qui-quadrado de Pearson
- Modelo saturado (inicial) tem ajuste perfeito: \(ln(O)\) é perfeitamente explicado por todos os efeitos principais e de interação de todas as ordens; \(\text{#efeitos}=2^{\text{#VI}}-1\).
Análise hierárquica de dependências
- Efeitos de interação inferior só podem ser testados se os de nível superior não são significativos: a interação de ordem mais alta não deve ser significativa
- Interpretar apenas os efeitos de interação de ordem mais alta significativos
- Classe geradora: conjunto de efeitos de interação de maior ordem candidatos a serem removidos no próximo passo para garantir que o modelo seja hierárquico; no próximo passo os efeitos da classe geradora com valor p maior que o nível de significância estabelecido pelo pesquisador são removidos do modelo; o modelo resultante é adequado se o valor p associado à razão de verossimilhança é maior que o nível de significância estabelecido pelo pesquisador (e.g., 5%)
Modelo saturado
Teste LR de qualidade de ajuste: \(X^2(0)=0,\; p=1\) (ajuste perfeito!)
Lógica da análise loglinear
Quando nossos dados são categóricos e incluímos todos os termos disponíveis (efeitos principais e interações), não erramos: nossos previsores podem perfeitamente prever nossa saída (os valores esperados).
Assim, se começarmos com o modelo mais complexo possível, não erramos.
A tarefa da análise loglinear é tentar ajustar um modelo mais simples aos dados sem uma perda substancial do poder preditivo.
Portanto, a análise loglinear tipicamente trabalha num princípio de eliminação de trás para frente.
Assim, começamos com o modelo saturado e, depois, removemos um previsor do modelo e usando esse novo modelo prevemos nossos dados (calculando as frequências esperadas como no teste do qui-quadrado).
Então, verificamos quão bem o modelo se ajusta aos dados (i.e., se as frequências esperadas estão próximas das frequências observadas).
Se a aderência do novo modelo não for muito diferente do modelo complexo, abandonamos o modelo complexo em favor do novo.
Colocando de outra maneira, presumimos que o termo que removemos não estava tendo um impacto significativo na habilidade do nosso modelo para prever os dados observados.
Entretanto, a análise não apenas remove os termos ao acaso, ela faz isso de forma hierárquica. Assim, começamos com o modelo saturado e, depois, removemos a interação da ordem mais alta e avaliamos o efeito que ela tinha.
Se a remoção desse termo de interação não afetou o modelo, ele obviamente não estava sendo muito útil.
Portanto, nos livramos dele e continuamos removendo outras interações de ordem mais baixa agora.
Se a remoção dessas interações não tem impacto, continuamos com os efeitos principais até encontrarmos um efeito que realmente afete a adequação do modelo se for removido.
Análise de concordância
Concordância é diferente de associação.
Há testes estatísticos de concordância entre dois métodos que avaliam dicotomicamente as mesmas unidades observacionais.
Portanto, testes baseados da estatística qui-quadrado são testes de associação.
Silveira & Siqueira (2022) avaliaram e compararam 11 medidas estatísticas de (supostamente!) concordância. Mostramos que das 11 medidas apenas duas são verdadeiramente medidas de concordância: G de Holley-Guilford (1964) e AC1 de Gwet (2008). Portanto, \(\kappa\) de Cohen não é uma medida de concordância; é uma medida de associação.
O resumo do artigo é o seguinte:
“Avaliamos vários coeficientes de concordância aplicados em tabelas de contingência 2×2, que são comumente utilizados em pesquisas devido à dicotomização. Aqui, não apenas estudamos alguns estimadores específicos, mas também desenvolvemos um método geral para o estudo de qualquer candidato a estimador para ser uma medida de concordância. Este método foi desenvolvido em códigos R e está disponível para os pesquisadores. Testamos este método verificando o desempenho de vários estimadores tradicionais em todas as configurações possíveis com tamanhos variando de 1 a 68 (total de 1.028.789 tabelas). Kappa de Cohen mostrou um comportamento prejudicado similar a r de Pearson, Q de Yule e Y de Yule. \(\pi\) de Scott e B de Shankar e Bangdiwala parecem avaliar melhor as situações de discordância do que de concordância entre avaliadores. Alfa de Krippendorff emula, sem qualquer vantagem, \(\pi\) de Scott em casos com variáveis nominais e dois avaliadores. F1 de Dice e qui-quadrado de McNemar avaliam de forma incompleta as informações da tabela de contingência, mostrando o pior desempenho entre todos. Concluímos que kappa de Cohen é uma medida de associação e qui-quadrado de McNemar não avalia nem associação nem concordância; os únicos dois estimadores de concordância autênticos são G de Holley e Guilford e AC1 de Gwet. Estes dois últimos estimadores também mostraram o melhor desempenho ao longo da gama de tamanhos de tabelas e devem ser considerados como as primeiras escolhas para medição de concordância em tabelas de contingência 2×2.”
AC1 de Gwet
Agreement Coefficient do tipo 1, AC1, assume que as duas variáveis são nominais dicotômicas (tabela 2 \(\times\) 2).
Em Sleigh et al. (1982), na análise de 315 amostras de fezes, verifica-se a concordância para a detecção Positiva ou Negativa de ovos de Schistosoma mansoni usando os métodos Bell e Kato-Katz.
A hipótese nula é
\[ \begin{cases} H_0: \text{ausência de concordância}\\ H_1: \text{presença de concordância} \end{cases}\\ \alpha=5\% \]
O exemplo está implementado em TC2x2_Gwet_concord.R:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.matrix2dataframe.R")
# GWET, KL (2008) Computing inter-rater reliability and its variance
# in the presence of high agreement. British Journal of Mathematical and
# Statistical Psychology 61: 29-48.
Tabela <- ("
BellxKK P N
P 184 54
N 14 63
")
cat(Tabela)
# Foram analisadas 315 amostras usando os métodos Bell e
# Kato-Katz para detecção ovos de Schistosoma mansoni nas fezes.
# Sleigh A, Hoff R, Mott K, Barreto M, de Paiva TM, Pedrosa Jde S,
# Sherlock I. (1982)
# Comparison of filtration staining (Bell) and thick smear (Kato)
# for the detection of quantitation of Schistosoma mansoni eggs
# in faeces. Transactions of the Royal Society of
# Tropical Medicine and Hygiene 76(3):403-6.
# doi: 10.1016/0035-9203(82)90201-2. PMID: 7112662.
TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE,
row.names=1))
print(gplots::balloonplot(t(as.table(TC)),
main ="Detecção ovos de Schistosoma mansoni",
xlab ="Kato-Katz", ylab="Bell",
label = TRUE,
show.margins = TRUE))
# AC1 de Gwet (2008)
cat("\nAC1 de Gwet",sep="")
AC1 <- irrCAC::gwet.ac1.table(ratings=TC)
AC1.f <- AC1$coeff.val
AC1.ic <- AC1$coeff.ci
AC1.p <- AC1$coeff.pval
cat("\n\tAC1 = ",round(AC1.f,2), ", IC95 = ",AC1.ic,", p = ",AC1.p,"\n",sep="")
cat("\nDuas variáveis nominais:\n")
alfa <- 0.05
seed <- 123
B <- 1e4
print(confintr::ci_cramersv(x=as.table(TC),
probs=c(alfa/2, 1-alfa/2),
type="bootstrap",
R=B,
seed=seed),
digits=4)
BellxKK P N
P 184 54
N 14 63
NULL
AC1 de Gwet
AC1 = 0.62, IC95 = (0.536,0.712), p = 0e+00
Duas variáveis nominais:
Two-sided 95% bootstrap confidence interval for the population Cramer's
V based on 10000 bootstrap replications and the bca method
Sample estimate: 0.5259
Confidence interval:
2.5% 97.5%
0.4231 0.6156 Neste exemplo, \(H_0\) é rejeitada para \(\alpha=5\%\) e, portanto, conclui-se que os dois métodos para a detecção de ovos podem ser concordantes.
Uma forma de expressar este coeficiente (adotando a notação para tabelas 2×2 que usamos acima) é:
\[ \text{AC}_1 = { {a^2+d^2- {\dfrac{(b+c)^2}{2}} }\over{ a^2 + d^2 + {\dfrac{(b+c)^2}{2}}} + (a+d)(b+c)} \]
na qual podemos ver que existe tensão entre as diagonais, pois aparece a soma dos quadrados de \(a\) e \(d\) da qual é subtraída a soma de \(b\) e \(c\) ao quadrado. O denominador da fórmula faz com que esta medida de concordância, convenientemente, fique entre 0 e 1. AC1 é uma medida de tamanho de efeito.
Fórmulas adaptadas de Gwet (2008)
As fórmulas a seguir estão implementadas na função
eiras2x2::AC1.
O pacote eiras2x2 está disponível em Harvard Dataverse.
Sample size: \(n = a + b + c + d\)
Probabilities:
\[ p_{11} = \dfrac{a}{n}, \quad p_{12} = \dfrac{b}{n}, \quad p_{21} = \dfrac{c}{n}, \quad p_{22} = \dfrac{d}{n} \]
- Overall agreement probability:
\[ p_a = p_{11} + p_{22} = \dfrac{a + d}{n} \]
- Marginal probabilities:
\[ p_{R1} = \dfrac{a + b}{n}, \quad p_{R2} = \dfrac{c + d}{n}, \quad p_{C1} = \dfrac{a + c}{n}, \quad p_{C2} = \dfrac{b + d}{n} \]
- Average marginal probabilities:
\[ \bar{p}_{1} =\dfrac{p_{R1}+p_{C1}}{2} = \dfrac{2a + b + c}{2n}, \quad \bar{p}_{2} =\dfrac{p_{R2}+p_{C2}}{2} = \dfrac{b + c + 2d}{2n} \]
- Chance-agreement probabilities:
\[ p_{c}^{\kappa} = p_{R1}p_{C1}+p_{R2}p_{C2} = \dfrac{(a + b)(a + c) + (c + d)(b + d)}{n^2} \]
\[ p_{c}^{\text{AC}_{1}} = \bar{p}_{1}(1 - \bar{p}_{1}) + \bar{p}_{2}(1 - \bar{p}_{2}) \]
\[ p_{c}^{G}=\dfrac{1}{2} \]
- Agreement measures:
\[ \kappa = \dfrac{p_{a}-p_{c}^{\kappa}}{1-p_{c}^{\kappa}}=\dfrac{2(ad - bc)}{(a + c)(c + d) + (a + b)(b + d)} \]
\[ \text{AC}_{1} = \dfrac{p_{a}-p_{c}^{\text{AC}_{1}}}{1-p_{c}^{\text{AC}_{1}}} = \dfrac{a^2 + d^2 - \dfrac{(b + c)^2}{2}}{a^2 + d^2 + \dfrac{(b + c)^2}{2} + (a + d)(b + c)} \]
\[ G = \dfrac{p_{a}-p_{c}^{G}}{1-p_{c}^{G}}=\dfrac{(a + d) - (b + c)}{n} \]
- Standard error:
\[ \text{se}_{\kappa}^2 = \dfrac{1}{n(1 - p_{c}^{\kappa})^2} \left(p_a(1 - p_a) \\ - 4(1 - \kappa)(p_{11}\bar{p}_{1} + p_{22}\bar{p}_{2} - p_a p_{c}^{\kappa}) \\ + 4(1 - \kappa)^2\left(p_{11}\left(\dfrac{p_{R1}+p_{C1}}{2}\right)^2 + p_{12}\left(\dfrac{p_{R1}+p_{C2}}{2}\right)^2 \\ + p_{21}\left(\dfrac{p_{R2}+p_{C1}}{2}\right)^2 + p_{22}\left(\dfrac{p_{R2}+p_{C2}}{2}\right)^2 - (p_{c}^{\kappa})^2\right)\right) \]
\[ \text{se}_{\text{AC}_{1}}^2 = \dfrac{1}{n\left(1 - p_{c}^{\text{AC}_{1}}\right)^2} \left(p_a(1 - p_a) \\ - 4\left(1 - \text{AC}_{1}\right)\left(p_{11}\left(1 - \bar{p}_{1}\right) + p_{22}\left(1 - \bar{p}_{2}\right) - p_a p_{c}^{\text{AC}_{1}}\right) \\ + 4(1 - \text{AC}_{1})^2\left(p_{11}\left(1 - \dfrac{\bar{p}_{1}+\bar{p}_{1}}{2}\right)^2 + p_{12}\left(1 - \dfrac{\bar{p}_{1}+\bar{p}_{2}}{2}\right)^2 \\ + p_{21}\left(1 - \dfrac{\bar{p}_{2}+\bar{p}_{1}}{2}\right)^2 + p_{22}\left(1 - \dfrac{\bar{p}_{2}+\bar{p}_{2}}{2}\right)^2 - \left(p_{c}^{\text{AC}_{1}}\right)^2\right)\right) \]
\[ \text{se}_{G}^2 = \dfrac{4}{n}p_a(1 - p_a) = \dfrac{4(a + d)(b + c)}{n^3} \]
Cálculo dos Testes t para Medidas de Concordância
Para conduzir testes t para cada medida de concordância (kappa, AC1 e G), precisamos calcular os valores t e seus respectivos valores p. O valor t é baseado na medida de concordância e seu erro-padrão. Podemos então usar uma distribuição t para calcular os valores p.
\[ t = \dfrac{\text{medida de concordância}}{\text{erro-padrão}} \]
\[ p = 2 \times \left(1 - \Phi(|t|)\right) \]
em que \(\Phi\) é a função de distribuição cumulativa da distribuição t com \(n - 1\) graus de liberdade.
# Define the input variables
a <- 118
b <- 5 # Number of disagreements in the first category
c <- 2 # Number of disagreements in the second category
d <- 0 # Number of agreements in the second category
# Sample size
n <- a + b + c + d
df <- n - 1
# Probabilities
p11 <- a / n
p12 <- b / n
p21 <- c / n
p22 <- d / n
# Overall agreement probability
pa <- p11 + p22
# Marginal probabilities
pR1 <- (a + b) / n
pR2 <- (c + d) / n
pC1 <- (a + c) / n
pC2 <- (b + d) / n
# Average marginal probabilities
p_bar1 <- (pR1 + pC1) / 2
p_bar2 <- (pR2 + pC2) / 2
# Chance-agreement probabilities
pc_kappa <- (pR1 * pC1) + (pR2 * pC2)
pc_AC1 <- p_bar1 * (1 - p_bar1) + p_bar2 * (1 - p_bar2)
pc_G <- 0.5
# Agreement measures
kappa <- (pa - pc_kappa) / (1 - pc_kappa)
AC1 <- (a^2 + d^2 - (b + c)^2 / 2) / (a^2 + d^2 + (b + c)^2 / 2 + (a + d) * (b + c))
G <- ((a + d) - (b + c)) / n
# Standard errors
se_kappa2 <- (1 / (n * (1 - pc_kappa)^2)) * (pa * (1 - pa)
- 4 * (1 - kappa) * (p11 * p_bar1 +
p22 * p_bar2 -
pa * pc_kappa) +
+ 4 * (1 - kappa)^2 * (p11 * (pR1 + pC1)^2 / 4
+ p12 * (pR1 + pC2)^2 / 4
+ p21 * (pR2 + pC1)^2 / 4
+ p22 * (pR2 + pC2)^2 / 4
- pc_kappa^2))
se_kappa <- sqrt(se_kappa2)
se_AC12 <- (1 / (n * (1 - pc_AC1)^2)) * (pa * (1 - pa)
- 4 * (1 - AC1) * (p11 * (1 - p_bar1) +
p22 * (1 - p_bar2) -
pa * pc_AC1)
+ 4 * (1 - AC1)^2 * (p11 * (1 - (p_bar1 + p_bar1) / 2)^2
+ p12 * (1 - (p_bar1 + p_bar2) / 2)^2
+ p21 * (1 - (p_bar2 + p_bar1) / 2)^2
+ p22 * (1 - (p_bar2 + p_bar2) / 2)^2
- pc_AC1^2))
se_AC1 <- sqrt(se_AC12)
se_G2 <- (4 / n) * pa * (1 - pa)
se_G <- sqrt(se_G2)
# t-values
t_kappa <- kappa / se_kappa
t_AC1 <- AC1 / se_AC1
t_G <- G / se_G
# p-values (two-tailed)
p_kappa <- 2 * (1 - pt(abs(t_kappa), df = n - 1))
p_AC1 <- 2 * (1 - pt(abs(t_AC1), df = n - 1))
p_G <- 2 * (1 - pt(abs(t_G), df = n - 1))
# Print the results
# cat("Kappa: ", kappa, "\n")
# cat("AC1: ", AC1, "\n")
# cat("G: ", G, "\n")
# cat("Standard Error (Kappa): ", se_kappa, "\n")
# cat("Standard Error (AC1): ", se_AC1, "\n")
# cat("Standard Error (G): ", se_G, "\n")
# cat("Degrees of Freedom: ", df, "\n")
# cat("t-value (Kappa): ", t_kappa, "\n")
# cat("p-value (Kappa): ", p_kappa, "\n")
# cat("t-value (AC1): ", t_AC1, "\n")
# cat("p-value (AC1): ", p_AC1, "\n")
# cat("t-value (G): ", t_G, "\n")
# cat("p-value (G): ", p_G, "\n")| Measure | Value | SE | t | p | df |
|---|---|---|---|---|---|
| Kappa | -0.0234 | 0.0812 | -0.2880 | 0.7738 | 124 |
| AC1 | 0.9408 | 0.0230 | 40.9665 | 0.0000 | 124 |
| G | 0.8880 | 0.0411 | 21.5903 | 0.0000 | 124 |
Teste AC2 de Gwet em tabela \(k\times k\) para dois avaliadores
O teste AC2 de Gwet em tabela \(k\times k\) é aplicável para a comparação de dois métodos de avaliação com desfecho ordinal politômico. AC2 é medida de tamanho de efeito.
Em Agresti (1990), apud Mehta & Patel (1996, p. 72), há o exemplo da avaliação de estadiamento tumoral submetido à apreciação de dois patologistas.
O objetivo é analisar a concordância de diagnóstico entre dois patologistas que classificaram conforme a sereridade de uma determinada lesão uterina de 118 lâminas de diferentes mulheres:
- N: Negativo
- HEA: Hiperplasia escamosa atipica
- CIS: carcinoma in situ
- CE: Carcinoma escamoso
- CI: Carcinoma invasivo
Formulamos:
\[ \begin{cases} H_0: \text{Ausência de concordância}\\ H_1: \text{Presença de concordância} \end{cases}\\ \alpha=5\% \]
Este teste é realizado pelo coeficiente de concordância \(\text{AC}_2\) de Gwet.
Ele permite também avaliar a concordância entre três ou mais avaliadores.
O conjunto está implementado em AC2Ordinal.R
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
# Agresti (1990), apud Mehta, C. R. &
# Patel, N. R. (1996) SPSS Exact Tests 7 for Windows. IL: SPSS, p. 72.
# Pap-Smear Classification by Two Pathologists
# O objetivo é analisar a concordância de diagnóstico entre 2 patologistas
# que classificaram conforme a sereridade de uma determinada
# lesão uterina de 118 lâminas de diferentes mulheres.
# N=Negativo, HEA=Hiperplasia escamosa atipica, CIS=carcinoma in situ
# CE=Carcinoma escamoso, CI=Carcinoma invasivo
alfa <- 0.05
Tabela <- ("
Patologistas N HEA CIS CE CI
N 22 2 2 0 0
HEA 5 7 14 0 0
CIS 0 2 36 0 0
CE 0 1 14 7 0
CI 0 0 3 0 3
")
print(TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE, row.names=1)))
print(gplots::balloonplot(t(as.table(TC)),
main ="Lâmina de lesão uterina",
xlab ="Patologista 1",
ylab="Patologista 2",
label = TRUE,
show.margins = TRUE))
cat("Teste AC2 de Gwet para duas variáveis ordinais\n")
cat("H0: AC2 = 0 vs H1: AC2 != 0\n",sep="")
# AC2 de Gwet
print(irrCAC::gwet.ac1.table(TC,
irrCAC::linear.weights(1:nrow(TC)),
conflev = 1-alfa))
print(vcdExtra::GKgamma(TC,
level = 1-alfa)) N HEA CIS CE CI
N 22 2 2 0 0
HEA 5 7 14 0 0
CIS 0 2 36 0 0
CE 0 1 14 7 0
CI 0 0 3 0 3
NULL
Teste AC2 de Gwet para duas variáveis ordinais
H0: AC2 = 0 vs H1: AC2 != 0
coeff.name coeff.val coeff.se coeff.ci coeff.pval
1 Gwet's AC2 0.7809185 0.03055234 (0.72,0.841) 0e+00
gamma : 0.923
std. error : 0.033
CI : 0.858 0.988 No teste estatístico, rejeitamos \(H_0\).
Kappa de Fleiss para três ou mais avaliadores
Suponha três ou mais avaliadores para os mesmos sujeitos com desfecho nominal politômico.
Por exemplo, seis psiquiatras avaliam os pacientes identificando:
- Depressão
- Desordem de personalidade
- Esquizofrenia
- Neurose
- Outros diagnósticos
Formulamos:
\[ \begin{cases} H_0: \text{Presença de concordância entre os avaliadores}\\ H_1: \text{Ausência de concordância entre os avaliadores} \end{cases}\\ \alpha=5\% \]
Neste código irr::kappam.fleiss verifica a concordância
global e, depois, verifica a concordância para cada um dos
diagnósticos.
Esta função precisa de dados completos. Qualquer valor faltante faz com que todos as observações do indivíduo sejam desconsideradas.
Os dados fazem parte do pacote irr e são usados em Fleiss Kappa for m raters.R
data(diagnoses, package = "irr")
print(head(diagnoses))
cat("\n...\n")
print(tail(diagnoses))
summary(diagnoses)
cat("\n")
# Fleiss' and category-wise Kappa
print(irr::kappam.fleiss(diagnoses, detail=TRUE)) rater1 rater2 rater3
1 4. Neurosis 4. Neurosis 4. Neurosis
2 2. Personality Disorder 2. Personality Disorder 2. Personality Disorder
3 2. Personality Disorder 3. Schizophrenia 3. Schizophrenia
4 5. Other 5. Other 5. Other
5 2. Personality Disorder 2. Personality Disorder 2. Personality Disorder
6 1. Depression 1. Depression 3. Schizophrenia
rater4 rater5 rater6
1 4. Neurosis 4. Neurosis 4. Neurosis
2 5. Other 5. Other 5. Other
3 3. Schizophrenia 3. Schizophrenia 5. Other
4 5. Other 5. Other 5. Other
5 4. Neurosis 4. Neurosis 4. Neurosis
6 3. Schizophrenia 3. Schizophrenia 3. Schizophrenia
...
rater1 rater2 rater3
25 1. Depression 4. Neurosis 4. Neurosis
26 2. Personality Disorder 2. Personality Disorder 2. Personality Disorder
27 1. Depression 1. Depression 1. Depression
28 2. Personality Disorder 2. Personality Disorder 4. Neurosis
29 1. Depression 3. Schizophrenia 3. Schizophrenia
30 5. Other 5. Other 5. Other
rater4 rater5 rater6
25 4. Neurosis 4. Neurosis 5. Other
26 2. Personality Disorder 2. Personality Disorder 4. Neurosis
27 1. Depression 5. Other 5. Other
28 4. Neurosis 4. Neurosis 4. Neurosis
29 3. Schizophrenia 3. Schizophrenia 3. Schizophrenia
30 5. Other 5. Other 5. Other
Fleiss' Kappa for m Raters
Subjects = 30
Raters = 6
Kappa = 0.43
z = 17.7
p-value = 0
Kappa z p.value
1. Depression 0.245 5.192 0.000
2. Personality Disorder 0.245 5.192 0.000
3. Schizophrenia 0.520 11.031 0.000
4. Neurosis 0.471 9.994 0.000
5. Other 0.566 12.009 0.000Sendo \(H_0\) rejeitada, conclui-se que não há concordância entre os avaliadores, globalmente e para cada diagnóstico em particular.
ICC: IntraClass Correlation
Este coeficiente de correlação intraclasse é aplicável para variável de desfecho intervalar e precisamos verificar a concordância entre dois ou mais avaliadores.
Em SPSS Base 10: Applications Guide (1999), há um exemplo sobre treinamento que o IOC (International Olympic Commitee) utiliza para que os juizes das competições de ginástica sejam concordantes.
Optar pela análise de concordância absoluta significa que está interessado na igualdade das avaliações dos juízes/métodos, i.e., se eles usam a mesma escala para avaliar as mesmas unidades experimentais. As hipóteses são:
\[ \begin{cases} H_0: r_0=0~\text{(ou ICC = 0, ausência de concordância entre os juízes)}\\ H_1: r_0>0~\text{(ou ICC > 0, concordância positiva entre os juízes)} \end{cases}\\ \alpha=5\% \]
Neste exemplo, 7 juízes são verificados com a implementação em ICC_IOC.R.
O oitavo juíz é removido porque é um entusiasta:
# O International Olympic Commitee (IOC) treina juizes para competicoes
# de ginastica.
# SPSS Base 10: Applications Guide, 1999.
Dados <- haven::read_sav("judges.sav")
Dados <- Dados[,1:7]
# Good absolute agreement between the raters, ICC=0.723, was observed while using
# the two-way random effect models and single rater with a p-value of 3.24e-11.
print(irr::icc(Dados,
model="twoway",
type="agreement",
unit="single"))
cat("\n\nMatriz de correlacao:\n")
print(mc <- round(cor(Dados),2))
for (rc in 1:7)
{
mc[rc,rc] <- NA
}
cat("\ncorrelacao media = ",round(mean(mc,na.rm=TRUE),3)) Single Score Intraclass Correlation
Model: twoway
Type : agreement
Subjects = 300
Raters = 7
ICC(A,1) = 0.723
F-Test, H0: r0 = 0 ; H1: r0 > 0
F(299,15.1) = 52.1 , p = 3.24e-11
95%-Confidence Interval for ICC Population Values:
0.513 < ICC < 0.832
Matriz de correlacao:
judge1 judge2 judge3 judge4 judge5 judge6 judge7
judge1 1.00 0.91 0.92 0.91 0.93 0.91 0.90
judge2 0.91 1.00 0.87 0.94 0.88 0.93 0.85
judge3 0.92 0.87 1.00 0.88 0.94 0.86 0.92
judge4 0.91 0.94 0.88 1.00 0.88 0.92 0.85
judge5 0.93 0.88 0.94 0.88 1.00 0.87 0.92
judge6 0.91 0.93 0.86 0.92 0.87 1.00 0.85
judge7 0.90 0.85 0.92 0.85 0.92 0.85 1.00
correlacao media = 0.897O modelo two-way random (parâmetro
model="twoway") pode ser utilizado e significa que os
juízes (avaliadores) e unidades experimentais (avaliandos) no estudo
constituem amostras aleatórias de suas respectivas populações de juízes
e de unidades observacionais. Portanto, nesse modelo, as conclusões são
extensíveis para a população dos juízes.
Neste exemplo, ICC(A,1) = 0.723 representa o grau de
concordância geral entre os sete juízes. Esta é uma medida de tamanho de
efeito associada à correlação alta entre as notas dadas pelos juízes
(para comparação, exibimos a matriz de correlações e sua média). Note
que \(H_0\) foi rejeitada e, portanto,
existe concordância entre as notas dadas pelos juízes.
Teste de simetria
O teste de simetria em tabela de contingência \(k\times k\) com duas variáveis intervalares (McNemar, Evans-Hoenig, Bowker) não depende de delineamento intraparticipantes no sentido longitudinal (repetição no mesmo sujeito ao longo do tempo), mas de observações emparelhadas entre dois métodos ou avaliadores aplicados aos mesmos sujeitos ou unidades.
É, portanto, um delineamento entre participantes quanto ao fenômeno medido, mas com medidas emparelhadas - cada unidade (ou paciente, ou lâmina) é avaliada duas vezes por métodos distintos (referência e alternativo).
O teste avalia a hipótese nula de que há ausência de assimetria sistemática entre os métodos. Assim, é um delineamento entre participantes no sentido de comparar dois métodos distintos aplicados às mesmas unidades, mas com observações emparelhadas (cross-classified).
O objetivo é detectar viés direcional, que é uma forma de concordância.
A hipótese nula do teste de simetria (McNemar ou Bowker) é:
\[ H_0: p_{ij} = p_{ji} \quad \forall i \neq j \]
Isto é, a probabilidade de uma unidade ser classificada como categoria i pelo método A e j pelo método B é igual à probabilidade reversa.
Em termos práticos: não há viés sistemático entre os dois métodos; as discordâncias se distribuem simetricamente em torno da diagonal principal da tabela.
A rejeição de \(H_0\) indica assimetria, isto é, que um método tende a superestimar ou subestimar as categorias em relação ao outro.
Exemplo: A verdadeira idade do peixe American Shad
Conforme Ogle (2016, p. 76), a verdadeira idade de cada peixe deve ser conhecida para que se possa determinar a exatidão das estimativas de idade. O viés entre estruturas ou entre técnicos é determinado a partir de estimativas diversas de idade para o mesmo peixe. Por fim, a precisão entre estruturas ou técnicos é examinada com duas ou mais estimativas de idade para o mesmo peixe.
Como exemplo, McBride et al. (2005) examinaram a validade do uso de escamas para estimar a idade do American Shad. As idades verdadeiras dos peixes em sua amostra eram conhecidas porque os peixes haviam sido marcados antes do repovoamento. Além disso, 13 biólogos estimaram duas vezes (de forma independente) a idade a partir das escamas de cada peixe.
A idade verdadeira de cada peixe (trueAge) e as
estimativas de idade de três dos 13 biólogos estão disponíveis no
arquivo ShadCR.txt. As variáveis de idade estimada são
rotuladas com o prefixo ager, seguido por uma letra que
identifica o biólogo (A, B ou C) e um número que indica a ocasião em que
a escama foi interpretada (1 ou 2). Alguns biólogos optaram por não
atribuir idade a certas escamas e, portanto, esses dados estão ausentes
(registrados como valores NA).
Esses dados do Shad serão usados para examinar: (1) a exatidão, comparando as idades estimadas por um biólogo com as idades verdadeiras (usando, por simplicidade, apenas a primeira estimativa de cada biólogo); (2) o viés entre biólogos, comparando as estimativas de idade de dois biólogos (novamente, apenas as primeiras estimativas); e (3) a precisão, comparando as duas estimativas do mesmo biólogo e também as primeiras estimativas entre os três biólogos.
library(fishmethods)
# Introductory fisheries analyses with R - Ogle - 2016
# https://fishr-core-team.github.io/fishR/
shad <- read.csv("ShadCR.txt")
head(shad)
tail(shad)
ab.tA1 <- FSA::ageBias(agerA1~trueAge,
data=shad,
ref.lab="True Age",
nref.lab="Ager A")
plot(ab.tA1,
col.CIsig="black")
summary(ab.tA1,
what="bias")
summary(ab.tA1,
what="table")
summary(ab.tA1,
what="symmetry")
ab.tC1 <- FSA::ageBias(agerC1~trueAge,
data=shad,
ref.lab="True Age",
nref.lab="Ager C")
summary(ab.tC1,
what="symmetry")
plot(ab.tC1,
col.CIsig="black")
ap.A <- FSA::agePrecision(~agerA1+agerA2,
data=shad)
summary(ap.A,
what="difference")
summary(ap.A,
what="absolute difference")
summary(ap.A,
what="precision")
ap.ABC <- FSA::agePrecision(~agerA1+agerB1+agerC1,
data=shad)
summary(ap.ABC,
what="difference")
summary(ap.ABC,
what="precision")
Data <- na.omit(subset(shad, select=c(2,3)))
fishmethods::compare2(readings=Data,
usecols=c(1,2),
twovsone=TRUE,
plot.summary=TRUE,
barplot=TRUE,
chi=TRUE,
pool.criterion=1,
cont.cor=FALSE,
#correct="Yates",
first.name="trueAge",
second.name="agerA1")
trueAge n min max mean SE t adj.p sig LCI UCI
3 9 2 4 3.11 0.200 0.555 1.000 FALSE 2.65 3.57
4 11 3 6 4.33 0.302 1.106 1.000 FALSE 3.66 5.01
5 10 4 7 5.00 0.298 0.000 1.000 FALSE 4.33 5.67
6 10 3 8 5.80 0.467 -0.429 1.000 FALSE 4.74 6.86
7 10 5 8 6.30 0.367 -1.909 0.443 FALSE 5.47 7.13
8 3 6 7 6.67 0.333 -4.000 0.343 FALSE 5.23 8.10
True Age
Ager A 2 3 4 5 6 7 8
2 - 1 - - - - -
3 - 6 2 - 1 - -
4 - 2 3 3 - - -
5 - - 3 5 3 3 -
6 - - 1 1 4 3 1
7 - - - 1 - 2 2
8 - - - - 2 2 -
symTest df chi.sq p
1 McNemar 1 1.580645 0.2086678
2 EvansHoenig 3 2.636364 0.4511502
3 Bowker 10 8.333333 0.5963127
symTest df chi.sq p
1 McNemar 1 30.42222 3.475241e-08
2 EvansHoenig 4 31.78261 2.119129e-06
3 Bowker 16 35.66667 3.218018e-03
-2 -1 0 1 2 3
3.922 15.686 50.980 19.608 7.843 1.961
0 1 2 3
50.980 35.294 11.765 1.961
n validn R PercAgree ASD ACV AAD APE
53 51 2 50.98 0.4575 9.739 0.3235 6.887
-2 -1 0 1 2 3 4
agerA1 - agerB1 18.182 15.152 45.455 12.121 9.091 0.000 0.000
agerA1 - agerC1 0.000 5.882 19.608 41.176 19.608 11.765 1.961
agerB1 - agerC1 0.000 6.061 6.061 36.364 36.364 15.152 0.000
n validn R PercAgree ASD ACV ACV2 AAD APE APE2 AD
53 33 3 6.061 1.001 22.98 24.94 0.7273 16.7 14.27 13.27
Delineamento intraparticipantes
Teste qui-quadrado de mudança de McNemar
O teste de McNemar (1947) é utilizado para testar a igualdade das proporções de resposta binária de duas populações nas quais os dados consistem em respostas pareadas e dependentes, uma de cada população. É tipicamente usado em situações de medidas repetidas, onde a resposta de cada sujeito é solicitada duas vezes, uma antes e outra depois de um evento especificado (tratamento) ocorrer. O teste então determina se a proporção de resposta inicial (antes do evento) é igual à proporção de resposta final (após o evento).
O teste de McNemar para a significância das mudanças é particularmente aplicável aos projetos “antes-depois“ de intervenção, nos quais cada sujeito é usado como seu próprio controle e em que as medidas são feitas em escala nominal dicotômica.
As condições podem ser usadas para testar a eficácia de um tratamento específico (reuniões, editoriais de jornais, discursos, visitas pessoais etc.) sobre as preferências dos eleitores sobre candidatos a cargos públicos ou para testar o efeito de migração do campo para a cidade pela afiliação política de pessoas.
Observe que nesses estudos as pessoas podem servir como seu próprio controle e que a escala nominal (ou categorização) é usada adequadamente para avaliar a mudança “antes-depois”.
Assim, \(b + c\) é o total de pessoas cujas respostas foram alteradas.
A hipótese nula é que o número de mudanças em cada direção é o mesmo na população.
Portanto, dos indivíduos com \(b + c\) que mudaram, esperamos que \((b + c) / 2\) indivíduos mudem de \(+\) para \(-\) e \((b + c) / 2\) pessoas mudem de \(-\) para \(+\).
Em outras palavras, quando \(H_0\) é verdadeira, a frequência esperada em cada uma das duas caselas é \((b + c) / 2\).
No teste de McNemar para a significância das mudanças, estamos interessados apenas nas caselas nas quais as mudanças podem ocorrer.
Assim, se \(b\) é o número de casos observados cujas respostas mudaram de \(+\) para \(-\), \(c\) é o número observado de casos que mudaram de \(-\) para \(+\) e \((b + c) / 2\) é o número esperado de casos nas caselas \(b\) e \(c\).
\[ \begin{align} X^2 &= \sum_{i=1}^{2}{{(o_i - e_i)^2} \over{e_i}} \\ &= {{\left(b-{{b+c}\over{2}}\right)^2} \over {{b+c}\over{2}}} + {{\left(c-{{b+c}\over{2}}\right)^2} \over {{b+c}\over{2}}}\\ X^2&={{(b - c)^2} \over{b+c}} \end{align} \]
A distribuição amostral de \(X^2\) quando \(H_1\) é verdadeira, é distribuída assintoticamente como qui-quadrado com graus de liberdade iguais a 1.
\[ \begin{cases} H_0: P(+ \to \, – ) = P(– \, \to +) \\ H_1: P(+ \to \, – ) \ne P(– \, \to +) \end{cases}\\ \alpha=5\% \] ### Exemplo: Intenção de voto antes e depois de um debate eleitoral
Um exemplo é dado por Siegel & Castellan (1988, p. 77). Neste exemplo verificou-se a intenção de voto antes e depois de um debate eleitoral. Quando há um debate eleitoral, se houver manutenção do carregamento da diagnonal principal, não se rejeitando a hipótese nula, significa que as intenções de voto não mudaram.
Este teste está implementado em coin::mh_test que já foi
utilizado antes, aqui com a opção
distribution = "exact".
No exemplo do debate eleitoral:
library(exact2x2)
source("eiras.matrix2dataframe.R")
# Siegel & Castellan (1988), p. 77
Tabela <- ("
TVDebate Carter Reagan
Carter 28 13
Reagan 7 27
")
cat(Tabela)
TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE, row.names=1))
print(mcnemar.test(TC,correct=FALSE)) # classico
print(coin::mh_test(as.table(TC), distribution = "exact"))
print(exact2x2::mcnemar.exact(TC))
print(gplots::balloonplot(t(as.table(TC)),
main ="Intenção de voto pre-pós debate",
xlab ="Depois", ylab="Antes",
label = TRUE,
show.margins = TRUE))
pre <- c(rep(1,sum(TC[1,1:2])),rep(0,sum(TC[2,1:2])))
pos <- c(rep(1,sum(TC[1:2,1])),rep(0,sum(TC[1:2,2])))
data <- data.frame(matrix(nrow=length(pre),ncol=3))
names(data) <- c("ID","pre","pos")
data$ID <- 1:length(pre)
data$pre <- pre
data$pos <- pos
data$change <- data$pos-data$pre
summary(data[-1])
print(t <- t.test(data$change,mu=0))
misty::ci.prop.diff(data$pre, data$pos, paired = TRUE, digits=4)
# RR+
print(ep <- epiR::epi.2by2(TC))
# RR-
print(ep <- epiR::epi.2by2(DescTools::Rev(TC, margin=1)))
TVDebate Carter Reagan
Carter 28 13
Reagan 7 27
McNemar's Chi-squared test
data: TC
McNemar's chi-squared = 1.8, df = 1, p-value = 0.1797
Exact Marginal Homogeneity Test
data: response by
conditions (Var1, Var2)
stratified by block
chi-squared = 1.8, p-value = 0.2632
Exact McNemar test (with central confidence intervals)
data: TC
b = 13, c = 7, p-value = 0.2632
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.6886514 5.4973372
sample estimates:
odds ratio
1.857143 
NULL
One Sample t-test
data: data$change
t = -2.5367, df = 74, p-value = 0.0133
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-0.14283923 -0.01716077
sample estimates:
mean of x
-0.08
Two-Sided 95% Confidence Interval: Difference in Proportions from Paired Samples
n nNA p1 p2 p.Diff Low Upp
75 0 0.5467 0.4667 -0.0800 -0.1388 -0.0190
Outcome+ Outcome- Total Inc risk *
Exposure+ 28 13 41 68.29 (51.91 to 81.92)
Exposure- 7 27 34 20.59 (8.70 to 37.90)
Total 35 40 75 46.67 (35.05 to 58.55)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 3.32 (1.66, 6.63)
Inc odds ratio 8.31 (2.88, 23.98)
Attrib risk in the exposed * 47.70 (28.02, 67.39)
Attrib fraction in the exposed (%) 69.85 (43.36, 85.23)
Attrib risk in the population * 26.08 (8.41, 43.75)
Attrib fraction in the population (%) 55.88 (35.28, 75.17)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 16.995 Pr>chi2 = <0.001
Fisher exact test that OR = 1: Pr>chi2 = <0.001
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units
Outcome+ Outcome- Total Inc risk *
Exposure+ 7 27 34 20.59 (8.70 to 37.90)
Exposure- 28 13 41 68.29 (51.91 to 81.92)
Total 35 40 75 46.67 (35.05 to 58.55)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 0.30 (0.15, 0.60)
Inc odds ratio 0.12 (0.04, 0.35)
Attrib risk in the exposed * -47.70 (-67.39, -28.02)
Attrib fraction in the exposed (%) -231.71 (-577.04, -76.54)
Attrib risk in the population * -21.63 (-39.80, -3.45)
Attrib fraction in the population (%) -46.34 (-48.11, -39.90)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 16.995 Pr>chi2 = <0.001
Fisher exact test that OR = 1: Pr>chi2 = <0.001
Wald confidence limits
CI: confidence interval
* Outcomes per 100 population units Como \(H_0\) não foi rejeitada, conclui-se que o debate não provocou desigualdade nas mudanças de intenções de voto.
Apêndices
Apêndice 1: distribuição qui-quadrado
O teste qui-quadrado usa as propriedades da distribuição de mesmo nome. Esta distribuição, diferente da distribuição normal que têm domínio de \([- \infty, + \infty]\), estas são distribuições assimétricas, iniciadas em zero, que tendem assintoticamente a zero com o valor crescente de \(\chi^2\), \([0, + \infty]\).
O comando para calcular a densidade de probabilidade para cada valor de qui-quadrado é:
dchisq(x, df, ncp = 0, log = FALSE)
em que x é o valor de \(\chi^2\), df são os graus de
liberdade (degrees of freedom) e ncp é o parâmetro
de não centralidade (non-centrality parameter = 0, por
default, correspondente à distribuição centrada de \(H_0\); voltaremos ao ncp
adiante). O parâmetro log, desligado por default é
para a possibilidade de devolver o logarítmo do valor (não o
exploraremos neste texto).
Pode-se gerar curvas desta distribuição com o seguinte código, mostrando seu aspecto entre 1 e 6 e 50 graus de liberdade:
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 1:6
for(g in c(1:length(graus.liberdade)))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
plot(chi2.valor, p, type="l",
main=paste("Distribuição Qui-quadrado, df=",
graus.liberdade[g],sep=""),
xlab="qui-quadrado", ylab="Densidade")
}
ou, caso prefira ver todas as curvas no mesmo gráfico (o eixo das
ordenadas foi truncado em 0.8 por causa da curva com df=1
para melhor visualização), sugere-se:
source("friendlycolor.R")
cor <- 3
padrao <- 1
leg.txt <- c()
leg.cor <- c()
chi2.valor <- seq(from=0, to=40, by=0.01)
graus.liberdade <- 1:8
for (g in 1:length(graus.liberdade))
{
p <- dchisq(chi2.valor, df=graus.liberdade[g])
if (g==1)
{
plot(chi2.valor, p, type="l", lty=padrao,
col=friendlycolor(cor), lwd=2,
ylim=c(0, 0.8),
main="Distribuição Qui-quadrado",
xlab="qui-quadrado", ylab="Densidade")
} else
{
lines(chi2.valor, p, lty=padrao,
col=friendlycolor(cor), lwd=2)
}
# guarda para a legenda
leg.txt <- c(leg.txt, paste("df = ",graus.liberdade[g],sep=""))
leg.cor <- c(leg.cor, friendlycolor(cor))
cor <- cor+6
padrao <- padrao+1
}
legend("topright",
leg.txt,
col=leg.cor,
lwd=2,
lty=1:6,
bty="n") 
A função qchisq faz parte do conjunto:
pchisq: dado um valor \(X^2\) fornece a probabilidade acumulada, i.e., a área sob a curva. Por default (comlower.tail=TRUE), computa a área que vai de zero ao valor \(X^2\) solicitado.qchisq: faz o contrário depchisq, dado um valor de probabilidade fornece o valor de \(X^2\) correspondente.dchisq: dado um valor \(X^2\) fornece a densidade de probabilidade para aquele ponto; é a que usamos para a construção dos gráficos da distribuição \(\chi^2\), fornecendo a altura, no eixo das ordenadas, para cada valor de \(X^2\) solicitado.rchisq: fornece números pseudo-aleatórios sorteados com distribuição \(\chi^2\).
O seguinte código gera 4 gráficos, mostrando o comportamento das 4 variantes para 3 graus de liberdade:
graus.liberdade = 3
# densidade de probabilidade
# mostra a distribuicao X^2
X2 <- seq(from=0, to=15, by=0.01)
dX2 <- dchisq(X2,df=graus.liberdade)
plot(X2,dX2,main="Distribuição X^2\nobtida com dchisq",type="l")
# da probabilidade para X^2
pX2 <- pchisq(X2,df=graus.liberdade)
plot(X2,pX2,main="Probabilidade acumulada\nobtida com pchisq",type="l")
# de X^2 para probabilidade
p <- seq(from=0, to=1, length.out = 1000)
X2 <- qchisq(p,df=graus.liberdade)
plot(X2,p,main="Valores de X^2\nobtidos com qchisq",type="l")
# 1000 valores rand com distribuicao X^2
r <- rchisq(1000,df=graus.liberdade)
density.r <- density(r)
plot(density.r,main="Distribuição de números aleatórios\nobtidos com rchisq")
Observe que as distribuições disponíveis no R-base, por convenção,
sempre estão em conjuntos iniciados pelas letras p,
q, d e r. Por exemplo, a
distribuição normal tem pnorm, qnorm,
dnorm e rnorm, a distribuição binomial tem
pbinom, qbinom, dbinom e
rbinom, a distribuição t tem pt,
qt, dt e rt e assim por
diante.
Apêndice 2: objeto de chisq.test
Quando executamos
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
print(tabela) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508a variável chi2 armazena tudo que a função
chisq.test retorna. É comum obtermos conjuntos de variáveis
como este a partir de muitas funções. São objetos dos quais podemos
obter dados para uso nos passos seguintes de um código R.
Neste exemplo, podemos ver seu conteúdo e os tipos das variáveis
contidas em chi2 com:
List of 9
$ statistic: Named num 7.31
..- attr(*, "names")= chr "X-squared"
$ parameter: Named int 1
..- attr(*, "names")= chr "df"
$ p.value : num 0.00686
$ method : chr "Pearson's Chi-squared test"
$ data.name: chr "tabela"
$ observed : 'table' num [1:2, 1:2] 17 130 138 508
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ expected : num [1:2, 1:2] 28.7 118.3 126.3 519.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ residuals: 'table' num [1:2, 1:2] -2.189 1.079 1.044 -0.515
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
$ stdres : 'table' num [1:2, 1:2] -2.7 2.7 2.7 -2.7
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "Capacete +" "Capacete -"
.. ..$ : chr [1:2] "Trauma +" "Trauma -"
- attr(*, "class")= chr "htest" statistic parameter p.value method data.name observed
"double" "integer" "double" "character" "character" "double"
expected residuals stdres
"double" "double" "double" Assim, o resultado principal do teste é exibido por:
Pearson's Chi-squared test
data: tabela
X-squared = 7.3099, df = 1, p-value = 0.006858Mas podemos acessar mais detalhes. Por exemplo em
chi2\$expected está a tabela dos valores esperados sob
\(H_0\).
Trauma + Trauma -
Capacete + 28.73266 126.2673
Capacete - 118.26734 519.7327O valor p isolado do teste está em
[1] 0.006857643Localize estas variáveis na saída de str e explore o
objeto para encontrar mais de seus detalhes.
Apêndice 3: transformando planilhas em matrizes
Há várias formas de fornecer uma tabela para uma função em R.
Para a entrada de dados na forma de uma matriz 2×2 o seguinte código funciona:
Tabela <- ("
TC2x2 Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25
")
TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE, row.names=1))
print(TC)Obtendo-se:
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25O problema está na entrada de dados. Somente espaços em branco podem ser usados para a montagem da tabela. O alinhamento é visual, mas não necessário para o R. Funciona igualmente
Tabela <- ("
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25
")
TC <- as.matrix(read.table(textConnection(Tabela),
header=TRUE,
row.names=1))
print(TC) Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25Portanto, pode ser ruim a visualização de sua entrada.
Uma opção mais “limpa” é sempre colocar os dados em planilhas no
formato .xls ou .xlsx e usar a função
readxl::read_excel para utilizá-los em R:
Esta planilha é tabagismo_e_etilismo_2x2.xlsx, lida por
New names:
• `` -> `...1`# A tibble: 2 × 3
...1 Tabagista NaoTabagista
<chr> <dbl> <dbl>
1 Etilista 50 15
2 NaoEtilista 20 25A função read_excel avisa que teve que suprir um nome
(por causa da célula A1 que está vazia), mas esta célula será desprezada
adiante.
No entanto, temos um problema a resolver. A função
read_excel produz um data frame e a função
chisq.test precisa de uma matriz para sua entrada:
[1] TRUE[1] FALSETemos, portanto, que converter TCnew para a
representação em matrix. Existe uma função para isto:
[1] TRUE ...1 Tabagista NaoTabagista
[1,] 1 50 15
[2,] 2 20 25Ainda não é o desejado. TCmatrix tem 2 linhas mas 3
colunas, além de ter perdido os nomes das linhas que TCnew
trazia em sua primeira coluna:
[1] 2[1] 3NULL# A tibble: 2 × 1
...1
<chr>
1 Etilista
2 NaoEtilistaA solução, portanto, é transferir o conteúdo de TCnew
como nomes das linhas de TCmatrix e, então, deletar a
coluna 1 de TCmatrix que contém valores
NA.
Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25Colocando tudo junto, para terminar com uma matriz chamada
TC (como era originalmente), podemos ler a planilha em uma
variável (e.g., TCtmp), reorganizar seu conteúdo em
TC e remover TCtmp da memória com a função
remove. O código enxuto é:
New names:
• `` -> `...1`TC <- data.matrix(TCtmp)
rownames(TC) <- as.character(unlist(TCtmp[1:2,1]))
TC <- TC[,-1]
remove(TCtmp)
print(TC) Tabagista NaoTabagista
Etilista 50 15
NaoEtilista 20 25Mensagens como:
New names: * `` -> …1
e
Warning in data.matrix(TCnew): NAs introduced by coercion
podem aparecer aqui. Warnings são avisos do R que não impedem o funcionamento do programa (avalie, caso a caso, se precisa corrigir algo; aqui não é necessário). Neste caso, o motivo é a célula A1 na planilha estar vazia. Preencha com algo (um ponto, uma palavra, qualquer coisa - esta célula será desprezada) e ambas as mensagens desaparecem.
A primeira acontece porque read_excel lê uma planilha
contando com a primeira linha para nomear as colunas de um data
frame. Como a célula está em branco, avisa que fez a atribuição de
um novo nome. A segunda mensagem, lidando com uma célula vazia, avisa
que utiliza NA, uma constante do R que assinala valor
faltante (do inglês, Not Available, não
disponível).
Referências
- Agresti, A (2002) Categorical Data Analysis. 2nd ed. NJ: Wiley.
- Appleton, DR et al. (1996) Ignoring a covariate: an example of Simpson’s paradox. The American Statistician 50(4): 340-1.
- Beckles et al. (1985) International Journal of Obesity 9:127-35, apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2×2 and larger contingency tables.
- Bilder, CR & Loughin,TM (2024) Analysis of categorical data with R. 2nd ed. USA: CRC. https://www.chrisbilder.com/categorical/
- Catania, AC (1999) Aprendizagem: Comportamento, Linguagem, Cognição. 4a ed. Porto Alegre: Artes Médicas.
- Cochran, WG (1954) Some Methods for Strengthening the Common \(\chi^2\) Tests. Biometrics 10(4): 417-51.
- Dancey C & Reidy J (2019) Estatística sem Matemática para Psicologia. 7a ed. Porto Alegre: PENSO.
- Duncan, OD (1979) How Destination Depends on Origin in the Occupational Mobility Table. American Journal of Sociology 84: 793-803.
- Edwardes, MD & Baltzan, M (2000) The generalization of the odds ratio, risk ratio and risk difference to r×k tables. Statistics in Medicine 19: 1901-14.
- Field, A et al. (2012) Discovering statistics using R. London: Sage.
- García-Pérez, MA & Núñez-Antón, V (2003) Cellwise Residual Analysis in Two-Way Contingency Tables. Educational and Psychological Measurement 63(5): 825–39. https://doi.org/10.1177/0013164403251280
- Goodman, LA (1979) Simple models for the analysis of association in cross-classifications having ordered categories. Journal of the American Statistical Association 74(367): 537–52. https://doi.org/10.1080/01621459.1979.10481639
- Gudziunaite, S & Moshammer, H (2022) Temporal patterns of weekly births and conceptions predicted by meteorology, seasonal variation, and lunar phases. Wien Klin Wochenschr 134(13-14):538-545. doi: 10.1007/s00508-022-02038-7.
- Gwet, KL (2008) Computing inter-rater reliability and its variance in the presence of high agreement. The British journal of mathematical and statistical psychology 61(1): 29–48. https://doi.org/10.1348/000711006X126600
- Khamis, H (2008) Measures of association: How to choose? Journal of Diagnostic Medical Sonography 24(3) 155–162. https://doi.org/10.1177/8756479308317006
- Kirkwood, BR & Sterne, JAC (2003) Essential Medical Statistics. 2nd Edition, Blackwell Science, Oxford.
- Ludbrook, J (2011) Is there still a place for Pearson’s chi-squared test and Fisher’s exact test in surgical research? ANZ Journal of Surgery 81(12): 923-6. https://doi.org/10.1111/j.1445-2197.2011.05906.x
- Mehta, CR & Patel, NR (1996) SPSS Exact Tests 7 for Windows. IL: SPSS.
- Ogle, DH (2016) Introductory fisheries analyses with R. USA: CRC. https://doi.org/10.1201/9781482235227
- Siegel, S & Castellan Jr., NJ (1988) Nonparametric statistics for the behavioral sciences. 2ª ed. NY: McGraw-Hill.
- Silveira, PSP & Siqueira, JO (2023) Better to be in agreement than in bad company: A critical analysis of many kappa-like tests. Behavior research methods 55(7): 3326–47. https://doi.org/10.3758/s13428-022-01950-0
- Sharpe, D. (2015) Your chi-square test is statistically significant: now what? Practical Assessment: Research & Evaluation 20(8).
- Sleigh, A et al. (1982) Comparison of filtration staining (Bell) and thick smear (Kato) for the detection of quantitation of Schistosoma mansoni eggs in faeces. Transactions of the Royal Society of Tropical Medicine and Hygiene 76(3):403-6. doi: 10.1016/0035-9203(82)90201-2.
- SPSS Base 10: Applications Guide (1999) IL: SPSS.
- Staboulidou, I et al. (2008) The influence of lunar cycle on frequency of birth, birth complications, neonatal outcome and the gender: A retrospective analysis. Acta Obstetricia et Gynecologica Scandinavica, 87: 875-879. https://doi.org/10.1080/00016340802233090
- Streiner, DL & Lin, E (1998) Life after chi-squared: An introduction to log-linear analysis. Canadian Journal of Psychiatry, 43(8) https://doi.org/10.1177/070674379804300809
- Wake R et al. (2010) The effect of the gravitation of the moon on frequency of births. Environ Health Insights. 4:65-9. doi: 10.4137/ehi.s5525.
- Wongpakaran N, Wongpakaran T, Wedding D, Gwet KL (2013) A comparison of Cohen’s Kappa and Gwet’s AC1 when calculating inter-rater reliability coefficients: a study conducted with personality disorder samples. BMC Med Res Methodol 13:61. https://doi.org/10.1186/1471-2288-13-61