Reação Alérgica

Neste exemplo são utilizados os seguintes pacotes: magrittr, knitr, kableExtra, lmtest e vcdExtra. Para a confecção do documento, é necessário o pacote rmarkdown.

Exemplo de aplicação de modelo logito com inserção de variável “dummy”.

Problema: Um antibiótico para pneumonia foi injetado em 100 pacientes com problemas renais e 100 pacientes normais. Alguns pacientes desenvolveram reações alérgicas, sendo 38 dos doentes e 21 dos normais.

freq <- c(38,62,100,21,79,100)
db2 <- matrix(freq,nrow = 2,byrow = T)
dimnames(db2) <- list(c("RENAL","NORMAL"),c("SIM","NÃO","TOTAL"))


kable(db2, booktabs = T) %>% 
  kable_styling(full_width = F) %>%
  add_header_above(c(" " = 1, "Reação alérgica" = 2, " " = 1))

	SIM	NÃO	TOTAL
RENAL	38	62	100
NORMAL	21	79	100

A seguir temos algumas estatísticas calculadas a partir da tabela acima.

CMHtest(db2[,-3]) # Mantel Haenzel Statistics

## Cochran-Mantel-Haenszel Statistics 
## 
##                  AltHypothesis  Chisq Df      Prob
## cor        Nonzero correlation 6.9132  1 0.0085561
## rmeans  Row mean scores differ 6.9132  1 0.0085561
## cmeans  Col mean scores differ 6.9132  1 0.0085561
## general    General association 6.9132  1 0.0085561

fisher.test(db2[,-3], alternative = "two.sided") # Fisher's Exact Test

## 
##  Fisher's Exact Test for Count Data
## 
## data:  db2[, -3]
## p-value = 0.01275
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.178639 4.565483
## sample estimates:
## odds ratio 
##   2.295953

Risco relativo e razão de chances:

Obteve-se as seguintes estimativas com os respectivos intervalos de confiança, considerando um nível de confiança de 95%.

freq <- c(38,62,21,79)
grupo <- rep(c("RENAL","NORMAL"),each = 2)
reage <- rep(c("SIM","NÃO"),2)
db1 <- data.frame(freq,grupo,reage)
db2 <- matrix(freq,nrow = 2,byrow = T)
rownames(db2) <- c("RENAL","NORMAL")
colnames(db2) <- c("SIM","NÃO")

### Função para calcular razão de chances e risco relativo

tab2x2 <- function(dados,conf.level = 0.95){
  # Soma em cada linha
  n1 <- sum(dados[1,])
  n2 <- sum(dados[2,])
  
  # Frequência em cada casela
  n11 <- dados[1,1]
  n12 <- dados[1,2]
  n21 <- dados[2,1]
  n22 <- dados[2,2]
  
  # Proporções
  p11 <- dados[1,1]/n1
  p12 <- dados[1,2]/n1
  p21 <- dados[2,1]/n2
  p22 <- dados[2,1]/n2
  
  # Auxílio IC
  conf.level = conf.level
  cfl <- 1-(1-conf.level)/2
  
  # Risco relativo
  RR <- p11/p21
  
  # intervalo de confiança para o risco relativo
  sigmaRR <- sqrt(((1-p11)/(p11*n1))+((1-p21)/(p21*n2)))
  ICLR <- exp(log(RR) - qnorm(cfl)*sigmaRR)
  ICUR <- exp(log(RR) + qnorm(cfl)*sigmaRR)
  
  ## Razão de chances
  OR <- (n11*n22)/(n12*n21)
  
  ## intervalo de confiança para razão de chances
  sigmaOR <- sqrt((1/n11)+(1/n21)+(1/n12)+(1/n22))
  ICLO <- exp(log(OR) - qnorm(cfl)*sigmaOR)
  ICUO <- exp(log(OR) + qnorm(cfl)*sigmaOR)
  
  ## Organização dos dados em tabela
  mat1 <- matrix(round(c(OR,ICLO,ICUO,RR,ICLR,ICUR),4), nrow = 2, byrow = T)
  rownames(mat1) <- c("Razão de chances","Risco relativo")
  colnames(mat1) <- c("Estimativa", "LI", "LS")
  
  lista <- list('Estimativas' = mat1,'Nível de confiança' = conf.level)
  return(lista)
}

