Ejemplo de análisis de un estudio de casos y controles

Investigando una epidemia: Oswego, NY

El 19 de abril de 1940, el funcionario de salud local de la aldea de Lycoming, condado de Oswego, Nueva York, informó al funcionario de salud del distrito de Siracusa de la aparición de un brote de enfermedad gastrointestinal aguda.  El Dr. A. M. Rubin, epidemiólogo en formación, fue asignado para llevar a cabo una investigación.

Cuando el Dr. Rubin llegó al campo, se enteró por el oficial de salud que todas las personas que se sabía que estaban enfermas habían asistido a una cena de la iglesia la noche anterior, el 18 de abril.  Los familiares que no habían asistido a la cena de la iglesia no se habían enfermado.  Por consiguiente, la investigación se centró en las circunstancias relacionadas con la cena.

Las entrevistas sobre la presencia de los síntomas, incluyendo el día y la hora de inicio, y la comida consumida en la cena de la iglesia, se completaron en 75 de las 80 personas que se sabía que estaban presentes.  Se identificó un total de 46 personas que habían experimentado enfermedades gastrointestinales.

Preguntas:

• ¿Esto es una epidemia? De ser así. ¿cuál es la causa?
• ¿Endémica para la región?
  
• ¿Debido a la variación estacional?
  
• ¿Debido a la variación aleatoria?

Ahora vamos a investigar esto mediante el uso de R y un paquete que resumen algunas medidas epidemiológicas de utilidad. El objetivo es identificar algún alimento que pueda ser asociado con el brote de la enfermedad.

Primero instalo los paquetes necesarios # Paquetes

# install.packages("epiDisplay") # si no están instalados
library(epiDisplay) # resumen de medidas estadísticas

Loading required package: foreign
Loading required package: survival
Loading required package: MASS
Loading required package: nnet

library(tidyverse) # para manejo de datos

── Attaching packages ────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.1          ✔ purrr   0.3.2     
✔ tibble  2.1.3          ✔ dplyr   0.8.1     
✔ tidyr   0.8.3.9000     ✔ stringr 1.4.0     
✔ readr   1.3.1          ✔ forcats 0.3.0     
── Conflicts ───────────────────────────────── tidyverse_conflicts() ──
✖ ggplot2::alpha() masks epiDisplay::alpha()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ dplyr::select()  masks MASS::select()

library(ggpubr) # para graficos listos para publicación

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

Dataset

data(Oswego)

Veo los datos

glimpse(Oswego)

Observations: 75
Variables: 20
$ age        <dbl> 11, 52, 65, 59, 13,…
$ sex        <I<chr>> M, F, M, F, F, F…
$ timesupper <dbl> NA, 2000, 1830, 183…
$ ill        <lgl> FALSE, TRUE, TRUE, …
$ onsetdate  <I<chr>> NA, 04/19, 04/19…
$ onsettime  <dbl> NA, 30, 30, 30, NA,…
$ bakedham   <lgl> FALSE, TRUE, TRUE, …
$ spinach    <lgl> FALSE, TRUE, TRUE, …
$ mashedpota <lgl> FALSE, TRUE, TRUE, …
$ cabbagesal <lgl> FALSE, FALSE, TRUE,…
$ jello      <lgl> FALSE, FALSE, FALSE…
$ rolls      <lgl> FALSE, TRUE, FALSE,…
$ brownbread <lgl> FALSE, FALSE, FALSE…
$ milk       <lgl> FALSE, FALSE, FALSE…
$ coffee     <lgl> FALSE, TRUE, TRUE, …
$ water      <lgl> FALSE, FALSE, FALSE…
$ cakes      <lgl> FALSE, FALSE, FALSE…
$ vanilla    <lgl> FALSE, TRUE, TRUE, …
$ chocolate  <lgl> TRUE, FALSE, TRUE, …
$ fruitsalad <lgl> FALSE, FALSE, FALSE…

Hay 20 columnas y 75 observaciones. Esto lo puedo ver también con

dim(Oswego)

[1] 75 20

Todas las instrucciones en R van como filas, columnas

Examino el dataset

knitr::kable(head(Oswego), digits = 2, align = c(rep("l", 4), rep("c", 4), rep("r", 4)))

