Contenido

Introducción

Este corto texto fue escrito acompaña la serie de videos sobre estadística aplicada a la investigación que está disponible en mi canal de YouTube DaniMedi. https://youtube.com/playlist?list=PLiR4mMxzSHWhy6hjLaYe0tDqOIr-RhikY

Muchas veces la parte estadística y el manejo de datos es una limitación para realizar trabajos de investigación, sin embargo, esto puede ser realizado de forma sencilla y no debe ser una limitante para investigar. En esta serie de videos se presentarán conceptos básicos de epidemiología, estadística y manejo de datos para los tipos de investigación más comunes. Vengo del área médica, por lo que los tipos de estudio y los ejemplos irán orientados hacia esta área, pero conceptos generales del manejo de datos, estadística descriptiva y muchas de las funciones que se presentarán son generalizables para diferentes campos.

Lógicamente, no cubriré todos los temas sobre estadística, epidemiología, ni manejo de datos, pero creo que pueden ser de utilidad para comenzar. En realidad, procuraré no tocar temas teóricos que ya están muy bien explicados en distintos libros y diversas otras fuentes, incluso las clases de universidad. Incluiré en lo posible algunas referencias en donde se explique la teoría detrás de algunos conceptos que no se expliquen en los videos y al final de la serie brindaré también algunos recursos con los cuales continuar aprendiendo sobre el tema.

Los materiales usados estarán referenciados en este texto y tanto este texto como el data set sobre cigarro y cáncer (que se presentará más adelante) estarán disponibles en la carpeta Drive. Dejaré el link de la carpeta en la descripción de todos los videos.

Yendo a aspectos más técnicos, los temas serán desarrollados en el programa R, debido a que es el programa en el que me siento más cómodo, tiene una gran comunidad, es muy completo, relativamente sencillo de aprender y, lo más importante, es completamente gratis. También se pueden realizar los mismos conceptos en programas como SAS, STATA, Excel, Python, o incluso a mano. Pero quizá sea más sencillo seguir el proceso si se usa el mismo programa. Además, este texto contiene el código en R que usaré en todos los videos, así que pueden acompañar los videos con este texto o pueden manejar el texto de forma independiente.

Para los que usarán el programa junto con los videos, dedicaré un video/subtítulo para descargar R y otro para explicar aspectos básicos del uso de R y RStudio, pueden saltarse si ya se tiene el programa o si no usarán el programa.

Descargando R, RStudio, Rtools

Descargar primero R y Rtools, RStudio es “opcional” pero es bastante bueno, también existen otros editores, pero la mayoría concuerda en que RStudio es el mejor, y comparto esa postura. Es necesario “activar” Rtools antes de comenzar a descargar paquetes.

Para activar Rtools:

writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")
Sys.which("make")
## "C:\\rtools40\\usr\\bin\\make.exe"

Aspectos básicos de R y RStudio

Podemos realizar operaciones matemáticas usando los mismos símbolos que una calculadora.

1 + 4 / (3 - 1)
## [1] 3

Además, podemos usar <- o = para asignar variables, para esta serie de videos podrán usarse de forma indistinta cuando se utilizan para asignar valores a variables, así que pueden elegir el que más les guste.

x <- 6
y = 10
z <- (x + y) / 2 ^ 3
z
## [1] 2

Podemos entender a los vectores como una serie de números, y podemos crearlos usando c(), como se ve:

x <- c(2, 5, 1, 3, 3)
x
## [1] 2 5 1 3 3

Podemos hacer operaciones matemáticas usando vectores. Aquí sumamos 1 a todos los valores del vector y guardamos el resultado en la variable y.

y <- x + 1
y
## [1] 3 6 2 4 4

Al sumar vectores de la misma longitud, la suma será elemento por elemento.

x + y
## [1]  5 11  3  7  7

Los comandos son funciones que realizan una acción específica. Para realizar una función simplemente se coloca el nombre de la función y seguido de paréntesis, en el interior de los paréntesis se colocan los parámetros de la función (dependiendo de la función), separados por una coma. Ya vimos a la función c() que crea vectores.

Aquí tenemos a la función sum() que permite obtener la suma.

sum(x)
## [1] 14

Pueden obtenerse información de ayuda para todas las funciones e incluso objetos utilizando ? o help(). Aquí estamos obteniendo ayuda para la función sum(), se puede usar cualquiera de las dos formas de obtener ayuda:

help(sum)

Los paquetes tienen gran importancia y agregan mucha funcionalidad al programa. Aquí instalaremos un paquete llamado “epiR” que agrega funciones relacionadas a epidemiología y usaremos una de sus funciones en videos futuros. Para instalar paquetes se usa la función install.packages()

install.packages("epiR")

R ya incluye algunos data sets y aquí observamos la información y las primeras 6 filas del data set infert.

# help(infert)
head(infert)
##   education age parity induced case spontaneous stratum pooled.stratum
## 1    0-5yrs  26      6       1    1           2       1              3
## 2    0-5yrs  42      1       1    1           0       2              1
## 3    0-5yrs  39      6       2    1           0       3              4
## 4    0-5yrs  34      4       2    1           0       4              2
## 5   6-11yrs  35      3       1    1           1       5             32
## 6   6-11yrs  36      4       2    1           1       6             36

Para visualizar el data set completo en RStudio podemos usar la función View(), la cual nos mostrará el data set en otra ventana.

View(infert)

Para trabajar con data sets es más sencillo utilizar la función attach() para tener acceso a las columnas (variables) de una forma más sencilla. Aquí, asignamos el data set infert a la variable dat y usamos attach() para trabajar con este data set.

dat <- infert
attach(dat)

Ahora, podemos acceder a las columnas como si fueran variables normales, en forma de vectores.

parity
##   [1] 6 1 6 4 3 4 1 2 1 2 2 4 1 3 2 2 5 1 3 1 1 2 2 1 2 2 2 3 4 4 1 2 1 2 1 2 1
##  [38] 3 3 1 3 1 1 2 1 1 2 4 2 3 3 5 3 1 2 2 2 2 1 2 2 1 1 2 2 1 1 3 3 1 1 1 1 6
##  [75] 2 1 2 1 1 1 2 1 1 6 1 6 4 3 4 1 2 1 2 2 4 1 3 2 2 5 1 3 1 1 2 2 1 2 2 2 3
## [112] 4 4 1 2 1 2 1 2 1 3 3 1 3 1 1 2 1 1 2 4 2 3 3 5 3 1 2 2 2 2 1 2 2 1 1 2 2
## [149] 1 1 3 3 1 1 1 1 2 1 2 1 1 1 2 1 1 6 1 6 4 3 4 1 2 1 2 2 4 1 3 2 2 5 1 3 1
## [186] 1 2 2 1 2 2 2 3 4 4 1 2 1 2 1 2 1 3 3 1 3 1 1 2 1 1 2 4 2 3 3 5 3 1 2 2 2
## [223] 2 1 2 2 1 1 2 2 1 1 3 3 1 1 1 1 6 2 1 2 1 1 1 2 1 1

Y del mismo modo, podemos realizar funciones con estas variables, como si fueran vectores normales. Aquí estamos sumando los valores en la columna (variable) parity.

sum(parity)
## [1] 519

Estas funciones y unas pocas más que explicaré más adelante serán las necesarias para realizar y comprender los contenidos de los videos, si se quiere aprender de una forma más completa, hay diferentes tutoriales gratis en internet, con el que aprendí yo es con ‘swirl’ (un paquete), si se desea puede usarse para aprender lo básico, pero no es necesario para comprender y realizar lo que se desarrollará en adelante.

install.packages("swirl")
library(swirl)

Cómo ingresar datos a un programa estadístico

Al ingresar los datos a un programa estadístico lo hacemos a través de una base de datos, hay diversas formas de crear esto, para mí, la más sencilla es crear CSV files usando Excel.

Correcta forma de agregar información a una base de datos:

Para leer un archivo CSV usamos el código read.csv() y en el interior colocamos entre paréntesis el nombre del archivo CSV que queremos leer, este archivo debe estar en nuestra carpeta de trabajo para poder leerlo, esto ya lo veremos más adelante; si se quiere buscar el archivo manualmente podemos usar dentro el código file.choose() como se muestra, pero es más cómodo simplemente tenerlo en la carpeta de trabajo.