tab2x2(db2)

## $Estimativas
##                  Estimativa     LI     LS
## Razão de chances     2.3057 1.2302 4.3213
## Risco relativo       1.8095 1.1478 2.8526
## 
## $`Nível de confiança`
## [1] 0.95

Preparação dos dados para aplicação dos modelos de regressão logística:

GRUPO <- c(rep("RENAL",100),rep("NORMAL",100))
#REAGE <- c(rep("SIM",38),rep("NÃO",62),rep("SIM",21),rep("NÃO",79))
REAGE <- c(rep(1,38),rep(0,62),rep(1,21),rep(0,79))
DUMMY <- c(rep(1,100),rep(0,100))
CODE <- c(rep(2,100),rep(1,100))
DEV <- c(rep(-1,100),rep(1,100))

db3 <- data.frame(GRUPO,REAGE,DUMMY,CODE,DEV)

db3

Codificações utilizadas para as variáveis:

GRUPO	DUMMY	CODE	DEV
RENAL	1	2	-1
NORMAL	0	1	1

A variável resposta, sob a denominação REAGE, assume o valor 1 para resposta SIM e 0 para NÃO.

1) Modelo com Codificação Deviation from the mean

Ajuste do modelo:

m3 <- glm(REAGE ~ DEV, family = binomial, data = db3)
summary(m3)

## 
## Call:
## glm(formula = REAGE ~ DEV, family = binomial, data = db3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9778  -0.9778  -0.6866   1.3911   1.7667  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9072     0.1603  -5.661  1.5e-08 ***
## DEV          -0.4177     0.1603  -2.606  0.00915 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 242.63  on 199  degrees of freedom
## Residual deviance: 235.60  on 198  degrees of freedom
## AIC: 239.6
## 
## Number of Fisher Scoring iterations: 4

-2logLik:

# -2logLik
-2*logLik(m3)

## 'log Lik.' 235.6042 (df=2)

Testes de hipóteses para \(\beta\) = 0:

# Razão de verossimilhanças
lmtest::lrtest(m3)

# Teste de Wald
lmtest::waldtest(m3)

# Score
anova(m3,test = "Rao")

Razão de chances (OR) e intervalo de confiança:

O parâmetro \(\beta\) representa o log da razão de chances, logo, para obter a estimativa da OR, deve-se exponenciar \(\beta\).

# estimativa de Beta
beta3 <- m3$coefficients[2]

# ASE
sigma3 <- summary(m3)[[12]][2,2]

# Razão de chances
OR3 <- round(exp(2*beta3),3)

# Intervalo de confiança
LI3 <- round(exp(2*beta3 - 2*1.96*sigma3),3)
LS3 <- round(exp(2*beta3 + 2*1.96*sigma3),3)

or3 <- data.frame(Estimativa = OR3, LI = LI3, LS = LS3)
row.names(or3) <- "OR"

or3

2) Modelo com codificação X = 1 para RENAL e X = 0 para NORMAL

O modelo será estimado utilizando a codificação “DUMMY”, X = 1 para pacientes com problemas renais (RENAL) e X = 0 para pacientes normais (NORMAL).

m1 <- glm(REAGE ~ DUMMY, family = binomial, data = db3)
summary(m1)

## 
## Call:
## glm(formula = REAGE ~ DUMMY, family = binomial, data = db3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9778  -0.9778  -0.6866   1.3911   1.7667  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3249     0.2455  -5.397 6.79e-08 ***
## DUMMY         0.8354     0.3205   2.606  0.00915 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 242.63  on 199  degrees of freedom
## Residual deviance: 235.60  on 198  degrees of freedom
## AIC: 239.6
## 
## Number of Fisher Scoring iterations: 4

-2logLik:

# -2logLik
-2*logLik(m1)

## 'log Lik.' 235.6042 (df=2)