age	sex	timesupper	ill	onsetdate	onsettime	bakedham	spinach	mashedpota	cabbagesal	jello	rolls	brownbread	milk	coffee	water	cakes	vanilla	chocolate	fruitsalad
11	M	NA	FALSE	NA	NA	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE
52	F	2000	TRUE	04/19	30	TRUE	TRUE	TRUE	FALSE	FALSE	TRUE	FALSE	FALSE	TRUE	FALSE	FALSE	TRUE	FALSE	FALSE
65	M	1830	TRUE	04/19	30	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE	FALSE	TRUE	TRUE	FALSE
59	F	1830	TRUE	04/19	30	TRUE	TRUE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE	TRUE	TRUE	TRUE	FALSE
13	F	NA	FALSE	NA	NA	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE
63	F	1930	TRUE	04/18	2230	TRUE	TRUE	FALSE	TRUE	TRUE	FALSE	FALSE	FALSE	FALSE	TRUE	FALSE	TRUE	FALSE	FALSE

Análisis exploratorio de datos o EDA

Primero veo que datos tengo. Puedo ver los nombres de las columnas

names(Oswego)

 [1] "age"        "sex"       
 [3] "timesupper" "ill"       
 [5] "onsetdate"  "onsettime" 
 [7] "bakedham"   "spinach"   
 [9] "mashedpota" "cabbagesal"
[11] "jello"      "rolls"     
[13] "brownbread" "milk"      
[15] "coffee"     "water"     
[17] "cakes"      "vanilla"   
[19] "chocolate"  "fruitsalad"

o la estructura de los datos

str(Oswego)

'data.frame':   75 obs. of  20 variables:
 $ age       : num  11 52 65 59 13 63 70 40 15 33 ...
 $ sex       : 'AsIs' chr  "M" "F" "M" "F" ...
 $ timesupper: num  NA 2000 1830 1830 NA 1930 1930 1930 2200 1900 ...
 $ ill       : logi  FALSE TRUE TRUE TRUE FALSE TRUE ...
 $ onsetdate : 'AsIs' chr  NA "04/19" "04/19" "04/19" ...
 $ onsettime : num  NA 30 30 30 NA 2230 2230 200 100 2300 ...
 $ bakedham  : logi  FALSE TRUE TRUE TRUE FALSE TRUE ...
 $ spinach   : logi  FALSE TRUE TRUE TRUE FALSE TRUE ...
 $ mashedpota: logi  FALSE TRUE TRUE FALSE FALSE FALSE ...
 $ cabbagesal: logi  FALSE FALSE TRUE FALSE FALSE TRUE ...
 $ jello     : logi  FALSE FALSE FALSE FALSE FALSE TRUE ...
 $ rolls     : logi  FALSE TRUE FALSE FALSE FALSE FALSE ...
 $ brownbread: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ milk      : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 $ coffee    : logi  FALSE TRUE TRUE TRUE FALSE FALSE ...
 $ water     : logi  FALSE FALSE FALSE FALSE FALSE TRUE ...
 $ cakes     : logi  FALSE FALSE FALSE TRUE FALSE FALSE ...
 $ vanilla   : logi  FALSE TRUE TRUE TRUE FALSE TRUE ...
 $ chocolate : logi  TRUE FALSE TRUE TRUE TRUE FALSE ...
 $ fruitsalad: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
 - attr(*, "prompts")= Named chr  "Age:" "Sex:" "Time of Supper(24 hour):" "Ill?" ...
  ..- attr(*, "names")= chr  "AGE" "SEX" "TIMESUPPER" "ILL" ...

o un resumen de los datos

Resumen de los datos

summary(Oswego)

      age            sex           
 Min.   : 3.00   Length:75         
 1st Qu.:16.50   Class :AsIs       
 Median :36.00   Mode  :character  
 Mean   :36.81                     
 3rd Qu.:57.50                     
 Max.   :77.00                     
                                   
   timesupper      ill         
 Min.   :   0   Mode :logical  
 1st Qu.:1930   FALSE:29       
 Median :1930   TRUE :46       
 Mean   :1931                  
 3rd Qu.:2200                  
 Max.   :2200                  
 NA's   :47                    
  onsetdate           onsettime   
 Length:75          Min.   :  30  
 Class :AsIs        1st Qu.: 100  
 Mode  :character   Median : 630  
                    Mean   :1132  
                    3rd Qu.:2230  
                    Max.   :2400  
                    NA's   :29    
  bakedham        spinach       
 Mode :logical   Mode :logical  
 FALSE:29        FALSE:32       
 TRUE :46        TRUE :43       
                                
                                
                                
                                
 mashedpota      cabbagesal     
 Mode :logical   Mode :logical  
 FALSE:37        FALSE:47       
 TRUE :37        TRUE :28       
 NA's :1                        
                                
                                
                                
   jello           rolls        
 Mode :logical   Mode :logical  
 FALSE:52        FALSE:38       
 TRUE :23        TRUE :37       
                                
                                
                                
                                
 brownbread         milk        
 Mode :logical   Mode :logical  
 FALSE:48        FALSE:71       
 TRUE :27        TRUE :4        
                                
                                
                                
                                
   coffee          water        
 Mode :logical   Mode :logical  
 FALSE:44        FALSE:51       
 TRUE :31        TRUE :24       
                                
                                
                                
                                
   cakes          vanilla       
 Mode :logical   Mode :logical  
 FALSE:35        FALSE:21       
 TRUE :40        TRUE :54       
                                
                                
                                
                                
 chocolate       fruitsalad     
 Mode :logical   Mode :logical  
 FALSE:27        FALSE:69       
 TRUE :47        TRUE :6        
 NA's :1                        
                                
                                
                                