dat <- read.csv("cigarro-cancer.csv")
dat <- read.csv(file.choose())

Estadística descriptiva básica

Variables nominales: frecuencia, proporción, razón

Variables continuas:

Aquí observamos que la distribución de este data set es normal.

# help(lh)
dat <- lh
hist(dat)

Para valores con distribución normal, es de gran utilidad la media aritmética (mean()) y la desviación estándar (sd()) para describir la información.

mean(dat)
## [1] 2.4
sd(dat) # desviación estándar de la muestra
## [1] 0.5515934

La distribución de este otro data set no es normal.

# help(infert)
dat <- infert
hist(dat$parity)

En estos casos es mejor usar diagramas de cajas en lugar de histogramas para visualizar la información (aunque es más útil cuando se compara entre diferentes grupos).

boxplot(dat$parity)

Para los valores que no tengan distribución normal (presencia de outliers), es preferible utilizar otras medidas, como la media (median()) y el rango intercuartil (IQR()), que son menos sensibles a los valores extremos (outliers).

median(dat$parity)
## [1] 2
IQR(dat$parity)
## [1] 2

p-value

La mayoría de análisis estadísticos se basa en la comparación entre diferentes grupos, pero diferencias entre los grupos pueden encontrarse simplemente por el azar, por lo que debe haber alguna manera de filtrar a las diferencias “sospechosas” que tengan una gran probabilidad de ser obtenidas por azar y no porque de verdad haya una diferencia por aquello que estemos estudiando.

El valor de p (p-value) es usado para determinar significancia estadística en una diferencia encontrada. Nos indica qué tan probable es que la diferencia encontrada o una diferencia más grande sea producida simplemente por el azar. Para encontrarla hay diversas formas como pruebas estadísticas o usando fórmulas y conceptos del teorema del límite central (central limit theorem).

Un valor de p de 0.05 nos indica que hay 5% de probabilidades de que la diferencia encontrada o una diferencia mayor sea producida por el azar, en la investigación en el campo biológico y de ciencias de la salud un valor de p < 0.05 indica que la probabilidad de que los hallazgos sean obtenidos por el azar es lo suficientemente baja y se considera que los hallazgos son estadísticamente significativos. Cabe resaltar que el valor de p de un resultado binario: lo encontrado es estadísticamente significativo o lo encontrado no es estadísticamente significativo, además, si se entendió el concepto podemos darnos cuenta de que el valor de p solo indica probabilidades y no indica nada sobre lo encontrado. Por otro lado, los intervalos de confianza, el tema del siguiente video, sí brindan una estimación de los resultados.

Por este motivo es muy importante diferenciar la significancia estadística de la significancia clínica. Por ejemplo, en un estudio en el que buscamos encontrar la diferencia entre dos grupos obtenemos una diferencia con un valor de p sumamente pequeño, pero la diferencia que encontramos es insignificante entre ambos grupos. La situación mencionada puede ser perfectamente posible, es más, el valor de p disminuye mientras aumenta el tamaño de la muestra (al disminuir la variabilidad en los grupos), por lo que se puede obtener significancia estadística con diferencias mínimas entre grupos comparados solo por contar con un gran tamaño de muestra. Por otro lado, a pesar de tener una gran diferencia entre grupos comparados, muchos estudios no alcanzan la significancia estadística debido a un pequeño tamaño de muestra.

Por este motivo es necesario diferenciar la significancia estadística y la significancia clínica y se debe considerar también qué tan grande es la diferencia encontrada y qué tanto impacto tiene esta diferencia en el mundo real.

Intervalos de confianza

Un intervalo de confianza es un rango en el que un determinado valor se encontrará con un determinado nivel de confianza, en el campo biológico y de ciencias de la salud se usa un nivel de confianza de 95% (que es equivalente al p-value de 0.05). Los intervalos de confianza pueden reemplazar en muchos casos al valor de p y tienen la ventaja de incluir a la estimación de la diferencia encontrada. Por ejemplo, si se utiliza los intervalos de confianza para hayar la diferencia entre los promedios de peso de dos grupos y el intervalo de confianza del 95% para esta diferencia no incluye al 0, esto nos indica que podemos decir con un 95% de certeza de que existe diferencia real entre los grupos (que la diferencia no sea 0), esto nos da significancia estadística.

Por lo mencionado, actualmente existe una tendencia a preferir a los intervalos de confianza sobre los valores de p, pero en realidad son conceptos complementarios y aún hay algunas revistas científicas que requieren el valor de p para los resultados. Cabe mencionar que hay muchos errores en la comprensión y uso de los conceptos de los intervalos de confianza y el valor de p, hay diferentes artículos que desarrollan este tema, encontré bastante informativo el siguiente: Confidence Interval or P-Value?

En los siguientes videos se usará a los intervalos de confianza para determinar la significancia estadística de los hallazgos y también dedicaré un video para mostrar como obtener el valor de p usando la prueba de chi cuadrado y la prueba exacta de Fisher.

Tablas de contingencia

Tablas de contingencia son tablas 2x2 que se utilizan para comparar datos, aquí hay un ejemplo:

Enfermos No Enfermos
Expuestos a b
No Expuestos c d

Usaremos tablas similares para realizar métodos estadísticos en videos posteriores.

En este ejemplo se usará un pequeño archivo CSV que se usará en los siguientes videos para ejemplificar los conceptos. Pueden descargarlo desde la carpeta Drive del curso. Asegurarse de que el archivo esté en la misma carpeta en la que esté usando el programa (generalmente en Documentos si no se ha cambiado).

Otra opción es usar el siguiente código que descargará el archivo y lo pondrá en la carpeta de trabajo:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Abrimos el archivo y observamos las 6 primeras filas para comprobar que el archivo se leyó de forma adecuada.

dat <- read.csv("cigarro-cancer.csv")
head(dat)
##   cigarro cancer
## 1       0      0
## 2       0      0
## 3       0      0
## 4       0      0
## 5       0      0
## 6       0      0

Ahora usamos attach() para comenzar a trabajar con este data set (ignorar los mensajes de advertencia):

attach(dat)

A esto se le conoce como data set y este contiene 2 columnas (variables): “cigarro” y “cancer”, con valores de 0 que representan “No” y valores de 1 que representan “Sí”. Estas variables adquirirán diferentes significados en los ejemplos que se utilizarán en videos posteriores. Por ahora, nos servirá para mostrar cómo crear una tabla de contingencia.

La tabla de contingencia puede crearse con la función table(), y en este ejemplo la guardamos como variable tab.

tab <- table(cigarro, cancer)
tab
##        cancer
## cigarro   0   1
##       0 330  59
##       1  84  27

Es necesario acomodar la tabla de contingencia para calcular OR, RR, etc., para esto tenemos que recordar la estructura de las tablas de contingencia. En una tabla de contingencia, la primera fila debería ser las personas expuestas al cigarro y la primera columna debería ser las personas con cáncer, si nos fijamos, en nuestra tabla de contingencia tenemos los valores al revés.

Comparar:

Enfermos No Enfermos
Expuestos a b
No Expuestos c d 
tab
##        cancer
## cigarro   0   1
##       0 330  59
##       1  84  27

Podemos invertir estos valores de diferentes maneras y cada persona es libre de elegir la forma que desee, si conoce otra. En mi opinión la más sencilla, y la que usaré en los videos, es la siguiente:

tab[c(2,1), c(2,1)]
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Aclaración

Nuestro data set solo tiene 2 columnas, por lo que no es necesario especificar las columnas para generar nuestra tabla de contingencia, puede hacerse de esta forma simplificada, sin especificar los nombres de las columnas (en los siguientes videos usaré esa forma, pero recuerden que en sí estamos generando una tabla de contingencia a partir de dos columnas).

tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Al tener la tabla de contingencia en la forma que queremos ya estamos listos para calcular valores como el odds ratio, riesgo relativo, y otros conceptos que explicaré en videos posteriores.

Odds ratio y risk ratio (riesgo relativo)

