Análise de Tabela de Contingência
MSP4112: Estatística Aplicada à Medicina
15 agosto 2025 20:53h
Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
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(effectsize, warn.conflicts=FALSE))
suppressMessages(library(epiR, warn.conflicts=FALSE))
suppressMessages(library(exact2x2, 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(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")
options(warn=0)Material
- HTML de Rmarkdown em
RPubs
Objetivos
- Identificar e sumariar variáveis nominais e ordinais em tabela de contingência.
- Identificar distribuição amostral \(\chi^2\) (qui-quadrado).
- Formular hipótese estatística na análise de dados qualitativos.
- Identificar, aplicar e interpretar os testes 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 tabelas de contingencia com variáveis ordinais.
- Analisar tabelas tridimensionais \(l\times c \times k\), com a inclusão de variável de estratificação.
- Realizar os procedimentos estatísticos em R.
Roteiro em R
Este roteiro serve para facilitar a localização das principais funções e procedimentos envolvidos neste texto, de acordo com o tipo de delineamento dos estudos.
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 exato de Fisher de OR bruto:
exact2x2::fisher.exact,coin::chisq_test(...,distribution="exact"),confintr::ci_cramersv - Teste qui-quadrado de Pearson:
chisq.test(simulate.p.value=TRUE,B=1e6,...),confintr::ci_cramersv - Teste kappa de Cohen:
epiR::epi.kappa
- Teste exato de Fisher de OR bruto:
- Duas variáveis nominais com pelo menos uma politômica
- Teste qui-quadrado de Pearson:
chisq.test(simulate.p.value=TRUE,B=1e6,...),coin::chisq_test - 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
- 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
- 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/2 de Gwet:
rrCAC::gwet.ac1.table,irrCAC::gwet.ac1.raw,eiras2x2::AC1
- Testes AC1/2 de Gwet:
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)
Testes qui-quadrado
Os testes baseados na estatística qui-quadrado de Pearson podem ser divididos em dois tipos:
Teste de aderência de uma variável nominal com duas ou mais categorias: teste de proporções;
Teste 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:
Independência das observações: cada observação da unidade observacional deve ser alocada em apenas uma categoria da variável nominal.
Siegel & Castellan Jr., 1988, p. 111
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:
Esta é 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 número de pessoas que escolheu apenas uma entre três marcas
diferentes está em Chocolate.rds, uma planilha organizada
da seguinte forma:
marca prefreq probH0
Dream Brown 12 0.3333333
Best Cocoa 23 0.3333333
Wonka 44 0.3333333O teste está implementado em TesteQuiQuadradoAderencia_Chocolate.R,
resultando:
# TesteQuiQuadradoAderencia_Chocolate.R
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# carrega a função chisqgof()
source("eiras.chisqgof.R")
# arquivo de dados no formato _rds_ com os dados e hipótese.
# primeira coluna é o fator
# segunda coluna com frequências observadas
# coluna **probH0** contém 1/3 em todas as linhas (frequências esperadas para a hipótese nula)
Dados <- readRDS("Chocolate.rds")
# executa o teste
fit <- chisqgof(Dados)
--------------------------------------------
Goodness of Fit, robust test (bootstrapping)
--------------------------------------------
----
Data
----
marca prefreq probH0
Dream Brown 12 0.3333333
Best Cocoa 23 0.3333333
Wonka 44 0.3333333
------------------------------------
Analysis of Statistical Significance
------------------------------------
X^2 = 20.07595 df = 2 p.value = 9.9999e-06
Chi-squared test for given probabilities with
simulated p-value (based on 1e+05 replicates)
Heuristic for significance
(rejection of H0 is expected if X^2/df > 2)
computed X^2/df = 10.03797
----------------------
Standardized residuals
----------------------
Dream Brown -3.421
Best Cocoa -0.796
Wonka 4.216
-----------------
Post hoc analysis
-----------------
Pairwise comparisons using exact binomial tests
data: preference
Dream Brown Best Cocoa
Best Cocoa 0.269 -
Wonka 6.3e-05 0.042
P value adjustment method: bonferroni
Populational inference, between:
- Best Cocoa e Dream Brown: there is no difference of proportion for marca
- Wonka e Dream Brown: there is difference of proportion for marca
- Wonka e Best Cocoa: there is difference of proportion for marca
----------------------------------
Analysis of Practical Significance
----------------------------------
Cramer V = 0.3564589
intermediate degree of non adherenceCom \(p \approx\) 0.001000%, rejeita-se \(H_0\) e, 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: genética
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?
A diferença com o exemplo anterior é que o arquivo de dados Flores.rds tem 4
colunas:
CorPetalas Frequencia Codominancia Epistasis
laranja 182 0.50 0.5625
vermelha 77 0.25 0.2500
amarela 61 0.25 0.1875 As duas primeiras são as cores das flores e os números observados de
plantas, respectivamente similares às marcas e número de pessoas em Chocolate.rds. 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 \(\frac{9}{16}, \frac{4}{16}, \frac{3}{16}\),
a \(H_0\) da epistasis recessiva. Por
isso, chamamos chisqgof 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:
# TesteQuiQuadradoAderencia_Flores.R
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")
Sys.setlocale("LC_ALL", "pt_BR.UTF-8")
source("eiras.chisqgof.R")
# Lê o arquivo de dados
Dados <- readRDS("Flores.rds")
# chama o teste estatístico duas vezes
cat("\n### Hipótese de codominância: ###\n")
chisqgof(Dados[,c(1,2,3)])
cat("\n### Hipótese de epistasis recessiva: ###\n")
chisqgof(Dados[,c(1,2,4)])
### Hipótese de codominância: ###
--------------------------------------------
Goodness of Fit, robust test (bootstrapping)
--------------------------------------------
----
Data
----
CorPetalas Frequencia Codominancia
laranja 182 0.50
vermelha 77 0.25
amarela 61 0.25
------------------------------------
Analysis of Statistical Significance
------------------------------------
X^2 = 7.65 df = 2 p.value = 0.02244978
Chi-squared test for given probabilities with
simulated p-value (based on 1e+05 replicates)
Heuristic for significance
(rejection of H0 is expected if X^2/df > 2)
computed X^2/df = 3.825
----------------------
Standardized residuals
----------------------
laranja 2.46
vermelha -0.387
amarela -2.453
-----------------
Post hoc analysis
-----------------
Pairwise comparisons using exact binomial tests
data: preference
laranja vermelha
vermelha 0.6 -
amarela 1.1e-14 1.7e-10
P value adjustment method: bonferroni
Populational inference, between:
- vermelha e laranja: there is no difference of proportion for CorPetalas
- amarela e laranja: there is difference of proportion for CorPetalas
- amarela e vermelha: there is difference of proportion for CorPetalas
----------------------------------
Analysis of Practical Significance
----------------------------------
Cramer V = 0.1093303
small degree of non adherence
### Hipótese de epistasis recessiva: ###
--------------------------------------------
Goodness of Fit, robust test (bootstrapping)
--------------------------------------------
----
Data
----
CorPetalas Frequencia Epistasis
laranja 182 0.5625
vermelha 77 0.2500
amarela 61 0.1875
------------------------------------
Analysis of Statistical Significance
------------------------------------
X^2 = 0.1513889 df = 2 p.value = 0.9328107
Chi-squared test for given probabilities with
simulated p-value (based on 1e+05 replicates)
Heuristic for significance
(rejection of H0 is expected if X^2/df > 2)
computed X^2/df = 0.07569444
----------------------
Standardized residuals
----------------------
laranja 0.225
vermelha -0.387
amarela 0.143
-----------------
Post hoc analysis
-----------------
This analysis is not required
(omnibus test has p=0.9328107).
----------------------------------
Analysis of Practical Significance
----------------------------------
Cramer V = 0.01538002
minimal degree of non adherenceNeste exemplo, buscamos não rejeitar \(H_0\). Estes resultados sugerem que a herança das cores das pétalas desta espécie é explicada por epistasis 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 epistasis 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. |
Relações entre 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: Gatos e Alergia
- gatos causam alergia …
- estresse
- mecanismo imune
- desinfetantes
- alergia e estresse
- alergia é a causa dos gatos
Exemplo 2: Vinho e Doença Coronariana
“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.”
Pharma Hoje (2017) Vinte alimentos nutritivos que ajudam a manter o bom funcionamento do organismo
- ambiente inseguro
- moradia precária
- má nutrição
- boa moradia
- segurança
- boa alimentação
- outros cuidados
- é o vinho?
- Coronary heart disease: what is the evidence?
“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.”
Traduzido de Pandit K (2017) (grifos nossos)
Mediterranean-style diet for the prevention of cardiovascular disease
“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.”
Traduzido de Rees K, Takeda A, Martin N, Ellis L, Wijesekara D, Vepa A, Das A, Hartley L, Stranges S (2019) (grifos nossos)
Exemplo 3: Tabagismo e Câncer de Pulmão
- The Cigarette Century
- How could one attribute lung cancer to cigarette smoking?
“The public impression of scientific and medical uncertainty that resulted was a crucial element in maintaining the market of current smokers as well as recruiting new ones. Industry literature, for example, frequently pointed to the fact that nonsmokers could also develop lung cancer. Therefore, they argued, how could one attribute lung cancer to cigarette smoking?”
- Medical knowledge is provisional
medical knowledge was always provisional and contingent. Just as drugs deemed “effective” do not work in every case, so too a cause of disease does not always result in disease. As Richard Doll would later explain, the epidemiologists had identified a cause of lung cancer (and other diseases), not the cause.
- Fisher was against the association between smoking and lung cancer
Another vocal critic of the lung cancer findings was Sir Ronald Fisher, the leading biometrician and geneticist in Great Britain during the first half of the twentieth century and a man deeply committed to bringing statistical analysis to genetics and agricultural experimentation. His 1925 book, Statistical Methods for Research Workers, quickly became a classic, leading to academic appointments at University College London and Cambridge University. Fisher’s critiques were similar to Berkson’s. The ethical impossibility of conducting a randomized experiment led him to question the results of the epidemiological studies. As a believer in genetic* notions of cancer causality, Fisher speculated that there was some constitutional factor that led individuals both to become smokers and to get lung cancer, even though smoking and lung cancer might not be causally related. Doll and Hill repeatedly rebutted this theory, returning to the critical question of how to account for the rise in lung cancers during the twentieth century if the disease was simply “constitutional.”
* 1923: One of the [A. Bradford] Hill’s innovations was the first randomized, double-blind clinical trial, designed to reduce investigator bias in the evaluation of clinical outcomes. […] This method, which drew on Fisher’s agricultural experimentation in genetics, became a critical new tool for evaluating medical treatments.
By the 1930s, the relationship of smoking to cancer was a topic of uresolved debate. It was the life insurance industry, which like the tobacco corporations grew by leaps and bounds in the first half of the twentieth century, that took the lead in understanding the effects of smoking on health.
By the early 1950s, however, it was abundantly clear that the evidence implicating cigarette smoking as a risk to health was now of a different order. First, the link between smoking and disease was categorical, outside the realm of individual clinical judgment.
They [tobacco industry] responded with a new and unprecedent public relations strategy. Its goals was to produce and sustain scientific skepticism and controversy in order to disrupt the emerging consensus on the harms of cigarette smoking. This strategy required intrusions into scientific process and procedure. […] So long as there appeared to be doubt, so long as the industry could assert “not proven”. […] The future of the cigarette would now depend on the successful production of a scientific controversy.
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:
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.
Siegel & Castellan Jr., 1988, p. 111
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{align} H_0&: \text{não existe associação entre uso do capacete e trauma craniano}\\ H_1&: \text{existe associação entre uso do capacete e trauma craniano}\\ \alpha &= 5\% \end{align} \]
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., se não existe associaçã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 do 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.006858
Two-sided 95% chi-squared confidence interval for the population
Cramer's V
Sample estimate: 0.09601
Confidence interval:
2.5% 97.5%
0.04415 0.16938 print(confintr::ci_cramersv(chi2,
probs = c(alpha/2, 1-alpha/2),
type = "bootstrap",
R = 1e4,
seed = 123,),
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 É mais fácil entender a interpretação deste teste, representando-o
graficamente (capacetetrauma_testeX2.R):
source("friendlycolor.R")
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Capacete +","Capacete -")
alpha <- 0.05
chi2 <- chisq.test(tabela, correct=FALSE)
# vetor com valores de 0 a 8
chi2.valor <- seq(from=0, to=9, by=0.01)
graus.liberdade <- as.numeric(chi2$parameter)
# vetor com probabilidades correspondentes
p <- dchisq(chi2.valor, df=graus.liberdade)
# exibe o gráfico
plot(chi2.valor, p, type="l", ylim=c(0, 1.2),
main="Teste X^2 (capacete e trauma), df=1",
xlab="qui-quadrado", ylab="Densidade")
# vetor com os valores acima de qui critico
chi2.critico <- qchisq(p=1-alpha, df=graus.liberdade)
chi2maioralfa <- chi2.valor[chi2.valor>chi2.critico]
# probabilidades correspondentes
pmaioralfa <- dchisq(chi2maioralfa, df=graus.liberdade)
# hachura sob a curva
chi2maioralfa <- c(min(chi2maioralfa), chi2maioralfa, max(chi2maioralfa))
pmaioralfa <- c(0, pmaioralfa, 0)
polygon(chi2maioralfa, pmaioralfa, col=friendlycolor(10), border = NA)
# reforca a linha
lines(chi2maioralfa[2:(length(chi2maioralfa)-1)],
pmaioralfa[2:(length(pmaioralfa)-1)],
col=friendlycolor(10), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2.critico, lty=2)
# ########## Area a direita de p ############
# vetor com os valores acima de qui calculado
chi2maiorp <- chi2.valor[chi2.valor>chi2$statistic]
# probabilidades correspondentes
pmaiorp <- dchisq(chi2maiorp, df=graus.liberdade)
# hachura sob a curva
chi2maiorp <- c(min(chi2maiorp), chi2maiorp, max(chi2maiorp))
pmaiorp <- c(0, pmaiorp, 0)
polygon(chi2maiorp, pmaiorp, col=friendlycolor(20), border = NA)
# reforca a linha
lines(chi2maiorp[2:(length(chi2maiorp)-1)],
pmaiorp[2:(length(pmaiorp)-1)],
col=friendlycolor(20), lwd=2)
# linha pontilhada, valor critico
abline(v=chi2$statistic, lty=3)
# legenda
legend("topright",
c(paste("Qui^2 crítico=", round(chi2.critico,5), sep=""),
paste("Qui^2 obs.=", round(chi2$statistic,5), sep=""),
paste("alfa=", alpha, sep=""),
paste("p=", round(chi2$p.value,5),sep="")),
col=c("black", "black", friendlycolor(10), friendlycolor(20)),
lty=c(2, 3, 1, 1),
lwd=c(1, 1, 5, 5),
bty="n",
cex=0.8) 
(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
set.seed(123)
chi2.robusto <- chisq.test(tabela,
simulate.p.value=TRUE,
B=1e5)
print(chi2.robusto)
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: tabela
X-squared = 7.3099, df = NA, p-value = 0.00763
Two-sided 95% chi-squared confidence interval for the population
Cramer's V
Sample estimate: 0.09601047
Confidence interval:
2.5% 97.5%
0.04415044 0.16937524 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 \times d}{b \times 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 função para o cálculo de \(\text{OR}\) implementado como parte do Teste Exato de Fisher. No exemplo de capacete e trauma, temos:
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 <- fisher.test(tabela) # Teste de OR robusto
print(ft) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508
Fisher's Exact Test for Count Data
data: tabela
p-value = 0.005688
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.2630993 0.8357342
sample estimates:
odds ratio
0.4817645
|
A função nativa Existe uma alternativa, implementada em
Para ver qual é o problema, usando outra tabela de 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 |
Decisão pelo intervalo de confiança
A estimativa pontual de odds ratio 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: O odds ratio Aliás, observe que são os mesmos valores obtidos anteriormente e guardados no objeto
|
|
Lembre-se de nunca decidir a respeito de odds ratio sem considerar o intervalo de confiança. Por isso, evite fazer o cálculo manualmente: use R! |
Risco relativo em epiR
Os cálculos feitos acima são implementados no pacote
epiR com:
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,outcome="as.columns")) Trauma + Trauma -
Capacete + 17 138
Capacete - 130 508
Outcome+ Outcome- Total Inc risk *
Exposure+ 17 138 155 10.97 (6.52 to 16.98)
Exposure- 130 508 638 20.38 (17.32 to 23.71)
Total 147 646 793 18.54 (15.89 to 21.42)
Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio 0.54 (0.34, 0.86)
Inc odds ratio 0.48 (0.28, 0.83)
Attrib risk in the exposed * -9.41 (-15.24, -3.58)
Attrib fraction in the exposed (%) -85.78 (-199.77, -17.56)
Attrib risk in the population * -1.84 (-5.97, 2.29)
Attrib fraction in the population (%) -9.92 (-10.71, -8.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 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 alteracas 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 (demo_epiR_2.R).
tabela <- as.table(matrix(c(17, 138, 130, 508), nrow = 2, byrow = TRUE))
colnames(tabela) <- c("Trauma +","Trauma -")
rownames(tabela) <- c("Uso de Capacete","Falta de Capacete")
tabela <- DescTools::Rev(tabela,margin=c(1))
print(tabela)
print(tb <- epiR::epi.2by2(dat=tabela,method="cohort.count",
units=100,outcome="as.columns")) Trauma + Trauma -
Falta de Capacete 130 508
Uso de Capacete 17 138
Outcome+ Outcome- Total Inc 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:
-------------------------------------------------------------------
Inc risk ratio 1.86 (1.16, 2.98)
Inc odds ratio 2.08 (1.21, 3.56)
Attrib risk in the exposed * 9.41 (3.58, 15.24)
Attrib fraction in the exposed (%) 46.17 (14.94, 66.64)
Attrib risk 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 (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}} \]
Portanto, caso epiR::epi.2by2 calculasse \(\text{RR}_-\), obteria (demo_epiR_3.R):
RRneg <- as.numeric(tb$massoc.summary$est[1]/tb$massoc.summary$est[2])
RR.ci <- sort(c(as.numeric(tb$massoc.summary$lower[1]/tb$massoc.summary$lower[2]),
as.numeric(tb$massoc.summary$upper[1]/tb$massoc.summary$upper[2])))
cat("RR- = ",round(RRneg,4),
", IC95(RR-) = [",round(RR.ci[1],4),",",round(RR.ci[2],4),"]", sep="")RR- = 0.8943, IC95(RR-) = [0.8378,0.9547]Riscos 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.
Compare com os resultados obtidos com epiR::epi.2by2: as
estimativas pontuais coincidem e os intervalos de confiança são muito
similares.
epiR::epi.2by2
Trauma + Trauma -
Capacete - 130 508
Capacete + 17 138
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
Codificando (conferindo):
1=Capacete -
0=Capacete +
1=Trauma +
0=Trauma - Trauma
Capacete 1 0
1 130 508
0 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 Modelo linear generalizado (GzLM): stats::glm
- 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
Regressão logística: OR
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]Regressão de Poisson: RR+ e RR-
Poisson Regression: RR+ e 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]RR- = 0.89, IC95(RR-) = [0.87, 0.93]Regressão log-binomial: RR+ e RR-
log-binomial Regression: RR+ e 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- = RR+/OR = 0.89IC95(RR-) = [0.81, 0.93]Saí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
O V de Cramer é uma medida do grau de correlação de Pearson absoluta 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 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))
print(try(rcompanion::cramerV(tcrxc,ci=TRUE,R=1e3)))
# 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.2758
V de Cramer = 0.3428
X-squared
"large"
(Rules: funder2019)
Cramer.V lower.ci upper.ci
1 0.3428 0.09266 0.7078
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.
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
source("eiras.show.MCSTAR.R")
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- 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)
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)
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="")
print(effectsize::interpret_cramers_v(V))
print(try(rcompanion::cramerV(TC,ci=TRUE,R=1e3)))
cat("\nTeste de razao de chances (OR) robusto:\n")
resft <- exact2x2::fisher.exact(TC,conf.level = 1-alfa) # Teste de OR robusto
print(resft)
# reabilita warnings
options(warn=0)New names:
• `` -> `...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.00069
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
X-squared
"large"
(Rules: funder2019)
Cramer.V lower.ci upper.ci
1 0.332 0.1535 0.5183
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{observado}} } \]
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{observado}} } \]
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):
o 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.”
Ludbrook, 2011
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
Em determinado estudo mulheres foram avaliadas quanto à sobrevida e o tabagismo, obtendo-se.
| 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")
# desabilita warnings
options(warn=-1)
TCtmp <- 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)
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="")
print(effectsize::interpret_cramers_v(V))
print(try(rcompanion::cramerV(TC,ci=TRUE,R=1e3)))
cat("\nTeste de razao de chances (OR):\n")
resft <- exact2x2::fisher.exact(TC,conf.level=1-alfa) # Teste de OR
print(resft)
# reabilita warnings
options(warn=0)
------------------------------------------
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.00288
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
X-squared
"very small"
(Rules: funder2019)
Cramer.V lower.ci upper.ci
1 0.08331 0.03232 0.1373
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).
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Dados
dados <- matrix(c(
"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
), ncol = 9, byrow = TRUE)
# Converter em data.frame
colunas <- c("Tabagista", "Desfecho", "18-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75-")
df <- as.data.frame(dados)
colnames(df) <- colunas
# Imprimir tabela formatada
tabela <- knitr::kable(df, align = "c", format = "html", table.attr = 'class="center"')
kableExtra::kable_styling(tabela, full_width = FALSE, position = "center", bootstrap_options = c("striped", "hover"))| 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 |
Appleton et al., 1996
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
# desabilita warnings
options(warn=-1)
# 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)")
# reabilita warnings
options(warn=0), , 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.01664
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, o 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} \]
Além disto, existe o \(\text{OR}\) generalizado.
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.
Os dados desta tabela são, originalmente, de Duncan (1979).
As categorias são as seguintes:
- Professional and high administrative
- Managerial and executive
- Inspectional, supervisory, and other nonmanuar (high grade)
- Inspectional, supervisory, and other nonmanual (lower grade)
- (V)(a) Routine grades of nonmanual
- (V)(b) Skilled manual
- Semiskilled manual
- Unskilled manual
Portanto, a variável de status ocupacional é, claramente, qualitativa 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
é:
# TCrXc_StatusOcupacional.R
library(datasets)
library(DescTools)
source ("eiras.cramerV.R")
source ("eiras.show.MCSTAR.R")
# desabilita warnings
options(warn=-1)
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
prmatrix(tcrxc)
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 = ", gama$gamma,
", IC",round((1-alfa)*100,1)," [",gama$CI[1],",",gama$CI[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 + ",gama$gamma,")/(1 - ",gama$gamma,") = ",ORg,sep="")
cat(", IC",round((1-alfa)*100,1)," [",ORg.LI, ",", ORg.LS,"]\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))
print(try(rcompanion::cramerV(tcrxc,ci=TRUE,R=1e3)))
# reabilita warnings
options(warn=0)
------------------------------------------
Dados:
------------------------------------------
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.4209053, IC95 [0.3912413,0.4505693]
OR bruto generalizado = (1 + 0.4209053)/(1 - 0.4209053) = 2.453667, IC95 [2.285374,2.640132]
------------------------------------------
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)
Cramer.V lower.ci upper.ci
1 0.2405 0.223 0.2668É 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 variáveis ordinais
A seguir é apresentado um exemplo com o uso de funções do 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 é:
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")))
ftable(Job)
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 = 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("(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))
print(try(rcompanion::cramerV(Job,ci=TRUE,R=1e3)))
options(warn=0)
Teste qui-quadrado de Pearson exato
(assumindo as variáveis como nominais):Warning in chisq.test(Job): Aproximação do qui-quadrado pode estar incorreta
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.7707
------------------------------------------
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.09326
alternative hypothesis: two.sided
------------------------------------------
Analise de significancia pratica:
------------------------------------------
(assumindo as variáveis como nominais):
V de Cramer = 0.1439
X-squared
"small"
(Rules: funder2019)
Cramer.V lower.ci upper.ci
1 0.1439 NA NAObserve 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
Para situações em que uma variável é ordinal e a outra é nominal, 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,
usando a planilha menarca.xlsx, que produz:
# Teste qui-quadrado de Cochran-Armitage.R
# Cochran-Armitage test for trend
# Menarca precoce & espessura da dobra cutânea do tríceps
# Fonte: Beckles et al. (1985) International Journal of Obesity 9:127-35,
# apud Kirkwood & Sterne (2003): Chapter 17: Chi-squared tests for 2x2 and
# larger contingency tables.
library(readxl)
library(DescTools)
# desabilita warnings
options(warn=-1)
# Tabela de Dancey&Reidy (2019, p. 275)
TCtmp <- read_excel("menarca.xlsx")
TC <- data.matrix(TCtmp)
rownames(TC) <- as.character(unlist(TCtmp[,1]))
TC <- TC[,-1]
remove(TCtmp)
cat("\n------------------------------------------\n")
cat("Dados:")
cat("\n------------------------------------------\n")
print(TC)
cat("\nTeste qui-quadrado de Pearson exato:\n")
res <- chisq.test(TC,simulate.p.value=TRUE,B=1e5)
res$parameter <- (nrow(TC)-1)*(ncol(TC)-1)
print(res)
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("Heuristica de significancia (rej. H0 se X^2/gl > 2): ", X2/df, "\n", sep="")
cat("\n------------------------------------------\n")
cat("Analise de significancia estatistica:")
cat("\n------------------------------------------\n")
alfa <- 0.05
cat("\nTeste qui-quadrado de Cochran-Armitage:\n")
res <- DescTools::CochranArmitageTest(TC)
print(res)
x2ca <- res$statistic^2
p <- res$p.value
cat("\nX-squared trend =",x2ca, ", df = 1, ", "p-value =", p,"\n")
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="")
print(effectsize::interpret_cramers_v(V))
print(try(rcompanion::cramerV(TC,ci=TRUE,R=1e3)))
# reabilita warnings
options(warn=0)
------------------------------------------
Dados:
------------------------------------------
menor12 maiorigual12
Pequena 15 156
Intermediaria 29 197
Grande 36 186
Teste qui-quadrado de Pearson exato:
Pearson's Chi-squared test with simulated p-value (based on 1e+05
replicates)
data: TC
X-squared = 4.7594, = 2, p-value = 0.09431
X^2 critico de 95% = 3.841459
Heuristica de significancia (rej. H0 se X^2/gl > 2): 2.379692
------------------------------------------
Analise de significancia estatistica:
------------------------------------------
Teste qui-quadrado de Cochran-Armitage:
Cochran-Armitage test for trend
data: TC
Z = 2.1783, dim = 3, p-value = 0.02938
alternative hypothesis: two.sided
X-squared trend = 4.744926 , df = 1, p-value = 0.02938481
------------------------------------------
Analise de significancia pratica:
------------------------------------------
V de Cramer = 0.0877
X-squared
"very small"
(Rules: funder2019)
Cramer.V lower.ci upper.ci
1 0.08769 0.02831 0.1643Compare 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
Para desfecho com três ou mais categorias (os desfechos são tratados
como variáveis nominais) 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:
library(coin)
library(rcompanion)
# 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)
# 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"= c(1, 2, 3, 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)
}
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.7701
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 variáveis ordinais estratificado por variável nominal
Esta é uma generalização do teste de Mantel-Haenszel estratificado, aqui aproveitando a informação de ordem das duas variáveis principais, analisadas em estratos (apenas esta variável é nominal).
São usadas as funções coin::lbl_test e
coin::cmh_test.
Agresti, 2002.
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"))))
ftable(Satisfaction)
print(coin::lbl_test(Satisfaction,
distribution = coin::approximate(nresample = 1e6)))
Approximative 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.01104
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 é regressão!
“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.”
Field, 2012, p. 833.
Streiner & Lin, 1998
Por mais de quatro décadas, a General Social Survey (GSS) 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.
The General Social Survey (https://gss.norc.org/)
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
Estes procedimentos avaliam concordância entre dois métodos que avaliam dicotomicamente as mesmas unidades observacionais.
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 \(AC_1\) 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 de R de código aberto 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). O kappa de Cohen mostrou um comportamento prejudicado similar ao r de Pearson, Q de Yule e Y de Yule. O \(\pi\) de Scott e o B de Shankar e Bangdiwala parecem avaliar melhor as situações de discordância do que de concordância entre avaliadores. O alfa de Krippendorff emula, sem qualquer vantagem, o \(\pi\) de Scott em casos com variáveis nominais e dois avaliadores. O F1 de Dice e o 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 o kappa de Cohen é uma medida de associação e o qui-quadrado de McNemar não avalia nem associação nem concordância; os únicos dois estimadores de concordância autênticos são o G de Holley e Guilford e o 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.”
\(AC_1\) de Gwet
O Agreement Coefficient do tipo 1 assume que as duas variáveis são nominais dicotômicas (tabela 2 \(\times\) 2).
Por exemplo, 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.
Sleigh et al., 1982
A hipótese nula é
\[ \begin{align} H_0: \;&\text{ausência de concordância}\\ H_1: \;&\text{presença de concordância} \end{align} \]
O exemplo está implementado em TC2x2_Gwet_concord.R:
library(irrCAC)
source("eiras.matrix2dataframe.R")
# sink("TC2x2_Gwet.txt")
# pdf("TC2x2_Gwet.pdf")
# 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 ="", xlab ="", ylab="",
label = FALSE, show.margins = FALSE))
# 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 = ",AC1.f, ", IC95 ",AC1.ic,", p=",AC1.p,"\n",sep="")
# dev.off()
# sink()
BellxKK P N
P 184 54
N 14 63
NULL
AC1 de Gwet
AC1 = 0.6237683, IC95 (0.536,0.712), p=0e+00Neste exemplo, \(H_0\) é rejeitada e, portanto, conclui-se que os dois métodos para a detecção de ovos são 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- {{(b+c)^2}\over{2}} }\over{ a^2 + d^2 + {{(b+c)^2}\over{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, convenientemente, fique entre 0 e 1.
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, \(\text{AC}_1\) 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 \(AC_2\) de Gwet em tabela \(k\times k\) para dois avaliadores
Aplicável para a comparação de dois métodos de avaliação com desfecho ordinal politômico, como por exemplo na avaliação de estadiamento tumoral submetido à apreciação de dois patologistas.
Agresti (1990), apud Mehta & Patel, 1996, 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
Formulamos:
\[ \begin{cases} H_0: &\text{Discordância}\\ H_1: &\text{Concordância} \end{cases} \]
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 MHOrdinal.R
# 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 slides de diferentes mulheres.
# N=Negativo, HEA=Hiperplasia escamosa atipica, CIS=carcinoma in situ
# CE=Carcinoma escamoso, CI=Carcinoma invasivo
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 ="", xlab ="", ylab="",
label = FALSE,
show.margins = FALSE))
cat("Teste AC2 de Gwet para duas variáveis ordinais\n")
cat("H0: AC2 = 0 vs H1: AC2 != 0\n",sep="")
print(irrCAC::gwet.ac1.table(TC,
irrCAC::linear.weights(1:nrow(TC)),
conflev = 0.95))
print(vcdExtra::GKgamma(TC,
level = 0.95)) 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} \]
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))
cat("\n")
print(irr::kappam.fleiss(diagnoses, detail=TRUE)) # Fleiss' and category-wise Kappa 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 ou para cada diagnóstico em particular.
ICC: intraclass correlation
Este coeficiente de correlação intraclasse é aplicável quando a variável de desfecho é intervalar e precisamos verificar a concordância entre vários avaliadores.
Um exemplo bastante utilizado é o treinamento que o IOC (International Olympic Commitee) utiliza para que os juizes das competições de ginástica sejam concordantes.
SPSS Base 10: Applications Guide, 1999.
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 observacionais. 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} \]
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.
library(psych)
library(haven)
library(irr)
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 = ",mean(mc,na.rm=TRUE)) 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.8971429O modelo Two-Way Random (parâmetro
model="twoway") precisa ser utilizado e significa que os
juízes e unidades observacionais (avaliados) 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.
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 um.
\[ \begin{cases} H_0&: P(+ \to \, – ) = P(– \, \to +) \\ H_1&: P(+ \to \, – ) \ne P(– \, \to +) \end{cases} \]
Um exemplo é dado por Siegel & Castellan (1988, p. 77). Neste exemplo verificou-se a intenção de votos antes e depois de um debate eleitoral.
Siegel & Castellan Jr., 1988
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 (TC2x2_MN_concord2.R):
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(coin::mh_test(as.table(TC), distribution = "exact")) # robusto
print(mcnemar.test(TC,correct=FALSE)) # classico
print(exact2x2::mcnemar.exact(TC))
print(gplots::balloonplot(t(as.table(TC)), main ="",
xlab ="", ylab="",
label = FALSE,
show.margins = FALSE))
TVDebate Carter Reagan
Carter 28 13
Reagan 7 27
Exact Marginal Homogeneity Test
data: response by
conditions (Var1, Var2)
stratified by block
chi-squared = 1.8, p-value = 0.2632
McNemar's Chi-squared test
data: TC
McNemar's chi-squared = 1.8, df = 1, p-value = 0.1797
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 
NULLComo \(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
- 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
- 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.
- 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