Resumen con skim

Otra forma de hacer un resumen de lso datos es con el paquete skim

skimr::skim(Oswego)

Skim summary statistics
 n obs: 75 
 n variables: 20 

── Variable type:AsIs ───────────────────
  variable missing complete  n n_unique
 onsetdate      29       46 75        2
       sex       0       75 75        2
 min_length max_length
          1          1
          1          1

── Variable type:logical ────────────────
   variable missing complete  n  mean
   bakedham       0       75 75 0.61 
 brownbread       0       75 75 0.36 
 cabbagesal       0       75 75 0.37 
      cakes       0       75 75 0.53 
  chocolate       1       74 75 0.64 
     coffee       0       75 75 0.41 
 fruitsalad       0       75 75 0.08 
        ill       0       75 75 0.61 
      jello       0       75 75 0.31 
 mashedpota       1       74 75 0.5  
       milk       0       75 75 0.053
      rolls       0       75 75 0.49 
    spinach       0       75 75 0.57 
    vanilla       0       75 75 0.72 
      water       0       75 "print" 0.32 
                   count
 TRU: 46, FAL: 29, NA: 0
 FAL: 48, TRU: 27, NA: 0
 FAL: 47, TRU: 28, NA: 0
 TRU: 40, FAL: 35, NA: 0
 TRU: 47, FAL: 27, NA: 1
 FAL: 44, TRU: 31, NA: 0
  FAL: 69, TRU: 6, NA: 0
 TRU: 46, FAL: 29, NA: 0
 FAL: 52, TRU: 23, NA: 0
 FAL: 37, TRU: 37, NA: 1
  FAL: 71, TRU: 4, NA: 0
 FAL: 38, TRU: 37, NA: 0
 TRU: 43, FAL: 32, NA: 0
 TRU: 54, FAL: 21, NA: 0
 FAL: 51, TRU: 24, NA: 0

── Variable type:numeric ────────────────
   variable missing complete  n    mean
        age       0       75 75   36.81
  onsettime      29       46 75 1132.07
 timesupper      47       28 75 1930.71
      sd p0    p25  p50    p75 p100
   21.45  3   16.5   36   57.5   77
 1057.9  30  100    630 2230   2400
  440.64  0 1930   1930 2200   2200
     hist
 ▆▇▂▇▂▅▆▃
 ▇▁▁▁▁▁▁▇
 ▁▁▁▁▁▁▂▇

Resumen epidemiológico

Es conveniente tener a la vista un resumen de datos

epiDisplay::summ(Oswego)

No. of observations = 75

   Var. name  obs. mean    median 
1  age        75   36.81   36     
2  sex                            
3  timesupper 28   1930.71 1930   
4  ill        75   0.61    1      
5  onsetdate                      
6  onsettime  46   1132.07 630    
7  bakedham   75   0.61    1      
8  spinach    75   0.57    1      
9  mashedpota 74   0.5     0.5    
10 cabbagesal 75   0.37    0      
11 jello      75   0.31    0      
12 rolls      75   0.49    0      
13 brownbread 75   0.36    0      
14 milk       75   0.05    0      
15 coffee     75   0.41    0      
16 water      75   0.32    0      
17 cakes      75   0.53    1      
18 vanilla    75   0.72    1      
19 chocolate  74   0.64    1      
20 fruitsalad 75   0.08    0      
   s.d.   min.   max.  