Odds ratio (OR) y risk ratio (RR), también conocido como riesgo relativo, son medidas de fuerza y de dirección de una asociación, por lo que son de gran utilidad en estudios sobre factores de riesgo o de protección. Esta dirección de la asociación es algo que no te brindan de forma directa las pruebas estadísticas de chi-cuadrado o de Fisher, las que también podrían aplicarse de forma similar; además, los valores de OR y RR pueden interpretarse de forma relativamente sencilla, por este motivo, en videos posteriores se mostrará cómo podemos calcular estos valores en los estudios epidemiológicos más comunes.

Las fórmulas para calcular OR y RR pueden mostrarse en base a la tabla de contingencia:

Enfermos No Enfermos
Expuestos a b
No Expuestos c d

El riesgo relativo es más sencillo de interpretar, indica la probabilidad de que una persona expuesta a algo presente o no presente un determinado resultado (generalmente enfermedad), en comparación con los no expuestos.

\[ RR = \frac{a/(a+b)}{c/(c+d)} \]

El odds ratio nos indica la probabilidad de que una persona con un determinado resultado (generalmente enfermedad) haya tenido exposición o no a un determinado factor, en comparación con los que no presentan la condición (sanos).

\[ OR = \frac{a/c}{b/d} = \frac{a*d}{b*c} \]

Si analizamos ambas ecuaciones obtenidas de la tabla de contingencia, un resultado de 1 nos indica que no existe diferencia en el riesgo entre dos grupos (sea por exposición o por resultado/enfermedad). Valores mayores que 1 nos indican un aumento en el riesgo y valores menores que 1 nos indican una disminución en el riesgo.

Los conceptos de RR y OR son muy similares, pero si se analizan, tienen un diferente sentido. El RR nos permite determinar el riesgo a partir de la causa, por lo que el RR nos permite obtener una relación de causalidad. Es mejor obtener el RR, sin embargo, este solo puede ser obtenido mediante estudios cuyo diseño vaya de la causa al efecto, de lo contrario solo se podrá obtener el OR, se mencionarán estos conceptos en los siguientes videos de estudios.

Como se puede observar, las fórmulas son muy sencillas y pueden obtenerse de forma simple, sin embargo, es de suma utilidad poder calcular también los intervalos de confianza que nos permitan determinar la significancia estadística de estos valores, como ya se explicó en el video sobre los intervalos de confianza.

Intervalos de confianza para OR y RR

Es importante obtener los intervalos de confianza para saber si estos resultados son estadísticamente significativos, si el intervalo de confianza incluye al valor de 1, esto nos indica que el resultado no es estadísticamente significativo, por otro lado, si el intervalo de confianza no incluye al valor de 1, esto nos indica que el resultado sí es estadísticamente significativo.

Como veremos en otros videos podremos calcular los intervalos de confianza con simplemente una función, pero se pueden calcular estos intervalos de confianza mediante una fórmula, la presento aquí para que puedan calcular estos intervalos si es que no se está usando el programa.

\[ Intervalo = OR ^ {1 \pm Z/X} \]

\[ Intervalo = RR ^ {1 \pm Z/X} \]

\(Z\) = z-score (1.96)
\(X\) = raíz cuadrada del chi-cuadrado (prueba estadística)

Nótese que ambas ecuaciones son iguales y que nos brindan un límite inferior y otro superior, que son los límites de nuestro intervalo de confianza. Otra cosa a tener en cuenta es que para el cálculo necesitamos de ‘X’ que proviene de la prueba estadística de chi-cuadrado, esta es una prueba sencilla de realizar usando nuestra tabla de contingencia, esta prueba se explicará en un video futuro.

Clasificaciones de los estudios científicos

El objetivo de esta serie de videos es brindar herramientas para poder realizar estudios científicos sencillos, por este motivo es importante tener en cuenta aspectos básicos sobre estudios científicos.

Hay diversas formas de clasificar los estudios científicos, pero la forma que encontré más sencilla es la que explicaré, pero si no les gusta pueden usar la clasificación de su mayor agrado.

Cabe resaltar que en esta clasificación no entra en ningún momento los términos de “retrospectivo” o “prospectivo”. Estos términos son muy usados y generalmente hacen referencia al tiempo, los estudios retrospectivos estudian datos que ya existen, datos que fueron obtenidos en el pasado, mientras que los estudios prospectivos son estudios en los que no contamos con los datos aún, obtenemos los datos en el futuro, es decir, generamos datos nuevos.

Estudios de casos y controles

Un estudio de casos y controles es un estudio observacional que se caracteriza porque se comienza de un determinado resultado y se examinan los factores que han podido producir este resultado, en otras palabras se parte de un efecto (generalmente enfermedad) y se busca encontrar una causa a este efecto. Un estudio de casos y controles generalmente comienza con 2 grupos: un grupo de personas con el efecto (casos) y un grupo de personas sin el efecto (controles), la exposición es lo que se examinará en el estudio. Es importante hacer notar que estos estudios se definen por el diseño del estudio y no por el tiempo en el que se realizan, a pesar de que generalmente se asocian a estudios retrospectivos, también pueden haber estudios de casos y controles prospectivos. Aquí hay un ejemplo de un estudio de casos y controles prospectivo: Increased risk for an atypical autism diagnosis following Thimerosal-containing vaccine exposure in the United States: A prospective longitudinal case-control study in the Vaccine Safety Datalink.

Para entender en qué consiste este tipo de estudios usaremos un ejemplo, imaginemos que comenzamos un estudio con 2 grupos, un grupo de personas con cáncer (casos) y un grupo de personas sin cáncer (controles); entonces, obtenemos cuántas personas fuman o fumaban en ambos grupos. Estos datos se encuentran en nuestro data set de cigarro y cáncer, si no lo hemos descargado, podemos hacerlo con el siguiente código:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Leemos el data set, lo guardamos en una variable y usamos attach() para trabajar con los datos. Además, debemos cargar el paquete “epiR” usando la función library() (ya descargamos este paquete en un video previo, para descargarlo nuevamente se usa el comando: install.packages("epiR")).

dat <- read.csv("cigarro-cancer.csv")
attach(dat)
library(epiR)
# es normal que aparezcan mensajes de advertencia

Ahora podemos crear la tabla de contingencia y ponerla en el formato adecuado, como lo presentamos en el video de tablas de contingencia, y observamos los datos obtenidos.

tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Si recordamos del video sobre el odds ratio y el risk ratio (riesgo relativo), lo ideal es poder obtener el RR, pero en este caso no podemos, ¿por qué?, porque el diseño de este estudio no va de la causa al efecto, ya que comenzamos el estudio con el efecto. No se desarrollará la teoría detrás de este concepto.

Para obtener el OR de esta asociación en nuestro estudio (no podemos obtener el RR), usaremos la función epi.2by2() en la que usamos el argumento "case.control" para precisar de que se trata de un estudio de casos y controles. Podemos usar ? para obtener ayuda sobre esta función.

# help(epi.2by2)
epi.2by2(tab, "case.control")
##              Outcome +    Outcome -      Total        Prevalence *        Odds
## Exposed +           27           84        111                24.3       0.321
## Exposed -           59          330        389                15.2       0.179
## Total               86          414        500                17.2       0.208
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Odds ratio (W)                               1.80 (1.07, 3.01)
## Attrib prevalence *                          9.16 (0.42, 17.90)
## Attrib prevalence in population *            2.03 (-2.83, 6.90)
## Attrib fraction (est) in exposed  (%)        44.31 (2.81, 67.54)
## Attrib fraction (est) in population (%)      13.93 (-0.10, 26.00)
## -------------------------------------------------------------------
##  Test that OR = 1: chi2(1) = 5.085 Pr>chi2 = 0.02
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Como resultado tenemos una tabla con los expuestos y no expuestos y otros valores, y también obtenemos una lista con resultados. En esta lista podemos observar un valor de OR de 1.80, lo que nos indica que es más probable que las personas con cáncer fumen o hayan fumado. Ahora, evaluamos el intervalo de confianza para el OR y encontramos que es de 1.07 a 3.01, por lo que no incluye al 1, esto nos dice que nuestros hallazgos son estadísticamente significativos, genial!, como recordamos del video de OR y RR.