Testes de hipóteses para \(\beta\) = 0:

# Razão de verossimilhanças
lmtest::lrtest(m1)

# Teste de Wald
lmtest::waldtest(m1)

# Score
anova(m1,test = "Rao")

Razão de chances (OR) e intervalo de confiança:

O parâmetro \(\beta\) representa o log da razão de chances, logo, para obter a estimativa da OR, deve-se exponenciar \(\beta\).

# estimativa de Beta
beta1 <- m1$coefficients[2]

# ASE
sigma1 <- summary(m1)[[12]][2,2]

# Razão de chances
OR1 <- round(exp(beta1),3)

# Intervalo de confiança
LI1 <- round(exp(beta1 - 1.96*sigma1),3)
LS1 <- round(exp(beta1 + 1.96*sigma1),3)

or1 <- data.frame(Estimativa = OR1, LI = LI1, LS = LS1)
row.names(or1) <- "OR"

or1

Verifica-se que o intercepto, \(\hat{\alpha}\), é o log da chance para X = 0 (Grupo de referência);

\(\hat{\beta}\) é o log da razão de chances para X = 1 relativo a X = 0.

3) Modelo com codificação X = 2 para RENAL e X = 1 para NORMAL

O modelo seguinte considera uma codificação alternativa para as variáveis, X = 2 para RENAL e X = 1 para NORMAL.

Ajuste do modelo:

m2 <- glm(REAGE ~ CODE, family = binomial, data = db3)
summary(m2)

## 
## Call:
## glm(formula = REAGE ~ CODE, family = binomial, data = db3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9778  -0.9778  -0.6866   1.3911   1.7667  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1603     0.5325  -4.057 4.97e-05 ***
## CODE          0.8354     0.3205   2.606  0.00915 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 242.63  on 199  degrees of freedom
## Residual deviance: 235.60  on 198  degrees of freedom
## AIC: 239.6
## 
## Number of Fisher Scoring iterations: 4

-2logLik:

# -2logLik
-2*logLik(m2)

## 'log Lik.' 235.6042 (df=2)

Testes de hipóteses para \(\beta\) = 0:

# Razão de verossimilhanças
lmtest::lrtest(m2)

# Teste de Wald
lmtest::waldtest(m2)

# Score
anova(m2,test = "Rao")

Razão de chances (OR) e intervalo de confiança:

# estimativa de Beta
beta2 <- m2$coefficients[2]

# ASE
sigma2 <- summary(m2)[[12]][2,2]

# Razão de chances
OR2 <- round(exp(beta2),3)

# Intervalo de confiança
LI2 <- round(exp(beta2 - 1.96*sigma2),3)
LS2 <- round(exp(beta2 + 1.96*sigma2),3)

or2 <- data.frame(Estimativa = OR2, LI = LI2, LS = LS2)
row.names(or2) <- "OR"

or2

Em relação à codificação anterior, as novas estimativas, \(\hat\alpha^{*}\) e \(\hat\beta^{*}\), possuem a seguinte relação com as estimativas anteriores:

\(\hat\beta^{*}\) = \(\hat\beta\)
\(\hat\alpha^{*}\) = \(\hat\alpha\) - \(\hat{\beta}\)

Referências Bibliográficas

AGRESTI, ALAN. An Introduction to Categorical Data Analysis. John Wiley & Sons, second edition, 2007.

NOTAS DE AULA. Análise de Dados Categorizados. Curso de Graduação em Estatística, UnB, 2018.

R CORE TEAM. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Reação Alérgica

Regressão Logística - Análise de Dados Categorizados - Prova 3

Allan Vieira - 14/0128492

Caio Balena - 13/0104566

Frederico de Lucca - 15/0072139

José Cezário - 15/0072414

Junho de 2018

Exemplo de aplicação de modelo logito com inserção de variável “dummy”.

1) Modelo com Codificação Deviation from the mean

2) Modelo com codificação X = 1 para RENAL e X = 0 para NORMAL

3) Modelo com codificação X = 2 para RENAL e X = 1 para NORMAL

Referências Bibliográficas