1  21.45  3      77    
2                      
3  440.64 0      2200  
4  0.49   0      1     
5                      
6  1057.9 30     2400  
7  0.49   0      1     
8  0.5    0      1     
9  0.5    0      1     
10 0.49   0      1     
11 0.46   0      1     
12 0.5    0      1     
13 0.48   0      1     
14 0.23   0      1     
15 0.5    0      1     
16 0.47   0      1     
17 0.5    0      1     
18 0.45   0      1     
19 0.48   0      1     
20 0.27   0      1     

Visualización de datos faltantes

Y ver los datos faltantes con visdat

visdat::vis_dat(Oswego)

Resúmenes por variable

Veo la distribución de la edad

Puedo ver una piramide con la distribución de la edad por sexo

epiDisplay::pyramid(
  Oswego$age,
  Oswego$sex,
  binwidth = 10,
  percent = "each",
  decimal = 2,
  col.gender = c("pink", "lightblue")
  )

o con un histograma dividido por edad

Oswego %>%
  ggplot(aes(x = age, fill = sex)) +
  geom_histogram(bins = 10) +
  facet_grid(sex ~ .) + 
  labs(title = "Distribución de la edad por sexo", 
       x = "Edad", y = "n") + 
  ggpubr::theme_pubclean()

Error in grDevices::col2rgb(colour, TRUE) : invalid color name 'F'

Puedo comenzar a examinar algunas distribuciones, por ejemplo, edad y enfermo

Exploro la edad

Oswego$sex <- as.factor(Oswego$sex)

Hay diferencia en la edad? Ejemplo general para comparar diferencias de una variabla continua entre dos grupos

statStack (
  cont.var = age,
  by = sex,
  dataFrame = Oswego,
  iqr = "auto",
  var.labels = TRUE,
  decimal = 1,
  assumption.p.value = .01
  )

  Total  age   median(IQR)     Test         P value
  sex                          Ranksum test 0.987  
F 44         35.5 (16.8,53.2)                      
M 31         36 (16,61)                            
                                                   

Ok, no hay diferencia de edad por sexo Veo ahora si hay distinta edad por enfermedad

epiDisplay::statStack (
  cont.var = age,
  by = ill,
  dataFrame = Oswego,
  iqr = "auto",
  var.labels = TRUE,
  decimal = 1,
  assumption.p.value = .01
  )

      Total  age   median(IQR)     Test         P value
      ill                          Ranksum test 0.284  
FALSE 29         35 (14,50)                            
TRUE  46         38.5 (17.2,58.8)                      
                                                       

Tampoco. Lo interesante es que el valor por defecto para ’iqr’ es “auto”, lo que significa que hace una comparación asumiendo una distribución no paramétrica y el test elige automáticamente si el valor de p para el test de Bartlett está bajo el ’assumption.p.value’ argument (0.01).

Esto también lo podría hacer mediante un t-test

t.test(age~ill, data = Oswego)

    Welch Two Sample t-test

data:  age by ill
t = -1.2663, df = 62.329, p-value = 0.2101
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -16.320837   3.661167
sample estimates:
mean in group FALSE  mean in group TRUE 
           32.93103            39.26087 

pero lo más importante es verlo

Este gráfico es esencialmente el mismo que este

Puedo además verlo por sexo

Intoxicación

Ahora vamos a ver que podría haber causado la intoxicación

Veo cuantos enfermos hay

table(Oswego$ill)

FALSE  TRUE 
   29    46 

Veo los enfermos por sexo

table(Oswego$ill, Oswego$sex)

       
         F  M
  FALSE 14 15
  TRUE  30 16

Siempre es una buena idea comenzar a ver los datos mediante gráficos, tratando de entender la historia que quieren contar Por ejemplo: se enferman igual los hombres y mujeres? Grafico. Como son datos categóricos, el gráfico de elecciónes mosaicplot, que es un tipo de gráfico de barras apilado. La orden shade marca diferencias significativas

Esto lo confirmo mediante un chi2

chisq.test(table(Oswego$ill, Oswego$sex))

    Pearson's Chi-squared test with Yates' continuity correction

data:  table(Oswego$ill, Oswego$sex)
X-squared = 1.4646, df = 1, p-value = 0.2262

A continuación, comienzo a examinar los alimentos

Otra manera, mediante el paquete epiDisplay

epiDisplay::tabpct(Oswego$ill, Oswego$spinach)

Original table 
          Oswego$spinach
Oswego$ill  FALSE  TRUE  Total
     FALSE     12    17     29
     TRUE      20    26     46
     Total     32    43     75