Recordar nuevamente que resultados estadísticamente significativos no son lo mismo que resultados clínicamente significativos y se tiene que evaluar la magnitud de este OR en la práctica clínica, además, tenemos que tener en cuenta que con el OR no encontramos una relación de causalidad, por lo que con los datos brindados, el tipo de estudio y nuestros resultados, no podemos concluir que fumar aumenta el riesgo de sufrir cáncer.

¿Entonces por qué haríamos este tipo de estudio?

Por último, mostraré que las fórmulas (del video de OR y RR) dan resultados similares (el código es solo para ejemplificar, no es necesario aprenderlo).

OR <- (tab[1,1] * tab[2,2]) / (tab[1,2] * tab[2,1])
Z <- 1.96
chi <- as.numeric(chisq.test(tab)$statistic)
ivl1 <- OR^(1 + Z / sqrt(chi))
ivl2 <- OR^(1 - Z / sqrt(chi))

conf_int <- ifelse(ivl1 < ivl2, paste(ivl1,"-",ivl2), paste(ivl2,"-",ivl1))
data.frame(c("odds ratio" = OR, "intervalo de confianza" = conf_int))
##                        c..odds.ratio....OR...intervalo.de.confianza....conf_int.
## odds ratio                                                      1.79782082324455
## intervalo de confianza                        1.0432192595209 - 3.09825540795335

Observamos que el OR calculado y los intervalos de confianza son muy similares a los obtenidos anteriormente, y en ambos casos los resultados son estadísticamente significativos, ya que el intervalo de confianza no incluye al 1.

Estudios de cohortes

Un estudio de cohortes es un estudio observacional que se caracteriza porque se comienza de la exposición y se estudian los resultados que se presentan por esta exposición. Esto quiere decir que comenzamos con grupos de sujetos de estudio expuestos a un factor y no expuestos a un factor, entonces, examinamos la aparición de resultados. Es importante hacer notar que estos estudios se definen por el diseño del estudio y no por el tiempo en el que se realizan, a pesar de que generalmente se asocian a estudios prospectivos, también pueden haber estudios de cohortes retrospectivos. Aquí hay un ejemplo de un estudio de cohortes retrospectivo: Predicting the success of vaginal birth after caesarean delivery: a retrospective cohort study in China.

Ahora, vamos con el ejemplo, imaginemos que hacemos un estudio en el que tengamos dos grupos, un grupo de fumadores (expuestos) y un grupo de no fumadores (no expuestos); entonces, obtenemos cuántas personas presentan cáncer después de un periodo de tiempo. Estos datos se encuentran en el data set de cigarro y cáncer, podemos descargarlo, si no lo hemos hecho aún, usando el siguiente código:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Leemos este data set, lo guardamos en una variable y usamos attach() para trabajar con los datos. Además, debemos cargar el paquete “epiR” usando la función library() (ya descargamos este paquete en un video previo, para descargarlo nuevamente se usa el comando: install.packages("epiR")).

dat <- read.csv("cigarro-cancer.csv")
attach(dat)
library(epiR)
# es normal que aparezcan mensajes de advertencia

Ahora podemos crear la tabla de contingencia y ponerla en el formato adecuado, como lo presentamos en el video de tablas de contingencia, y observamos los datos obtenidos.

tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Recordemos que es preferible obtener el risk ratio (riesgo relativo) en un estudio y este tipo de estudio sí nos permite obtener el RR, ¿por qué?, porque el diseño de este estudio sí va de la causa al efecto.

Para obtener el RR de esta asociación en nuestro estudio usaremos la función epi.2by2() acompañado del argumento "cohort.count" con lo que precisamos que se trata de un estudio de cohortes. Podemos usar ? para obtener ayuda sobre esta función.

# help(epi.2by2)
epi.2by2(tab, "cohort.count")
##              Outcome +    Outcome -      Total        Inc risk *        Odds
## Exposed +           27           84        111              24.3       0.321
## Exposed -           59          330        389              15.2       0.179
## Total               86          414        500              17.2       0.208
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                               1.60 (1.07, 2.40)
## Odds ratio                                   1.80 (1.07, 3.01)
## Attrib risk *                                9.16 (0.42, 17.90)
## Attrib risk in population *                  2.03 (-2.83, 6.90)
## Attrib fraction in exposed (%)               37.65 (6.64, 58.35)
## Attrib fraction in population (%)            11.82 (-0.24, 22.43)
## -------------------------------------------------------------------
##  Test that OR = 1: chi2(1) = 5.085 Pr>chi2 = 0.02
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Como resultado tenemos una tabla con los expuestos y no expuestos y otros valores, y también obtenemos una lista con resultados. En esta lista podemos observar que tenemos un RR de 1.60, lo que nos indica que es más probable que las personas que fuman presenten cáncer. Ahora, evaluamos el intervalo de confianza para el RR y encontramos que es de 1.07 a 2.40, por lo que no incluye al 1, esto nos dice que nuestros hallazgos son estadísticamente significativos, genial!, como recordamos del video de OR y RR.

Recordar nuevamente que resultados estadísticamente significativos no son lo mismo que resultados clínicamente significativos y se tiene que evaluar la magnitud de este RR en la práctica clínica; además, tenemos que tener en cuenta otros factores para determinar una relación de causalidad, por ejemplo, eliminar la presencia de variables que podrían ocasionar confusión en los resultados y otras limitaciones en nuestro estudio.

¿Ventajas de este tipo de estudio (en comparación a un estudio de casos y controles)?

¿Por qué no hacemos siempre este estudio (en lugar de estudios de casos y controles)?

Por último, mostraré que las fórmulas (del video de OR y RR) dan resultados similares (el código es solo para ejemplificar, no es necesario aprenderlo).

RR <- (tab[1,1] / sum(tab[1,])) / (tab[2,1] / sum(tab[2,]))
Z <- 1.96
chi <- as.numeric(chisq.test(tab)$statistic)
ivl1 <- RR^(1 + Z/ sqrt(chi))
ivl2 <- RR^(1 - Z/ sqrt(chi))

conf_int <- ifelse(ivl1 < ivl2, paste(ivl1,"-",ivl2), paste(ivl2,"-",ivl1))
data.frame(c("risk ratio" = RR, "intervalo de confianza" = conf_int))
##                        c..risk.ratio....RR...intervalo.de.confianza....conf_int.
## risk ratio                                                      1.60375629867155
## intervalo de confianza                       1.03465896526348 - 2.48587636301369

Observamos que el RR calculado y los intervalos de confianza son similares a los obtenidos anteriormente, y en ambos casos los resultados son estadísticamente significativos, ya que el intervalo de confianza no incluye al 1.

Estudios transversales

Los estudios transversales, son estudios observacionales que se realizan en un solo punto en el tiempo. Estos estudios son muy útiles para describir características de poblaciones, por ejemplo, son muy convenientes para obtener la prevalencia de enfermedades. Estos estudios pueden limitarse a la descripción, pero también podría incluirse en el estudio la presencia de posibles factores de riesgo o de protección, para poder realizar un análisis de los datos obtenidos, a esto se le conoce como estudio transversal analítico.

Esto podrá entenderse mejor con un ejemplo en el que usaremos el data set de cigarro y cáncer. Para este ejemplo pensemos que hicimos un estudio en una población, en donde, realizamos una encuesta casa por casa y preguntamos si en la casa hay algún integrante que fuma y si en la casa alguien padeció de cáncer. Entonces podríamos proceder a hacer un análisis de la asociación entre fumar y el cáncer. Podemos descargar nuestro data set de cigarro y cáncer, si no lo hemos hecho, con el siguiente código:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Para esto, leemos este data set, lo guardaremos en una variable y usamos attach() para trabajar con los datos. Además, debemos cargar el paquete “epiR” usando la función library() (ya descargamos este paquete en un video previo, para descargarlo nuevamente se usa el comando: install.packages("epiR")).

dat <- read.csv("cigarro-cancer.csv")
attach(dat)
library(epiR)
# es normal que aparezcan mensajes de advertencia

Ahora podemos crear la tabla de contingencia y ponerla en el formato adecuado, como lo presentamos en el video de tablas de contingencia.

tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Con los datos que tenemos podemos obtener el odds ratio (OR), para esto usamos al argumento "cross.sectional" que precisa que se trata de un estudio transversal.

# help(epi.2by2)
epi.2by2(tab, "cross.sectional")
##              Outcome +    Outcome -      Total        Prevalence *        Odds
## Exposed +           27           84        111                24.3       0.321
## Exposed -           59          330        389                15.2       0.179
## Total               86          414        500                17.2       0.208
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Prevalence ratio                             1.60 (1.07, 2.40)
## Odds ratio                                   1.80 (1.07, 3.01)
## Attrib prevalence *                          9.16 (0.42, 17.90)
## Attrib prevalence in population *            2.03 (-2.83, 6.90)
## Attrib fraction in exposed (%)              37.65 (6.64, 58.35)
## Attrib fraction in population (%)           11.82 (-0.24, 22.43)
## -------------------------------------------------------------------
##  Test that OR = 1: chi2(1) = 5.085 Pr>chi2 = 0.02
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Vemos que nuestros resultados son estadísticamente significativos, pero, si lo pensamos, podemos darnos cuenta de los problemas con este diseño. Recordemos que en este tipo de estudios no contamos con la variable tiempo, si recordamos, en ningún momento determinamos si la presencia de un integrante que fuma fue antes o después de la aparición de cáncer en esa casa. Esto hace que en este tipo de estudios sea muy difícil encontrar una relación de causalidad, debido a lo mencionado y debido a que este tipo de estudio es muy propenso a muchas otras variables que pueden influir en la asociación a analizar. Sin embargo, es posible hablar de asociación o potencialmente causalidad cuando sabemos que la exposición estudiada es estable con el tiempo (ej. genes).

Entonces, un estudio transversal puede brindarnos datos de una población en poco tiempo, ya que solo recolectamos la información una vez, y es muy útil para realizar descripciones y también para generar hipótesis. Además, podemos incluir un análisis en el estudio también e incluso llegar a conclusiones similares a estudios de casos y controles. Sin embargo, tenemos que ser cuidadosos con el diseño del estudio y prudentes al momento de llegar a conclusiones, debido a los factores ya mencionados.

Un artículo donde se expone lo mencionado en mayor detalle es el siguiente: “Cross‐sectional studies – what are they good for?”.

Estudios experimentales

Los estudios descritos en los anteriores videos corresponden a estudios observacionales, esto significa que, como investigadores, no separamos los grupos de acuerdo a una intervención que nosotros mismos generamos, simplemente observamos sucesos que ocurrirían aunque no hicieramos el estudio. En los estudios experimentales sí realizamos una intervención y creamos los grupos de estudio en base a esta intervención. En realidad, podríamos pensar en estos estudios como si fueran de tipo cohortes y prospectivo, la intervención que generamos sería la exposición y analizamos los cambios que ocurran en el futuro por nuestra intervención; al igual que un estudio de cohortes, necesitamos un grupo control que sería el grupo al que no realizamos una intervención o, en todo caso, usamos un placebo (este sería el grupo de no expuestos en el estudio de cohortes). Es importante hacer notar también que la intervención puede ser de cualquier tipo y no solo llevada a cabo en un laboratorio, por ejemplo, en el estudio del impacto de una charla en un grupo de estudiantes la charla sería una intervención, por lo que nuestro estudio sería un estudio experimental.

Ahora, vamos con un ejemplo. En este caso imaginemos que hacemos un estudio con ratas de laboratorio con las que formamos dos grupos: grupo de ratas expuestas a humo de cigarro las 24 horas del día y grupo de ratas no expuestas al humo de cigarro en ningún momento. Los resultados los tenemos en el mismo data set que usamos en los videos pasados. Podemos descargarlo, si no lo hemos hecho, con el siguiente código:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Leemos el data set, lo guardaremos en una variable y usaremos attach() para trabajar con los datos. Además, cargamos el paquete “epiR” usando la función library() (ya descargamos este paquete en un video previo, para descargarlo nuevamente se usa el comando: install.packages("epiR")).

dat <- read.csv("cigarro-cancer.csv")
attach(dat)
library(epiR)
# es normal que aparezcan mensajes de advertencia

Ahora podemos crear la tabla de contingencia y ponerla en el formato adecuado, como lo presentamos en el video de tablas de contingencia, y observamos los resultados.

tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Ahora, dijimos que podíamos pensar en este estudio como un estudio de cohortes prospectivo, así que podemos obtener el riesgo relativo (RR) del mismo modo, entonces usamos la función epi.2by2() con el argumento "cohort.count". Podemos usar ? para obtener ayuda sobre esta función.

# help(epi.2by2)
epi.2by2(tab, "cohort.count")
##              Outcome +    Outcome -      Total        Inc risk *        Odds
## Exposed +           27           84        111              24.3       0.321
## Exposed -           59          330        389              15.2       0.179
## Total               86          414        500              17.2       0.208
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                               1.60 (1.07, 2.40)
## Odds ratio                                   1.80 (1.07, 3.01)
## Attrib risk *                                9.16 (0.42, 17.90)
## Attrib risk in population *                  2.03 (-2.83, 6.90)
## Attrib fraction in exposed (%)               37.65 (6.64, 58.35)
## Attrib fraction in population (%)            11.82 (-0.24, 22.43)
## -------------------------------------------------------------------
##  Test that OR = 1: chi2(1) = 5.085 Pr>chi2 = 0.02
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

También se pueden obtener estos resultados mediante las fórmulas, como ya se demostró en el video sobre el estudio de cohortes.

Este diseño de estudio es el que tiene más nivel de evidencia entre los estudios epidemiológicos mencionados en el video de clasificaciones de los estudios científicos y es entendible, ya que permite controlar de mejor manera variables externas a lo que queremos estudiar y generamos datos nuevos con el propósito de contestar nuestra pregunta de investigación, por lo que se dice que son datos con mayor calidad.

Otras consideraciones a tener en cuenta en estos estudios es que se prefiere que los grupos sean formados de forma aleatoria, que los sujetos de investigación y los investigadores no puedan distinguir los grupos (doble ciego), y otras. Por ejemplo, las investigaciones experimentales realizadas en el área clínica reciben el nombre de ensayos clínicos y tienen una estructura ya establecida.

Es evidente que los estudios experimentales brindan una información de calidad superior a los estudios observacionales, pero también tienen diferentes desventajas a tener en cuenta. Los estudios experimentales suelen requerir de más infraestructura para ser llevados a cabo, por ejemplo, pueden requerir de algún laboratorio con un equipamiento necesario, además, para realizar estudios experimentales hay que tener muy presente la parte ética debido a que nosotros, como investigadores, somos los que estamos generando un cambio en los sujetos de estudio. Lógicamente hay otras limitaciones dependiendo del estudio que queramos realizar, así como puede no haber muchas limitaciones y podría ser realizado de forma sencilla, todo depende del estudio.

Prueba de chi cuadrado y prueba exacta de Fisher

Ya se mencionó las ventajas del odds ratio (OR) y el risk ratio (RR) sobre la prueba de chi-cuadrado o la prueba exacta de Fisher en el video de odds ratio y risk ratio. La diferencia está en que las pruebas estadísticas de chi-cuadrado o prueba exacta de Fisher no brindan dirección de la asociación, solo brindan el valor de p, que nos sirve para determinar la significancia estadística, algo que ya se puede obtener con los intervalos de confianza para el OR o RR. Muestro estas pruebas porque a veces se requiere el cálculo del valor de p, pero es preferible usar intervalos de confianza, de todas formas pueden usarse ambos, un artículo donde explica esto es el siguiente: Confidence Interval or P-Value?

Usaremos nuevamente nuestro data set de cigarro y cáncer, que podemos descargarlo, si aún no lo hemos hecho, con el siguiente código:

url <- "https://drive.google.com/uc?export=download&id=1lIa55z3lsS91VWat9jaDzXCmMKuBkgzk"
filename <- "cigarro-cancer.csv"
download.file(url, filename)

Hacemos lo mismo que los videos pasados y generamos la tabla de contingencia:

