Capítulo 18 Estatística em Epidemiologia
18.1 Pacotes necessários
pacman::p_load(BiocManager,
car,
cowplot,
DescTools,
dplyr,
epiR,
epitools,
ggplot2,
gmodels,
kableExtra,
MASS,
MKmisc,
mlbench,
pROC,
readxl,
survival,
survminer,
vcd)O pacote limma também vai ser usado indiretamente neste capítulo e a sua instalação é feita através do BiocManager da seguinte maneira:
## 'getOption("repos")' replaces Bioconductor standard repositories, see
## 'help("repositories", package = "BiocManager")' for details.
## Replacement repositories:
## CRAN: http://cran.rstudio.com/## Bioconductor version 3.19 (BiocManager 1.30.23), R 4.4.1 (2024-06-14 ucrt)## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'limma'## Installation paths not writeable, unable to update packages
## path: C:/Program Files/R/R-4.4.1/library
## packages:
## boot, foreign, MASS, nlme, survival## Old packages: 'BiocManager', 'broom.helpers', 'corrplot', 'cpp11', 'curl',
## 'data.table', 'datawizard', 'DescTools', 'digest', 'expm', 'gdtools',
## 'ggrepel', 'HSAUR3', 'insight', 'matrixStats', 'microbenchmark', 'minqa',
## 'mvtnorm', 'openssl', 'parameters', 'pcaPP', 'performance', 'ragg',
## 'RcppEigen', 'RcppParallel', 'rjags', 'rmarkdown', 'robust', 'robustbase',
## 'rrcov', 'servr', 'sf', 'wk', 'xfun'18.2 Raciocínio bayesiano no diagnóstico médico
O processo diagnóstico é o centro da atenção da atividade médica na busca de reduzir as incertezas e reconhecer a que classe pertence determinado paciente. Portanto, é extremamente importante saber quão bem os testes diagnósticos podem prever que um indivíduo é portador de certa condição ou doença. Entende-se aqui como teste diagnóstico todo o processo diagnótico, desde o exame clínico até o mais sofisticado exame de imagem ou laboratorial. A ideia é saber como o teste diagnóstico se comporta para separar um “doente” e um “não doente”; qual a sua validade neste processo?
Deve-se sempre ter em mente que o estabelecimento do diagnóstico é um processo imperfeito que resulta em uma probabilidade ao invés de uma certeza de estar correto. Ou seja, cada vez mais os médicos têm que aplicar as leis da probabilidade na avaliação de testes diagnósticos e sinais clínicos.
A abordagem bayesiana denomina de probabilidade a priori a probabilidade estabelecida inicialmente, baseada apenas na experiência do médico, em seu conhecimento em relação a doença suspeitada. Diante de uma evidência de doença, pode ser solictado um teste diagnóstico. Quando ele recebe um teste positivo para uma doença, a probabilidade muda, passa a ser uma probabilidade condicional, probabilidade da doença dado que o teste é positivo, denominada probabilidade a posteriori.
Um teste que define corretamente quem é doente e quem não é doente é denominado de padrão-ouro ou padrão de referência. Algumas vezes, o teste padrão de referência é simples e barato. Outras vezes, é caro, difícil de obter, tecnicamente complexo, arriscado ou pouco prático. Inclusive, pode não hver padrão-ouro. Em função dessas limitações, outros testes são usados e, como consequência, podem ocorrer erros. Em outras palavras, no processo diagnóstico podem ocorrer falsos positivos e falsos negativos.
Esta incerteza, na utilização de testes diagnósticos, gera a necessidade de o médico conferir a probabilidade de falsos positivos e falsos negativos na elaboração de um diagnóstico ao receber o resultado positivo ou negativo de um exame. Uma maneira simples de mostrar as relações de um teste diagnóstico e o verdadeiro diagnóstico, é mostrada na tabela de contingência \(2\times2\) (Figura 18.1).
Figura18.1: Falsos positivos e falsos negativos
18.2.1 Sensibilidade e Especificidade
As estatísticas mais utilizadas para descrever a validade dos testes de diagnóstico em contextos clínicos são sensibilidade e a especificidade.
Sensibilidade é a habilidade do teste em identificar corretamente quem tem a doença. É a taxa de verdadeiros positivos (VP) de um teste e corresponde a probabilidade de um indivíduo com a doença ter um teste positivo.
Um teste sensível raramente deixará passar pessoas que tenham a doença. Testes com sensibilidade alta são úteis para excluir a presença de uma doença. Isto é, um teste negativo exclui virtualmente a possibilidade de o paciente ter a doença de interesse, pois tem pouca probabilidade de produzir resultados falsos negativos. Isto pode ser lembrado pelo mnemônico SnNout, do inglês: High Sensivity, a Negative result rules out the diagnosis (154).
Especificidade é a habilidade do teste em identificar corretamente quem não tem a doença. É a taxa de verdadeiros negativos (VN) de um teste e corresponde a probabilidade de um indivíduo sem a doença ter um teste negativo. Um teste específico raramente classificará de forma errônea indivíduos sendo portadores da doença quando eles não são. Os testes muito específicos são usados para confirmar a presença da doença. Se o teste é altamente específico, um teste positivo sugere fortemente a presença da doença de interesse.
De forma similar que a sensibilidade pode-se usar o mnemônico SpPin, do inglês: High Specificity, a Positive result rules in the diagnosis (154).
Estas estatísticas de diagnóstico podem ser calculadas a partir das equações, cujas letras representam as caselas da tabela \(2 \times 2\), acima;
\[ Sensibilidade = \frac {a}{\left (a + c\right )} \quad \quad Especificidade = \frac {d}{\left (b + d\right )} \]
A taxa de falsos negativos (TFN) é a proporção de indivíduos que têm a doença e que têm um resultado de teste negativo e a taxa de falsos positivos (TFP) é a proporção de pacientes que não possuem a doença e que apresentam resultados positivos. Podem ser expressas pelas equações:
\[ TFN= \frac {c}{\left (a + c\right )} \quad ou \quad \left (1 - sensibilidade\right) \]
\[ TFP= \frac {b}{\left (b + d\right )} \quad ou \quad \left (1 - especifcidade\right) \]
Idealmente, um teste de diagnóstico deveria ter altos níveis de sensibilidade e especificidade. No entanto, isso não é possível, pois existe um balanço entre sensibilidade e especificidade. À medida que a especificidade aumenta, a sensibilidade diminui e vice-versa. As curvas ROC, que serão discutidas mais adiante neste capítulo, podem ser usadas para identificar um ponto de corte em uma medição contínua que maximize a sensibilidade e a especificidade.
Quando um clínico tem um paciente cujo teste apresentou resultado positivo, a pergunta mais importante é a seguinte: dado que o teste é positivo, qual é a probabilidade de o paciente ter a doença? A sensibilidade do teste não responde a este questionamento, mas sim a probabilidade de um resultado positivo, dado que o paciente tem a doença (155).
18.2.1.1 Exemplo
O conjunto de dados dadosApendicite.xlsx contém informações de 156 pacientes que realizaram ultrassonografia abdominal para o diagnóstico de apendicite aguda.Para obter arquivo, clique aqui e salve o mesmo em seu diretório de trabalho.
Foram avaliados pacientes com diagnóstico clínico de apendicite aguda, submetidos à ultrassonografia abdominal e apendicectomia laparoscópica, acompanhado de estudo anatomopatológico dos apêndices extirpados (156). Será avaliado o teste diagnóstico usado.
Leitura e observação do conjunto de dados
Será usado a função read_excel()do pacote readxl e a função glimpse() do pacote dplyr:
## Rows: 156
## Columns: 3
## $ id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ apendicite <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ eco <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…As variáveis apendicite e eco foram exibidas como variáveis numéricas e serão transformadas em fatores:
dados$apendicite <- factor(dados$apendicite,
levels = c(1,2),
labels = c("Presente",
"Ausente"))
dados$eco <- factor(dados$eco,
levels = c(1,2),
labels = c("Positivo",
"Negativo"))Construção de uma tabela de contingência \(2\times2\)
## Apendicite
## Eco Presente Ausente
## Positivo 85 7
## Negativo 46 18Cálculo da sensibilidade e da especificidade
Pode-se usar a função epi.tests() do pacote epiR (157) que calcula, junto com os intervalos de confiança, a prevalência aparente e verdadeira, sensibilidade, especificidade, valores preditivos positivos e negativos e razões de probabilidade positivas e negativas a partir de dados de contagem fornecidos em uma tabela \(2\times2\). Utiliza os argumentos
- dat dados sob a forma de vetor ou matriz
- conf.level magnitude do intervalode confiança, entre 0 e 1.
Os resultados serão atribuídos a um objeto de nome diag:
## Outcome + Outcome - Total
## Test + 85 7 92
## Test - 46 18 64
## Total 131 25 156
##
## Point estimates and 95% CIs:
## --------------------------------------------------------------
## Apparent prevalence * 0.59 (0.51, 0.67)
## True prevalence * 0.84 (0.77, 0.89)
## Sensitivity * 0.65 (0.56, 0.73)
## Specificity * 0.72 (0.51, 0.88)
## Positive predictive value * 0.92 (0.85, 0.97)
## Negative predictive value * 0.28 (0.18, 0.41)
## Positive likelihood ratio 2.32 (1.22, 4.40)
## Negative likelihood ratio 0.49 (0.35, 0.68)
## False T+ proportion for true D- * 0.28 (0.12, 0.49)
## False T- proportion for true D+ * 0.35 (0.27, 0.44)
## False T+ proportion for T+ * 0.08 (0.03, 0.15)
## False T- proportion for T- * 0.72 (0.59, 0.82)
## Correctly classified proportion * 0.66 (0.58, 0.73)
## --------------------------------------------------------------
## * Exact CIsAssim, a sensibilidade é igual a 65% (IC95%: 56 – 73%) e a especificidade é igual a 72% (IC95%: 51 – 88%). Isto significa que um indivíduo com apendicite aguda tem 65% de probabilidade de ter uma ecografia alterada; um indivíduo sem apendicite aguda tem 72% de probabilidade de ter uma ecografia normal. O objetivo do teste de diagnóstico é usá-lo para fazer um diagnóstico, então há necessidade de saber a probabilidade que o teste fornece para um diagnóstico correto. A sensibilidade e a especificidade não fornecem esta informação. Para atingir esse objetivo, usa-se o valor preditivo (158).
18.2.2 Valor Preditivo
O propósito de um teste diagnóstico é usar seus resultados para fazer um diagnóstico, portanto, é necessário conhecer a probabilidade de que o resultado do teste forneça o diagnóstico correto (158).
Os valores preditivos positivo e negativo descrevem a probabilidade de um paciente ter doença, uma vez que os resultados de seus testes são conhecidos.
O valor preditivo positivo (VPP) de um teste é definido como a proporção de pessoas com um resultado de teste positivo que realmente têm a doença.
O valor preditivo negativo (VPN) é a proporção de pacientes com resultados de teste negativos que não têm doença.
Como a sensibilidade e a especificidade, estas estatísticas de diagnóstico também podem ser calculadas a partir da tabela \(2\times2\), mostrada no início:
\[ VPP = \frac {a}{\left (a + b\right )} \quad \quad VPN = \frac {d}{\left (c + d\right )} \]
Observando os resultados anteriores da função epi.tests(), verifica-se que 92% (85/92) dos indivíduos que tiveram teste positivo (ultrassonografia alterada) tinham doença (apendicite aguda).
Isso significa que seu VPP é igual a 92% (IC95%: 18 – 41%), ou dito de outra forma, uma pessoa com ultrassonografia positiva tem 92% de probabilidade de ter a apendicite aguda. O VPP é também conhecido como probabilidade pós-teste de doença dado um teste positivo.
Dos 64 pacientes que tiveram ultrassonografia sem alterações, 18 não apresentaram apendicite aguda, portanto, um VPN de 28% (IC95%: 56 – 73%). Isso significa que uma pessoa quem tem um teste negativo tem 28,1% de probabilidade de não ter apendicite aguda.
Entretanto, essas proporções são de validade limitada. Os valores preditivos de um teste, na prática clínica, dependem criticamente da prevalência da anormalidade nos pacientes testados. No estudo, a prevalência de apendicite aguda é igual a
\[ \frac {total\ de\ casos\ de \ apendicite \ aguda}{total\ de\ casos\ no\ estudo} = \frac {131}{156} = 0,84\ ou\ 84\% \left(IC_{95\%}:77\ a\ 89\%\right) \]
Levando-se em consideração que a prevalência de apendicite aguda na população é de 7% (159), mantendo a sensibilidade (64%) e a especificidade (72%) da ultrassonografia, entre 156 pacientes, selecionados aleatoriamente, se esperaria encontrar aproximadamente 11 casos (7% de 156) de apendicite aguda. Para facilitar a compreensão, observe a a tabela \(2\times2\) (Figura 18.2):
Figura18.2: Prevalencia e valor preditivo
O VPP e o VPN são iguais a:
## [1] 14.6## [1] 96.3Ao se comparar o VPP obtido, agora, com o VPP do estudo, observa-se que o mesmo diminuiu bastante, de 92% para 14,6%. O contrário ocorre com a VPN que aumenta substancialmente de 28% para 96,3%, mostrando claramente a influência da prevalência.
Se a prevalência diminui, o VPP diminui e o VPN aumenta. Portanto, será errado aplicar diretamente os valores preditivos publicados de um teste ao seu pacciente, quando a prevalência da doença em sua população for diferente da prevalência da doença na população em que o estudo publicado foi realizado. Um teste pode ser útil em um lugar e não ter validade em outro onde a prevalência é muito baixa.
Pode-se chegar aos mesmos resultados, usando as equações:
\[ VPP =\frac{sens \times prev}{\left(sens \times prev\right) + \left [\left (1- espec\right) \times \left (1- prev\right)\right ]} \]
\[ VPN =\frac{espec\times \left (1- prev\right)}{\left[\left (1 - sens \right)\times prev\right]+\left[espec\times \left (1 - prev\right)\right]} \]
A prevalência pode ser interpretada como a probabilidade antes da realização do teste, conhecida como probabilidade pré-teste. A diferença entre as probabilidades pré e pós-teste é uma forma de avaliar a utilidade do teste. Esta diferença pode ser mensurada pela razão de probabilidade (likelihood ratio).
18.2.3 Razão de Probabilidade
A Razão de Probabilidades (likelihood ratio) é uma forma alternativa de descrever o desempenho de um teste diagnóstico. Alguns autores a denominam de razão de verossimilhança 59.
A razão de probabilidades para um resultado de teste é definida como a razão entre a probabilidade de observar aquele resultado em indivíduos com a doença em questão e a probabilidade desse resultado em indivíduos sem a doença (160).
Razões de probabilidade são, clinicamente, mais úteis do que sensibilidade e especificidade. Fornecem um resumo de quantas vezes mais (ou menos) a probabilidade de os indivíduos com a doença apresentarem aquele resultado específico do que os indivíduos sem a doença, e também podem ser usados para calcular a probabilidade de doença para pacientes individuais (161). Cada vez mais as razões de probabilidade estão se tornando populares para relatar a utilidade dos testes de diagnóstico.
Quando os resultados do teste são relatados como sendo positivos ou negativos, dois tipos de razões de probabilidades podem ser descritos, a razão de probabilidades para um teste positivo (denotada LR +) e a razão de probabilidades para um teste negativo (denotada LR−).
A razão de probabilidades para um teste positivo é definida como a probabilidade de um indivíduo com doença ter um teste positivo dividida pela probabilidade de um indivíduo sem doença ter um teste positivo. A fórmula para calcular LR + é
Ou seja,
\[ LR(+)=\frac{sensibilidade}{1 - especificidade} \]
Razão de probabilidades positiva maior que 1 significa que um teste positivo tem mais probabilidade de ocorrer em pessoas com a doença do que em pessoas sem a doença. De um modo geral, para os indivíduos que apresentam um resultado positivo, LR (+) > 10 aumenta significativamente a probabilidade de doença (“confirma” a doença), enquanto LR (+) < 0,1, virtualmente, exclui a probabilidade de uma pessoa ter a doença (162).
Usando os dados da do objeto diag, obtido com a função epi.tests() do pacote epiR, tem-se que a LR (+) da ultrassonografia para o diagnóstico de apendicite aguda é igual 2.32 (IC95%: 1,22 – 4,40). Significa que uma pessoa com apendicite aguda tem cerca de 2,32 vezes mais probabilidade de ter um teste positivo do que uma pessoa que não tem a doença.
A razão de probabilidade negativa é definida como a probabilidade de um indivíduo com doença ter um teste negativo dividido pela probabilidade de um indivíduo sem doença ter um teste negativo. A fórmula para calcular a LR− é:
Ou seja,
\[ LR(-)=\frac{sensibilidade}{1-especificidade} \]
Razão de probabilidade negativa menor que 1 significa que um teste negativo é menos provável de ocorrer em pessoas com a doença do que em pessoas sem a doença. Um LR muito baixo (abaixo de 0,1) praticamente exclui a chance de que uma pessoa tenha a doença (162).
Voltando aos dados anteriores, a LR (-) para a ultrassonografia é igual a 0.49 (IC95%: 0.35 - 0.68). Significa que a probabilidade de ter um teste negativo para indivíduos com doença A é 0,49 vezes ou cerca de metade daqueles sem a doença. Dito de outra forma, os indivíduos sem a doença têm cerca o dobro probabilidade de ter um teste negativo do que os indivíduos com a doença.
18.2.3.1 Estimando a probabilidade de doença
Uma grande vantagem das razões de probabilidade é que elas podem ser usadas para ajudar o médico a adaptar a sensibilidade e a especificidade dos testes aos pacientes individuais. Ao se atender um paciente em uma clínica, pode-se decidir realizar um teste específico, após uma anamnese e um exame físico. A decisão de fazer o teste baseia-se nos sintomas e sinais do paciente e na experiência pessoal. Existe suspeita de um determinado diagnóstico e o objetivo é excluir ou confirmar esse diagnóstico. Antes de solicitar o teste, geralmente existe uma estimativa aproximada da probabilidade do paciente de ter essa doença, conhecida como probabilidade pré-teste ou a priori, que geralmente é estimada com base na experiência pessoal do médico, dados de prevalência local e publicações científicas.
A razão mais importante pela qual um teste é realizado é tentar modificar a probabilidade de doença. Um teste positivo pode aumentar a probabilidade pós-teste e um teste negativo pode reduzir essa probabilidade. A probabilidade pós-teste de doença é o que mais interessa aos médicos e pacientes, pois isso pode ajudar a decidir se devem confirmar, descartar um diagnóstico ou realizar outros testes.
Os resultados dos testes clínicos são geralmente usados não para fazer ou excluir categoricamente um diagnóstico, mas para modificar a probabilidade do pré-teste a fim de gerar a probabilidade do pós-teste. O teorema de Bayes é uma relação matemática que permite estimar a probabilidade pós-teste.
Para se compreender este conceito, é importante entender a diferença entre probabilidade e odds (163).
Probabilidade é a proporção de pessoas que apresentam uma determinada característica (teste positivo, sinal clínico). Odds (chance) representa a razão entre duas características complementares, ou seja, a probabilidade de um evento dividido pela probabilidade do não evento (1 – evento). Ambos contêm as mesmas informações de maneiras diferentes. Por exemplo, usando os dados da tabela tab_ap, , verifica-se que a probabilidade (p) de uma ultrassonografia positiva para apendicite aguda é igual
## [1] 0.5897436e que o odds da ultrassonografia positiva60 é
## [1] 1.4375Para transformar a odds em probabilidades e vice-versa, procede-se da seguinte maneira:
\[ p=\frac{odds}{1+odds} \] Voltando ao exemplo:
## [1] 0.5897436e
\[ odds=\frac{p}{1-p} \]
## [1] 1.4375Pelo teorema de Bayes, sabendo-se a probabilidade a priori ou probabilidade pré-teste, é possível obter a probabilidade pós-teste ou a posteriori, usando a razão de probabilidades.
Para atingir este objetivo, basta, inicialmente, multiplicar o odds pré-teste pela razão de probabilidades:
\[ odds_{pos} = odds_{pre \times LR} \]
Após, para encontrar a probabilidade pós-teste, basta converter o odds pós-teste em probabilidade:
\[ p_{pos} = \frac{odds_{pos}}{1-odds_{pos}} \]
No exemplo da tab_ap, foi verificado que o LR (+) é igual a 2,32 e a prevalência de apendicite aguda é em torno de 7% pode-se prever a probabilidade de haver apendicite aguda, diante de uma ultrassonografia alterada:
prev <- 0.07
LR <- 2.32
odds_pre <- 0.07/(1 -0.07)
odds_pos <- odds_pre * LR
p_pos <- odds_pos/(odds_pos +1)
round(p_pos, 3)## [1] 0.149Ou, em outras palavras, diante de um teste positivo, a probabilidade de o paciente ter apendicite aguda passa de 7% antes do teste para praticamente 15%!
Estes cálculos podem ser simplificados, utilizando o nomograma de Fagan (164), extremamente fácil de se usar (165), pois basta unir a probabilidade pré-teste ao LR que a reta apontará para a probabilidade pós-teste (Figura 18.3).
Figura18.3: Nomograma de Fagan
18.2.4 Curva ROC
Nem sempre o resultado de um teste é dicotômico (positivo/negativo). Com frequência, trabalha-se com variáveis contínuas (pressão arterial, glicemia, dosagem do sódio, dosagens hormonais, etc.). Neste caso, não há um resultado “positivo” ou “negativo”. Um “ponto de corte” precisa ser criado, para definir quem será considerado positivo ou negativo.
A escolha do ponto de corte depende das consequências de um resultado falso positivo ou de um falso negativo. Falsos positivos estão associados com custos (emocional ou financeiro) e com a dificuldade de “desrotular” alguém que recebeu o rótulo de “positivo”. Resultados falsos negativos podem “tranquilizar” pessoas doentes que não são seguidas ou tratadas precocemente.
A distribuição dos níveis glicêmicos em diabéticos e não diabéticos não tem um ponto de corte bem nítido. As duas populações se sobrepõem (Figura 18.4), gerando falso positivos ou falso negativos, dependendo do ponto de corte escolhido (163).
Figura18.4: Populações de indivíduos normais (curva em azul) e diabétticos (curva em vermelho
Suponha que ao se examinar uma população fosse escolhido o ponto de corte de 80mg/dL, haveria um aumento no número de indivíduos com teste positivo com uma taxa de falsos positivos elevada, diminuindo a especificidade do teste. Se, por outro lado, o ponto de corte fosse elevado para 200mg/dL, o número de falsos negativos teria um grande aumento, reduzindo a sensibilidade. Esta oscilação entre a sensibilidade e a especificidade ocorre pelo fato de a localização do ponto de corte ser uma decisão arbitrária num contínuo entre o normal e anormal.
Ao se escolher um ponto de corte deve-se fazer um balanço entre a sensibilidade e a especificidade, levando em conta as consequências da escolha. Por exemplo, a triagem para fenilcetonúria em recém-nascidos valoriza a sensibilidade em vez de especificidade; o custo da perda de um caso é alto, pois existe tratamento eficaz. Uma desvantagem é que ocorre um grande número de testes falso positivos que causam angústia e a realização de mais testes.
Em contraste, a triagem para o câncer de mama deve favorecer a especificidade sobre a sensibilidade, uma vez que uma avaliação mais aprofundada daquelas com teste positivo, implica em biopsias dispendiosas e invasivas.
As curvas ROC (Receiver Operating Characteristic) são uma ferramenta inestimável para encontrar o ponto de corte em uma medida com distribuição contínua que melhor prediz se uma condição está presente, por exemplo, se pacientes são positivos ou negativos para a presença de uma doença (166). As curvas ROC são usadas para encontrar um ponto de corte que separa um resultado de teste “normal” de um “anormal” quando o resultado do teste é uma medida contínua. As curvas ROC são traçadas calculando a sensibilidade e a especificidade do teste na predição do diagnóstico para cada valor da medida. A curva permite determinar um ponto de corte para a medição que maximiza a taxa de verdadeiros positivos (sensibilidade) e minimiza a taxa de falsos positivos (1 – especificidade) e, portanto, maximiza a razão de probabilidades (likelihood ratio).
18.2.4.1 Exemplo
O conjunto de dados dadosTestes.xlsx contém informações para os resultados hipotéticos de três testes bioquímicos diferentes e uma variável (doença) que indica se foi confirmada a doença (padrão-ouro). Para obter arquivo, clique aqui e salve o mesmo em seu diretório de trabalho.
Leitura e observação dos dados
Como é um arquivo em Excel, a leitura será realizada pela função read_excel() do pacote readxl:
## tibble [145 × 5] (S3: tbl_df/tbl/data.frame)
## $ id : num [1:145] 1 2 3 4 5 6 7 8 9 10 ...
## $ teste1: num [1:145] 25 2.2 46.2 9.9 46.5 36.1 34.8 44.9 36.9 7.1 ...
## $ teste2: num [1:145] 25 2.2 15.6 20.4 15.7 35.7 34.8 55.4 36.9 7.1 ...
## $ teste3: num [1:145] 15 2.2 25 20.4 15.7 36.1 24 55.4 36.9 7.1 ...
## $ doenca: num [1:145] 2 2 1 1 2 2 2 1 2 2 ...A variável doença será transformada em fator:
As curvas ROC são usadas para avaliar qual teste é mais útil para prever quais pacientes serão positivos para a doença. A hipótese nula é que a área sob a curva ROC é igual a 0,5, ou seja, a habilidade do teste para identificar casos positivos e negativos é a esperada por acaso.
A Figura 18.5 mostra a quantidade de sobreposição na distribuição da medição dos testes bioquímicos contínuos em ambos os grupos doença positiva e doença negativa. No Teste 1, a sobreposição é completa e não haverá um ponto de corte que separe efetivamente os dois grupos. Nos Testes 2 e 3, há uma maior separação das medidas de teste entre os grupos, particularmente para Teste 3.
Figura18.5: Resultado do teste vs doença.
18.2.4.2 Construção da curva ROC
A validade dos testes, na distinção entre os grupos doença-positivo e doença-negativo, pode ser quantificada pelas curvas ROC, usando a função roc() do pacote pROC (167). Este pacote tem várias funções:
- auc: calcula a área da curva ROC;
- ci: calcula o intervalo de confiança da curva ROC;
- ci.auc: calcula o intervalo de confiança da AUC;
- ci.se: calcula o intervalo de confiança de sensibilidades em determinadas especificidades;
- ci.sp: calcula o intervalo de confiança de especificidades em determinadas sensibilidades;
- ci.thresholds: calcula o intervalo de confiança dos limites;
- coords: Retorna as coordenadas (sensibilidades, especificidades, pontos de corte) de uma curva ROC;
- roc: Constroi uma curva ROC;
- roc.test: Compara a AUC de duas curvas ROC correlacionadas;
- smooth: suaviza a curva ROC
Usar a função com os argumentos variável resposta (doenca), variável preditora (teste3, teste2 e teste1), indicação de que o gráfico deve ser desenhado (plot = TRUE). Como por padrão o gráfico é plotado com a sensibilidade no eixo x e a especificidade no eixo y; deve-se acrescentar o argumento legacy.axes = TRUE para aparecer o seu complemento, os falsos positivos (\(1 – especificidade\)).
Além desses, pode-se usar vários outros argumentos como: print.auc = TRUE, que imprime no gráfico a AUC e ci que é o intervalo de confiança da AUC. Para que a sensibilidade e especificidade apareçam como uma percentagem, deve-se usar o argumento percent = TRUE, pois o padrão é FALSE. Os demais argumentos são os rótulos dos eixos, cor da curva, largura da curva (lwd).
roc3 <- roc (dados$doenca,
dados$teste3,
plot=TRUE,
quiet = TRUE,
legacy.axes=TRUE,
print.auc=TRUE,
print.auc.y = 0.2,
ci = TRUE,
ylab="Sensibilidade",
xlab="1 - Especificdade",
col="steelblue",
smooth = TRUE,
lwd=2)
roc2 <- roc (dados$doenca,
dados$teste2,
plot=TRUE,
quiet = TRUE,
legacy.axes=TRUE,
print.auc=TRUE,
ci = TRUE,
print.auc.y=0.13,
col="chartreuse4",
lwd=2,
smooth = TRUE,
add=TRUE)
roc1 <- roc (dados$doenca,
dados$teste1,
plot=TRUE,
quiet = TRUE,
legacy.axes=TRUE,
print.auc=TRUE,
ci = TRUE,
print.auc.y=0.06,
col="tomato",
lwd=2,
smooth = TRUE,
add=TRUE)
# Legendas das curvas ROC
text (0.73,0.80,"Teste 3", col="steelblue", cex = 1)
text (0.53,0.73,"Teste 2", col="chartreuse4", cex = 1)
text (0.35,0.65,"Teste 1", col="tomato", cex = 1) 
Figura18.6: Curvas ROC para os Testes 1, 2 e 3.
Interpretação do resultado
A Figura 18.6, que exibe o desempenho diagnóstico dos três testes. Em uma curva ROC, a sensibilidade é calculada usando cada valor do teste no conjunto de dados como um ponto de corte e é plotada em relação à (1 – especificidade) correspondente nesse ponto, como mostrado na Figura 18.6.
Assim, a curva são os Verdadeiros Positivos (VP) plotados em relação aos Falsos Positivos (FP), calculados usando cada valor do teste como ponto de corte. A reta diagonal indica onde o teste cairia se os resultados não fossem melhores do que o acaso para predizer a presença de uma doença. O Teste 1 está próximo desta reta, confirmando que ele tem pouca capacidade de discriminar os pacientes doentes e não doentes.
A área abaixo da reta diagonal é equivalente a 0,5 da área total. Quanto maior a área sob a curva ROC, mais útil é o teste para predizer os pacientes que têm a doença. Uma curva que cai substancialmente abaixo da linha diagonal indica que o teste tem pouca capacidade de diagnosticar a doença. Quando há uma separação perfeita dos valores dos dois grupos, isto é, sem sobreposição das distribuições, a área sob a curva ROC é igual a 1 (a curva ROC alcançará o canto superior esquerdo do gráfico).
A área sob a curva (Area Under the Curve – AUC) e seu intervalo de confiança de 95% podem ser obtidos com os comandos usados na construção da Figura 18.6 ou separadamente usando as funções auc() e ci.auc() do pacote pROC.
## Area under the curve: 0.5891## 95% CI: 0.4856-0.681 (2000 stratified bootstrap replicates)## Area under the curve: 0.7616## 95% CI: 0.6743-0.8385 (2000 stratified bootstrap replicates)## Area under the curve: 0.898## 95% CI: 0.8337-0.9409 (2000 stratified bootstrap replicates)A acurácia geral de um teste pode ser descrita como a área sob a curva; quanto maior for a área, melhor será o teste. Na Figura 18.6, o Teste 3 tem uma AUC maior que os outros dois testes.
Usa-se a seguinte estimativa (Tabela 18.1 para avaliar a acurácia de um teste ou da capacidade de identificar corretamente uma condição usando curva ROC (168):
| AUC | Qualidade do Teste |
|---|---|
| >0,90 | excelente |
| 0,80 a 0,90 | muito bom |
| 0,70 a 0,80 | bom |
| 0,60 a 0,70 | suficiente |
| 0,50 a 0,60 | ruim |
| <0,50 | ignorar teste |
Desta forma, o Teste 3 pode ser considerado um bom teste e o Teste 1 é um teste ruim.
Comparando duas curvas
Pode-se comparar duas curvas ROC com a função roc.test(), por exemplo, comparando as curvas dos Teste 3 e 2 (169):
##
## Bootstrap test for two correlated ROC curves
##
## data: roc3 and roc2
## D = 4.6643, boot.n = 2000, boot.stratified = 1, p-value = 3.097e-06
## alternative hypothesis: true difference in AUC is not equal to 0
## sample estimates:
## Smoothed AUC of roc1 Smoothed AUC of roc2
## 0.8980454 0.7616201O Teste 3 tem uma AUC que o caracteriza como um bom teste e o teste de DeLong, entregue na saída do roc.test(), resultou que a diferença entre ele o Teste 2 é estatisticamente significativa (P < 0,0001).
18.2.4.3 Melhor ponto de corte
O melhor ponto de corte (Best Critical Value), que às vezes é chamado de ponto de diagnóstico ótimo ou de Youden, é o ponto da curva mais próximo da parte superior do eixo y (Figura 18.6, Teste 3). Este é o ponto em que a taxa de verdadeiros positivos é otimizada e a taxa de falsos positivos é minimizada. O melhor ponto de corte para o Teste 3 é mostrado na Figura 18.7. Este melhor ponto de corte pode ser identificado a partir dos pontos de coordenadas da curva, usando a função roc() com os seguintes argumentos:
best <- roc (dados$doenca,
dados$teste3,
plot = TRUE,
ci=TRUE,
thresholds="best",
print.thres="best",
legacy.axes=TRUE,
main="",
ylab="Sensibilidade",
xlab="1 - Especificidade",
col="steelblue",
lwd=2)
Figura18.7: Curvas ROC para os Testes 1, 2 e 3.
##
## Call:
## roc.default(response = dados$doenca, predictor = dados$teste3, ci = TRUE, plot = TRUE, thresholds = "best", print.thres = "best", legacy.axes = TRUE, main = "", ylab = "Sensibilidade", xlab = "1 - Especificidade", col = "steelblue", lwd = 2)
##
## Data: dados$teste3 in 48 controls (dados$doenca 1) > 97 cases (dados$doenca 2).
## Area under the curve: 0.8973
## 95% CI: 0.8444-0.9502 (DeLong)Assim, para o Teste 3, o ponto de corte ideal é 24,8, onde a especificidade é igual a 0,854 e a sensibilidade é igual 0,845. Estes dados, fornecem um LR para um resultado positivo igual a:
\[ LR \left(+\right) = \frac{0.845}{\left (1-0.854\right)} = 5,79 \]
As coordenadas da curva ROC podem ser obtida com a seguinte programação, a partir de uma sensibilidade e especificidade acima de 0 (zero):
## threshold specificity sensitivity
## 1 Inf 0.00000000 1
## 2 57.50 0.02083333 1
## 3 54.65 0.04166667 1
## 4 53.45 0.06250000 1
## 5 52.80 0.08333333 1
## 6 51.30 0.10416667 1
## 7 49.65 0.12500000 1
## 8 48.65 0.16666667 1
## 9 47.50 0.18750000 1
## 10 46.50 0.35416667 1A estatística J de Youden (170) é calculada deduzindo 1 a partir da soma de sensibilidade e especificidade do teste e não é expressa como porcentagem, mas como parte de um número inteiro: \(\left (sensibilidade + especificidade\right) - 1\). A estatística J de Youden no melhor ponto de corte do Teste 3 é igual a \(\left (0,845+ 0,854\right) - 1 = 0,699\).
Este é o maior valor de todos os valores das coordenadas (91 valores) usadas.
## [1] 0.6995275A Figura 18.7 mostra o ponto de corte ideal. Ele também pode ser obtido com a função coords() do pacote pRoc:
roc3 <- dados %>% roc(doenca, teste3)
coords(roc3, x = "best", ret="threshold", transpose = FALSE,
best.method="youden")## threshold
## 1 24.8O método para obter o melhor ponto de corte (best.method) pode ser pelo método de youden ou closest.topleft. No exemplo, o resultado é o mesmo. Para maiores detalhes consulte a ajuda da função (?coord).
18.3 Estatística kappa
A estatística de concordância kappa (k) de Cohen é utilizada para descrever a concordância entre dois ou mais avaliadores quando realizam uma avaliação nominal ou ordinal de uma mesma amostra (171). A estatística kappa corrige a chance do acaso nas avaliações e é obtida pela fórmula igual a:
\[ k= \frac{p_{o} - p_{e}}{1 - p_{e}} \]
Onde \(p_{o}\) = proporção observada de concordância e \(p_{e}\) = proporção esperada de concordância apenas pelo acaso.
Por exemplo, dois radiologistas podem revisar independentemente uma série de radiografias do tórax de pacientes para determinar a presença ou ausência de pneumonia. Para avaliar o grau de concordância entre as classificações dos dois médicos, pode ser relatado o percentual de concordância entre os avaliadores (por exemplo, 50% dos avaliadores responderam “sim” nas duas ocasiões). No entanto, esse percentual pode ser enganoso, pois não leva em conta o nível de concordância entre os dois avaliadores que pode ocorrer por acaso. A estatística kappa pode ser usada para avaliar a concordância das respostas para dois ou mais avaliadores após considerar a concordância casual. Portanto, a estatística kappa é uma estimativa da proporção de concordância entre avaliadores que excede a concordância que ocorreria por acaso.
A interpretação dos valores de kappa é mostrada na Tabela 18.3 (172). Quando a proporção observada de concordância é menor que a esperada por acaso, o kappa terá um valor negativo indicando não concordância. Um valor de kappa igual a 0 indica que a concordância observada é igual à concordância casual.
O teste de hipóteses testa a hipótese de que a concordância entre os dois avaliadores seja puramente aleatória. Quando o valor P é menor que 0,05, rejeitamos a hipótese de que a concordância foi puramente aleatória. As premissas para o kappa de Cohen são que os participantes ou itens a serem classificados são independentes e também que os avaliadores e categorias são independentes.
| Valor kappa | Concordância |
|---|---|
| <0,00 | pobre |
| 0,00 - 0,20 | leve |
| 0,21 - 0,40 | razoável |
| 0,41 - 0,60 | moderada |
| 0,61 - 0,80 | substancial |
| 0,81 - 1,00 | quase perfeita |
Existem diferentes tipos de estatísticas kappa. Para dados com três ou mais categorias possíveis (por exemplo, concordo, concordo parcialmente, discordo) ou para dados categóricos ordenados, o kappa ponderado deve ser usado para que as respostas que estão mais distantes da concordância tenham maior peso do que aquelas próximas à concordância. No exemplo usado, as categorias possíveis são dicotômicas (sim e não), portanto, o kappa não ponderado (unweighted) e o ponderado (weighted) retornam o mesmo resultado.
18.3.1 Exemplo
O arquivo dadosPneumonia.xlsx contém os dados de 54 crianças com suspeita de pneumonia, cujas radiografias foram avaliadas por dois radiologistas. O objetivo foi medir a concordância diagnóstica dos dois profissionais. Para o cálculo do coeficiente kappa será usada a função Kappa() do pacote vcd (173). Essa função tem os seguintes argumentos:
- x \(\longrightarrow\) matriz ou tabela
- weights \(\longrightarrow\) matriz especificada pelo usuário com as mesmas dimensões de x, desnecessário para kappa não ponderado.
Na impressão do kappa pode-se usar print (k, digits = 3, CI = TRUE, level = 0.95). Onde k é o coeficiente de kappa, calculado pela função Kappa(), CI é o intervalo de confiança e o nível de confiança padrão é 95%.
18.3.1.1 Leitura e exploração dos dados
O conjunto de dados dadosPneumonia.xlsx pode ser obtido aqui. Após salvar o arquivo em seu diretório, ele pode ser carregado com a função read_excel() do pacote readxl:
18.3.1.2 Construção da tabela
O cálculo do kappa com a função Kappa() exige uma tabela, onde os dados dos dois radiologistas são cruzados. As variáveis a serem cruzadas são rx1 e rx2:
dados$rx1 <- factor(dados$rx1,
ordered=TRUE,
levels = c("sim", "não"))
dados$rx2 <- factor(dados$rx2,
ordered=TRUE,
levels = c("sim", "não"))
tabk <- table (dados$rx1,
dados$rx2,
dnn = c ("Radiologista 1", "Radiologista 2"))
tabk## Radiologista 2
## Radiologista 1 sim não
## sim 32 5
## não 3 1418.3.1.3 Cálculo do kappa
O kappa é dado pela execução da função:
## value ASE z Pr(>|z|) lower upper
## Unweighted 0.667 0.107 6.21 5.42e-10 0.456 0.878
## Weighted 0.667 0.107 6.21 5.42e-10 0.456 0.878A saída exibe o kappa pontual e os intervalos de confiança de 95%, podendo-se concluir, desses resultados, que existe uma boa confiabilidade nos diagnósticos dos radiologistas (k = 0,67, concordância substancial,de acordo com a Tabela 18.3).
18.4 Medidas de frequência
18.4.1 Prevalência
A prevalência, ou mais adequadamente, a prevalência pontual de uma doença é a proporção da população portadora da doença em um determinado ponto do tempo. É uma medida instantânea por excelência e fornece uma medida estática da frequência da doença. É também conhecida como taxa de prevalência e é expressa em percentagem ou por \(10^{n}\) habitantes. As medidas de prevalência geram informações úteis para o planejamento e administração de serviços de saúde.
A prevalência por período descreve os casos que estavam presentes em qualquer momento durante um determinado período de tempo. Diz o número total de casos de uma doença que se sabe haver existido durante um período de tempo.
Um tipo especial de prevalência de período é a prevalência ao longo da vida, que mede a frequência cumulativa ao longo da vida de um resultado até o momento presente (ou seja, a proporção de pessoas que tiveram o evento em qualquer momento no passado).
As doenças, quanto a sua duração, podem ser agudas e de longa duração ou crônicas. A prevalência é proporcional ao tempo de duração da doença. Hipoteticamente, se o surgimento de novos casos de doença ocorre em ritmo constante e igual para doenças agudas e crônicas, estas últimas acumularão casos, aumentando a prevalência. As doenças agudas tenderão a manter uma prevalência constante. A terapêutica, diminuindo o tempo de duração das doenças, também reduz a prevalência. A prevalência é dada pela razão:
\[ prevalência = \frac{número \ de \ casos \ conhecidos \ da \ doença}{total \ da \ População} \times 10^{n} \]
18.4.1.1 Exemplo
Como exemplo, será verificada a frequência de tabagismo entre as puérperas da maternidade do HGCS. O banco de dados dadosMater.xlsx contém informação de 1368 nascimentos e pode ser consultado na Seção 5.3. Depois de salvo em seu diretório de trabalho, ele pode ser carregado com a função read_excel() do pacote readxl.
Inicialmente, será verificado quantas fumantes existem. O conjunto de dados contém uma variável fumo, onde 1 = fumante e 2 = não fumante. Portanto, há necessidade de transformar a variável numérica em um fator:
dados$fumo <- factor (dados$fumo,
ordered = TRUE,
levels = c(1,2),
labels = c("fumante", "não fumante"))
tabFumo <- table(dados$fumo)
tabFumo##
## fumante não fumante
## 301 1067Além de relatar a estimativa pontual da frequência da doença, é importante fornecer uma indicação da incerteza em torno dessa estimativa pontual. A função epi.conf(), do pacote epiR (157), permite calcular intervalos de confiança para prevalência, motivo da escolha dessa função.
A função epi.conf() usa os seguintes argumentos:
- dat \(\longrightarrow\) matriz ou tabela;
- ctype \(\longrightarrow\) tipo de intervalo de confiança a ser calculado. Opções: mean.single, mean.unpair, mean.pair, prop.single, prop.unpaired, prevalence, inc.risk, inc.rate, odds e smr (standardized mortality rate);
- method \(\longrightarrow\) método a ser usado. Quando
ctype = "inc.risk"ouctype = "prevalence", as opções sãoexact,wilsonefleissQuandoctype = "inc.rate"as opções sãoexactebyar; - N \(\longrightarrow\) tamanho da população;
- conf.level \(\longrightarrow\) magnitude do intervalo de confiança retornado. Deve ser um único número entre 0 e 1.
Construção da matriz
Com os dados da tabFumo, constrói-se uma matriz de duas colunas:
## n1 N1
## [1,] 301 1368Cálculo da prevalência
Usando a função epiR(), tem-se:
## est lower upper
## 1 0.2200292 0.1983313 0.2429365A saída mostra que a prevalência de fumantes entre as puérperas do HGCS é igual a 22,0% (IC95%: 19,8 – 24,3%).
18.4.2 Incidência
A incidência fornece uma medida da frequência com que os indivíduos suscetíveis se tornam casos de doenças, à medida que são observados ao longo do tempo.
Um caso de incidente ocorre quando um indivíduo deixa de ser suscetível e passa a ser doente. A contagem de casos de incidentes é o número de tais eventos que ocorrem em uma população durante um período de acompanhamento definido. Existem duas maneiras de expressar a incidência:
A incidência cumulativa (risco) é a proporção de indivíduos inicialmente suscetíveis em uma população que se tornam novos casos durante um período de acompanhamento definido.
Para calcular a incidência cumulativa, é necessário primeiro identificar os doentes e após acompanhar por um determinado tempo os não doentes (Figura 18.8).
Figura18.8: Incidência
A taxa de incidência (densidade de incidência ou taxa de incidência) é o número de novos casos da doença que ocorrem por unidade de tempo em risco durante um período de acompanhamento definido. Este período é expresso como pessoas-tempo (pessoas-ano, por exemplo).
O conceito de pessoas-tempo pode ser ilustrado com o seguinte exemplo: a Figura 18.9 representa um estudo epidemiológico hipotético com duração de cinco anos, onde D é o desfecho e C representa os sujeitos que deixaram o estudo por migração ou morte (censurados) por causa não relacionada ao desfecho
Figura18.9: Pessoas-tempo (estudo epidemiológico hipotético).
Nesse estudo hipotético, o indivíduo 1 permaneceu no estudo 3,5 anos; o indivíduo 2,5 anos; o indivíduo 3, 4,5 anos e, assim por diante, totalizando 32,5 pessoas-anos. Em outras palavras, ocorreram 4 desfechos durante os 5 anos do estudo, consequentemente, a taxa de incidência (TI) foi de
\[ TI = \frac{4}{32,5} \times 1000 = \frac{123}{1000\ pessoas-ano} \]
Isto significa que se fossem acompanhadas 1000 pessoas por um ano, 123 delas apresentariam o desfecho D.
18.4.2.1 Exemplo
Aparentemente, pessoas cegas tem uma menor incidência de câncer e esse efeito parece ser mais pronunciado em pessoas totalmente cegas do que em pessoas com deficiência visual grave.
Para testar essa hipótese, foi identificada uma coorte de 1.567 pessoas totalmente cegas e 13.292 sujeitos com deficiência visual grave. As informações sobre a incidência de câncer foram obtidas do Registro Sueco de Câncer (174). Foram diagnosticados de 136 casos de câncer em 22050 pessoas-ano em risco totalmente cegas e 1709 casos de câncer em 127650 pessoas-anos em risco com deficiência visual grave.
A taxa de incidência pode ser calculada, usando-se a mesma função epi.conf(), usada para o cálculo da prevalência, mudando o argumento ctype = “prevalence” para ctype = “inc.rate”, conforme recomendado:
Pessoas totalmente cegas
Inicialmente, contrói-se a matriz:
## n2 N2
## [1,] 136 22050Logo, a incidência de câncer nos totalmente cegos é:
## est lower upper
## n2 6.1678 5.174806 7.295817Pessoas com grave deficiência visual
Inicialmente, contrói-se a matriz:
## n3 N3
## [1,] 1709 127650Logo, a incidência de câncer nos com grave deficiência visual é:
## est lower upper
## n3 13.38817 12.76088 14.03832As saídas mostram que para cada 1000 pessoas cegas (a função foi multiplicada por 1000) acompanhadas por um ano, ocorreu 6,2 ((IC95%: 5,2 – 7,3) casos de câncer. Uma taxa de incidência, praticamente, metade da taxa de incidências das pessoas com deficiência visual grave. Os IC95% não são coincidentes, o que significa que essa diferença é significativa. Houve, na amostra, uma incidência menor de câncer entre os indivíduos totalmente cegos, sugerindo que a melatonina possa ser um fator protetor contra o câncer.
18.4.3 Relação entre prevalência e incidência
A incidência é uma medida de risco. A prevalência, por não levar em consideração o tempo de duração da doença (t), não tem esta capacidade. Em uma população onde a situação da doença encontra-se em estado estacionário (ou seja, sem grandes migrações ou mudanças ao longo do tempo na incidência/prevalência), a relação entre prevalência e incidência e duração da doença pode ser expressa pela seguinte fórmula (175):
\[ prevalência \ pontual = incidência \times prevalência \]
Por exemplo, se a incidência da doença for de 0,8% ao ano e sua duração média (sobrevida após o diagnóstico) for de 10 anos, a prevalência pontual será de aproximadamente 8%.
18.5 Medidas de associação
18.5.1 Odds Ratio
Odds Ratio (OR) é a razão entre dois odds. A Odds Ratio, traduzida como Razão de Chances, está associada, usualmente, com estudos retrospectivos tipo caso-controle com desfechos dicotômicos.
A odds ratio (OR) expressa a odds de exposição entre os que têm o desfecho (casos) pela odds de exposição nos livres de desfecho (controles).
Figura18.10: Tabela de contingência 2 x 2
Usando a Figura 18.10, a fórmula \(odds =\frac{p}{1 -p}\) e que
\[ p_{exp \ doentes} = \frac{a}{a+c} \] \[ p_{exp \ não \ doentes} = \frac{b}{b+d} \]
tem-se:
\[ odds_{exp} \ {casos} = \frac{\frac{a}{a+c}}{1- \frac{a}{a+c}}=\frac{a}{c} \] \[ odds_{exp} \ {controles} = \frac{\frac{b}{b+d}}{1- \frac{b}{b+d}}=\frac{b}{d} \]
Portanto, a OR é igual a:
\[ OR = \frac{odds_{exp}\ {casos}}{odds_{exp}\ {controles}}=\frac{\frac{a}{c}}{\frac{b}{d}}=\frac{a \times d}{c \times b} \] Em decorrência da última fórmula, a OR é definida como a razão dos produtos cruzados em uma tabela de contingência 2×2.
18.5.1.1 Exemplo
Em um estudo de caso-controle hipotético, a distribuição das exposições entre os casos e um grupo de pessoas saudáveis (“controles”) é comparada entre si. Os casos correspondem a um tipo raro de câncer, onde se suspeita que exista uma associação à exposição a um determinado fator de risco.
Os dados desse estudo hipotético estão no arquivo dadosCasoControle.xlsx. O conjunto de dados pode ser obtido aqui. Depois de salvo em seu diretório de trabalho, ele pode ser carregado com a função read_excel() do pacote readxl.
As variáveis cc$exposto e cc$desfecho devem ser transformadas em fatores e na ordem sim, não, uma vez que o R coloca em ordem alfabética (não, sim):
cc$exposto <- factor (cc$exposto,
levels = c("sim", "não"))
cc$desfecho <- factor (cc$desfecho,
levels = c("sim", "não"))Após essa etapa, construir uma tabela \(2 \times 2\):
## Desfecho
## Exposição sim não Sum
## sim 48 20 68
## não 12 40 52
## Sum 60 60 120A OR será obtida utilizando a função epi.2by2() do pacote epiR (157). Esta função tem os seguintes argumentos:
- dat \(\longrightarrow\) tabela de contingência \(2 \times 2\);
- method \(\longrightarrow\) as opções são “cohort.count”, “cohort.time”, “case.control” ou “cross.sectional”.;
- conf.level \(\longrightarrow\) padrão = 0.95;
- units \(\longrightarrow\) multiplicador para incidência e prevalência;
- outcome \(\longrightarrow\) indicação de como a variável desfecho é representada na tabela de contingência (“as.columns” ou “as.rows”).
epiR::epi.2by2(tab_cc,
method = "case.control",
conf.level = 0.95,
units = 100,
outcome = "as.columns")## Outcome + Outcome - Total Odds
## Exposed + 48 20 68 2.40 (1.43 to 4.23)
## Exposed - 12 40 52 0.30 (0.13 to 0.53)
## Total 60 60 120 1.00 (0.69 to 1.45)
##
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Exposure odds ratio 8.00 (3.49, 18.34)
## Attrib fraction (est) in the exposed (%) 87.24 (69.26, 95.03)
## Attrib fraction (est) in the population (%) 70.00 (48.69, 82.46)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 26.606 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
## Wald confidence limits
## CI: confidence intervalA saída exibe os dados em uma tabela \(2 \times 2\), mostrando as odds e os IC95% e outras estatísticas epidemiológicas relacionadas.
A OR varia de zero ao infinito. Quando o valor da OR se aproxima de 1, a doença e o fator de risco não estão associados. Acima de 1 significa que existe associação e valores menores de 1 indicam uma associação negativa (efeito protetor).
No exemplo hipotético, os indivíduos que se expuseram ao fator de risco têm uma chance 8 vezes maior de apresentar este tipo de câncer. O valor P do qui-quadrado é altamente significativo (P < 0,001).
18.5.2 Risco Relativo
O Risco relativo (RR) é a razão entre a incidência de desfecho em indivíduos expostos e a incidência de desfecho em indivíduos não expostos. O RR estima a magnitude da associação entre a exposição e o desfecho (doença). Em outras palavras, compara a probabilidade de ocorrência do desfecho entre os indivíduos expostos com a probabilidade de ocorrência do desfecho nos indivíduos não expostos.
A partir da tabela de contingência \(2 \times 2\) (Figura 18.10), tem-se que o estimador do RR é dado por:
\[ RR = \frac{incidência_{exp}}{incidência_{não \ exp}}=\frac{\frac{a}{a + b}}{\frac{c}{c + d}} \]
18.5.2.1 Exemplo
Em 1940, ocorreu um surto de gastroenterite, após um jantar, em uma igreja, na cidade de Lycoming, Condado de Oswego, Nova York. Das 80 pessoas presentes, 75 foram entrevistadas. Quarenta e seis relataram doença gastrointestinal, atendendo à definição de caso.
As taxas de ataque (incidência) foram calculadas para aqueles que comeram e não comeram cada um dos 14 itens alimentares consumidos na ceia (176). O pacote epitools (177) contém os dados desta investigação no arquivo oswego.
## Rows: 75
## Columns: 21
## $ id <int> 2, 3, 4, 6, 7, 8, 9, 10, 14, 16, 17, 18, 20, 21, 2…
## $ age <int> 52, 65, 59, 63, 70, 40, 15, 33, 10, 32, 62, 36, 33…
## $ sex <chr> "F", "M", "F", "F", "M", "F", "F", "F", "M", "F", …
## $ meal.time <chr> "8:00 PM", "6:30 PM", "6:30 PM", "7:30 PM", "7:30 …
## $ ill <chr> "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", …
## $ onset.date <chr> "4/19", "4/19", "4/19", "4/18", "4/18", "4/19", "4…
## $ onset.time <chr> "12:30 AM", "12:30 AM", "12:30 AM", "10:30 PM", "1…
## $ baked.ham <chr> "Y", "Y", "Y", "Y", "Y", "N", "N", "Y", "N", "Y", …
## $ spinach <chr> "Y", "Y", "Y", "Y", "Y", "N", "N", "Y", "N", "Y", …
## $ mashed.potato <chr> "Y", "Y", "N", "N", "Y", "N", "N", "Y", "N", "N", …
## $ cabbage.salad <chr> "N", "Y", "N", "Y", "N", "N", "N", "N", "N", "N", …
## $ jello <chr> "N", "N", "N", "Y", "Y", "N", "N", "N", "N", "N", …
## $ rolls <chr> "Y", "N", "N", "N", "Y", "N", "N", "Y", "N", "Y", …
## $ brown.bread <chr> "N", "N", "N", "N", "Y", "N", "N", "Y", "N", "N", …
## $ milk <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N", …
## $ coffee <chr> "Y", "Y", "Y", "N", "Y", "N", "N", "N", "N", "Y", …
## $ water <chr> "N", "N", "N", "Y", "Y", "N", "N", "Y", "N", "N", …
## $ cakes <chr> "N", "N", "Y", "N", "N", "N", "Y", "N", "N", "Y", …
## $ vanilla.ice.cream <chr> "Y", "Y", "Y", "Y", "Y", "Y", "N", "Y", "Y", "Y", …
## $ chocolate.ice.cream <chr> "N", "Y", "Y", "N", "N", "Y", "Y", "Y", "Y", "Y", …
## $ fruit.salad <chr> "N", "N", "N", "N", "N", "N", "N", "N", "N", "N", …Existem 75 observações de 21 variáveis, algumas características dos indivíduos como idade, sexo, etc. Importante para a análise é a variável ill (Y – sim, doente; N – não doente) e a variáveis relacionadas aos alimentos ingeridos durante o jantar na igreja. O sorvete de baunilha foi considerado o principal responsável pelo surto.
A seguir, as variáveis oswego$vanilla.ice.cream e oswego$ill 61 serão transformadas em fator e os níveis colocados na ordem Y, N, uma vez que o R coloca em ordem alfabética (N, Y) :
oswego$ill <- factor (oswego$ill,
levels = c ("Y", "N"))
oswego$vanilla.ice.cream <- factor (oswego$vanilla.ice.cream,
levels = c ("Y", "N"))Realizada essa etapa, será construída uma tabela para o cálculo do RR:
## Ill
## Vanilla Y N
## Y 43 11
## N 3 18O RR será obtido, utilizando a função epi.2by2() do pacote epiR, cujos argumentos foram mostrados no cálculo da OR, mudando a tabela para tab_vanilla e method = “cohort.count”:
epiR::epi.2by2(tab_vanilla,
method = "cohort.count",
conf.level = 0.95,
units = 100,
outcome = "as.columns")## Outcome + Outcome - Total Inc risk *
## Exposed + 43 11 54 79.63 (66.47 to 89.37)
## Exposed - 3 18 21 14.29 (3.05 to 36.34)
## Total 46 29 75 61.33 (49.38 to 72.36)
##
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio 5.57 (1.94, 16.03)
## Inc odds ratio 23.45 (5.84, 94.18)
## Attrib risk in the exposed * 65.34 (46.92, 83.77)
## Attrib fraction in the exposed (%) 82.06 (48.41, 93.76)
## Attrib risk in the population * 47.05 (28.46, 65.63)
## Attrib fraction in the population (%) 76.71 (37.11, 91.37)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 27.223 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
## Wald confidence limits
## CI: confidence interval
## * Outcomes per 100 population unitsOs resultados da saída indicam que os indivíduos que ingeriram sorvete de baunilha (n = 54) tiveram um risco maior de desenvolver gastrenterite aguda quando comparado aos que não ingeriram (n = 21). Dividindo o risco dos indivíduos expostos (incidência = 79,6) pelo risco dos não expostos (incidência = 14,3), encontra-se o RR = 5,57. Isso confirma que o sorvete de baunilha foi o principal responsável.
Quanto maior o RR mais forte é a associação entre a doença em questão e a exposição ao fator de risco. Um RR = 1 indica que a doença e a exposição ao fator de risco não estão associadas. Valores < 1 indicam uma associação negativa entre o fator de risco e a doença (efeito protetor).
18.5.3 Odds Ratio vs Risco Relativo
A OR não deve ser entendida como uma medida aproximada do RR, exceto para doenças raras (doenças, em geral com prevalência menor do que 10%). Caso contrário, a OR tenderá a superestimar a magnitude da associação e o OR afasta-se da hipótese nula da não associação (OR =1), independentemente de ser um fator de risco ou de proteção. A discrepância (d)62 entre as estimativas do RR e OR pode ser definido como a razão entre o OR e o RR estimados (178). Em outras palavras, a discrepância corresponde a uma proporção do RR (179).
\[ d = \frac {1- p_{não \ exp}}{1- p_{exp}}= \frac{\frac{c}{c + d}}{\frac{a}{a + b}} \]
Logo,
\[ OR = RR \times d \]
Para finalizar, uma comparação entre OR e RR é mostrada na Tabela 18.5 (180).
| OR | RR | Magnitude |
|---|---|---|
| 1,0 | 1,0 | insignificante |
| 1,5 | 1,2 | pequena |
| 3,5 | 1,9 | moderada |
| 9,0 | 3,0 | grande |
| 32 | 5,7 | muito grande |
| 360 | 19 | quase perfeita |
| infinito | infinito | perfeita |
18.5.4 Razão de Prevalência
Quando dados transversais estão disponíveis, muitas vezes as associações são avaliadas, usando a razão de prevalência pontual (RPP).
Tendo o mesmo princípio das duas medidas anteriores, a razão de prevalência (RPP) compara a prevalência do desfecho entre os expostos com a prevalência do desfecho entre os não expostos.
Matematicamente, a RPP é calculada de maneira semelhante ao RR. Apenas, deve-se ter em mente que o desfecho e a exposição foram medidos no mesmo momento, enquanto para o cálculo do RR há necessidade de calcular a incidência.
Usando uma tabela de contingência 2 x 2 (Figura 18.10), tem-se:
\[ RPP = \frac{prevalência \ de \ doença_{exp}}{prevalência \ de \ doença_{não \ exp}}=\frac{\frac{a}{a + b}}{\frac{c}{c + d}} \]
Também é possível verificar a prevalência de exposição entre doentes e não doentes:
\[ RPP = \frac{prevalência \ de \ exposição_{doentes}}{prevalência \ de \ exposição_{não \ doentes}}=\frac{\frac{a}{a + c}}{\frac{b}{b + d}} \]
18.5.4.1 Exemplo
Em um estudo transversal (181), foi verificada a prevalência de infecções congênitas entre as puérperas com idade igual ou acima de 20 anos comparadas às mulheres com menos de 20 anos (adolescentes). A hipótese foi de que as adolescentes tinham uma prevalência maior de infecções.
Parte dos dados estão no arquivo dadosMater.xlsx, que contém, como já mencionado, informações de 1368 nascimentos. Entre essas, tem-se a idade das mães (idadeMae) e se foi diagnosticada infecção congênita (infCong).
O arquivo pode ser obtido aqui. Depois de salvo em seu diretório de trabalho, ele pode ser carregado com a função read_excel() do pacote readxl.
A partir da variável idadeMae, criar a variável faixaEtaria, dividindo as parturientes em menores de 20 anos (adolescentes) e ≥ 20 anos. Para isso, usou-se a função cut() do pacote base. Revise os argumentos desta função.
dados$faixaEtaria <- cut (dados$idadeMae,
breaks=c(13,20,46),
labels = c("<20a","=>20a"),
right = FALSE,
include.lowest = TRUE)A variável ìnfCong encontra-se como uma variável numérica e deve ser transformada em fator:
dados$infCong <- factor (dados$infCong,
ordered = TRUE,
levels = c (1,2),
labels = c ("sim", "não"))Após estes procedimentos, constroi-se uma tabela \(2 \times 2\):
tab_infCong <- table(dados$faixaEtaria,
dados$infCong,
dnn = c("Faixa Etária", "Inf. Cong."))
addmargins(tab_infCong) ## Inf. Cong.
## Faixa Etária sim não Sum
## <20a 7 212 219
## =>20a 119 1030 1149
## Sum 126 1242 1368Cálculo da RPP
Usando a tabela tab_infCong com a função epi.2by2() do pacote epiR, cujos argumentos foram mostrados no cálculo da OR e RR, e mudando a tabela para tab_infCong e method = “cross.sectional”, obtem-se:
epiR::epi.2by2(tab_infCong,
method = "cross.sectional",
conf.level = 0.95,
units = 100,
outcome = "as.columns")## Outcome + Outcome - Total Prev risk *
## Exposed + 7 212 219 3.20 (1.29 to 6.47)
## Exposed - 119 1030 1149 10.36 (8.65 to 12.26)
## Total 126 1242 1368 9.21 (7.73 to 10.87)
##
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Prev risk ratio 0.31 (0.15, 0.65)
## Prev odds ratio 0.29 (0.13, 0.62)
## Attrib prev in the exposed * -7.16 (-10.08, -4.24)
## Attrib fraction in the exposed (%) -224.02 (-584.89, -53.29)
## Attrib prev in the population * -1.15 (-3.48, 1.19)
## Attrib fraction in the population (%) -12.45 (-17.53, -7.58)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 11.278 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
## Wald confidence limits
## CI: confidence interval
## * Outcomes per 100 population unitsA saída exibe várias informações. Foi feita a hipótese de uma maior prevalência entre as mulheres com menos de 20 anos. Por este motivo, elas aparecem como as expostas (Exposed +) e tem uma prevalência de 3,20/100, enquanto as mulheres com mais de 20 anos tiveram uma prevalência de 10,36/100. Isto mostra que a razão de prevalência é igual a 0,31 (IC95%: 0,15-0,65)63, ou seja, < 1, sugerindo que ao contrário da hipótese inicial, as adolescentes têm, neste estudo, uma menor prevalência de infecções congênitas.
18.6 Medidas de impacto
18.6.1 Risco Atribuível
O Risco Atribuível (RA) possui características de medida de impacto. O RA, ao invés de concentrar-se na associação em si, refere-se mais às consequências e às repercussões da exposição sobre a ocorrência do desfecho.
O RA é a medida do excesso ou acréscimo absoluto de risco que pode ser atribuído à exposição (182). Com o RA é possível estimar o número de casos que podem ser prevenidos se a exposição for eliminada e assim estimar a magnitude do impacto, em termos de saúde pública, imposto por esta exposição.
O risco de desenvolver o desfecho (incidência) está aumentado em RA nos indivíduos expostos em comparação com os que não estão expostos. Nos estudos de coorte, costuma-se usar mais a expressão Risco Atribuível ou Diferença de Risco. Nos ensaios clínicos, usa-se mais a expressão Redução Absoluta do Risco (RAR), pois se espera que a intervenção reduza o risco.
Calcula-se o RA ou a RAR pela diferença absoluta entre as incidências dos expostos e não expostos:
\[ RA = \left|I_{expostos} - I_{não \ expostos}\right| \]
Utilizando a tabela de contingência \(2 \times 2\) (Figura 18.10), o RA fica expresso da seguinte maneira:
\[ RA = \left|\frac{a}{a + b} -\frac{c}{c + d}\right| \]
No exemplo do Risco Relativo, o RA pode ser calculado usando a mesma tabela de contingência, repetida aqui para facilitar a leitura (Figura 18.11):
Figura18.11: Taxa de ataque de gastrenterite com sorvete de baunilha - Oswego
Logo,
\[ RA = \left|\frac{43}{43 + 11} -\frac{3}{3 + 18}\right| = \left|0,796 - 0,143\right| = 0,653 \]
O risco atribuível na exposição mede o excesso de risco associado a uma determinada categoria de exposição. Por exemplo, com base no exemplo, a incidência cumulativa de gastrenterite aguda entre os indivíduos que comeram o sorvete de baunilha é de 79,6% e para os que não ingeriram o sorvete (categoria de referência ou não exposta) foi de 14,3%. Desta forma, o risco excessivo associado à exposição 79,6 – 14,3 = 65,3%. Ou seja, assumindo uma associação causal (sem confusão ou viés), a não ocorrência da festa diminuiria o risco no grupo exposto de 79,6% para 14,3%.
O RA expresso em relação à incidência nos expostos e apresentado em percentual é denominado de Risco Atribuível Proporcional (RAP) ou Fração Atribuível nos Expostos.
O RAP informa qual a proporção de desfecho, expresso em percentagem, entre os expostos que poderia ter sido prevenida se a exposição fosse eliminada. É dado pela fórmula:
\[ RAP = \left(\frac{I_{expostos} - I_{não \ expostos}}{I_{expostos}}\right) \times 100 \]
No exemplo do surto de gastrenterite aguda no jantar da igreja de Oswego, tem-se:
\[ RAP = \left(\frac{0,796 - 0,143}{0,796}\right) \times 100 = 82,06 \% \]
Se a causalidade foi estabelecida, essa medida pode ser interpretada como a porcentagem do risco total de gastrenterite aguda que é atribuível à ingesta de sorvete de baunilha.
Outra maneira de se chegar a este mesmo resultado é através do RR, usando a seguinte fórmula
\[ RAP = \left(\frac{I_{expostos} - I_{não \ expostos}}{I_{expostos}}\right) \times 100 \]
\[ RAP = \left(\frac{I_{expostos}}{I_{expostos}} - \frac{I_{não \ expostos }}{I_{expostos}}\right) \times 100 \]
\[ RAP = \left(1 - \frac{1}{\frac{I_{expostos }}{I_{não \ expostos}}}\right) \times 100 \]
\[ RAP = \left(1 - \frac{1}{RR}\right) \times 100 \]
\[ RAP = \left(\frac{RR - 1}{RR}\right) \times 100 \]
No exemplo, o RR é igual a 5,57, logo:
\[ RAP = \left(\frac{5,57 - 1}{5,57}\right) \times 100 = 82,05\% \]
18.6.2 Redução Relativa do Risco
Quando se avalia um tratamento ou alguma intervenção onde se supõe que haja uma redução do risco, por exemplo, uso da aspirina para reduzir infarto agudo de miocárdio, o termo Risco Atribuível é substituído por Redução do Risco Atribuível e é calculado da mesma forma visto na equação do Risco Atribuível.
Neste caso, ao invés de usar o Risco Atribuível Proporcional (RAP), onde se pressupõe que a exposição é um fator de risco para a doença e o RR > 1, usa-se a Redução Relativa do Risco, pois a exposição é supostamente um fator protetor e o RR < 1, como se espera que ocorra nos ensaios clínicos.
Esta medida análoga ao RAP é também chamada de Eficácia, definida como a proporção da incidência nos indivíduos não tratados (por exemplo, o grupo de controle) que é reduzida pela intervenção (183).
O cálculo da Redução Relativa do Risco (RRR) é semelhante ao Risco Atribuível Proporcional (RAP), onde a incidência nos expostos é a incidência no grupo que recebeu a intervenção (ou taxa de eventos no grupo tratamento) e a incidência nos não expostos é incidência nos controles (ou taxa de eventos nos controles – TEC). Como se supõe que a incidência nos controles seja maior que a incidência no grupo de tratamento, a equação fica:
\[ RRR = \left(\frac{I_{controle} - I_{tratamento}}{I_{controle}}\right) \times 100 \]
Alternativamente, a RRR pode ser estimada pela equação:
\[ RRR = \left(1 - RR\right) \times 100 \]
O Physicians’ Health Study (50) é um ensaio clinico randomizado controlado, duplo cego, desenhado com o objetivo de determinar se uma dose baixa de aspirina (325 mg a cada 48 horas) diminui a mortalidade cardiovascular e se o betacaroteno reduz a incidência de câncer. Participaram deste estudo 22071 indivíduos por uma média de 60,2 meses.
O estudo do componente aspirina mostrou os seguintes resultados (Figura 18.12):
Figura18.12: Physicians’ Health Study, componente aspirina e IAM.
A incidência cumulativa de Infarto Agudo de Miocárdio (IAM) em ambos os grupos foi:
\[ Incidencia_{aspirina} = \frac{139}{11037} = 0,0126 \]
\[ Incidencia_{placebo} = \frac{239}{11034} = 0,0217 \]
\[ RR = \frac{0,0126}{0,0217} = 0,58 \]
Logo, a RRR é igual a:
\[ RRR = \left(1 - 0,58\right) \times 100 = 42\% \]
Ou seja, houve uma redução de 42% no risco de IAM no grupo que usou aspirina e a conclusão dos autores foi que este ensaio clínico demonstrou, em relação à prevenção primária de doença cardiovascular, uma diminuição no risco de IAM.
Estes cálculos podem ser realizados com a função risks() do pacote MKmisc (184). Esta função calcula o risco relativo (RR), odds ratio (OR), redução relativa do risco (RRR) e outras estatísticas epidemiológicas, como RAR, NNT.
A função risks() usa como argumento:
- p0 \(\longrightarrow\) incidência do desfecho de interesse no grupo não exposto;
- p1 \(\longrightarrow\) incidência do desfecho de interesse no grupo exposto.
Além disso, para o seu funcionamento, deve-se ter instalado o pacote BiocManager para poder instalar o pacote limma, necessário para a execução do pacote MKmisc. Veja início do capítulo em pacotes usados neste capítulo.
A função risks() será usada dentro da função round() para reduzir o número de dígitos decimais:
## p0 p1 RR OR RRR ARR NNT
## 0.0217 0.0126 0.5806 0.5753 0.4194 0.0091 109.890118.6.3 Número Necessário para Tratar
Os resultados da função risks() entrega junto o Número Necessário para Tratar (NNT) que deve ser arredondado para o número inteiro mais próximo (no caso, 110) e significa a estimativa do número de indivíduos que devem receber uma intervenção terapêutica, durante um período específico de tempo, para evitar um efeito adverso ou produzir um desfecho positivo.
O NNT equivale à recíproca do RAR (Redução Absoluta do Risco ou Diferença de Risco):
\[ NNT = \frac{1}{RAR} = \frac{1}{I_{não \ expostos} - I_{expostos}} \]
No exemplo do Physicians’ Health Study, o RAR igual a:
\[ RA = \left|I_{expostos} - I_{não \ expostos}\right| = \left|0,0126 - 0,0217\right| = 0,0091 \]
\[ NNT = \frac{1}{0,0091} = 109,89 \simeq 110 \]
Pode-se calcular os IC95%, calculando o NNT para os limites do RAR usando a seguinte equação (185):
\[ IC_{95\%} \longrightarrow RAR \pm z_{\left({1 - \frac{\alpha}{2}}\right)} \times EP_{RAR} \]
Onde,
\[ EP_{RAR} = \sqrt{\frac{p0\left(1 - p0\right)}{n_{1}}+\frac{p1\left(1 - p1\right)}{n_{2}}} \]
Usando os dados do Physicians’ Health Study, pode-se criar um script no RStudio para os cálculos:
Vetor dos dados
Matriz dos dados64
mat_iam <- matrix (dados, byrow = TRUE, nrow = 2)
tratamento <- c ("aspirina", "placebo")
desfecho <- c ("IAM", "s/IAM")
rownames (mat_iam) <- tratamento
colnames (mat_iam) <- desfecho
mat_iam## IAM s/IAM
## aspirina 139 10898
## placebo 239 10795Cálculo das incidências no grupo tratamento e no grupo placebo
Na matriz o que está entre colchetes [1,1] significa: linha 1 e coluna 1, ou seja, o valor 139.
## [1] 11037## [1] 0.0126## [1] 11034## [1] 0.0217Os resultados da matriz de dados e o cálculo das incidências p0 (incidência no grupo placebo) e p1 (incidência no grupo de tratamento) já eram conhecidos e foram repetidos apenas para entrar na programação do cálculo do IC95%.
Cálculo do erro padrão da RAR
## [1] 1.96EP_RAR <- sqrt((((p0*(1-p0)) / n0)) + (((p1*(1-p1)) / n1)))
# Limite inferior
li_RAR <- RAR - (z * EP_RAR)
round (li_RAR, 4)## [1] 0.0056## [1] 0.0125## [1] 0.005645 0.009066 0.012488## [1] 0 0 0Portando, ao Redução Absoluta do Risco foi igual a 0,0091 (IC95%: 0,0056-0,0125). A partir destes resultados, pode-se calcular o intervalo de confiança para o NNT:
## [1] 80.07881## [1] 177.1497Concluindo, o uso da aspirina no Physicians’ Health Study reduziu o risco de infarto agudo do miocárdio em 42% (RRR), ou seja, foi eficaz. Por outro lado, para ter este impacto será necessário tratar 110 (IC95%: 80-177) pacientes para que um tenha benefício. Este NNT é grande; o ideal é um NNT < 10. Apesar disso, como a aspirina tem baixo custo e seus benefícios suplantam os efeitos adversos, seu uso pode estar justificado.
18.6.4 Número Necessário para Causar Dano
Deve-se comparar o NNT com o Número Necessário para causar Dano (NND), em inglês, Number Needed to Harm (NNH). Deve ser interpretado como o número de pacientes tratados para que um deles apresente um efeito adverso.
O NND é calculado pela recíproca do aumento absoluto do risco (ARA), equivalente a diferença de risco ou redução absoluta do risco:
\[ NND = \frac{1}{ARA} = \frac{1}{I_{expostos} - I_{não \ expostos}} \]
18.6.4.1 Exemplo
No Physicians’ Health Study sobre o uso de aspirina na prevenção de IAM, foi verificado também os efeitos colaterais da aspirina, como acidentes vasculares cerebrais (AVC), Figura 18.13.
Figura18.13: Physicians’ Health Study, componente aspirina e AVC.
Cálculo das incidências
## [1] 0.0089## [1] 0.0108Para o cálculo do NND, usa-se a função risk(), como mencionado antes:
## p0 p1 RR OR RRI ARI NNH
## 0.0089 0.0108 1.2135 1.2158 0.2135 0.0019 526.3158Os resultados mostram que o NND65 é igual a 526. Ou seja, para evitar um IAM há necessidade de tratar 110 pacientes e a cada 526 tratados espera-se um caso de AVC, havendo um benefício bem maior quando comparado ao risco de AVC.
18.7 Análise de sobrevida
A análise de sobrevida é utilizada quando se pretende investigar o tempo entre o início de um estudo e a ocorrência subsequente de um evento que modifica o estado de saúde do indivíduo. É bastante usada em estudos sobre câncer, por exemplo, analisando o tempo desde a cirurgia até a morte, o tempo desde o início do tratamento até a progressão da doença, o tempo desde a resposta até a recorrência da doença. Ela também é usada para medir a ocorrência de outros eventos como o tempo desde a infecção pelo vírus da imunodeficiência humana (HIV) até o desenvolvimento da Síndrome de Imunodeficiência Adquirida (SIDA), o tempo de hospitalização, tempo de amamentação, etc.
O interesse está centrado na verificação do efeito dos fatores de risco ou de prognóstico sobre o tempo de sobrevida de um indivíduo ou de um grupo, bem como definir as probabilidades de sobrevida em diversos momentos no seguimento do grupo. Considera-se tempo de sobrevida, ou simplesmente sobrevida, o tempo a entre a entrada do indivíduo no estudo e a ocorrência do evento de interesse. Com relação aos dados relacionados ao tempo, podem ocorrer problemas. O tempo para um evento geralmente não tem distribuição normal. Além disso, nem sempre se pode esperar até que o evento ocorra em todos os pacientes e alguns pacientes abandonam o estudo mais cedo. Todos devem ser considerados e as análises de sobrevida contornam esses problemas.
Em estudos de sobrevida, os indivíduos são observados até a ocorrência de um evento final que, geralmente, corresponde à morte, ou à variação de um parâmetro biológico ou outro evento que indique a modificação do estado inicial (cura, recorrência, retorno ao trabalho, etc.) O evento final é denominado de falha, por referir-se, em geral, a algo indesejável.
18.7.1 Dados Censurados
Quando, em um estudo de sobrevida, os pacientes que saem do estudo ou que não vivenciam o evento são chamados de observações censuradas.
Esses tempos de sobrevida censurados subestimam o verdadeiro (mas desconhecido) tempo para o evento. Quando o evento (supondo que ocorreria) está além do final do período de acompanhamento, a censura costuma ser chamada de censura à direita.
A censura também pode ocorrer quando se observa a presença de um evento, mas não se sabe onde começou. Por exemplo, considere um estudo que investigue o tempo para a recorrência de um câncer após a remoção cirúrgica do tumor primário. Se os pacientes forem examinados 3 meses após a cirurgia e já tinham recorrência, então o tempo de sobrevida será censurado a esquerda, porque o tempo real (desconhecido) de recorrência ocorreu menos de 3 meses após a cirurgia.
Os dados de tempo do evento também podem ser censurados em intervalos, o que significa que os indivíduos entram e saem da observação. Se considerarmos o exemplo anterior e os pacientes também forem examinados aos 6 meses, aqueles que estão livres da doença aos 3 meses e perdem o acompanhamento entre 3 e 6 meses são considerados censurados no intervalo. A maioria dos dados de sobrevivência incluem observações censuradas à direita (186).
18.7.2 Método de Kaplan-Meier
O método de Kaplan-Meier (KM) é um método não paramétrico usado para estimar a probabilidade de sobrevivência a partir dos tempos de sobrevivência observados (187).
A função de sobrevida é a probabilidade de sobreviver a pelo menos um determinado ponto no tempo e o gráfico desta probabilidade é a curva de sobrevida. O método de sobrevida de Kaplan-Meier pode ser usado para comparar as curvas de sobrevida de dois ou mais grupos, como comparar um grupo tratado a um grupo não tratado (placebo), ou homens comparados a mulheres.
A curva de sobrevida KM, um gráfico da probabilidade de sobrevida de Kaplan-Meier em relação ao tempo, fornece um resumo útil dos dados que podem ser usados para estimar medidas como a mediana de sobrevida.
18.7.2.1 Pressupostos do método de Kaplan-Meier
Os pressupostos para o uso da análise de sobrevida são as seguintes (188):
- os participantes devem ser independentes, ou seja, cada participante aparece apenas uma vez no grupo;
- os grupos devem ser independentes, ou seja, cada participante está apenas em um grupo;
- todos os participantes são livres de eventos quando se inscrevem no estudo;
- a medição do tempo até o evento deve ser precisa;
- o ponto inicial e o evento são claramente definidos;
- as perspectivas de sobrevida dos participantes permanecem constantes, ou seja, os participantes inscritos no início ou no final do estudo devem ter as mesmas perspectivas de sobrevida;
- a probabilidade de censura não está relacionada à probabilidade do evento.
Como em todas as análises, se o número total de pacientes em qualquer grupo for pequeno, digamos menos de 30 participantes em cada grupo, os erros padrão em torno das estatísticas resumidas serão grandes e, portanto, as estimativas de sobrevida serão imprecisas. Para estudos de sobrevida, recomenda-se fazer o cálculo do tamanho amostral previamente. O R dispõe de um pacote que possibilita este cálculo, o powerSurvEpi (189).
18.7.2.2 Exemplo
O arquivo dadosSobrevida.xlsx contém as informações de 60 pacientes selecionados para um ensaio clínico randomizado hipotético de dois tratamentos nos quais 32 pacientes receberam o novo tratamento e 28 pacientes receberam o tratamento padrão. Para obter o arquivo, clique aqui e salve o mesmo em seu diretório de trabalho.
Destes pacientes, 33 eram mulheres e 27 homens. Durante o estudo (65 meses), um total de 21 pacientes morreram (7 mulheres e 14 homens).
Carregar o conjunto de dados
A partir do diretório de trabalho, carregue para um objeto que será denominado de sobrevida, usando a função read_excel() do pacote readxl e observe os dados com a função head().
## # A tibble: 6 × 5
## id evento tempo sexo grupo
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 22 0 5 fem novo
## 2 21 0 7 masc novo
## 3 19 0 8 fem novo
## 4 13 0 9 fem novo
## 5 50 1 9 masc novo
## 6 20 1 12 masc novoA Saída exibe um banco de dados com cinco variáveis:
- id \(\longrightarrow\) Identificação do indivíduo
- evento \(\longrightarrow\) Desfecho. 0 = censurado; 1 = morte
- tempo \(\longrightarrow\) Sobrevida em meses
- sexo \(\longrightarrow\) 1 = masculino; 2 = feminino
- grupo \(\longrightarrow\) Grupo de tratamento: 1 = nova droga; 2 = padrão
Construir uma tabela tratamento vs evento
## Evento
## Tratamento 0 1
## novo 24 8
## padrão 15 13A saída mostra o número em cada grupo, o número de eventos e o número censurados. Houve menos eventos, mas mais pacientes censurados no grupo do tratamento novo.
Calcular as estimativas de sobrevida de Kaplan-Meier para a construção da Curva de Sobrevida de cada tratamento
Para isso, usa-se a função survfit() do pacote survival(190). Seus principais argumentos incluem:
- objeto de sobrevida, criado usando a função
Surv()aninhada na funçãosurvfit() - e o conjunto de dados contendo as variáveis.
Para a construção da tabela e da curva de sobrevida, digite e execute o seguinte:
## Call: survfit(formula = Surv(tempo, evento) ~ grupo, data = sobrevida)
##
## grupo=novo
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 9 29 1 0.966 0.0339 0.9013 1.000
## 12 27 1 0.930 0.0479 0.8404 1.000
## 15 26 1 0.894 0.0579 0.7874 1.000
## 16 25 1 0.858 0.0657 0.7387 0.997
## 32 15 1 0.801 0.0826 0.6545 0.980
## 36 13 1 0.739 0.0965 0.5725 0.955
## 40 11 1 0.672 0.1086 0.4897 0.923
## 58 2 1 0.336 0.2438 0.0811 1.000
##
## grupo=padrão
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 28 3 0.893 0.0585 0.7853 1.000
## 2 25 1 0.857 0.0661 0.7369 0.997
## 3 24 1 0.821 0.0724 0.6911 0.976
## 4 23 2 0.750 0.0818 0.6056 0.929
## 7 20 1 0.712 0.0859 0.5625 0.902
## 17 19 1 0.675 0.0892 0.5210 0.875
## 21 17 2 0.596 0.0947 0.4361 0.813
## 38 9 1 0.529 0.1048 0.3592 0.780
## 52 2 1 0.265 0.1944 0.0628 1.000A Tabela de sobrevida é uma tabela descritiva com a coluna time, indicando o dia em que o evento ocorreu. A coluna n.risk indica o número de pacientes sob risco naquele momento. A coluna denominada n.event indica o número total de pacientes que sofreram o evento desde o início do estudo até o momento avaliado. A coluna survival indica a proporção de pacientes que sobreviveram desde o início do estudo até aquele momento. Por exemplo, a sobrevida cumulativa é de 0,801 aos 32 meses no grupo tratamento novo e de 0,529 aos 38 meses no grupo tratamento padrão.
O método Kaplan-Meier produz uma única estatística resumida do tempo de sobrevida, isto é, a média ou mediana. O tempo médio de sobrevida é estimado a partir dos tempos observados e é mostrado para cada grupo na tabela de médias e medianas para o tempo de sobrevida.
A sobrevida média é calculada como a soma do tempo dividido pelo número de pacientes que permanecem sem censura. Essa estatística pode ser usada para indicar o período de tempo em que um paciente pode sobreviver. O tempo mediano de sobrevida é o ponto em que metade dos pacientes experimentou o evento. Se a curva de sobrevida não cair para 0,5 (ou seja, probabilidade de sobrevida de 50%), o tempo mediano de sobrevida não poderá ser calculado.
Estes dados podem ser visualizados na Saída, obtida com o comando:
## records n.max n.start events rmean se(rmean) median 0.95LCL
## grupo=novo 32 32 32 8 49.92533 4.078218 58 40
## grupo=padrão 28 28 28 13 36.62437 5.403382 52 21
## 0.95UCL
## grupo=novo NA
## grupo=padrão NAVisualização da curva de sobrevida
Pode-se visualizar a curva (Figura 18.14) de uma maneira simples, utilizando a função plot() do pacote básico do R:
plot (tabsurv, col = c ("steelblue", "rosybrown"), lwd = 2)
legend (legend = c ("Tratamento novo", "Tratamento padrão"),
fill = c ("steelblue", "rosybrown"),
bty="n",
cex = 1,
y = 0.3,
x = 5)
Figura18.14: Curva de sobrevida comparando dois grupos de tratamento.
Outra maneira, mais sofisticada, de produzir a curva de KM é usando a função ggsurvplot(), incluída no pacote survminer (191) que utiliza o pacote ggplot2 (Figura 18.15)
Com essa função é possível mostrar:
- os intervalos de confiança de 95% da função de sobrevida, usando o argumento
conf.int = TRUE; - o número e/ou a porcentagem de indivíduos em risco por tempo, utilizando a opção
risk.table. Os valores permitidos para arisk.tableincluem: TRUEouFALSEespecificando se deve mostrar ou não a tabela de risco. O padrão é FALSE.absoluteoupercentage: para mostrar o número absoluto e o percentual de sujeitos em risco por tempo, respectivamente. Useabs_pctpara mostrar o número absoluto e a porcentagem.- o nrisk_cumcensor e nrisk_cumevents . Mostra o número em risco e o número acumulado de censura e eventos, respectivamente.
- o valor P do teste
Log-Rankcomparando os grupos usandopval = TRUE. - linha horizontal/vertical na sobrevida mediana usando o argumento
surv.median.line. Os valores permitidos incluem um dec(“nenhum”, “hv”, “h”, “v”). Onde v = vertical, h = horizontal.
ggsurvplot (tabsurv,
pval = TRUE,
conf.int = FALSE,
risk.table = "abs_pct",
risk.table.col = "strata",
surv.median.line = "hv",
ggtheme = theme_bw (),
legend.labs = c ("Tratamento Novo",
"Tratamento padrão"),
palette = c ("steelblue", "tomato"))## Warning in geom_segment(aes(x = 0, y = max(y2), xend = max(x1), yend = max(y2)), : All aesthetics have length 1, but the data has 2
## rows.
## ℹ Please consider using `annotate()` or provide
## this layer with data containing a single row.
## All aesthetics have length 1, but the data has 2
## rows.
## ℹ Please consider using `annotate()` or provide
## this layer with data containing a single row.
## All aesthetics have length 1, but the data has 2
## rows.
## ℹ Please consider using `annotate()` or provide
## this layer with data containing a single row.
## All aesthetics have length 1, but the data has 2
## rows.
## ℹ Please consider using `annotate()` or provide
## this layer with data containing a single row.
Figura18.15: Curva de sobrevida comparando dois grupos de tratamento, usando ggsurvplot().
O teste Log Rank pondera todos os pontos de tempo igualmente e é a estatística de sobrevida mais usada (192). O teste de log rank é um teste não paramétrico, que não faz suposições sobre as distribuições de sobrevivência. Os pressupostos deste teste são os mesmos do método de Kaplan-Meier. No exemplo, o valor P do teste é fornecido na Figura 18.15 e é igual a 0,083, ou seja, acima de 0,05, indicando não rejeição da \(H_{0}\). A hipótese nula diz que não há diferença na sobrevivência entre os dois grupos.
Essencialmente, o teste de log rank compara o número observado de eventos em cada grupo com o que seria esperado se a hipótese nula fosse verdadeira (ou seja, se as curvas de sobrevivência fossem idênticas). A estatística de log rank é aproximadamente distribuída como uma estatística de teste qui-quadrado.
A função survdiff(), também do pacote survival, pode ser usada para calcular o teste de log-rank comparando duas ou mais curvas de sobrevida e pode ser usado da seguinte forma:
## Call:
## survdiff(formula = Surv(tempo, evento) ~ grupo, data = sobrevida)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo=novo 32 8 11.91 1.28 3.01
## grupo=padrão 28 13 9.09 1.68 3.01
##
## Chisq= 3 on 1 degrees of freedom, p= 0.08A suposição de que o risco de um evento em um grupo em comparação com o outro grupo não muda ao longo do tempo é chamado de risco proporcional. Se as curvas de sobrevida se cruzam, isso sugere que os riscos não são proporcionais. Nessa situação, o teste log rank será menos poderoso e um teste alternativo deve ser considerado, como a Regressão de Cox ou Modelo de Riscos Proporcionais.
Regressão de Cox ou Modelo de Riscos Proporcionais
O modelo tem como objetivo a examinar simultaneamente como os fatores especificados influenciam a taxa de ocorrência de um determinado evento (por exemplo, infecção, morte) em um determinado ponto no tempo. Essa taxa é referida como hazard ratio.
Geralmente, as variáveis preditoras (ou fatores) são denominadas covariáveis.
O modelo de Cox é expresso pela função de risco denotada por h(t). Pode ser interpretada como o risco de morrer no tempo t e estimada da seguinte forma:
\[ h\left(t\right) = h_{0} \left(t\right) \times e^{\left( {b_{1}x_{1}+b_{2}x_{2}+...+b_{n}x_{n}} \right)} \]
Onde,
- t é o tempo de sobrevida, indica que o risco varia com o tempo;
- h(t) é a função de risco (hazard) determinada por um conjunto de n covariáveis (\(x_{1}, x_{2}, ..., x_{n}\));
- Os coeficientes (\(b_{1}, b_{2}, ..., b_{n}\)) medem o tamanho do efeito das covariáveis;
- h(0) é o risco basal, o valor do risco se todos os \(x_{i}\) fossem iguais a zero (\(exp(0) = 1\)).
As quantidades exp(\(b_{i}\)) são chamadas de hazard ratio (HR). Uma hazard ratio acima de 1 indica uma covariável que está positivamente associada à probabilidade do evento e, portanto, negativamente associada ao tempo de sobrevida.
Resumindo,
- HR = 1: Sem efeito
- HR <1: Redução do risco
- HR> 1: Aumento do risco
Para calcular o modelo de Cox no R serão utilizados os mesmos dados da do arquivo dadosSobrevida.xlsx.
O pacote survival tem uma função para calcular o modelo de Cox, coxph(), que usa os argumentos:
- formula \(\longrightarrow\) é o modelo linear com um objeto de sobrevivida como variável desfecho. O objeto de sobrevida é criado usando a função
Surv()como segue:Surv(tempo, evento). - data \(\longrightarrow\) um banco de dados contendo as variáveis.
## Call:
## coxph(formula = Surv(tempo, evento) ~ grupo, data = sobrevida)
##
## coef exp(coef) se(coef) z p
## grupopadrão 0.7698 2.1593 0.4505 1.709 0.0875
##
## Likelihood ratio test=3.03 on 1 df, p=0.08171
## n= 60, number of events= 21A função summary() fornece um relatório mais completo:
## Call:
## coxph(formula = Surv(tempo, evento) ~ grupo, data = sobrevida)
##
## n= 60, number of events= 21
##
## coef exp(coef) se(coef) z Pr(>|z|)
## grupopadrão 0.7698 2.1593 0.4505 1.709 0.0875 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## grupopadrão 2.159 0.4631 0.893 5.221
##
## Concordance= 0.637 (se = 0.054 )
## Likelihood ratio test= 3.03 on 1 df, p=0.08
## Wald test = 2.92 on 1 df, p=0.09
## Score (logrank) test = 3.06 on 1 df, p=0.08Os resultados da regressão de Cox, podem ser interpretados da seguinte forma:
Significância estatística. A coluna marcada com z fornece o valor da estatística Wald. Corresponde à razão de cada coeficiente de regressão para seu erro padrão (\(z = \frac{coef}{EP_{coef}}\)). A estatística Wald avalia se o coeficiente beta (\(\beta\)) de uma determinada variável é estatisticamente diferente de 0. A partir da saída, pode-se concluir que não há diferença estatisticamente significativa entre os grupos (P = 0,0875).
Coeficientes de regressão. A seguir deve-se observar, no modelo de Cox, o sinal dos coeficientes de regressão (coef). Um sinal positivo significa que o hazard (risco) é maior e, portanto, pior o prognóstico, para sujeitos com valores mais elevados dessa variável. No exemplo, a variável grupo é codificada como 1=novo, 2=padrão. O resumo do modelo de Cox fornece a hazard ratio (HR) para o segundo grupo em relação ao primeiro grupo, ou seja, tratamento padrão versus tratamento novo. O coeficiente beta para grupo = 0,7698 indica que os indivíduos do tratamento padrão têm maior risco de morte (taxas de sobrevivência mais baixas) do que os do grupo tratamento novo, nesses dados. Entretanto, esta diferença não é estatisticamente significativa.
Hazard ratios. Os coeficientes exponenciados (exp(coef) = exp(0,7698) = 2,1593), também conhecidos como hazard ratio, fornecem o tamanho do efeito das covariáveis. Por exemplo, ser do grupo padrão aumenta o risco por um fator de 2,1593. Se esta diferença fosse significativa (P < 0,05), pertencer ao grupo padrão estaria associado a um mau prognóstico.
Intervalos de confiança das taxas de risco. O resultado do resumo também fornece intervalos de confiança de 95% para a razão de risco (exp(coef)), limite inferior de 95% = 0,893, limite superior de 95% = 5,221, mostrando a não significância estatística, pois cruza o 1.
Significância estatística global do modelo. Finalmente, a saída fornece valores de P para três testes alternativos para significância geral do modelo: O teste de razão de verossimilhança (Likelihood ratio test), teste de Wald e a estatística logrank. Esses três métodos são equivalentes. Para um tamanho amostral grande, eles darão resultados semelhantes. Para n pequeno, eles podem diferir um pouco. O teste de razão de verossimilhança tem melhor comportamento para tamanhos de amostra pequenos, por isso é geralmente preferido.
18.8 Regressão logística binária
A regressão logística (também conhecida como regressão logit ou modelo logit) foi desenvolvida como um a extensão do modelo linear pelo estatístico David Cox em 1958. Pertence a uma família, denominada Modelo Linear Generalizado (GLM) e é um modelo de regressão em que o desfecho Y é categórico.
A regressão logística permite estimar a probabilidade de uma resposta categórica com base em uma ou mais variáveis preditoras (X). Possibilita informar se a presença de um preditor aumenta (ou diminui) a probabilidade de um determinado desfecho em uma porcentagem específica. No caso em que Y é binário - ou seja, assume apenas dois valores, 0 e 1, que representam desfechos como aprovação/reprovação, sim/não, vivo/morto ou saudável/doente, tem-se a regressão logística binária.
Na regressão logística binária, as variáveis que afetam a probabilidade do resultado são medidas como Odds Ratio, que são chamadas de Odds Ratios ajustadas (193).
Na regressão linear, os valores das variáveis desfecho são preditos a partir de uma ou mais variáveis explicativas. Na regressão logística, uma vez que o desfecho é binário, a probabilidade de o desfecho ocorrer é calculada com base nos valores das variáveis explicativas. A regressão logística é semelhante à regressão linear na medida em que uma equação de regressão pode ser usada para prever a probabilidade de ocorrência de um desfecho. No entanto, a equação de regressão logística é expressa em termos logarítmicos (ou logits) e, portanto, os coeficientes de regressão devem ser convertidos para serem interpretados.
Embora as variáveis explicativas ou preditores no modelo possam ser variáveis contínuas ou categóricas, a regressão logística é mais adequada para medir os efeitos das exposições ou variáveis explicativas que são variáveis binárias. Variáveis contínuas podem ser incluídas, mas a regressão logística produzirá uma estimativa de risco para cada unidade de medida. Assim, a suposição de que o efeito de risco é linear sobre cada unidade da variável deve ser atendida e a relação não deve ser curva ou ter um ponto de corte sobre o qual o efeito ocorre. Além disso, as interações entre variáveis explicativas podem ser incluídas (193). Os casos em que a variável dependente tem mais de duas categorias de resultados podem ser analisados com regressão logística multinomial, não mostrada neste livro.
18.8.1 Pressupostos da regressão logística
O método de regressão logística assume que:
- O desfecho é uma variável binária ou dicotômica como sim vs não, positivo vs negativo, 1 vs 0.
- Existe uma relação linear entre o logit do desfecho e cada variável preditora. A função logit é \(logit\left(P\right) = log \ P\left(\frac{P}{\left(1 - P\right)}\right)\) , onde P é a probabilidade do desfecho. Na regressão linear, é assumido que o desfecho tem uma correlação linear com os preditores. Na regressão logística, o desfecho é categórico e, portanto, essa suposição é violada. Por isso que se usa o log (ou logit) dos dados. A suposição de linearidade na regressão logística, portanto, é que existe uma correlação linear entre quaisquer preditores contínuos e o logit da variável de desfecho.
- Os casos devem ser independentes. Os dados dos casos não devem ser relacionados; por exemplo, não pode medir as mesmas pessoas em pontos diferentes no tempo (medidas repetidas)
- Não há valores influentes (valores extremos ou outliers) nos preditores contínuos
- Não há altas intercorrelações (ou seja, multicolinearidade) entre os preditores.
Para melhorar a precisão de seu modelo, você deve certificar-se de que essas suposições sejam verdadeiras para seus dados.
18.8.2 Dados para a regressão logística
Os dados foram obtidos do banco de dados PimaIndiansDiabetes2, incluído no pacote mlbench (194) que contém 768 observações sobre 9 variáveis. Este conjunto de dados é originalmente do Instituto Nacional de Diabetes e Doenças Digestivas e Renais. São dados de mulheres com 21 anos ou mais com herança indígena Pima.
As variáveis traduzidas são:
- gesta \(\longrightarrow\) Número de vezes que engravidou
- glicose \(\longrightarrow\) Concentração de glicose plasmática após 2 horas de um teste oral de tolerância à glicose (mg%)
- pd \(\longrightarrow\) Pressão arterial diastólica (mm Hg)
- tríceps \(\longrightarrow\) Espessura da dobra cutânea do tríceps (mm)
- insulina \(\longrightarrow\) Insulina sérica após 2 horas (mu U/ml)
- imc \(\longrightarrow\) Índice de massa corporal (peso em kg/(altura em m)2)
- pedigree \(\longrightarrow\) Função de pedigree de diabetes
- idade \(\longrightarrow\) Idade (anos)
- diabetes \(\longrightarrow\) 0 = neg, 1 = pos
Depois de removidos os dados omissos, o banco de dados ficou como encontrado em dadosPima.xlsx. Estes dados serão utilizados para executar uma regressão logística binária para verificar se alguma dessas variáveis podem predizer a presença de diabetes em mulheres com herança Pima.
Para obter arquivo, clique aqui e salve o mesmo em seu diretório de trabalho.
18.8.3 Preparação dos dados
Carregar os dados, usando a função read_excel () do pacote readxl:
## tibble [392 × 9] (S3: tbl_df/tbl/data.frame)
## $ gesta : num [1:392] 1 0 3 2 1 5 0 1 1 3 ...
## $ glicose : num [1:392] 89 137 78 197 189 166 118 103 115 126 ...
## $ pd : num [1:392] 66 40 50 70 60 72 84 30 70 88 ...
## $ triceps : num [1:392] 23 35 32 45 23 19 47 38 30 41 ...
## $ insulina: num [1:392] 94 168 88 543 846 175 230 83 96 235 ...
## $ imc : num [1:392] 28.1 43.1 31 30.5 30.1 25.8 45.8 43.3 34.6 39.3 ...
## $ pedigree: num [1:392] 0.167 2.288 0.248 0.158 0.398 ...
## $ idade : num [1:392] 21 33 26 53 59 51 31 33 32 27 ...
## $ diabetes: num [1:392] 0 1 1 1 1 1 1 0 1 0 ...Para realizar a regressão logística, os dados devem ser inseridos como na regressão linear: organizados em colunas (uma representando cada variável). Olhando a saída da função glimpse, nota-se que as variáveis categóricas foram carregadas como <dbl> (numéricas). A variável diabetes foi codificada como 1 (evento ocorreu) e 0 (evento não ocorreu); neste caso, 1 representa ter diabetes e 0 representa uma ausência de diabetes. Como o R codifica os fatores na ordem 0 e 1, mantém-se assim e transforma-se a variável diabetes em fator.
As variáveis idade e gesta (número de gravidezes) são altamente assimétricas (Figura: 18.16):
# Este comando coloca os gráficos em uma mesma linha, em duas colunas:
par(mfrow=c(1,2))
# Histogramas
hist(dados$idade, breaks = 8, main = "", ylab = "Frequência", xlab = "Idade (anos)")
hist(dados$gesta, breaks = 12, main = "", ylab = "Frequência", xlab = "Nº de gravidezes")
Figura18.16: Histogramas das variáveis idade e número de gravidezes.
Por isso, elas serão categorizadas como fatores com a função cut(), consultar a seção que trata da construção de tabelas de frequência:
dados$idadeCateg <- cut (dados$idade,
breaks= c (20,31,41,51,81),
labels = c ("20-30", "31-40","41-50", ">50"),
right = FALSE,
include.lowest = TRUE)
dados$gestaCateg <- cut (dados$gesta,
breaks= c (0,6,11,17),
labels = c ("0-5", "6-10", ">10"),
right = FALSE,
include.lowest = TRUE)Após a visualização e modificações das variáveis, serão retiradas as variáveis contínuas idade e gesta, mantendo as demais variáveis:`:
## Rows: 392
## Columns: 9
## $ glicose <dbl> 89, 137, 78, 197, 189, 166, 118, 103, 115, 126, 143, 125, 9…
## $ pd <dbl> 66, 40, 50, 70, 60, 72, 84, 30, 70, 88, 94, 70, 66, 82, 76,…
## $ triceps <dbl> 23, 35, 32, 45, 23, 19, 47, 38, 30, 41, 33, 26, 15, 19, 36,…
## $ insulina <dbl> 94, 168, 88, 543, 846, 175, 230, 83, 96, 235, 146, 115, 140…
## $ imc <dbl> 28.1, 43.1, 31.0, 30.5, 30.1, 25.8, 45.8, 43.3, 34.6, 39.3,…
## $ pedigree <dbl> 0.167, 2.288, 0.248, 0.158, 0.398, 0.587, 0.551, 0.183, 0.5…
## $ diabetes <fct> 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1,…
## $ idadeCateg <fct> 20-30, 31-40, 20-30, >50, >50, >50, 31-40, 31-40, 31-40, 20…
## $ gestaCateg <fct> 0-5, 0-5, 0-5, 0-5, 0-5, 0-5, 0-5, 0-5, 0-5, 0-5, >10, 6-10…18.8.4 Criação do modelo de regressão (enter)
Agora, pode-se prosseguir com a análise com uma regressão logística binária, usando a função glm() do pacote stats que pode usar vários argumentos:
- formula \(\longrightarrow\) objeto da classe “formula”. Um preditor típico tem a forma “resposta ~ preditor” em que resposta, na regressão logística binária, é uma variável dicotômica e o preditor pode ser uma série de variáveis numéricas ou categóricas;
- family \(\longrightarrow\) uma descrição da distribuição de erro e função de link a ser usada no modelo glm, pode ser uma string que nomeia uma função de family. O padrão é
family = gaussian(). No caso da regressão logística bináriafamily = binomial()oufamily = binomial (link = ”logit”). Para outras informações, usehelp(glm)ouhelp(family); - data \(\longrightarrow\) banco de dados.
- … \(\longrightarrow\) …
A função glm() diz ao R para executar um modelo linear generalizado. Dentro dos parênteses, são fornecidas informações importantes sobre o modelo. À esquerda do til (~) está a variável dependente. Deve ser codificada como 0 e 1 para a função ler como binária. Após o til (~), são listadas as variáveis preditoras. Quando depois do sinal til é colocado um ponto (~.), significa a presença de todas as variáveis preditoras. Quando é colocado o asterisco (*) entre duas variáveis preditoras, isto indica que não se deseja apenas o efeito principal, mas também um termo de interação entre elas. No exemplo, não foi pedido essa análise. Finalmente, após a vírgula, especifica-se que a distribuição é binomial. A função link padrão em glm para uma variável desfecho binomial é o logit, portanto, não há necessidade de especificar no modelo.
Inicialmente, será feita uma regressão logística do tipo entrada forçada (enter), método padrão de conduzir uma regressão que consiste em simplesmente colocar todos os preditores no modelo de regressão em um bloco e estimar parâmetros para cada um (195). O dataframe dados será usado no modelo com todos os preditores dentro da função. O objeto criado será denominado de mod1.
Para ver o modelo gerado, há necessidade de executar a função summary():
##
## Call:
## glm(formula = diabetes ~ ., family = binomial(), data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.8655465 1.1691937 -7.583 3.39e-14 ***
## glicose 0.0392471 0.0059152 6.635 3.25e-11 ***
## pd -0.0045595 0.0119562 -0.381 0.70294
## triceps 0.0147965 0.0173948 0.851 0.39497
## insulina -0.0006838 0.0013559 -0.504 0.61402
## imc 0.0630043 0.0273349 2.305 0.02117 *
## pedigree 1.0165275 0.4388864 2.316 0.02055 *
## idadeCateg31-40 0.8544642 0.3767035 2.268 0.02331 *
## idadeCateg41-50 1.5745534 0.5203173 3.026 0.00248 **
## idadeCateg>50 1.3840205 0.6367486 2.174 0.02974 *
## gestaCateg6-10 -0.2430053 0.4202225 -0.578 0.56308
## gestaCateg>10 0.8284550 0.7673085 1.080 0.28028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 337.98 on 380 degrees of freedom
## AIC: 361.98
##
## Number of Fisher Scoring iterations: 5Essas estatísticas resumidas ajudam a entender melhor o modelo, fornecendo-nos as seguintes informações:
- Distribuição dos desvios residuais;
- Estimativas do Intercepto e inclinação junto com o erro padrão, valor z (estatística de Wald) e valor P;
- Valor AIC e
- Desviância residual e desviância nula.
Em uma regressão logística, a resposta que está sendo modelada é o log(odds) ou logit de que o desfecho é igual a 1. Os coeficientes de regressão fornecem a mudança no log(odds) no desfecho para a mudança de uma unidade na variável preditora, mantendo todas as outras variáveis preditivas constantes (196).
Nas variáveis contínuas, para cada aumento de uma unidade na concentração da glicose, por exemplo, o log(odds) de ser diabético “1 = pos” (versus ser diabético “0 = neg”) aumenta em 0,039. Da mesma forma, para um aumento de unidade na pressão diastólica, a probabilidade de log(odds) de ser diabético “1 = pos” (versus ser diabético “0 = neg”) diminui em 0,0045, pois o coeficiente é negativo.
Para variáveis categóricas, o desempenho de cada categoria é avaliado em relação a uma categoria de base. A categoria de base para a variável idadeCateg é 20–30 (primeira categoria que aparece quando se observa os níveis 66) e para gestaCateg é 0–5. A interpretação de tais variáveis é a seguinte:
- Estar na faixa etária de 31-40 anos, versus faixa etária de 20-30, muda o
log(odds)de ser diabético “1 = pos” (versus ser diabético “0 = neg”) em 0,854. - * Estar na faixa de 6–10 gravidezes, versus 0–5 gravidezes, muda olog(odds)de ser diabético “pos” (versus ser diabético “neg”) em -0,24.
Como o log(odds) é difícil de interpretar, é possível exponenciar o log(odds) para colocar os resultados em uma escala de odds. O R calcula as odds ratios e seus IC95, usando a função exp() e, dentro dela, a função coef() – do coeficiente b para as variáveis preditoras – e a função confint(), que fornece os intervalos de confiança. A função cbind() combinará as duas. A função round() é dispensável, foi usada apenas para que o resultado tenha 5 dígitos:
## OR 2.5 % 97.5 %
## (Intercept) 0.00014 0.00001 0.00127
## glicose 1.04003 1.02846 1.05267
## pd 0.99545 0.97246 1.01935
## triceps 1.01491 0.98076 1.05019
## insulina 0.99932 0.99667 1.00200
## imc 1.06503 1.00999 1.12487
## pedigree 2.76358 1.18886 6.64511
## idadeCateg31-40 2.35011 1.11963 4.92734
## idadeCateg41-50 4.82858 1.75113 13.60151
## idadeCateg>50 3.99091 1.17784 14.52935
## gestaCateg6-10 0.78427 0.33953 1.77546
## gestaCateg>10 2.28978 0.53772 11.29542A função exp() transformou o log(odds) em odds. A Saída mostra que a chance de ser diabético aumenta em um fator de 1,04 para um aumento de uma unidade na concentração de glicose, mantendo todas as outras variáveis constantes. Todo o IC95% encontra-se acima de 1. Isto ocorre com a triceps (prega cutânea tricipital), com o imc, pedigree (herança Pima), idade acima de 30 anos, indicando haver significância estatística, para essas variáveis, até este momento da análise.
18.8.5 Criação do modelo passo a passo (stepwise)
Existem vários métodos para seleção de variáveis, aqui será usado o método passo a passo (stepwise). O procedimento de seleção é realizado automaticamente por pacotes estatísticos como o pacote MASS (197) através da função stepAIC() que utiliza o AIC (Critério de Informação de Akaike 67) para a seleção (198). Esta função usa os seguintes argumentos:
- object \(\longrightarrow\) um objeto que representa um modelo de uma classe apropriada;
- direction \(\longrightarrow\) o modo de pesquisa passo a passo pode ser “backward” (para trás), “forward” (para frente) ou “both” (ambos), que é o padrão . Se o argumento
scopeestiver faltando, o padrão para a direção é “backward”; - scope \(\longrightarrow\) define a gama de modelos examinados na pesquisa passo a passo
- … \(\longrightarrow\) para maiores informações, use o help.
Execute o comando abaixo para realizar a seleção do modelo, que será recebido pelo objeto mod2:
## Start: AIC=361.98
## diabetes ~ glicose + pd + triceps + insulina + imc + pedigree +
## idadeCateg + gestaCateg
##
## Df Deviance AIC
## - pd 1 338.13 360.13
## - insulina 1 338.24 360.24
## - gestaCateg 2 340.36 360.36
## - triceps 1 338.71 360.71
## <none> 337.98 361.98
## - imc 1 343.41 365.41
## - pedigree 1 343.61 365.61
## - idadeCateg 3 349.12 367.12
## - glicose 1 391.66 413.66
##
## Step: AIC=360.13
## diabetes ~ glicose + triceps + insulina + imc + pedigree + idadeCateg +
## gestaCateg
##
## Df Deviance AIC
## - insulina 1 338.36 358.36
## - gestaCateg 2 340.43 358.43
## - triceps 1 338.84 358.84
## <none> 338.13 360.13
## - imc 1 343.49 363.49
## - pedigree 1 343.89 363.89
## - idadeCateg 3 349.25 365.25
## - glicose 1 391.88 411.88
##
## Step: AIC=358.36
## diabetes ~ glicose + triceps + imc + pedigree + idadeCateg +
## gestaCateg
##
## Df Deviance AIC
## - gestaCateg 2 340.79 356.79
## - triceps 1 339.10 357.10
## <none> 338.36 358.36
## - imc 1 343.49 361.49
## - pedigree 1 344.00 362.00
## - idadeCateg 3 349.53 363.53
## - glicose 1 406.62 424.62
##
## Step: AIC=356.79
## diabetes ~ glicose + triceps + imc + pedigree + idadeCateg
##
## Df Deviance AIC
## - triceps 1 341.42 355.42
## <none> 340.79 356.79
## - imc 1 346.32 360.32
## - pedigree 1 346.82 360.82
## - idadeCateg 3 361.77 371.77
## - glicose 1 408.50 422.50
##
## Step: AIC=355.42
## diabetes ~ glicose + imc + pedigree + idadeCateg
##
## Df Deviance AIC
## <none> 341.42 355.42
## - pedigree 1 347.78 359.78
## - imc 1 355.10 367.10
## - idadeCateg 3 363.70 371.70
## - glicose 1 409.44 421.44##
## Call:
## glm(formula = diabetes ~ glicose + imc + pedigree + idadeCateg,
## family = binomial(), data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.995423 1.008994 -8.915 < 2e-16 ***
## glicose 0.037429 0.005086 7.360 1.84e-13 ***
## imc 0.072975 0.020377 3.581 0.000342 ***
## pedigree 1.055869 0.430381 2.453 0.014154 *
## idadeCateg31-40 0.818513 0.334366 2.448 0.014367 *
## idadeCateg41-50 1.633945 0.399650 4.088 4.34e-05 ***
## idadeCateg>50 1.352995 0.528396 2.561 0.010450 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 341.42 on 385 degrees of freedom
## AIC: 355.42
##
## Number of Fisher Scoring iterations: 5Depois de implementar a função stepAIC(), restam agora quatro variáveis independentes - glicose, imc, pedigree e idadeCateg. De todos os modelos possíveis, este modelo (mod2) possui o valor mínimo de AIC. Além disso, as variáveis selecionadas têm valor P < 0,05.
Observe o cálculo das odds ratios com o mod2 e compare com o mod1
## OR 2.5 % 97.5 %
## (Intercept) 0.00012 0.00002 0.00082
## glicose 1.03814 1.02820 1.04897
## imc 1.07570 1.03450 1.12092
## pedigree 2.87447 1.25945 6.80648
## idadeCateg31-40 2.26713 1.17331 4.36901
## idadeCateg41-50 5.12405 2.35800 11.37105
## idadeCateg>50 3.86900 1.40027 11.30915Na Saída, observa-se que ter idade entre 41 a 50 anos tem uma chance 5,1 (IC95%: 2,4 – 11,4) maior de ser diabético comparado com a faixa etária 20-30 anos. As odds ratios ajustados da regressão logística binária fornecem uma estimativa que não é enviesada por confusão.
18.8.6 Uso do modelo para fazer predições
A regressão logística fornece uma curva de probabilidade em forma de S (Figura 18.17). Pode ser utilizada a função predict(), uma função R genérica para fazer previsões a partir de modelos de ajuste. Se nenhum conjunto de dados for fornecido à função predict(), as probabilidades são calculadas a partir dos dados que foram usados para ajustar o modelo de regressão logística, no caso o mod2. A função predict() toma como argumentos o modelo de regressão e os valores da variável preditora (199). Para a função mutate() e o operador pipe, consulte help().
- object \(\longrightarrow\) um objeto que representa um modelo de regressão
- data \(\longrightarrow\) banco de dados
- type \(\longrightarrow\) “reponse” ou “terms”
- … \(\longrightarrow\)…
Esta função verifica o impacto da variação dos níveis de uma variável preditora sobre a probabilidade do desfecho. Por exemplo, verifica o impacto da concentração sérica de glicose na probabilidade de ser diabético.
predito <- predict(mod2, data=dados, type = "response")
dados %>%
mutate(predito = ifelse(diabetes == "1", 1, 0)) %>%
ggplot(aes(glicose, predito)) +
geom_point() +
geom_smooth(method = "glm", method.args = list(family = "binomial")) +
labs(title = NULL,
x = "Glicemia (mg%)",
y = "Probabilidade de diabetes")## `geom_smooth()` using formula = 'y ~ x'
Figura18.17: Probabilidade de diabetes de acordo com a glicemia em mulheres de herança Pima.
Uma vez que os coeficientes para a glicose foram estimados para o mod2, é uma questão simples calcular a probabilidade de diabetes para qualquer concentração da glicose. Por exemplo, qual a probabilidade de uma mulher da etnia Pima ser diabética com uma glicose de 180 mg%?
Em primeiro lugar constrói-se um modelo logístico novo, usando apenas a glicose como preditor:
##
## Call:
## glm(formula = diabetes ~ glicose, family = binomial(link = "logit"),
## data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.095521 0.629787 -9.679 <2e-16 ***
## glicose 0.042421 0.004761 8.911 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 386.67 on 390 degrees of freedom
## AIC: 390.67
##
## Number of Fisher Scoring iterations: 4A seguir, utiliza-se a glicose = c(90, 120, 150, 180), criando assim um novo dataframe, chamado diabetes.
## glicose
## 1 90
## 2 120
## 3 150
## 4 180A partir daí calcula-se as predições, colocando, no dataframe diabetes, uma variável de probabilidade de predição (pred.prob):
## glicose pred.prob
## 1 90 0.09299244
## 2 120 0.26795891
## 3 150 0.56651015
## 4 180 0.82350197A Saída mostra que a probabilidade de uma mulher de herança Pima ter diabete com uma glicemia de 90 mg% é 9,2%. Esta probabilidade sobe para 82,3% quando a glicemia é de 180 mg%.
Esta ferramenta permite também que se faça a comparação das probabilidades usando outros preditores. Para isso basta acrescentar ao modelo de regressão logística outras variáveis, por exemplo, IMC igual a 25.
##
## Call:
## glm(formula = diabetes ~ glicose + imc, family = binomial(link = "logit"),
## data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.301460 0.927405 -8.951 < 2e-16 ***
## glicose 0.040713 0.004825 8.437 < 2e-16 ***
## imc 0.071794 0.019606 3.662 0.00025 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 372.12 on 389 degrees of freedom
## AIC: 378.12
##
## Number of Fisher Scoring iterations: 4Acrescentar o ICM = 25 aos dados:
## glicose imc
## 1 90 25
## 2 120 25
## 3 150 25
## 4 180 25A seguir, faz-se as predições:
## glicose imc pred.prob
## 1 90 25 0.05507376
## 2 120 25 0.16506181
## 3 150 25 0.40139839
## 4 180 25 0.69460853Ou seja, uma mulher com herança Pima e um ICM = 25 kg/m2 e uma glicemia de 180 mg% tem 69,4% de probabilidade de ser diabética.
Qual a probabilidade de ter diabete mudando o IMC para 30?
Tente responder, seguindo os passos mostrados!
18.8.7 Avaliação do modelo
18.8.7.1 Linearidade
Esta suposição pode ser testada examinando se o termo de interação entre o preditor e sua log transformação é significativo (200). Portanto há necessidade de criar os termos de interação de cada uma das variáveis independentes contínuas que estão no mod2 com seu log, usando a função log().
Glicose
IMC
Pedigree
Essas variáveis foram adicionadas ao conjunto de dados dados. podem ser observadas com a função str().
## tibble [392 × 12] (S3: tbl_df/tbl/data.frame)
## $ glicose : num [1:392] 89 137 78 197 189 166 118 103 115 126 ...
## $ pd : num [1:392] 66 40 50 70 60 72 84 30 70 88 ...
## $ triceps : num [1:392] 23 35 32 45 23 19 47 38 30 41 ...
## $ insulina : num [1:392] 94 168 88 543 846 175 230 83 96 235 ...
## $ imc : num [1:392] 28.1 43.1 31 30.5 30.1 25.8 45.8 43.3 34.6 39.3 ...
## $ pedigree : num [1:392] 0.167 2.288 0.248 0.158 0.398 ...
## $ diabetes : Factor w/ 2 levels "0","1": 1 2 2 2 2 2 2 1 2 1 ...
## $ idadeCateg : Factor w/ 4 levels "20-30","31-40",..: 1 2 1 4 4 4 2 2 2 1 ...
## $ gestaCateg : Factor w/ 3 levels "0-5","6-10",">10": 1 1 1 1 1 1 1 1 1 1 ...
## $ glicoseInt : num [1:392] 399 674 340 1041 991 ...
## $ imcInt : num [1:392] 93.7 162.2 106.5 104.2 102.5 ...
## $ pedigreeInt: num [1:392] -0.299 1.894 -0.346 -0.292 -0.367 ...As variáveis pd, insulina, triceps e gestaCateg serão removidas, pois mostraram-se não significativas quando se criou o mod2 no método passo a passo (stepwise)::
## [1] "glicose" "imc" "pedigree" "diabetes" "idadeCateg"
## [6] "glicoseInt" "imcInt" "pedigreeInt"Com as variáveis restantes, será criado um quinto modelo , mod5, para examinar o comportamento da interação entre o preditor e sua log transformação:
##
## Call:
## glm(formula = diabetes ~ ., family = binomial, data = dados)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.80900 7.05884 -2.523 0.01164 *
## glicose 0.09919 0.23220 0.427 0.66924
## imc 0.90065 0.57720 1.560 0.11867
## pedigree 2.08633 0.77570 2.690 0.00715 **
## idadeCateg31-40 0.81470 0.33792 2.411 0.01591 *
## idadeCateg41-50 1.62746 0.40570 4.011 6.04e-05 ***
## idadeCateg>50 1.37610 0.54080 2.545 0.01094 *
## glicoseInt -0.01047 0.03948 -0.265 0.79091
## imcInt -0.17971 0.12529 -1.434 0.15149
## pedigreeInt -1.62801 0.95970 -1.696 0.08981 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.10 on 391 degrees of freedom
## Residual deviance: 336.65 on 382 degrees of freedom
## AIC: 356.65
##
## Number of Fisher Scoring iterations: 5A Saída exibe a parte que testa a suposição de linearidade dos logit. O interesse está centrado apenas nos termos de interação. Qualquer interação que seja significativa indica que o efeito principal violou o pressuposto de linearidade do logit. Todas as interações têm valores de significância (valor P) maiores que 0,05, indicando que o pressuposto de linearidade do logit foi atendido.
18.8.7.2 Multicolinearidade
A regressão logística é propensa ao efeito do viés da colinearidade e é essencial testar a colinearidade, após uma análise de regressão logística (201).
Os problemas de multicolinearidade consistem em incluir, no modelo, diferentes variáveis que possuem relação preditiva semelhante com o desfecho. Isso pode ser avaliado para cada preditor calculando o valor VIF (Variance Inflation Factor), usando a função vif() do pacote car (107) para o mod2 que tem mais de dois termos.
## GVIF Df GVIF^(1/(2*Df))
## glicose 1.050569 1 1.024973
## imc 1.023676 1 1.011769
## pedigree 1.026294 1 1.013062
## idadeCateg 1.088364 3 1.014213Qualquer variável com um valor VIF alto (acima de 5 ou 10) deve ser removida do modelo. A Saída mostra que todos os valores estão abaixo de 10. Isso leva a um modelo mais simples sem comprometer a precisão do modelo, o que é bom.
18.8.7.3 Estatística \(R^2\)
R ao quadrado é uma métrica útil para regressão linear simples e múltipla, mas não tem o mesmo significado na regressão logística.
Os estatísticos descobriram uma variedade de análogos do R ao quadrado para regressão logística que eles referem, coletivamente, como pseudo R ao quadrado. Não têm a mesma interpretação, na medida em que não são simplesmente a proporção da variância explicada pelo modelo.
Infelizmente, existem muitas maneiras diferentes de calcular um \(R^2\) para regressão logística e nenhum consenso sobre qual é a melhor. Os dois métodos mais frequentemente relatados em softwares estatísticos são um proposto por McFadden (1974) e outro por Cox-Snell (1989). No entanto, também é bastante relatado o de Nagelkerke. Os valores mais altos indicam um melhor ajuste do modelo.
O \(R^2\) de McFadden é uma versão, baseada no log-likelihood para o modelo somente com o intercepto e o modelo estimado completo. O \(R^2\) de Cox e Snell é baseado no log-likelihood para o modelo em comparação com o log-likelihood para um modelo basal. No entanto, com resultados categóricos, tem um valor máximo teórico inferior a 1, mesmo para um modelo “perfeito”. O \(R^2\) de Nagelkerke é uma versão ajustada do \(R^2\) de Cox e Snell que ajusta a escala da estatística para cobrir todo o intervalo de 0 a 1.
Em seguida, serão calculados os vários valores de \(R^2\), usando a função PseudoR2() do pacote DescTools:
## Nagelkerke McFadden CoxSnell
## 0.4580199 0.3145606 0.3294780Como se verifica, na Saída, todos os valores de \(R^2\) diferem ligeiramente, mas podem ser usados como medidas de tamanho de efeito para o modelo.
Então, basicamente, o pseudo R quadrado pode ser interpretado como \(R^2\), mas não se espera que seja tão grande. Uma regra prática, bastante útil é que o pseudo R quadrado de McFadden variando de 0,2 a 0,4 indica um ajuste muito bom do modelo (202).
18.8.7.4 Pesquisa de valores atípicos e pontos de alavancagem
Um outlier é uma observação que não é bem prevista pelo modelo de regressão ajustado (ou seja, tem um grande resíduo positivo ou negativo). Uma observação com um alto valor de alavancagem possui uma combinação incomum de valores preditores. Ou seja, é um outlier no espaço do preditor. O valor da variável dependente não é usado para calcular a alavancagem de uma observação.
Uma observação influente é uma observação que tem um impacto desproporcional na determinação dos parâmetros do modelo. As observações influentes são identificadas usando uma estatística chamada distância de Cook ou D de Cook.
A identificação dos outliers é feita essencialmente através dos resíduos padronizados, com a função rstandard():
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.96350 -0.63838 -0.36367 -0.08147 0.63008 2.68035Em uma amostra normalmente distribuída, ao redor de 95% dos valores estão entre –1,96 e +1,96, 99% deve estar entre –2,58 e +2,58 e quase todos (99,9%) devem situar-se entre –3,09 e +3,09. Portanto, resíduos padronizados com um valor absoluto maior que 3 são motivo de preocupação porque em uma amostra média é improvável que aconteça um valor tão alto por acaso (203). Na Saída, observa-se que os valores estão dentro de –3 e +3.
Essa avaliação pode ser acompanhada de um gráfico diagnóstico (veja seção da regressão linear para maiores detalhes), usando a função plot() para o modelo mod2.
Figura18.18: Gráfico diagnóstico dos resíduos e pontos de alavancagem.
São produzidos vários gráficos com a função plot(), mas o foco é o gráfico 5 (Figura 18.18) que confirma os achados de não existirem valores atípicos que comprometam o ajuste do modelo.
Para a análise dos pontos influentes, pode-se verificar os pontos de alavancagem (leverage) com a função hatvalues() do pacote stats, incluído no R base.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001847 0.005862 0.013332 0.017857 0.023846 0.090653Os valores de alavancagem podem estar entre 0 (indicando que o caso não tem qualquer influência) e 1 (indicando que o caso tem grande influência). Se nenhum caso exercer influência indevida sobre o modelo, se espera que todos os valores de alavancagem estivessem próximos do valor médio. Alguns autores (204), recomendam investigar casos com valores superiores a duas vezes a média (2 x 0,0179 = 0,0358) como ponto de corte para identificar casos com influência indevida. Alguns valores estão acima de duas vezes a média. Entretanto, o maior valor está bem longe do valor igual a 1.
É interessante fazer a análise, observando junto a distância de Cook para ver se um ponto é um outlier significativo.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.220e-06 8.955e-05 5.703e-04 2.868e-03 2.908e-03 1.121e-01Se a distância de Cook é < 1, não há necessidade real de excluir esse ponto, uma vez que não tem um grande efeito na análise de regressão (205). Na Figura 18.18, verifica-se que os casos mais distantes não alcançam a distância de Cook.
18.8.7.5 Matriz de confusão
É uma representação tabular de valores observados versus valores previstos. Ajuda a quantificar precisão (ou acurácia) do modelo. Agora será realizada uma comparação dos valores observados de “diabetes” com os valores previstos.
Inicialmente, será criada uma variável correspondente aos valores previstos (mod2$fitted.values) classificando como “pos” se o valor ajustado exceder 0,5, caso contrário, “neg”.
Pode-se avaliar o desempenho do modelo (acurácia) com a comparação tabular entre os valores observados e os previstos:
tabDiabetes <- table(dados$diabetes, dados$predito)
rownames(tabDiabetes) <- c("Obs.neg", "Obs.pos")
colnames(tabDiabetes) <- c("Pred.neg", "Pred.pos")
tabDiabetes##
## Pred.neg Pred.pos
## Obs.neg 237 25
## Obs.pos 48 82A acurácia é obtida pela soma dos valores das caselas na diagonal dividido pelo total.
## [1] 0.8137755De acordo com a matriz de confusão, a acurácia do modelo é de 81,4%.
18.8.7.6 Avaliação do modelo com a Curva ROC
A curva ROC permite explicar o desempenho do modelo avaliando a sensibilidade versus especificidade.
Para se obter a curva ROC será usada a função roc() do pacote pROC (veja Estatísticas Diagnósticas).
Os comandos abaixo exibem em seus resultados (Figura 18.19) a curva ROC e a área sob a curva (AUC) e seus IC95%. Quanto maior a área sob a curva melhor é o poder de predição do modelo (206).
roc (dados$diabetes,
mod2$fitted.values,
plot=TRUE,
legacy.axes=FALSE,
print.auc=TRUE,
print.auc.y = 0.2,
ci = TRUE,
ylab="Sensibilidade",
xlab="1 - Especificdade",
col="steelblue",
lwd=2)## Setting levels: control = 0, case = 1## Setting direction: controls < cases
Figura18.19: Curva ROC do modelo de regressão.
##
## Call:
## roc.default(response = dados$diabetes, predictor = mod2$fitted.values, ci = TRUE, plot = TRUE, legacy.axes = FALSE, print.auc = TRUE, print.auc.y = 0.2, ylab = "Sensibilidade", xlab = "1 - Especificdade", col = "steelblue", lwd = 2)
##
## Data: mod2$fitted.values in 262 controls (dados$diabetes 0) < 130 cases (dados$diabetes 1).
## Area under the curve: 0.8631
## 95% CI: 0.8254-0.9007 (DeLong)A Figura 18.19 exibe uma a curva ROC do mod2 onde se observa que AUC = 0,863. Isto significa que a performance do modelo como preditor é muito boa (ver Tabela 18.1).
Verossimilhança no sentido de a qualidade de algo que parece verdadeiro ou provável, que não contraria a verdade↩︎
Existem duas maneiras de descrever uma estimativa de odds: ou como um número isolado, por exemplo, 0,25, subentendendo que expressa uma razão, 0,25: 1,0, ou de forma clara como uma razão 1:4. Ou seja, para cada indivíduo com o fator existem quatro sem o fator. Tradicional e comumente usados no mundo das apostas em corridas de cavalos.↩︎
Foi mantido o nome das variáveis em inglês, pois no banco de dados
oswegoelas estão nessa língua.↩︎em inglês, built-in bias↩︎
Observem que todo o intervalo de confiança de 95% encontra-se abaixo de 1, indicando que existe significância estatística.↩︎
Aproveite para revisar como construir matriz↩︎
Em inglês, NNH (number needed to harm).↩︎
Se houver uma justificativa, esta referência pode ser modificada, usando a função
relevel()do pacotestats. Por exemplo, para mudar a referência da variávelidadeCategpara 31-40, deve-se executar:dados$idadeCateg <- relevel (dados$idadeCateg, ref = "31-40").↩︎https://www.scribbr.com/statistics/akaike-information-criterion/↩︎