Row percent 
          Oswego$spinach
Oswego$ill   FALSE    TRUE  Total
     FALSE      12      17     29
            (41.4)  (58.6)  (100)
     TRUE       20      26     46
            (43.5)  (56.5)  (100)

Column percent 
          Oswego$spinach
Oswego$ill  FALSE       %  TRUE       %
     FALSE     12  (37.5)    17  (39.5)
     TRUE      20  (62.5)    26  (60.5)
     Total     32   (100)    43   (100)

otro alimento

Ajá! al parecer la vainilla está asociada de alguna manera con la intoxicación. Para aislar la exposición de la vainilla, hago una regresión

Riesgo

Luego de familiarizarme con todos y cada uno de los datos, comienzo a hacer análisis multivariados. Los términos multivariado y multivariable a menudo se utilizan indistintamente en la literatura de salud pública. Sin embargo, estos términos representan en realidad dos tipos de análisis muy distintos

La mayoría de los modelos de regresión se describen en términos de la forma en que se modela la variable de resultado: en la regresión lineal el resultado es continuo, la regresión logística tiene un resultado dicotómico, y el análisis de supervivencia implica un resultado de tiempo hasta el evento. Estadísticamente hablando, el análisis multivariado se refiere a modelos estadísticos que tienen 2 o más variables dependientes o de resultado, y el análisis multivariable se refiere a modelos estadísticos en los que hay múltiples variables independientes o de respuesta.

Un modelo multivariable puede ser considerado como un modelo en el que múltiples variables se encuentran en el lado derecho de la ecuación del modelo. Este tipo de modelo estadístico se puede utilizar para intentar evaluar la relación entre una serie de variables; se pueden evaluar relaciones independientes al tiempo que se ajustan a los posibles factores de confusión.

En este caso, la respuesta es una sola, si enferma o no, y las variables independientes son varias, por eso es un análisis multivariable.

El objetivo del análisis es identificar algúna variable independiente (alimento) asociada a la variable dependiente (enfermedad) y lograr aislar el efecto de esa variable, esto es, si todo lo demás está igual, es la presencia de esta variable un factor asociado a la enfermedad?

La medida de resultados de un estudio de casos y controles es el odds ratio que se interpretra como el riesgo asociado a una exposición. Por ejemplo, un OR de 20 para fumar y cáncer al pulmón quiere decir que una persona que fuma tiene 20 veces más riesgo de tener cáncer al pulmón comparado a una persona que no fuma. Nota que esto no quiere decir que una persona que fuma tendrá cáncer, o que una persona que no fuma estará sana. Solo quiere decir que hay un mayor riesgo asociado al fumar para presentar cáncer al pulmón.

Ahora vamos a hacer una regresión, esto es, ver de que manera cada variable se asocia de manera individual a la respyesta y luego vamos a calcular los OR para cada uno de los alimentos.

model1 <- glm(
  ill ~ bakedham +  spinach + mashedpota + cabbagesal + jello + rolls + brownbread +
  milk + coffee + water + cakes + vanilla + chocolate + fruitsalad,
  data = Oswego,
  family = "binomial"
  )

Veo el modelo

model1

Call:  glm(formula = ill ~ bakedham + spinach + mashedpota + cabbagesal + 
    jello + rolls + brownbread + milk + coffee + water + cakes + 
    vanilla + chocolate + fruitsalad, family = "binomial", data = Oswego)

Coefficients:
   (Intercept)    bakedhamTRUE     spinachTRUE  mashedpotaTRUE  cabbagesalTRUE       jelloTRUE       rollsTRUE  
       -2.1486          0.9605         -0.4917         -0.3472          0.1143         -0.3085         -0.2051  
brownbreadTRUE        milkTRUE      coffeeTRUE       waterTRUE       cakesTRUE     vanillaTRUE   chocolateTRUE  
        0.7349         -1.9674         -0.5125          0.1614          0.8512          3.2256         -0.1184  
fruitsaladTRUE  
       -0.3998  

Degrees of Freedom: 72 Total (i.e. Null);  58 Residual
  (2 observations deleted due to missingness)
Null Deviance:      97.2 
Residual Deviance: 66.36    AIC: 96.36

Los coeficientes sugieren que efectivamente la vainilla está asociada

summary(model1)