dat <- read.csv("cigarro-cancer.csv")
attach(dat)
library(epiR)
tab <- table(dat)
tab <- tab[c(2,1), c(2,1)]
tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Ahora, a partir de esta tabla de contingencia podemos ya realizar las pruebas de chi-cuadrado o de Fisher, con sus funciones respectivas (chisq.test() y fisher.test()), también podemos obtener ayuda para estas funciones mediante ? o help() como ya se explicó en el video de aspectos básicos de R y RStudio.

# H0: no hay diferencia entre los grupos
chisq.test(tab)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 4.4621, df = 1, p-value = 0.03465
fisher.test(tab)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab
## p-value = 0.03183
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.028886 3.080742
## sample estimates:
## odds ratio 
##   1.795568

Ahora, pueden revisar la teoría de esto por su cuenta, aquí solo diré que estas pruebas pueden usarse cuando queremos hallar la diferencia en una variable nominal entre dos grupos independientes (esto es lo que hemos tenido todo este tiempo en nuestro data set de cigarro y cáncer). La elección de una u otra es en función al tamaño de la muestra y a los números que encontremos en la tabla de contingencia de resultados esperados; podemos resumir esto de la siguiente manera: primero obtenemos la tabla de resultados esperados, luego observamos los valores de estas celdas, si tenemos celdas con valores de 5 o más usamos la prueba de chi-cuadrado, además, en ninguna celda debemos tener algún valor de 0.

Esta es la tabla de contingencia con los valores observados:

tab
##        cancer
## cigarro   1   0
##       1  27  84
##       0  59 330

Esta es la tabla de contingencia con los valores esperados (observemos que podemos extraer esta tabla con expected):

chisq.test(tab)$expected
##        cancer
## cigarro      1       0
##       1 19.092  91.908
##       0 66.908 322.092

En nuestro ejemplo es más apropiado usar la prueba de chi-cuadrado.

En realidad estas son pruebas muy sencillas y hay múltiples recursos en los que pueden explorar estos temas, vuelvo a recodar que estas pruebas estadísticas nos brindan el valor de p, pero esto solo nos indica significancia estadística y no nos indica la dirección de la asociación o una estimación de los resultados, a diferencia de los valores de RR y OR acompañados de sus intervalos de confianza.

Análisis de artículo científico publicado y manejo de variables continuas

Este video es algo especial, y esta vez estudiaremos un estudio científico publicado el año 2014. El estudio es el siguiente:

Comparison of Low-Dose Rosuvastatin with Atorvastatin in Lipid-Lowering Efficacy and Safety in a High-Risk Pakistani Cohort: An Open-Label Randomized Trial

Y podemos encontrar el data set con los datos del estudio en: https://doi.org/10.7910/DVN/U8NLQC

Podemos ir a ese link para descargar el data set como tab separated file (.tab) y luego colocar el archivo en el la misma carpeta en la que estemos trabajando. También pueden usar este código con el que descargarán el archivo en la carpeta de trabajo:

url <- "https://dataverse.harvard.edu/api/access/datafile/3992518?gbrecs=false"
filename <- "statins_comparison_raw.tab"
download.file(url, filename)

Para leer el archivo usamos esta vez read.table(), que nos permite leer archivos delimitados por tabulaciones (.tab), con los argumentos mostrados. Del mismo modo, usamos attach() para acceder de forma más sencilla a las columnas del data set. Podemos usar la función colnames() que nos dará los nombres de las columnas, también podemos visualizar el data set completo con la función View().

dat <- read.table("statins_comparison_raw.tab", header = TRUE, sep = "\t")
attach(dat)
colnames(dat)
##  [1] "StatinGROUP" "age"         "gender"      "Edn"         "HTN"        
##  [6] "Diabetes"    "Height"      "Weight"      "BMI"         "WC"         
## [11] "smoking"     "sysBP1"      "diaBP1"      "PGF1"        "CHOL1"      
## [16] "TG1"         "HDL1"        "LDL1"        "CoMorbid"    "CHOL2"      
## [21] "TG2"         "HDL2"        "LDL2"        "Bodyaches"   "CK"
View(dat)

Podemos observar que este data set tiene más variables (columnas) que nuestro data set de cigarro y cáncer que usamos de ejemplo y no solo tiene variables nominales, sino también tiene variables numéricas como peso, talla, índice de masa corporal, etc. Pero, a pesar de esta mayor complejidad, los fundamentos y la forma de trabajar con la información sigue los mismos principios. Lo que haremos a continuación será replicar algunos de los datos obtenidos en el estudio, como motivación y para demostrar los aspectos más básicos de cómo trabajar con variables continuas.

Ahora, accedemos a las columnas del mismo modo, como ya lo vimos en videos pasados, y podemos usar las mismas funciones que usamos para vectores.

mean(age)
## [1] 54.78295
mean(Height)
## [1] 161.3953
mean(Weight)
## [1] 69.46512

Ahora, en el paper vemos que los datos se obtienen de forma separada, y observamos que en el data set separan los grupos de forma adecuada con la primera columna StatinGROUP en la cual ‘1’ representa al grupo de atorvastatina y ‘2’ representa al grupo de rosuvastatina. Podemos especificar el grupo de forma muy sencilla usando los corchetes [], los observadores pudieron darse cuenta que ya estuvimos usando estos corchetes al momento de generar nuestra tabla de contingencia en videos pasados.

De este modo podemos obtener el promedio de las edades en el grupo de atorvastatina y en el grupo de rosuvastatina usando mean(), además, podemos obtener la desviación estándar usando sd(). Podemos comparar nuestros resultados con los resultados en la tabla 1 del paper y podemos hacer lo mismo para los otros valores, pueden comprobarlo.

# promedio en el grupo de atorvastatina
mean(age[StatinGROUP == 1])
## [1] 54.09524
# desviación estándar en el grupo de atorvastatina
sd(age[StatinGROUP == 1])
## [1] 12.65894
# promedio en el grupo de rosuvastatina
mean(age[StatinGROUP == 2])
## [1] 55.43939
# desviación estándar en el grupo de la rosuvastatina
sd(age[StatinGROUP == 2])
## [1] 9.42685

Podemos darnos cuenta de que estos valores están acompañados de un p-value. En realidad hay muchas maneras de obtener este p-value, en el artículo se menciona que se usa el t-test para las variables numéricas y el chi-cuadrado para las variables discretas. Podemos realizar esta prueba usando la función t.test(), podemos obtener ayuda de esta función como ya sabemos, usando ? o help(). Para realizar un t-test simple debemos agregar el argumento var.equal = TRUE (el artículo no menciona datos específicos sobre esto así que probablemente se realizó un t-test simple).