Call:
glm(formula = ill ~ bakedham + spinach + mashedpota + cabbagesal + 
    jello + rolls + brownbread + milk + coffee + water + cakes + 
    vanilla + chocolate + fruitsalad, family = "binomial", data = Oswego)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0451  -0.5118   0.4803   0.6850   1.8073  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.1486     1.0818  -1.986 0.047022 *  
bakedhamTRUE     0.9605     1.4919   0.644 0.519715    
spinachTRUE     -0.4917     1.3119  -0.375 0.707801    
mashedpotaTRUE  -0.3472     0.9166  -0.379 0.704872    
cabbagesalTRUE   0.1143     0.8956   0.128 0.898413    
jelloTRUE       -0.3085     1.0255  -0.301 0.763525    
rollsTRUE       -0.2051     1.1887  -0.173 0.863031    
brownbreadTRUE   0.7349     0.9627   0.763 0.445254    
milkTRUE        -1.9674     1.6379  -1.201 0.229681    
coffeeTRUE      -0.5125     1.1531  -0.444 0.656705    
waterTRUE        0.1614     0.9552   0.169 0.865835    
cakesTRUE        0.8512     0.7170   1.187 0.235159    
vanillaTRUE      3.2256     0.8517   3.787 0.000152 ***
chocolateTRUE   -0.1184     0.7407  -0.160 0.872973    
fruitsaladTRUE  -0.3998     1.1679  -0.342 0.732115    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 97.204  on 72  degrees of freedom
Residual deviance: 66.365  on 58  degrees of freedom
  (2 observations deleted due to missingness)
AIC: 96.365

Number of Fisher Scoring iterations: 4

la diferencia es significativa.

Odds ratio

Ahora voy a ver los odds ratio, los que los obtengo del exponente de los coeficientes

(coefficients(model1))

   (Intercept)   bakedhamTRUE    spinachTRUE mashedpotaTRUE cabbagesalTRUE      jelloTRUE      rollsTRUE brownbreadTRUE 
    -2.1486370      0.9604614     -0.4917336     -0.3471575      0.1143420     -0.3085216     -0.2050688      0.7349155 
      milkTRUE     coffeeTRUE      waterTRUE      cakesTRUE    vanillaTRUE  chocolateTRUE fruitsaladTRUE 
    -1.9673923     -0.5125336      0.1613828      0.8511587      3.2255750     -0.1184222     -0.3997926 

El problema es que al trabajar con muestras, el indicador puntual, el OR, puede ser engañoso. Por ejemplo, alguien podría preguntar por la profundidad promedio de una piscina y el hecho que tenga 1m promedio de profundidad no es un dato muy útil. Nos interesa saber dentro de qué rango está la profundida. Así, por ejemplo es más util que nos diga, que tiene 30cm en la parte más baja y 3m en la parte más profunda.

Ahora bien, el intervalo de confianza nos indica un rango dentro del cual, con una probabilidad conocida, se encuentra el valor que estamos buscando. Usualmente se utiliza el IC 95%, lo que quiere decir que con un 95% de confianza, el valor real a estimar está dentro de ese rango.

Así que para ver ese rango, del 2.5 al 97.5% del OR calculo

(confint(model1))

Waiting for profiling to be done...

                    2.5 %     97.5 %
(Intercept)    -4.4718504 -0.1456918
bakedhamTRUE   -1.9287861  4.0615325
spinachTRUE    -3.3418368  1.8972819
mashedpotaTRUE -2.2293234  1.4365270
cabbagesalTRUE -1.6822460  1.8930439
jelloTRUE      -2.3763566  1.7450059
rollsTRUE      -2.7455423  2.0866226
brownbreadTRUE -1.1094364  2.7515784
milkTRUE       -5.5774125  1.2752807
coffeeTRUE     -2.9114795  1.7274071
waterTRUE      -1.7398567  2.0996029
cakesTRUE      -0.5290465  2.3317015
vanillaTRUE     1.7085300  5.1209226
chocolateTRUE  -1.6011794  1.3505816
fruitsaladTRUE -2.7190521  2.0817094

Super!

Tabla odds ratio

Ahora los puedo unir para crear una tabla:

exp(cbind(coef(model1), confint(model1))) 

Waiting for profiling to be done...

                               2.5 %      97.5 %
(Intercept)     0.1166430 0.01142615   0.8644241
bakedhamTRUE    2.6129019 0.14532450  58.0632233
spinachTRUE     0.6115652 0.03537193   6.6677464
mashedpotaTRUE  0.7066940 0.10760120   4.2060629
cabbagesalTRUE  1.1211355 0.18595584   6.6395480
jelloTRUE       0.7345321 0.09288839   5.7259352
rollsTRUE       0.8145912 0.06421347   8.0576556
brownbreadTRUE  2.0853058 0.32974474  15.6673412
milkTRUE        0.1398210 0.00378234   3.5797062
coffeeTRUE      0.5989761 0.05439519   5.6260470
waterTRUE       1.1751348 0.17554556   8.1629277
cakesTRUE       2.3423595 0.58916649  10.2954442
vanillaTRUE    25.1680412 5.52083967 167.4898205
chocolateTRUE   0.8883209 0.20165854   3.8596697
fruitsaladTRUE  0.6704591 0.06593722   8.0181636

Ahora busco aquel donde el intervalo menor sea mayor a 1 Veo que el único alimento que tiene esa condición es la vainilla, que en el mejor de los casos expone a un paciente a un riesgo 5 veces mayor de tener la enfermedad, independientemente de que haya comido alguna otra cosa.

Modelo final de riesgo

logistic.display(model1, simplified = F)