#help(t.test)
# asignamos variables con los grupos a comparar
age1 <- age[StatinGROUP == 1]
age2 <- age[StatinGROUP == 2]
# two sample t-test
t.test(age1, age2, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  age1 and age2
## t = -0.6861, df = 127, p-value = 0.4939
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.220911  2.532599
## sample estimates:
## mean of x mean of y 
##  54.09524  55.43939

Observamos que obtenemos el mismo valor de p que en el estudio, en ambos casos no encontramos diferencias significativas. Podemos hacer lo mismo para las otras variables continuas, pueden comprobarlo.

Ahora, para las variables nominales el artículo menciona que se realizó la prueba de chi-cuadrado, para esto generamos la tabla de contingencia con las columnas StatinGROUP y gender (en el tema de tablas de contingencia se explica esto), observemos que los números de nuestra tabla concuerdan con los mostrados en la tabla.

gender_tab <- table(StatinGROUP, gender)
gender_tab
##            gender
## StatinGROUP  1  2
##           1 34 29
##           2 26 40

Posteriormente usamos chisq.test() como ya lo hicimos en el video de prueba de chi cuadrado y prueba exacta de Fisher. En el artículo no menciona si se aplicó corrección de continuidad a la prueba, así que probablemente se usó la prueba simple, esto hay que especificarlo con el argumento correct = FALSE. Podemos ver que obtenemos el mismo resultado que el mostrado en la tabla 1.

chisq.test(gender_tab, correct = FALSE)
## 
##  Pearson's Chi-squared test
## 
## data:  gender_tab
## X-squared = 2.752, df = 1, p-value = 0.09713

Ahora, ya yendo a la tabla 3, en la que se muestra información de la reducción de los niveles de diferentes fracciones de lípidos en los grupos de atorvastatina y rosuvastatina. Podemos obtener esta información de una forma similar a lo hecho anteriormente. La diferencia está en que ahora tenemos que restar el nivel de colesterol inicial CHOL1 con el nivel de colesterol final CHOL2.

# obtenemos la diferencia entre la primera medida y la segunda en ambos grupos
chol_red_ator <- CHOL1[StatinGROUP == 1] - CHOL2[StatinGROUP == 1] #atorvastatina
chol_red_rosu <- CHOL1[StatinGROUP == 2] - CHOL2[StatinGROUP == 2] #rosuvastatina
# obtenemos el promedio y la desviación estándar
mean(chol_red_ator)
## [1] 0.7704762
sd(chol_red_ator)
## [1] 0.8513849
mean(chol_red_rosu)
## [1] 1.095152
sd(chol_red_rosu)
## [1] 1.104344

Vemos que los valores obtenidos son similares a los del artículo, ahora solo faltaría obtener el valor de p, para esto lo hacemos también de forma muy similar, realizamos el t-test con las diferencias entre los grupos de fármacos.

t.test(chol_red_ator, chol_red_rosu, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  chol_red_ator and chol_red_rosu
## t = -1.8639, df = 127, p-value = 0.06465
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.66937542  0.02002477
## sample estimates:
## mean of x mean of y 
## 0.7704762 1.0951515

Observemos que obtenemos el mismo valor que en el artículo en el que, para este caso, no se encuentra diferencias significativas. Podemos hacer lo mismo para las otras fracciones lipídicas.

Este ejemplo sirvió para motivar al aprendizaje, para mostrar cómo trabajar con variables continuas y también para introducir al t-test. Con el fin de no hacer este video muy largo o aburrido no profundicé en algunos detalles o conceptos, como la desviación estándar, sobre el t-test, el uso de los corchetes y de operadores (==), o cómo realizar estas pruebas sin usar el programa; podrá ser de utilidad complementar lo aprendido con un poco de teoría si no está muy familiarizado con estos temas.

Manejo de información y correlación (ejemplo COVID-19)

En este video practicaremos un poco más el manejo de la información y de paso explicaré en qué consiste la correlación y cómo calcularla, haré esto usando información epidemiológica y demográfica sobre el COVID-19, que es un tema en verdad muy interesante y quizá haya interesados en trabajar con este tipo de información.

Hasta ahora, si lo pensamos, solo hemos desarrollado diferencias entre grupos de estudio. Ahora examinaremos la correlación, esto nos permite probar si los cambios en una variable guardan una relación con cambios en otra variable en un mismo grupo de estudio, para que esto tenga sentido, ambas variables deben ser variables continuas. Esto se entenderá mejor con un ejemplo.

Para este ejemplo, como mencionamos, usaremos información sobre el COVID-19 de los países, esta información podemos obtenerla del repositorio sobre este tema de Our World in Data. El data set se encuentra en el repositorio guardado como “owid-covid-data.csv”, podemos descargarlo y colocarlo en la carpeta en la que estemos trabajando. Otra opción, probablemente más sencilla, es utilizar el siguiente código:

url="https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv"
filename <- "owid-covid-data.csv"
download.file(url, filename)

Ahora podemos leer este documento de la misma manera que videos anteriores, con la función read.csv(), usamos attach() para acceder a las columnas de forma más sencilla y observamos los nombres de las columnas con colnames(), también podemos visualizar el data set completo con la función View().

dat <- read.csv("owid-covid-data.csv")
attach(dat)
colnames(dat)
##  [1] "iso_code"                        "continent"                      
##  [3] "location"                        "date"                           
##  [5] "total_cases"                     "new_cases"                      
##  [7] "total_deaths"                    "new_deaths"                     
##  [9] "total_cases_per_million"         "new_cases_per_million"          
## [11] "total_deaths_per_million"        "new_deaths_per_million"         
## [13] "new_tests"                       "total_tests"                    
## [15] "total_tests_per_thousand"        "new_tests_per_thousand"         
## [17] "new_tests_smoothed"              "new_tests_smoothed_per_thousand"
## [19] "tests_units"                     "stringency_index"               
## [21] "population"                      "population_density"             
## [23] "median_age"                      "aged_65_older"                  
## [25] "aged_70_older"                   "gdp_per_capita"                 
## [27] "extreme_poverty"                 "cardiovasc_death_rate"          
## [29] "diabetes_prevalence"             "female_smokers"                 
## [31] "male_smokers"                    "handwashing_facilities"         
## [33] "hospital_beds_per_thousand"      "life_expectancy"
View(dat)

Observemos que tiene muchas variables (columnas) sobre información epidemiológica y demográfica en los distintos países, y esta está ordenada de forma cronológica. Notemos también que este es un data set muy grande, más grande que los que usamos anteriormente, a la fecha en la que descargué este tiene estas dimensiones (se actualiza todos los días, así que cuando lo descarguen será incluso más grande):

dim(dat)
## [1] 33627    34

Hay muchas maneras de trabajar con toda esta información, los que mantienen este repositorio usan esta información para hacer sus gráficas, pueden revisar esto en su página web sobre el COVID-19: Coronavirus Pandemic (COVID-19).

En este video usaremos solo usaremos una pequeña parte de la información: las columnas de total_cases y total_deaths, y solo las filas que corresponden a la información del mundo en general, sin separación de países (esto se encuentra guardado como world en la columna location). Lo que analizaremos es la correlación del número de casos de COVID-19 con el número de muertes por COVID-19. Podemos hacer esto de la siguiente manera:

cases <- total_cases[location == "World"]
deaths <- total_deaths[location == "World"]

Ahora, podemos crear un scatter plot con los casos en el eje ‘x’ y las muertes en el eje ‘y’ usando plot(). Observamos que el resultado es el esperado, a medida que el número de casos por COVID-19 aumenta, el número de muertes por COVID-19 aumenta también.

plot(cases, deaths)

La correlación puede notarse a simple vista, pero necesitamos de una forma de cuantificar esta correlación, para esto usamos una prueba que nos permita esto. Para esto podemos usar la función cor.test(), con la que podremos aplicar los métodos estadísticos para hallar correlación más comúnmente usados (Pearson, Spearman, Kendall).

cor.test(cases, deaths, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  cases and deaths
## t = 65.741, df = 211, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9692447 0.9819795
## sample estimates:
##       cor 
## 0.9764479
cor.test(cases, deaths, method = "spearman")
## Warning in cor.test.default(cases, deaths, method = "spearman"): Cannot compute
## exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  cases and deaths
## S = 165.52, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9998972

No explicaré los detalles de cuándo usar un método u otro, pero cabe decir que el método más usado es el de Pearson, y que el de Spearman se usa cuando se tiene outliers (no afectan tanto la correlación). Más detalles sobre la correlación pueden ser revisados y también pueden consultar la manera de cómo realizarlos en otros programas o incluso a mano.

Aprovecho también para mencionar que los datos que son obtenidos fácilmente, como la información epidemiológica del COVID-19 mostrada, suelen ser más difíciles de analizar, hay que tener en cuenta que los países tienen diferentes características demográficas, cada país está en una etapa de la pandemia diferente, la prevalencia de enfermedades es distinta, la prevalencia de otras enfermedades también es distinta, realizan las pruebas de diagnóstico en diferentes proporciones, etc. Por otro lado, la información que se obtiene de forma más “difícil” (ej. estudios experimentales) puede ser analizada de una forma más sencilla, como se pudo observar en el ejemplo del video sobre análisis de un estudio científico y manejo de variables continuas, ya que podemos controlar las variables que harían más difícil el análisis de los datos obtenidos.

Poder estadístico y cálculo del tamaño de muestra

Ya nos acercamos al final de la serie y esta vez quiero desarrollar un tema muy importante el cual no suele profundizarse mucho, especialmente la parte práctica. El cálculo del tamaño de muestra y el poder de un estudio se realiza antes de comenzar el estudio, por lo que tenemos que tenerlo presente desde que planeamos hacer el estudio científico. Este es el último tema y lo dejé para el final porque en realidad no es terriblemente emocionante, pero es importante.

Sin entrar mucho en teoría cabe decir que el valor de p y los intervalos de confianza que desarrollamos en anteriores videos giran alrededor de lo que se conoce como error tipo 1 o \(\alpha\) (concluir falsamente que nuestros resultados son significativos), en el ámbito clínico y otros más solemos tolerar hasta un probabilidad de 5% de cometer este tipo de error, es por esto que usamos al p-value < 0.05 y a los intervalos de confianza del 95%. Sin embargo, existe otro tipo de error, que es el error tipo 2 o \(\beta\) (concluir que nuestros resultados no son significativos cuando sí lo son), para este tipo de error se suele permitir hasta un 20% de probabilidades de cometer este error.

El poder estadístico se relaciona al error de tipo 2, de hecho se define como \(1-\beta\) y este poder estadístico se incrementa mientras aumentamos el tamaño de muestra (de forma no lineal). Pueden profundizar estos conceptos por su cuenta, ahora mejor lo explicaré con un ejemplo.

Para calcular el tamaño de muestra podemos hacerlo de diferentes formas, las más sencilla son el uso de programas que permitan obtenerlas, para esto podemos usar programas como PS: Power and Sample Size Calculation, que es el que usaremos. Usaremos R para obtener los valores mencionados.

Para este ejemplo usaremos el artículo que vimos en el video sobre el análisis de artículo científico publicado y manejo de variables continuas. El artículo es el siguiente: Comparison of Low-Dose Rosuvastatin with Atorvastatin in Lipid-Lowering Efficacy and Safety in a High-Risk Pakistani Cohort: An Open-Label Randomized Trial y el data set se encuentra en: https://doi.org/10.7910/DVN/U8NLQC . Podemos usar el siguiente código para descargar el data set (si aún no lo hemos hecho):

url <- "https://dataverse.harvard.edu/api/access/datafile/3992518?gbrecs=false"
filename <- "statins_comparison_raw.tab"
download.file(url, filename)

Y leemos el archivo como lo hemos estado haciendo.

dat <- read.table("statins_comparison_raw.tab", header = TRUE, sep = "\t")
attach(dat)
View(dat)

Recordemos que obtuvimos los valores del artículo que corresponden a la reducción absoluta en los niveles de colesterol sérico que se muestran en la tabla 3. Lo hacemos de nuevo.

# diferencias entre la primera medida y la segunda en ambos grupos
chol_red_ator <- CHOL1[StatinGROUP == 1] - CHOL2[StatinGROUP == 1] #atorvastatina
chol_red_rosu <- CHOL1[StatinGROUP == 2] - CHOL2[StatinGROUP == 2] #rosuvastatina
# obtenemos el promedio y la desviación estándar
mean(chol_red_ator)
## [1] 0.7704762
sd(chol_red_ator)
## [1] 0.8513849
mean(chol_red_rosu)
## [1] 1.095152
sd(chol_red_rosu)
## [1] 1.104344
# valor de p
t.test(chol_red_ator, chol_red_rosu, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  chol_red_ator and chol_red_rosu
## t = -1.8639, df = 127, p-value = 0.06465
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.66937542  0.02002477
## sample estimates:
## mean of x mean of y 
## 0.7704762 1.0951515

Ahora, imaginemos que queremos realizar un estudio similar a este resultado, en el que queramos encontrar si existe una diferencia en la reducción de los niveles de colesterol entre la atorvastatina y la rosuvastatina, para esto tenemos que saber cuántas personas tenemos que reclutar para tener un estudio con un determinado poder estadístico. Para esto necesitamos los siguientes datos (el programa te ayuda):

Comenzamos determinando la prueba que usamos, en el caso del artículo se usó la prueba del t-test para la diferencia en la reducción de los niveles de colesterol entre la atorvastatina y la rosuvastatina, después los otros valores, menos el último son sencillos de obtener:

alfa <- 0.05
power <- 0.8
m <- 1
delta <- abs(mean(chol_red_ator) - mean(chol_red_rosu))

Para la desviación estándar, tenemos que calcularla con una fórmula, la siguiente:

\[ pooled \ SD = \sqrt{\frac{SD_1^2 + SD_2^2}{2}} \]

Así que simplemente la aplicamos.

SD1 <- sd(chol_red_ator)
SD2 <- sd(chol_red_rosu)
SD_pooled <- sqrt(abs((SD1^2 + SD2^2) / 2))

Ahora ya tenemos todos los valores, así que podemos ingresarlos al programa.

c(alfa=alfa, power=power, delta=delta, sigma=SD_pooled, m=m)
##      alfa     power     delta     sigma         m 
## 0.0500000 0.8000000 0.3246753 0.9860102 1.0000000

Observamos que para nuestro estudio necesitamos de alrededor de 146 personas en cada grupo para poder tener una potencia del 80%, en el artículo los grupos eran de los siguientes tamaños:

nrow(dat[dat$StatinGROUP==1,])
## [1] 63
nrow(dat[dat$StatinGROUP==2,])
## [1] 66

Vemos que en el estudio usaron muchos menos de los necesarios según nuestros cálculos realizados en ese mismo estudio, estudios previos brindan información valiosa. Un muy buen ejercicio para comprender es el obtener los tamaños de muestra necesarios para realizar estudios sobre la diferencia entre el uso atorvastatina y rosuvastatina de los valores de los triglicéridos, HDL y LDL. Para esto se puede volver a ver este video y el video sobre el análisis del artículo científico publicado y manejo de variables continuas.

Esta es la forma de hacerlo en R (debemos implementarnos con el paquete pwr), observemos que obtenemos el mismo valor.

# install.packages(pwr)
library(pwr)
# help(package = pwr) # elegimos la función
# help(pwr.t.test) # vemos qué debemos incluir en la función
# NOTA: necesitamos 3 de 4 valores en la función, para calcular el faltante
# effect size (d)
M1 <- mean(chol_red_ator)
M2 <- mean(chol_red_rosu)
SD1 <- sd(chol_red_ator)
SD2 <- sd(chol_red_rosu)
SD_pooled <- sqrt(abs((SD1^2 + SD2^2) / 2))
d <- abs(M1 - M2) / SD_pooled
# otros valores
sig.level <- 0.05
power <- 0.8
# usamos la función
pwr.t.test(d=d, sig.level=sig.level, power=power)
## 
##      Two-sample t test power calculation 
## 
##               n = 145.744
##               d = 0.3292819
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

También podemos replicar en R la curva de poder y tamaño de muestra, y además tenemos la opción de personalizarla a nuestro antojo.

pts <- 15
pwr <- 1:pts/pts
size <- c()
for (i in 1:pts) {
  j <- i/pts
  size[i] <- pwr.t.test(d=d, sig.level=sig.level, power=j)$n
}
plot(
  x=size, y=pwr, type="l", col="red", lwd="2", xlim=c(0,220),
  xlab="Tamaño de muestra", ylab="Poder",
  main="Curva de Poder vs. Tamaño de muestra"
)
grid()

Aquí dejo algunos recursos donde se explica como realizar estos cálculos en R con mayor profundidad:

Bien, ya habiendo explicado los cálculos es importante tener en cuenta de que realizar estos cálculos depende de los requerimientos de la revista o el tipo de investigación, sin embargo, es importante saber que contar con ellos da más validez a un trabajo de investigación, de hecho, para aplicar a subvenciones otorgadas por gobiernos u organizaciones se exige contar con estos cálculos en la aplicación, por ejemplo, múltiples guías sobre el reporte de estudios científicos tienen en sus checklists al cálculo del tamaño de muestra, como en la checklist de la declaración CONSORT.

Cómo continuar aprendiendo

Hay muchas formas de continuar aprendiendo, pero creo que la mejor forma es aplicando lo que vamos aprendiendo, así que mientras intentamos investigar descubrimos qué es lo que nos falta y podemos aprender de eso. En lo que he podido aprender de todo esto hasta ahora, podría hacer una lista de lo que podemos aprender (pero repito que no es tan bueno aprender sin aplicar lo aprendido):

Referencias