Logistic regression predicting ill : TRUE vs FALSE 
 
           crude OR(95%CI)     adj. OR(95%CI)      P(Wald's test) P(LR-test)
bakedham   1.24 (0.47,3.23)    2.61 (0.14,48.64)   0.52           0.516     
                                                                            
spinach    0.94 (0.36,2.43)    0.61 (0.05,8)       0.708          0.702     
                                                                            
mashedpota 0.96 (0.37,2.46)    0.71 (0.12,4.26)    0.705          0.704     
                                                                            
cabbagesal 1.41 (0.52,3.8)     1.12 (0.19,6.49)    0.898          0.898     
                                                                            
jello      1.66 (0.58,4.73)    0.73 (0.1,5.48)     0.764          0.763     
                                                                            
rolls      0.76 (0.29,1.95)    0.81 (0.08,8.37)    0.863          0.862     
                                                                            
brownbread 1.28 (0.47,3.47)    2.09 (0.32,13.76)   0.445          0.439     
                                                                            
milk       0.3 (0.03,3.42)     0.14 (0.01,3.47)    0.23           0.227     
                                                                            
coffee     0.89 (0.34,2.31)    0.6 (0.06,5.74)     0.657          0.654     
                                                                            
water      0.73 (0.27,2)       1.18 (0.18,7.64)    0.866          0.866     
                                                                            
cakes      1.73 (0.67,4.49)    2.34 (0.57,9.55)    0.235          0.228     
                                                                            
vanilla    21.64 (5.36,87.34)  25.17 (4.74,133.6)  < 0.001        < 0.001   
                                                                            
chocolate  0.42 (0.15,1.18)    0.89 (0.21,3.79)    0.873          0.873     
                                                                            
fruitsalad 0.93 (0.15,5.93)    0.67 (0.07,6.61)    0.732          0.734     
                                                                            
Log-likelihood = -33.1825
No. of observations = 73
AIC value = 96.3649

Y puedo hacer una linda tabla con el paquete stargazer

stargazer::stargazer(
  model1,
  title = "Intoxicación en Oswego",
  align = TRUE,
  type = "text",
  ci = TRUE,
  ci.level = 0.95,
  single.row = TRUE,
  apply.coef = exp, # coloca los OR
  apply.se   = exp # coloca los OR
  )

length of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changedlength of NULL cannot be changed

Intoxicación en Oswego
=============================================
                      Dependent variable:    
                  ---------------------------
                              ill            
---------------------------------------------
bakedham            2.613 (-6.100, 11.326)   
spinach              0.612 (-6.667, 7.890)   
mashedpota           0.707 (-4.195, 5.608)   
cabbagesal           1.121 (-3.679, 5.921)   
jello                0.735 (-4.731, 6.200)   
rolls                0.815 (-5.619, 7.249)   
brownbread           2.085 (-3.048, 7.218)   
milk                0.140 (-9.943, 10.222)   
coffee               0.599 (-5.610, 6.808)   
water                1.175 (-3.919, 6.269)   
cakes                2.342 (-1.672, 6.357)   
vanilla           25.168*** (20.575, 29.761) 
chocolate            0.888 (-3.222, 4.999)   
fruitsalad           0.670 (-5.631, 6.972)   
Constant             0.117 (-5.665, 5.899)   
---------------------------------------------
Observations                  73             
Log Likelihood              -33.182          
Akaike Inf. Crit.           96.365           
=============================================
Note:             *p<0.1; **p<0.05; ***p<0.01

Otros gráficos

Puedo además graficar el OR para cada alimento en particular

epiDisplay::cc(Oswego$ill, Oswego$chocolate, design="case-control", main = "OR para el chocolate")

          Oswego$chocolate
Oswego$ill FALSE TRUE Total
     FALSE     7   22    29
     TRUE     20   25    45
     Total    27   47    74

OR =  0.4 
95% CI =  0.14, 1.12  
Chi-squared = 3.14, 1 d.f., P value = 0.076
Fisher's exact test (2-sided) P value = 0.089 

epiDisplay::cc(Oswego$ill, Oswego$vanilla, design="case-control", main = "OR para el chocolate")

          Oswego$vanilla
Oswego$ill FALSE TRUE Total
     FALSE    18   11    29
     TRUE      3   43    46
     Total    21   54    75

OR =  23.45 
95% CI =  5.84, 94.18  
Chi-squared = 27.22, 1 d.f., P value = 0
Fisher's exact test (2-sided) P value = 0 

epiDisplay::cs(Oswego$ill, Oswego$chocolate)

          Exposure
Outcome    Non-exposed Exposed Total
  Negative 7           22      29   
  Positive 20          25      45   
  Total    27          47      74   
                                    
           Rne         Re      Rt   
  Risk     0.74        0.53    0.61 
                                         Estimate Lower95ci Upper95ci
 Risk difference (Re - Rne)              -0.21    -0.42     0.01     
 Risk ratio                              0.72     0.47      1.1      
 Protective efficacy =(Rne-Re)/Rne*100   28.2     -9.85     52.81    
   or percent of risk reduced                                        
 Number needed to treat (NNT)            4.79     -76.39    2.4      
   or -1/(risk difference)                                           

epiDisplay::cc(Oswego$ill, Oswego$vanilla, design = "case-control")

          Oswego$vanilla
Oswego$ill FALSE TRUE Total
     FALSE    18   11    29
     TRUE      3   43    46
     Total    21   54    75

OR =  23.45 
95% CI =  5.84, 94.18  
Chi-squared = 27.22, 1 d.f., P value = 0
Fisher's exact test (2-sided) P value = 0 

Codebook

epiDisplay tiene una función util para generar un libro de codigos

epiDisplay::codebook(Oswego)

 
 
age      :    
 obs. mean   median  s.d.   min.   max.  
 75   36.813 36      21.452 3      77    

 ================== 
sex      :    
  Frequency Percent
F        44    58.7
M        31    41.3

 ================== 
timesupper   :    
 obs. mean     median  s.d.    min.   max.  
 28   1930.714 1930    440.639 0      2200  

 ================== 
ill      :    
      Frequency Percent
FALSE        29    38.7
TRUE         46    61.3

 ================== 
onsetdate    :    
A character vector 

 ================== 
onsettime    :    
 obs. mean     median  s.d.     min.   max.  
 46   1132.065 630     1057.898 30     2400  

 ================== 
bakedham     :    
      Frequency Percent
FALSE        29    38.7
TRUE         46    61.3

 ================== 
spinach      :    
      Frequency Percent
FALSE        32    42.7
TRUE         43    57.3

 ================== 
mashedpota   :    
      Frequency Percent
FALSE        37      50
TRUE         37      50

 ================== 
cabbagesal   :    
      Frequency Percent
FALSE        47    62.7
TRUE         28    37.3

 ================== 
jello    :    
      Frequency Percent
FALSE        52    69.3
TRUE         23    30.7

 ================== 
rolls    :    
      Frequency Percent
FALSE        38    50.7
TRUE         37    49.3

 ================== 
brownbread   :    
      Frequency Percent
FALSE        48      64
TRUE         27      36

 ================== 
milk     :    
      Frequency Percent
FALSE        71   94.67
TRUE          4    5.33

 ================== 
coffee   :    
      Frequency Percent
FALSE        44    58.7
TRUE         31    41.3

 ================== 
water    :    
      Frequency Percent
FALSE        51      68
TRUE         24      32

 ================== 
cakes    :    
      Frequency Percent
FALSE        35    46.7
TRUE         40    53.3

 ================== 
vanilla      :    
      Frequency Percent
FALSE        21      28
TRUE         54      72

 ================== 
chocolate    :    
      Frequency Percent
FALSE        27    36.5
TRUE         47    63.5

 ================== 
fruitsalad   :    
      Frequency Percent
FALSE        69      92
TRUE          6       8

 ==================