Estadística con R

author: Javier Ballesteros, UPV/EHU date: Agosto - Septiembre 2015

Objetivos generales del curso

Introducir R como software estadístico
Presentar las pruebas de estadística clásica con R

Requisitos para optimizar un curso de bioestadística

Tener libre acceso a un ordenador, y a la aplicación estadística concreta
Tener, al menos, un problema real que resolver con dicha aplicación
Poseer un manual de la aplicación estadística
Romper la inercia nuevo aprendizaje - nuevo programa de ordenador, con un tercero que normalmente es el instructor
Contar con alguien a quien acudir en los momentos iniciales de duda

¿Qué es lo que no va a desarrollar este curso?

Probabilidad y distribuciones
distribuciones de probabilidad discretas o continuas
La estadística matemática que subyace a los tests de significación y a la inferencia
teorema central del límite
Contenidos especiales y avanzados:
análisis de supervivencia, psicometría
modelos lineales mixtos, estadística multivariada

Estadística con R: contenidos

Estadística descriptiva: tipos de variables y su visualización
Inferencia estadística: tests de significación, p-valores, intervalos de confianza, y tamaños del efecto
Pruebas estadísticas clásicas:
comparación de medias entre 2 grupos (t-tests)
comparación de medias entre más de 2 grupos (ANOVA)
Correlación
Regresión lineal
Análisis de frecuencias y proporciones
Librería lessR

¿Para qué sirven los conocimientos de bioestadística? Bradford Hill: Principles of Medical Statistics, 1937

Statistics are curious things…(They) tend to induce a strong emotional reaction in non-mathematical minds…It is exasperating, when we have studied a problem by methods that we have spent laborious years in mastering, to find our conclusions questioned, and perhaps refuted…There seems to be to me only one way of escape. The worker in medical problems, in the field of clinical as well as preventive medicine, must himself know something of statistical technique, both in experimental arrangements and in the interpretation of figures

Algunas ideas sobre la (bio)estadística

Es aburrida
Hace falta un profundo conocimiento matemático
Es como el bikini (enseña mucho, pero tapa lo más importante)
Es como una farola (apoyo y/o luz)
Es imprescindible
Ofrece información sobre la que basar decisiones
Es una herramienta más en la investigación y/o formación continuada

El proceso de investigación

Dominio teórico
Hipótesis
Operacionalización / medición de conceptos
Selección de sujetos
Diseño del estudio (observacional / experimental)
Recogida de datos
Análisis de resultados (univariados / bivariados / multivariados)
Presentación y diseminación de resultados

Tipos de datos / variables en cuanto a su medida

Cuantitativos: proporcionan medias
Continuos: toman valores continuos en un rango (altura, presión arterial)
Discretos: toman valores enteros en un rango (número de hijos, número de ataques de asma)
Categóricos: proporcionan porcentajes o proporciones
Binarios: presencia o no del acontecimiento en estudio (0/1)
Nominales: categorías no ordenadas (tipo sanguíneo, A, B, AB, O)
Ordinales: categorías ordenadas (estadíos cáncer, puntuaciones Likert)
Tiempo hasta acontecimiento : variable de supervivencia

Tipos de datos / variables en cuanto a su papel en el diseño

Variables de resultado
Variables explicativas, mediadoras, confusoras…

Estadística descriptiva: funciones generales

En toda la presentación se asume que las variables a analizar son vectores, o están incluidas en una base de datos rectangular (individuos por variables)

class(): tipo de variable
names(): devuelve los nombres de las variables de la base de datos
summary(): estadísticos de resumen de las variables de la base de datos
head(): visualización de datos, 6 primeros registros
tail(): visualización de datos, 6 últimos registros
which(): con condiciones lógicas localiza observaciones específicas

Estadística descriptiva: nombres y clases de las variables

data(iris)
names(iris)

[1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
[5] "Species"

class(iris$Species); class(iris$Sepal.Length) # estructura "datos$variable"

[1] "factor"

[1] "numeric"

Estadística descriptiva: visualización de los primeros registros

head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

Estadística descriptiva: visualización de los últimos registros

tail(iris)

    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
145          6.7         3.3          5.7         2.5 virginica
146          6.7         3.0          5.2         2.3 virginica
147          6.3         2.5          5.0         1.9 virginica
148          6.5         3.0          5.2         2.0 virginica
149          6.2         3.4          5.4         2.3 virginica
150          5.9         3.0          5.1         1.8 virginica

Estadística descriptiva: resumen de los estadísticos para cada variable

summary(iris)

  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50

Estadística descriptiva: identificación de observaciones

data(mtcars)
which(mtcars$mpg > median(mtcars$mpg)) # identificación de valores extremos

 [1]  1  2  3  4  8  9 18 19 20 21 26 27 28 30 32

Visualización gráfica de los datos

Objetivos:

¿Hay valores extremos aberrantes?
¿Hay valores sin sentido?
¿Cómo se distribuyen los datos?
Datos categóricos: frecuencias y porcentajes en cada categoría
Datos numéricos: forma de la distribución, tendencia central, y dispersión

Gráficos descriptivos: variables continuas

“Las reglas más importantes en estadística son utilizar el sentido común, y visualizar gráficamente los datos y asociaciones de interés”

stem(): distribución de los valores de la variable
hist(): histograma de los valores de la variable
boxplot(): gráfico de caja de la variable
qqnorm(): ajuste de la distribución de la variable a la normal teórica
qqline(): ajuste de la distribución de la variable a la normal teórica

Gráficos para variables continuas: función stem()

stem(mtcars$mpg) # distribución valores observados


  The decimal point is at the |

  10 | 44
  12 | 3
  14 | 3702258
  16 | 438
  18 | 17227
  20 | 00445
  22 | 88
  24 | 4
  26 | 03
  28 | 
  30 | 44
  32 | 49

Gráficos para variables continuas: función hist()

hist(mtcars$mpg) # forma de la distribución, tendencia central, y variabilidad

plot of chunk unnamed-chunk-7

Gráficos para variables continuas: función boxplot()

boxplot(mtcars$mpg) # mediana, IQR, máximos-mínimos según IQR, y observaciones extremas

plot of chunk unnamed-chunk-8

Gráficos para variables continuas: funciones qqnorm() y qqline()

qqnorm(mtcars$mpg); qqline(mtcars$mpg) # ajuste a distribución normal

plot of chunk unnamed-chunk-9

Descripción de variables continuas

summary(): estadístico de resumen de la variable (mínimo, primer cuartil, mediana, media, tercer cuartil, máximo)
fivenum(): estadísticos de resumen de la variable - Tukey (mínimo, “hinge” bajo, mediana, “hinge” alto, máximo)
mean(): media
median(): mediana
sd(): desviación estándar
var(): varianza
mad(): mediana de la desviación absoluta
IQR(): desviación intercuartil (Q75 - Q25), 50% de los datos

Descripción de variables continuas: funciones summary(), y fivenum()

summary(mtcars$mpg)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10.40   15.42   19.20   20.09   22.80   33.90

fivenum(mtcars$mpg)

[1] 10.40 15.35 19.20 22.80 33.90

Descripción de variables continuas: funciones mean(), y median()

mean(mtcars$mpg)

[1] 20.09062

median(mtcars$mpg)

[1] 19.2

Características de la media y la mediana

Ambas son medidas de la tendencia central de la distribución de una variable cuantitativa. La media se afecta por valores extremos, no así la mediana

\[\large{\bar {X} = \frac {\sum_{i=1}^n {x_i}}{n}}\]

Mientras que la \(\bar x\), está afectada por valores extremos de \(x_i\), la mediana al ser el punto medio de los valores ordenados por rangos no se afecta

Descripción de variables continuas: funciones sd(), y var()

var(mtcars$mpg); sd(mtcars$mpg)

[1] 36.3241

[1] 6.026948

\[\large{S^2 = \frac{\sum_{i=1}^n (x_i - \bar{X})^2}{n-1}}\]

\[\large{S = \sqrt{S^2}}\]

Características de la desviación estándar

La \(S^2\) presenta unidades al cuadrado lo que no es fácilmente interpretable. la \(S\) está en las mismas unidades que la \(\bar {X}\) lo que facilita su interpretación

La \(S\) al igual que la \(\bar {X}\) se afecta por valores extremos. Debido a la elevación al cuadrado de las distancias entre cada observación y la media de la distribución, los valores más alejados de la media contribuyen más a la desviación estándar que los valores cercanos a la media

Descripción de variables continuas: funciones mad() e IQR()

mad(mtcars$mpg)

[1] 5.41149

IQR(mtcars$mpg)

[1] 7.375

Características de los estadísticos de dispersión de datos basados en percentiles

Se basan en el ordenamiento de datos (de menor a mayor) y al igual que la mediana (percentil 50) no se afectan por valores extremos

El rango intercuartil (IQR) recoge los valores del primer al tercer cuartil de la distribución, el 50% de los datos (el box del gráfico boxplot), y no está afectado por valores extremos o outliers

Distribución normal en variables numéricas

Si una variable continua sigue una distribución normal, el 68% de los datos estarán comprendidos entre la media y 1 DE, el 95% entre la media y 2 DE, y el 99.7% entre la media y 3 DE

Un valor que presente más de 3 DE con respecto a la media puede ser un valor aberrante

Gráficos para variables categóricas: función genérica plot()

group <- rep(1:2, each=10)
groupf <- factor(group)
plot(groupf)

plot of chunk unnamed-chunk-14

Gráficos para variables categóricas: funciones table() y barplot()

table(groupf) # frecuencia de observaciones por grupo

groupf
 1  2 
10 10

Gráficos para variables categóricas: funciones table() y barplot()

barplot(table(groupf))

plot of chunk unnamed-chunk-16

Inferencia estadística

Inducción desde los valores de la muestra a los parámetros de la población. Tiene 2 posibles objetivos:

Test de hipótesis (p-valores): ¿hay algún efecto?
Estimación del efecto - intervalos de confianza: rango de valores plausibles para el efecto verdadero

Símbolos a utilizar

\(S^2\): varianza muestral
\(S\): desviación estándar muestral
\(\sigma^2\): varianza poblacional
\(\sigma\): desviación estándar poblacional
\(\bar{X}\): media muestral
\(\mu\): media poblacional
\(IQR\): rango intercuartil (50% medio de los valores de la distribución)

Desde la muestra a la población: SD a SE

Error estándar de la media: \(SE_{\bar{X}} = \frac {S}{\sqrt{n}}\)
Error estándar del coeficiente de correlación: \(SE_r = \sqrt{\frac{1-r^2}{n}}\)
Error estándar de una proporción: \(SE_p = \sqrt{\frac{p(1-p)}{n}}\)

El \(SE\) depende del estadístico, y disminuye al aumentar el tamaño muestral, y aumenta cuando se incrementa la variabilidad

Desde la muestra a la población: SD a SE

Error estándar de la diferencia de riesgos: \(SE_{RD} = \sqrt{\frac{p_1(1-p_1)}{n_1} + \frac{p_0(1-p_0)}{n_0}}\)
Error estándar del riesgo relativo: \(SE_{\ln{RR}} = \sqrt{\frac{1}{a} - \frac{1}{n_1} + \frac{1}{c} - \frac{1}{n_0}}\)
Error estándar de la odds-ratio: \(SE_{\ln{OR}} = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}}\)

El \(SE\) depende del estadístico, y disminuye al aumentar el tamaño muestral, y aumenta cuando se incrementa la variabilidad

Distribuciones de estadísticos ~ teorema central del límite

Distribución \(t \simeq {normal}\) en muestras grandes, \(n > 100\)
Medias / diferencia de medias
Coeficientes de correlación
Distribución normal
Proporciones / diferencia de proporciones, \(RD\)
\(\ln{RR}\), \(\ln{OR}\)

Distribución normal

\[\large{f(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}{e^{-{\frac{(x -\mu)^2}{2\sigma^2}}}}}\]

La distribución de los valores observados \(x\), depende sólo de la media de su distribución teórica \(\mu\), y de su desviación estándar \(\sigma\)

La ecuación contiene 3 clásicos números irracionales:

\(\sqrt 2 = 1.4142 \cdots\)
\(\pi = 3.1415 \cdots\)
\(e = 2.7182 \cdots\)

Tests de hipótesis

El p-valor informa sobre la probabilidad de que el efecto observado (u otro más extremo) pudiera haber aparecido por azar si la hipótesis nula fuera cierta

\(\large{p = \Pr{(datos|H_0)}}\)

Pasos en un test de hipótesis

Definir las hipótesis: \(H_0\), \(H_A\)
Especificar la distribución de \(H_0\) \[H_0 \sim{N(\mu, \sigma)}\]
Realizar el experimento para obtener datos
Calcular el p-valor observado
Rechazar o no \(H_0\)

Intervalos de confianza (IC)

Pretenden capturar el valor real del efecto en teóricas replicaciones
Un IC al 95% debería incluir el valor real del efecto al menos en 95 de 100 veces que se replicara el experimento
Un IC al 99% debería incluir el valor real del efecto al menos en 99 de 1000 veces que se replicara el experimento

IC asumiendo una distribución normal

IC 90%: \(statistic \pm{1.64} \times{SE_{statistic}}\)
IC 95%: \(statistic \pm{1.96} \times{SE_{statistic}}\)
IC 99%: \(statistic \pm{2.58} \times{SE_{statistic}}\)

Intervalos de confianza (IC) proporcionan:

Un rango plausible de valores para el parámetro poblacional que se está estimando
La precisión de la estimación según el IC sea amplio o estrecho
La significación estadística. Si un IC 95% no incluye el valor nulo, es significativo a nivel \(p < 0.05\)

Corespondencia entre p-valores e IC

\[P_{valor} = Z = \frac{Efecto - H_0}{SE}\]

\[IC_{1-\alpha{2}} = Efecto \mp{Z_{\alpha{2}}} \times({SE})\]

Los p-valores dependen del tamaño del efecto (señal), del tamaño muestral, y de la variabilidad (ruido)
¿Qué significa realmente que \(p < 0.05\) o que \(p \geq{0.05}\)?

Factores que afectan la significación estadística de un contraste de hipótesis

\[Significance \propto{\frac{Effect \times{Sample}}{Variability}}\]

Los p-valores distinguen entre efectos reales y la fluctuación aleatoria. Si la muestra es suficientemente grande, la fluctuación aleatoria no es un problema, y los p-valores son irrelevantes
Si el tamaño muestral es muy grande, deben ignorarse los p-valores
Prestar atención al tamaño del efecto y a su IC

Comparar los tamaños del efecto, no los p-valores

Expresar la reducción o incremento del riesgo en base a:
Los números observados
La reducción del riesgo absoluto: \(RD\)
El número de sujetos necesarios a tratar para beneficio o perjuicio: \[NNT = RD^{-1}\] \[NNH = RD^{-1}\]
Reducción o incremento en el \(RR\)

¿Qué es un valor estadísticamente significativo?

Karl Pearson
Razón señal-ruido: regla del 3 (según aparece en Biometrika)

Tome 3 DE como definitivamente significativo

Ronald Aylmar Fisher
p-valores

¿De donde procede que p < 0.05 sea estadísticamente significativo?

RA Fisher en su libro Statistical methods for research workers, con una primera edición en 1925, en el capítulo V titulado Tests of significance of means, differences of means, and regression coefficients, dice:

If the difference is many times greater than the standard error, it is certainly significant, and it is a convenient convention to take twice the standard error as the limit of significance; this is roughly equivalent to the corresponding limit \(\small{P = .05}\), already used for the \(\small{\chi^2}\) distribution

Más Fisher y p-valores

RA Fisher (1925): SMRW, sobre la tabla de la distribución \(\small{\chi^2}\):

In preparing this table we have borne in mind that in practice we do not want to know the exact value of \(\small{P}\) for any observed \(\small{\chi^2}\), but in the first place, whether or not the observed value is open to suspicion. If \(\small{P}\) is between 0.1 and 0.9 there is certainly no reason to suspect the hypothesis tested. If it is below 0.02 it is strongly indicated that the hypothesis fails to account for the whole of the facts. We shall not often be astray if we draw a conventional line at 0.05 and consider that higher values of \(\small{\chi^2}\) indicate a real discrepancy

Más Fisher y p-valores

RA Fisher (1929). The statistical method in psychical research. Proceedings of the Society for Psychical Research, xxxix:189-92

In the investigation of living beings by biological methods, statistical tests of significance are essential. Their function is to prevent us being deceived by accidental occurrences…It is a common practice to judge a result significant, if it is of such magnitude that it would have been produced by chance not more frequently than once in twenty trials

Interpretación usual de p-valores

Cuanto más pequeño sea el p-valor menos probable es que el valor observado en el experimento sea generado por la distribución de la hipótesis nula. Es un indicador de la sorpresa del resultado versus \(H_0\)

En base a su valor se rechaza o no se rechaza \(H_0\) en favor de \(H_1\)

\(p \simeq 0.1\), hay cierta evidencia de rechazo de \(H_0\)
\(p \simeq 0.05\), hay una evidencia modesta de rechazo de \(H_0\)
\(p \simeq 0.01\), hay una evidencia importante para rechazar \(H_0\)

Problemas con los p-valores: falacias en su interpretación

Falacia de la probabilidad inversa

Falsa creencia de que el p-valor indica la probabilidad de que \(\small{H_0}\) sea cierta dados los datos \(\small{Pr(H_0 | datos)}\)

Falacia del tamaño del efecto

Falsa creencia de que el p-valor puede interpretarse como un tamaño del efecto (cuanto más pequeño más importante el efecto)

Falacia de la significación clínica o práctica

Falsa creencia de que una significación estadística (habitualmente \(\small{p<{0.05}}\)), conlleva necesariamente una significación práctica

Problemas con los p-valores: significación vs relevancia

Significación versus relevancia del hallazgo

¿Hay algún efecto?: p-valor
¿Cuál es la magnitud del efecto?: tamaño del efecto, y en ciencia el tamaño importa
¿es importante el efecto?: aspectos de juicio personal (o social) en la toma de decisiones

Significación estadística versus significación clínica

plot of chunk unnamed-chunk-17

En resumen

Significación estadística: p-valor que dicotomiza unos resultados
Significación clínica: tamaño del efecto y su precisión (intervalo de confianza)

\[SMD=\frac{\widehat x_1 - \widehat x_2}{SD_{pooled}}\]

2 escuelas de pensamiento

Fisher-Neyman-Pearson: escuela frecuentista
NHST, intervalo de confianza, error tipo I (\(\alpha\)), error tipo II (\(\beta\)), poder del estudio (\(1-\beta\))

\[\large{P(datos|H_0)}\]

Bayes-Price-Laplace: escuela bayesiana
probabilidad subjetiva, distribución a priori, distribución posterior, intervalo de credibilidad

\[\large{P(H_0|datos)}\]

Falacia de la simetría condicional en la interpretación de los p-valores

Por supuesto no es cierto que:

\[\large{P(datos|H_0) \ne P(H_0|datos)}\]

Lo que es curioso ya que bastantes veces se interpreta el p-valor bajo un punto de vista bayesiano en vex de frecuentista (falacia de la probabilidad inversa)

Estadística clásica: asociación entre variables

Cuando la variable resultado es continua y la variable explicativa categórica: comparación de medias y varianzas entre grupos
Comparación de medias para datos emparejados (antes-después)
Comparación de medias en 2 grupos independientes
Comparación de medias en > 2 grupos independientes
Asociación entre 2 variables continuas: correlación
Cuando la variable resultado y la variable explicativa son continuas: regresión lineal
Asociación entre 2 variables categóricas: análisis de frecuencias y proporciones

Preliminares: ¿son mis datos “normales”?

Algunos de los tests estadísticos que vamos a ver en este curso asumen distribuciones normales, y esto es especialmente importante en el caso de muestras pequeñas.

No todos los datos continuos se distribuyen normalmente!

Evaluación de la distribución

¿Presenta el histograma forma de campana?: stem(); hist()
¿Es aproximadamente lineal el gráfico de probabilidad normal?: qqnorm(); qqline()
¿Son la media y la mediana aproximadamente iguales?: summary(); mean(); median()
¿Están 2/3 de las observaciones comprendidas entre la media y 1 DE?: summary()
¿Están el 95% de las observaciones comprendidas entre la media y 2 DE?: summary()
Los tests de normalidad están fuertemente influenciados por el tamaño de la muestra: ks.test(), shapiro.test()

Asunciones de los modelos lineales

t-test, ANOVA, r, y regresión lineal
La variable resultado se distribuye normalmente (es importante principalmente en muestras pequeñas, las muestras grandes son robustas a esta asunción por el Teorema Central del Límite)
Homogeneidad de varianzas

Comparación de medias y varianzas - 2 grupos

t.test(): comparación de medias, prueba paramétrica
wilcox.test(): comparación de medias, prueba no paramétrica
var.test(): comparación de varianzas entre 2 muestras
bartlett.test(): comparación de varianzas entre 2 o más muestras

Comparación de medias de datos emparejados (antes - después): datos

?sleep # mejora en horas de sueño tras tratamiento
pre <- c(0.7, -1.6, -0.2, -1.2, -0.1, 3.4, 3.7, 0.8, 0.0, 2.0)
post <- c(1.9, 0.8, 1.1, 0.1, -0.1, 4.4, 5.5, 1.6, 4.6, 3.4)

Comparación de medias de datos emparejados (antes - después): t-test

t.test(pre, post, paired = TRUE)


    Paired t-test

data:  pre and post
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of the differences 
                  -1.58

Comparación de medias de datos emparejados (antes - después): wilcoxon

wilcox.test(pre, post, paired = TRUE)


    Wilcoxon signed rank test with continuity correction

data:  pre and post
V = 0, p-value = 0.009091
alternative hypothesis: true location shift is not equal to 0

Comparación de medias de datos emparejados (antes - después): diferencia frente a 0

Se podría haber analizado de otra forma que reduce el problema a contrastar la diferencia frente a un valor teórico

diff <- pre - post
hist(diff)

plot of chunk unnamed-chunk-21

Comparación de medias de datos emparejados (antes - después): t-test

t.test(diff, mu = 0)


    One Sample t-test

data:  diff
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -2.4598858 -0.7001142
sample estimates:
mean of x 
    -1.58

Comparación de medias de datos emparejados (antes - después): wilcoxon

wilcox.test(diff, mu = 0)


    Wilcoxon signed rank test with continuity correction

data:  diff
V = 0, p-value = 0.009091
alternative hypothesis: true location is not equal to 0

Comparación de medias entre 2 grupos independientes: datos

data(sleep)
boxplot(extra ~ factor(group), data = sleep)

plot of chunk unnamed-chunk-24

Comparación de medias entre 2 grupos independientes: igualdad de varianzas

bartlett.test(extra ~ factor(group), data = sleep) # igualdad de varianzas


    Bartlett test of homogeneity of variances

data:  extra by factor(group)
Bartlett's K-squared = 0.10789, df = 1, p-value = 0.7426

Comparación de medias entre 2 grupos independientes: t-test asumiendo varianzas iguales

t.test(extra ~ factor(group), data = sleep, var.equal = TRUE)


    Two Sample t-test

data:  extra by factor(group)
t = -1.8608, df = 18, p-value = 0.07919
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.363874  0.203874
sample estimates:
mean in group 1 mean in group 2 
           0.75            2.33

Comparación de medias entre 2 grupos independientes: t-test sin asumir igualdad de varianzas

t.test(extra ~ factor(group), data = sleep)


    Welch Two Sample t-test

data:  extra by factor(group)
t = -1.8608, df = 17.776, p-value = 0.07939
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.3654832  0.2054832
sample estimates:
mean in group 1 mean in group 2 
           0.75            2.33

Comparación de medias entre 2 grupos independientes: wilcoxon

wilcox.test(extra ~ factor(group), data = sleep)


    Wilcoxon rank sum test with continuity correction

data:  extra by factor(group)
W = 25.5, p-value = 0.06933
alternative hypothesis: true location shift is not equal to 0

Fórmulas t-test

2 grupos emparejados:

\[\large{t_{n-1 gl} = \frac{\bar{D} - \mu_D}{S_{\bar{D}}}} \]

2 grupos independientes:

\[\large{t_{n_1 + n_2 -1 gl} = \frac{(\bar{X_1} - \bar{X_2}) - (\mu_1 - \mu_2)}{S_{(\bar{X_1}- \bar{X_2})}}} \]

Diferencia de medias entre más de 2 grupos: ANOVA

oneway.test(): análisis de la varianza de 1 vía, paramétrico
kruskal.test(): análisis de la varianza de 1 vía, no paramétrico
pairwise.t.test(): múltiples comparaciones (paramétrico)
pairwise.wilcox.test: múltiples comparaciones (no paramétrico)
aov(): anális de la varianza en general

ANOVA

Una prueba significativa de ANOVA (F-test) nos dice que al menos 2 de las medias difieren pero no nos dice cuáles. Es un omnibus test
Determinar entre qué grupos están las diferencias requieren de otros análisis que ya deben corregir el problema de las comparaciones múltiples

\[\large{F_{n, m} \sim \frac{\sigma^2_{between}}{\sigma^2_{witin}}}\]

ANOVA: datos

# Índice de pulsatilidad en la arteria cerebral media según exposición a cannabis
grupo <- c("control", "control", "control", "control", "control", "control", "activo", "activo", "activo", "activo", "activo", "activo", "ex", "ex", "ex", "ex", "ex", "ex")
index <- c(.9, .9, .9, .8, .8, NA, 1, 1, 1, 1, .9, NA, 1, 1, 1, 1, .9, .9)
pulsatilidad <- data.frame(grupo, index)
pulsatilidad

     grupo index
1  control   0.9
2  control   0.9
3  control   0.9
4  control   0.8
5  control   0.8
6  control    NA
7   activo   1.0
8   activo   1.0
9   activo   1.0
10  activo   1.0
11  activo   0.9
12  activo    NA
13      ex   1.0
14      ex   1.0
15      ex   1.0
16      ex   1.0
17      ex   0.9
18      ex   0.9

ANOVA: boxplot

boxplot(index ~ grupo, data = pulsatilidad)

plot of chunk unnamed-chunk-30

ANOVA: igualdad de varianzas

bartlett.test(index ~ grupo, data = pulsatilidad)


    Bartlett test of homogeneity of variances

data:  index by grupo
Bartlett's K-squared = 0.15376, df = 2, p-value = 0.926

ANOVA: omnibus test

oneway.test(index ~ grupo, data = pulsatilidad)


    One-way analysis of means (not assuming equal variances)

data:  index and grupo
F = 7.4891, num df = 2.0000, denom df = 8.4798, p-value = 0.01338

ANOVA: corrección para múltiples tests

pairwise.t.test(pulsatilidad$index, pulsatilidad$grupo)


    Pairwise comparisons using t tests with pooled SD 

data:  pulsatilidad$index and pulsatilidad$grupo 

        activo control
control 0.0073 -      
ex      0.6708 0.0082 

P value adjustment method: holm

ANOVA: omnibus test con la función aov()

aov.res <- aov(index ~ grupo, data = pulsatilidad) # ANOVA clásico
summary(aov.res)

            Df  Sum Sq  Mean Sq F value  Pr(>F)   
grupo        2 0.04417 0.022083   8.613 0.00415 **
Residuals   13 0.03333 0.002564                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 observations deleted due to missingness

ANOVA: corrección para múltiples tests

TukeyHSD(aov.res)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = index ~ grupo, data = pulsatilidad)

$grupo
                      diff         lwr         upr     p adj
control-activo -0.12000000 -0.20456166 -0.03543834 0.0064223
ex-activo      -0.01333333 -0.09429496  0.06762830 0.9017863
ex-control      0.10666667  0.02570504  0.18762830 0.0105962

ANOVA: gráfico para múltiples tests

plot(TukeyHSD(aov.res))

plot of chunk unnamed-chunk-36

ANOVA: prueba no paramétrica

kruskal.test(index ~ grupo, data = pulsatilidad)


    Kruskal-Wallis rank sum test

data:  index by grupo
Kruskal-Wallis chi-squared = 8.3714, df = 2, p-value = 0.01521

ANOVA no paramétrico: corrección para múltiples tests

pairwise.wilcox.test(index, grupo, data = pulsatilidad)


    Pairwise comparisons using Wilcoxon rank sum test 

data:  index and grupo 

        activo control
control 0.057  -      
ex      0.724  0.057  

P value adjustment method: holm

Correlación: una larga historia

Francis Galton (1822-1911): concepto de correlación entre 2 variables
Karl Pearson (1857-1936): descripción coeficiente correlación \(r\)
Charles Spearman (1863-1945): correlación no paramétrica - al igual que Pearson trabajó en el University College London, pero en el Departamento de Psicología, y tuvo malas relaciones con Pearson
Ronald Fisher (1890-1962): test de significación de \(r\), transformación \(z\) de \(r\), malas relaciones con Pearson
Maurice Kendall (1907-1983): London School of Economics, correlación no paramétrica

Covarianza y correlación

\[cov_{x,y} = \frac{\sum_{i=1} ^n{(x_i - \bar{x})(y_i - \bar{y})}}{n-1}\]

\[Var_x = cov_{x,x} = \frac{\sum_{i=1} ^n{(x_i - \bar{x})(x_i - \bar{x})}}{n-1}\]

\[r = \frac{cov_{x,y}}{S_x * S_y}\]

La covarianza tiene unidades que pueden ser difíciles de interpretar, mientras que la correlación es una covarianza estandarizada, y por lo tanto sin unidades

Correlación

\[\large{\hat r = \frac{cov_{x,y}}{S_x * S_y} = \frac{\frac{\sum_{i=1} ^n (x_i - \bar x)(y_i - \bar y)}{n-1}}{\sqrt{\frac{\sum_{i=1} ^n (x_i - \bar x)^2}{n-1}}{\sqrt{\frac{\sum_{i=1} ^n (y_i - \bar y)^2}{n-1}}}}}\]

\(R^2 = \hat r ^2\) es una medida del ajuste del modelo, y refleja la proporción de variabilidad en la variable resultado que es explicada por la variable o variables explicativas

Correlación: funciones en R

cor(): correlaciones paramétricas y no paramétricas
cor.test(): test de significación de la correlación
cov(): covarianza entre 2 variables

La obtención de correlaciones paramétrica/no paramétrica depende de sus argumentos:

method=“pearson” presenta la correlación paramétrica r de Pearson
method=“spearman” presenta la correlación no paramétrica rho de Spearman
method=“kendall” presenta la correlación no paramétrica tau de Kendall

Propiedades del coeficiente de correlación

El coeficiente de correlación es un número adimensional que presenta un rango entre -1 y +1, e indica la fortaleza de la relación lineal entre 2 variables numéricas. Gráficos de dispersión para valorar esa relación lineal antes de estimar su magnitud mediante los coeficientes de correlación aopropiados:

plot(xvar, yvar)
abline(lm(y ~ x))

Correlación lineal: Anscombe 1

?anscombe # base de datos histórica en **R**
data(anscombe)
cor.test(anscombe$x1, anscombe$y1)


    Pearson's product-moment correlation

data:  anscombe$x1 and anscombe$y1
t = 4.2415, df = 9, p-value = 0.00217
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4243912 0.9506933
sample estimates:
      cor 
0.8164205

Correlación lineal: Anscombe 1

plot(anscombe$x1, anscombe$y1)
abline(lm(y1 ~ x1, data=anscombe))

plot of chunk unnamed-chunk-40

Correlación lineal: Anscombe 2

cor.test(anscombe$x2, anscombe$y2)


    Pearson's product-moment correlation

data:  anscombe$x2 and anscombe$y2
t = 4.2386, df = 9, p-value = 0.002179
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4239389 0.9506402
sample estimates:
      cor 
0.8162365

Correlación lineal: Anscombe 2

plot(anscombe$x2, anscombe$y2)
abline(lm(y2 ~ x2, data=anscombe))

plot of chunk unnamed-chunk-42

Correlación lineal: Anscombe 3

cor.test(anscombe$x3, anscombe$y3)


    Pearson's product-moment correlation

data:  anscombe$x3 and anscombe$y3
t = 4.2394, df = 9, p-value = 0.002176
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4240623 0.9506547
sample estimates:
      cor 
0.8162867

Correlación lineal: Anscombe 3

plot(anscombe$x3, anscombe$y3)
abline(lm(y3 ~ x3, data=anscombe))

plot of chunk unnamed-chunk-44

Correlación lineal: Anscombe 4

cor.test(anscombe$x4, anscombe$y4)


    Pearson's product-moment correlation

data:  anscombe$x4 and anscombe$y4
t = 4.243, df = 9, p-value = 0.002165
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4246394 0.9507224
sample estimates:
      cor 
0.8165214

Correlación lineal: Anscombe 4

plot(anscombe$x4, anscombe$y4)
abline(lm(y4 ~ x4, data=anscombe))

plot of chunk unnamed-chunk-46

Correlación no implica causación: nidos de cigüeñas y nacimientos

births <- c(4, 20, 60, 30, 7, 6, 14, 20, 18) # número de nacimientos
nests <- c(2, 6, 8, 7, 0, 1, 2, 2, 3) # número de nidos
plot(nests, births) # ¿hay asociación?

plot of chunk unnamed-chunk-47

Correlación no implica causación: r de Pearson

cor.test(nests, births) # correlación de Pearson por defecto


    Pearson's product-moment correlation

data:  nests and births
t = 4.1991, df = 7, p-value = 0.00404
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4152802 0.9668960
sample estimates:
      cor 
0.8460611

Correlación no implica causación: rho de Spearman

cor.test(nests, births, method = "spearman")


    Spearman's rank correlation rho

data:  nests and births
S = 20.928, p-value = 0.006122
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8255992

Correlación no implica causación: tau de Kendall

cor.test(nests, births, method = "kendall")


    Kendall's rank correlation tau

data:  nests and births
z = 2.5669, p-value = 0.01026
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.7061879

Correlación no implica causación

La correlación evalúa la asociación entre 2 variables, y esto no significa necesariamente causación. Ambas variables podrían estar influenciadas simultáneamente por una tercera
Sólo si la variable explicativa está controlada por el experimentador, la asociación entre 2 variables puede presentar evidencia de una relación causa-efecto
El coeficiente de correlación es anumérico y puede conducir a falsas interpretaciones sobre la relación entre variables. Es aconsejable presentar además de r, las medias, DEs, y el número de observaciones
Las correlaciones basadas en proporciones o medias (correlaciones ecológicas) pueden conducir a serios errores de interpretación

Interpretación del coeficiente de correlación

\(r^2\) es la proporción estimada de la varianza en \(y\) que puede ser explicada por su regresión lineal en \(x\)
\(1-r^2\) es la proporción de varianza no explicada

\(r\): 0.50 - 0.70 - 0.90

\(r^2\): 0.25 - 0.49 - 0.81

\(1-r^2\): 0.75 - 0.51 - 0.19

Guía general valores absolutos correlación

\(\hat r > 0.70\) correlación importante
\(\hat r > 0.30, \hat r < 0.70\) correlación moderada
\(\hat r > 0.10, \hat r < 0.30\) correlación débil
\(\hat r < 0.10\) correlación despreciable

Coeficiente de correlación

\[\large{r = \frac {\sum (x_i - \bar x)(y_i - \bar y)}{\sqrt{\sum(x_i - \bar x)^2 \sum(y_i - \bar y)^2}}}\]

\[\large{r = \frac {S_{xy}}{S_x S_y}}\]

Test de significación estadística de r

\[\large{t_{n-2} = \frac{r \sqrt{n-2}}{\sqrt{1-r^2}}}\]

Test de significación estadística de r e IC: transformación z de Fisher

\[\large{z = \frac{1}{2} [\ln{(1+r)} - \ln{(1-r)}]}\]

\[\large{\sigma_z = \frac{1}{\sqrt{n-3}}}\]

Regresión lineal

lm()

Un modelo de regresión implica la existencia de una variable dependiente o resultado \(y\), y otra variable independiente o explicativa \(x\), (o un conjunto de variables explicativas), en un modelo \(y \sim x\)

Un modelo de regresión no es más que un modelo de diferencia de medias. Los coeficientes de regresión se interpretan como (a) la diferencia de medias en la variable resultado cuando se incrementa en 1 punto el valor de la variable explicativa cuando esta es numérica, o (b) la media del grupo de referencia \(b_0\), o la diferencia de medias entre el grupo de referencia y el de contraste \(b_1\) cuando la variable explicativa es categórica

Ecuaciones de la regresión

\[\large{E(y_i|x_i) = \alpha + \beta x_i}\]

\(\alpha\) es el intercepto, el valor de \(y\), cuando \(x = 0\)

\(\beta\) es la pendiente que refleja el cambio en \(y\) por cada cambio de 1 unidad en \(x\)

\[\large{\hat y_i = \alpha + \beta x_i + \epsilon_i}\]

\(\hat y_i\) es el valor predicho para una observación en base al valor \(x_i\)

\(\epsilon_i\) es el error aleatorio (residual) que sigue una distribución normal

Estimación del intercepto y de la pendiente, y relación con la correlación

\[\large{\hat \beta = \frac{cov_{x,y}}{Var_x}}\]

\[\large{\hat \alpha = \bar y - \hat \beta \bar x}\]

La línea de regresión siempre cruza por el punto \((\bar x, \bar y)\)

\[\large{\hat r = \hat \beta \frac{SD_x}{SD_y}}\]

Regresión lineal y los tests para diferencia de medias

Las diferencias de medias para 2 grupos (t-test), o para más de 2 grupos (ANOVA), no son más que situaciones particulares del modelo general de regresión lineal
Los siguientes ejemplos lo aclararán mejor

Dos tratamientos analgésicos y tiempo de latencia (segundos) en placa caliente

ind <- c(1:20)
arm <- rep(0:1, each = 10)
jump.time <- c(6.7, 3.8, 5.2, 4.2, 5.3, 3.2, 7.6, 7.4, 3.1, 4.2, 5.9, 6, 5.2, 6.8, 14.7, 15.5, 4.2, 14, 4.7, 4.7)
rats <- data.frame(ind, arm, jump.time)
head(rats)

  ind arm jump.time
1   1   0       6.7
2   2   0       3.8
3   3   0       5.2
4   4   0       4.2
5   5   0       5.3
6   6   0       3.2

Resultados según el t-test

t.test(jump.time ~ factor(arm), data = rats, var.equal = T)


    Two Sample t-test

data:  jump.time by factor(arm)
t = -2.0019, df = 18, p-value = 0.0606
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.3533931  0.1533931
sample estimates:
mean in group 0 mean in group 1 
           5.07            8.17

Resultados según regresión lineal

print(lm(jump.time ~ factor(arm), data = rats))


Call:
lm(formula = jump.time ~ factor(arm), data = rats)

Coefficients:
 (Intercept)  factor(arm)1  
        5.07          3.10

Pulsatilidad arteria cerebral y consumo cannabis (3 grupos)

ind <- c(1:18)
arm <- rep(0:2, each = 6) # 0:control, 1:ex-usuarios, 2:usuarios activos
pulse <- c(.9, .9, .9, .8, .8, NA, 1, 1, 1, 1, .9, .9, 1, 1, 1, 1, .9, NA)
drug <- data.frame(ind, arm, pulse)
head(drug)

  ind arm pulse
1   1   0   0.9
2   2   0   0.9
3   3   0   0.9
4   4   0   0.8
5   5   0   0.8
6   6   0    NA

Resultados según ANOVA

summary(aov(pulse ~ factor(arm), data = drug))

            Df  Sum Sq  Mean Sq F value  Pr(>F)   
factor(arm)  2 0.04417 0.022083   8.612 0.00415 **
Residuals   13 0.03333 0.002564                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 observations deleted due to missingness

Resultados según regresión lineal - ANOVA

anova(lm(pulse ~ factor(arm), data = drug))

Analysis of Variance Table

Response: pulse
            Df   Sum Sq   Mean Sq F value   Pr(>F)   
factor(arm)  2 0.044167 0.0220833  8.6125 0.004152 **
Residuals   13 0.033333 0.0025641                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Resultados según regresión lineal

print(lm(pulse ~ factor(arm), data = drug))


Call:
lm(formula = pulse ~ factor(arm), data = drug)

Coefficients:
 (Intercept)  factor(arm)1  factor(arm)2  
      0.8600        0.1067        0.1200

Regresión y regresión a la media

Francis Galton: Law of universal regression

“Each peculiarity in a man is shared by his kinsman, but on the average in a less degree”

Un ejemplo de regresión ecológica sin falacia (Doll, 1955)

En Friedman, Pisani & Purves (1978). smoke es el número de cigarrillos per capita consumidos en 1930. lungcancer.deaths son las muertes por cáncer de pulmon (por 10⁶) en 1950

pais <- c("Australia", "Canada", "Dinamarca", "Finlandia", "Inglaterra",
          "Holanda", "Islandia", "Noruega", "Suecia", "Suiza", "USA")
smoke <- c(480, 500, 380, 1100, 1100, 490, 230, 250, 300, 510, 1300)
lungcancer.deaths <- c(180, 150, 170, 350, 460, 240, 60, 90, 110, 250, 200)
Doll1955 <- data.frame(pais, smoke, lungcancer.deaths, stringsAsFactors = FALSE)

Regresión ecológica

Doll1955

         pais smoke lungcancer.deaths
1   Australia   480               180
2      Canada   500               150
3   Dinamarca   380               170
4   Finlandia  1100               350
5  Inglaterra  1100               460
6     Holanda   490               240
7    Islandia   230                60
8     Noruega   250                90
9      Suecia   300               110
10      Suiza   510               250
11        USA  1300               200

Regresión ecológica

fit <- lm(lungcancer.deaths ~ smoke, data = Doll1955)
plot(Doll1955$smoke, Doll1955$lungcancer.deaths)
abline(coef(fit))

plot of chunk unnamed-chunk-60

Análisis de frecuencias y proporciones: datos categóricos

binom.test(), test exacto 1 grupo
prop.test(), test para varios grupos
table(), tabulaciones cruzadas
chisq.test(), test de la ji al cuadrado
fisher.test(), test exacto de Fisher cuando el valor esperado en alguna de las celdas de la tabla a analizar es < 5

Retención en tratamiento con heroína vs metadona: análisis de proporciones (binom.test)

binom.test(x=22, n=31, p=0.61) # 22/31 es el éxito en DAM; 0.61 es el exito en MET


    Exact binomial test

data:  22 and 31
number of successes = 22, number of trials = 31, p-value = 0.2761
alternative hypothesis: true probability of success is not equal to 0.61
95 percent confidence interval:
 0.5196393 0.8577715
sample estimates:
probability of success 
             0.7096774

Retención en tratamiento con heroína vs metadona: análisis de proporciones (prop.test)

r <- c(22, 19); tot <- c(31, 31) # DAM/MET
prop.test(r, tot)


    2-sample test for equality of proportions with continuity
    correction

data:  r out of tot
X-squared = 0.28804, df = 1, p-value = 0.5915
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.1698584  0.3634068
sample estimates:
   prop 1    prop 2 
0.7096774 0.6129032

Retención en tratamiento con heroína vs metadona: tabulación de datos

heroin <- matrix(c(22, 19, 9, 12), nrow = 2, byrow = T)
colnames(heroin) <- c("DAM", "MET")
rownames(heroin) <- c("success", "failure")
heroin

        DAM MET
success  22  19
failure   9  12

Retención en tratamiento con heroína vs metadona: gráfico

barplot(heroin)

plot of chunk unnamed-chunk-64

Retención en tratamiento con heroína vs metadona: test Ji al cuadrado

chisq.test(heroin)


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

data:  heroin
X-squared = 0.28804, df = 1, p-value = 0.5915

Retención en tratamiento con heroína vs metadona: test exacto de Fisher

fisher.test(heroin)


    Fisher's Exact Test for Count Data

data:  heroin
p-value = 0.5921
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.4723836 5.1300383
sample estimates:
odds ratio 
  1.532994

Retención en tratamiento con heroína vs metadona: RD

p1 <- 22/31; p0 <- 19/31
n1 <- 31; n0 <- 31
RD <- p1 - p0
SE.RD <- sqrt(((p1*(1-p1))/n1) + ((p0*(1-p0))/n0))
CI95.RD <- RD + c(-1, 1) * (1.96 * SE.RD)
RD; CI95.RD

[1] 0.09677419

[1] -0.1376046  0.3311530

Retención en tratamiento con heroína vs metadona: lnRR

p1 <- 22/31; p0 <- 19/31
lnRR <- log(p1 / p0)
SE.lnRR <- sqrt((1/22)-(1/31)+(1/19)-(1/31))
CI95.lnRR <- lnRR + c(-1, 1) * (1.96 * SE.lnRR)
lnRR; CI95.lnRR

[1] 0.1466035

[1] -0.212510  0.505717

Retención en tratamiento con heroína vs metadona: RR

RR <- exp(0.1466035)
RR95.low <- exp(-0.212510)
RR95.high <- exp(0.505717)
RR; RR95.low; RR95.high

[1] 1.157895

[1] 0.8085522

[1] 1.658174

Retención en tratamiento con heroína vs metadona: lnOR

odds1 <- 22/9; odds0 <- 19/12
lnOR <- log(odds1 / odds0)
SE.lnOR <- sqrt((1/22)+(1/9)+(1/19)+(1/12))
CI95.lnOR <- lnOR + c(-1, 1) * (1.96 * SE.lnOR)
lnOR; CI95.lnOR

[1] 0.4342855

[1] -0.6258019  1.4943730

Retención en tratamiento con heroína vs metadona: OR

OR <- exp(0.4342855)
OR95.low <- exp(-0.6258019)
OR95.high <- exp(1.4943730)
OR; OR95.low; OR95.high

[1] 1.54386

[1] 0.5348324

[1] 4.456541

Eficacia vacuna Salk contra la polio (Snedecor & Cochran, 1980)

paralisis <- c(33, 115) # casos de parálisis en vacuna / placebo respectivamente
N <- c(200745, 201229) # Número total de sujetos aleatorizados a vacuna / placebo
(risk.vacuna <- 33 / 200745 * 100000)

[1] 16.43877

(risk.placebo <- 115 / 201229 * 100000)

[1] 57.14882

Eficacia vacuna Salk contra la polio (Snedecor & Cochran, 1980)

(RD <- risk.placebo - risk.vacuna)

[1] 40.71005

(NNT <- ceiling(1/RD))

[1] 1

Eficacia vacuna Salk contra la polio (Snedecor & Cochran, 1980)

prop.test(paralisis, N)


    2-sample test for equality of proportions with continuity
    correction

data:  paralisis out of N
X-squared = 44.153, df = 1, p-value = 3.038e-11
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.0005306031 -0.0002835980
sample estimates:
      prop 1       prop 2 
0.0001643877 0.0005714882

Más allá del bien y del mal: de la regresión lineal simple a los modelos multivariados

Funciones de los análisis multivariados
Control de variables de confusión (estudios experimentales/observacionales)
Mejorar la predicción (posibles problemas de sobreajuste)
Valorar modificaciones del efecto por interacción entre variables explicativas (tamaños muestrales necesarios)
Extensiones de la regresión lineal: GLMs

Variables de confusión

Se considera a una variable como potencialmente confusora si modifica el coeficiente \(\hat \beta\), entre el análisis simple y el ajustado en mas de un 10%
No se deben evaluar las variables como potencialmente confusoras por su efecto en los p-valores

Problemas de sobreajuste en modelos multivariados

En modelos multivariados se pueden obtener valores altamente significativos, pero sin ningún sentido, simplemente por introducir demasiados predictores en el modelo
El modelo ajusta perfectamente los datos pero no tiene valor predictivo en una nueva muestra (diferenciación entre muestra de entrenamiento y muestra de test)
Como regla común se considera que se necesitan al menos 10 sujetos por cada variable predictora (más el intercepto) en un modelo multivariado de regresión

Modelos lineales generales (GLMs)

El modelo de regresión lineal para una variable numérica de resultado, el modelo de regresión logística para una variable binaria/dicotómica de resultado, y el modelo de regresión de Poisson para contajes, no son más que aspectos específicos del modelo lineal general(GLM). En R se pueden ajustar con la función glm() que tiene el argumento family que describe la distribución del error y el argumento link que define la unión entre el resultado y la parte explicativa del modelo

Modelos y sintáxis general de GLM´s

Regresión lineal: glm(y ~ x, family=gaussian(link=“identity”), data=my.data)
Regresión logística: glm(y ~ x, family=binomial(link=“logit”), data=my.data)
Regresión de Poisson: glm(y ~ x, family=poisson(link=“log”), data=my.data)

Regresiones binomiales

glm(y ~ x, family=binomial(link=“logit”), data=my.data) coeficientes lnOR
glm(y ~ x, family=binomial(link=“log”), data=my.data) coeficientes lnRR
glm(y ~ x, family=binomial(link=“identity”), data=my.data) coeficientes RD

Análisis de supervivencia

R tiene la librería survival para realizar análisis de supervivencia, con funciones para obtener estimaciones de Kaplan-Meier (survfit()), log-rank test (survdiff()), regresión de Cox (coxph()), y para modelos paramétricos de supervivencia con la función survreg()

Resumen: test de hipótesis en R

t.test(): t-tests para 1 y 2 muestras
wilcox.test(): Wilcoxon tests para 1 y 2 muestras
var.test(): Varianza F-tests para 1 y 2 muestras
cor.test(): Correlación y p-valor (Pearson’s, Spearman’s, o Kendall’s)
binom.test(): Test de signos para una muestra binomial
prop.test(): Test binomial para comparar 2 proporciones
chisq.test(): Test chi-cuadrado para contajes
fisher.test(): Test exacto de Fisher para contajes

Resumen: test de hipótesis en R

friedman.test(): Test de la suma de rangos de Friedman
kruskal.test(): Test de suma de rangos de Kruskal-Wallis
shapiro.test(): Test de Shapiro-Wilk de normalidad para 1 muestra
ks.test(): Test de Kolmogorov-Smirnov para 1 o 2 muestras (varias distribuciones)
levene.test(): Test de igualdad de varianzas
bartlett.test(): Test de igualdad de varianzas
oneway.test(): Análisis de la varianza de 1 vía

Resumen: modelos lineales en R

Los modelos se estructuran en base a una fórmula, que es una descripción simbólica del modelo, y representa algebraicamente los efectos aditivos de acuerdo a los objetos que almacenan los datos:

respuesta ~ predictor lineal
\(y \sim 1\), modelo nulo
\(y \sim x_1\), modelo lineal simple
\(y \sim x_1 + x_2 + ... + x_{n-1} + x_n\), modelo lineal múltiple
\(y \sim x_1 * x_2\), es lo mismo que \(y \sim x_1 + x_2 + x_1:x_2\)

Resumen modelos lineales en R

Funciones para el ajuste de modelos lineales, y para extraer información de ellos

lm(): regresión lineal
aov(): análisis de la varianza
glm(): modelos lineales generalizados (GLMs)
update(): actualizar y ajustar un modelo lm o glm
step(): búsqueda paso-a-paso de modelos minimizando AIC
anova(): anális varianza modelo/comparación entre modelos anidados
summary(): resumen del ajuste del modelo

Resumen modelos lineales en R

Funciones para el ajuste de modelos lineales, y para extraer información de ellos

confint(): intervalos de confianza para los coeficientes del modelo
predict(): predicciones en una nueva matriz de datos
plot(): gráficos específicos
resid(): residuales (valores observados - valores predichos)
fitted(): valores predichos por el modelo
vcov(): matriz de varianzas-covarianzas del modelo
AIC(): Akaike’s index criterion para el modelo ajustado

Libreria lessR() y funciones estándar de R

La librería lessR() tiene como lema menos código, más resultados. Presenta funciones útiles para el análisis descriptivo y para tests estadísticos clásicos. Es recomendable tanto como herramienta de análisis como herrramienta de docencia

Como para cualquier package en CRAN, primero habría que descargarlo (install.packages(“lessR”)), y luego cargarlo en memoria (library(lessR)), para poder trabajar con él

t-test en la función estándar de R

data(sleep)
t.test(extra ~ factor(group), data = sleep, var.equal = TRUE)


    Two Sample t-test

data:  extra by factor(group)
t = -1.8608, df = 18, p-value = 0.07919
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.363874  0.203874
sample estimates:
mean in group 1 mean in group 2 
           0.75            2.33

t-test en la función de la libreria lessR - versión breve

library(lessR)
ttest(extra ~ group, data=sleep, brief=TRUE)


------------------------------------------------------------
Compare extra across group levels 2 and 1 
------------------------------------------------------------

extra for group 2:  n.miss = 0,  n = 10,  mean = 2.33,  sd = 2.00
extra for group 1:  n.miss = 0,  n = 10,  mean = 0.75,  sd = 1.79

---
t-cutoff: tcut =  2.101 
Standard Error of Mean Difference: SE =  0.85 

Hypothesis Test of 0 Mean Diff:  t = 1.861,  df = 18,  p-value = 0.079

Margin of Error for 95% Confidence Level:  1.78
95% Confidence Interval for Mean Difference:  -0.20 to 3.36

Sample Mean Difference of extra:  1.58
Standardized Mean Difference of extra, Cohen's d:  0.83
95% Confidence Interval for smd:  -0.10 to 1.74
--------------------------------------------------

plot of chunk unnamed-chunk-76

Comparación entre los resultados del t-test en R y en la librería lessR

La función t.test() devuelve el valor t, sus grados de libertad y el p-valor asociado, las medias de cada grupo, y el IC al 95% de su diferencia
La función ttest() de lessR, versión breve, devuelve además la media, DE, y número de observaciones por grupo, la diferencia de medias bruta y estandarizada, y su IC al 95%. Incluye un gráfico que permite apreciar la separación entre grupos

Comparación entre los resultados del t-test en R y en la librería lessR

La función ttest() de lessR, versión extendida, añade además (1) un análisis de las asunciones subyacentes (distribución normal de la variable resultado en cada grupo, igualdad de varianzas entre los grupos en comparación), (2) los resultados asumiendo - y sin asumir - igualdad de varianzas, y (3) un análisis de la importancia práctica de la diferencia de medias, y de la diferencia estandarizada de medias, en caso de que se hubieran establecido previamente

Relación entre las funciones de R y de la libreria lessR

merge() —> Merge()
barplot() —> BarChart()
boxplot() —> BoxPlot()
plot() —> Plot()
hist() —> Histogram()
summary() —> SummaryStats()
power.t.test() —> ttestPower()

Relación entre las funciones de R y de la librería lessR

dnorm() —> Density()
dnorm() —> prob.znorm()
pnorm() —> prob.norm()
qt(); pnorm() —> prob.tcut()

Relación entre las funciones de R y de la librería lessR

t.test() —> ttest()
aov() —> ANOVA()
cor(), cor.test(), cov() —> Correlation()
lm() —> Regression()
lm(); aov(); t.test() —> Model()
glm(family=“binomial”) —> Logit()

Algunos consejos finales sobre el trabajo con R

A diferencia de otros programas, R tiene una curva de aprendizaje larga. Aprender R es utilizarlo, cometer errores y resolverlos, sea por uno mismo, o con el apoyo de la muy extensa comunidad de R. En mi experiencia, la mejor forma de aprender los rudimentos de un nuevo programa ha sido siempre el utilizarlo en un proyecto nuevo, y desarrollarlo enteramente con ese programa

Continuar con R

El objetivo de este curso ha sido introducir R como herramienta de análisis estadístico y como tal no ha desarrollado contenidos de estadística teórica (teorema central del límite, probabilidad y pruebas de hipótesis), ni contenidos aplicados de interés para estudiantes de ramas biomédicas (regresión logística, análisis de supervivencia…). Tampoco ha presentado los aspectos de programación y simulación que son uno de los principales recursos de R en el trabajo académico. Para quién quiera profundizar en esos aspectos, las referencias bibliográficas del curso son un buen camino a seguir…

Referencias

NJ Horton & K Kleinman (2015). Using R and RStudio for data management, statistical analysis, and graphics. CRC Press.
TW McFarland (2014). Introduction to data analysis and graphical presentation in biostatistics with R. Springer.
P Dalgaard (2008). Introductory statistics with R. Springer.
HM Kaltenbach (2012). A concise guide to statistics. Springer.

Referencias para usuarios de otros paquetes estadísticos

RA Muenchen (2011). R for SAS and SPSS users.Springer.
RA Muenchen & JM Hilbe (2010). R for Stata users. Springer.

Anexo: ¿Cómo “pasar” de la estadística al preparar un manuscrito (Med Clin 1999; 113:138-49)

Si tiene datos por ahí, no los desperdicie. Introdúzcalos en el ordenador y presente como significativa cualquier comparación en la que \(P<0.05\)
Si el análisis previsto de los datos no le proporciona el resultado esperado, pruebe con otro
Si sus resultados no son interesantes, fuerce a su ordenador a comparar grupos y subgrupos. Puede que consiga engañar a alguien
Ignore todas las perdidas y no respuestas, y analice exclusivamente a los individuos que tenga completos

Anexo: ¿Como “pasar” de la estadística al preparar un manuscrito (Cont.)

Si le molestan algunos valores extremos, elimínelos del análisis. Pero si le convienen, aunque parezcan sospechosos, déjelos dentro
Suponga que todos los datos numéricos de su muestra son cuantitativos y proceda en consecuencia. Puede que el ordenador le calcule la media de su grupo sanguíneo
No compruebe si sus datos siguen una distribución normal. Si lo hace, se le pueden complicar las cosas
Si encuentra diferencias basales entre grupos a favor del grupo de estudio, recuerde no tenerlas en cuenta

Anexo: ¿Como “pasar” de la estadística al preparar un manuscrito (Cont.)

Suponga que todo se puede correlacionar con todo, calcule el coeficiente de correlación de Pearson, y concluya que una \(r\) significativa demuestra causalidad
Si no obtiene resultados significativos, no los presente. Y es más, ni siquiera los mencione en la discusión
Si las diferencias entre grupos ya se vislumbran antes de terminar el estudio, deténgalo y empiece a escribir. Alternativamente, si todavía no ha tenido esa suerte, prolónguelo un poquito más
En definitiva, no se preocupe demasiado por el diseño del estudio y reclute como pueda el máximo número de individuos. Analícelos, asuma que su muestra es representativa de la población y concluya

Anexo: Aspectos a tener en cuenta en los análisis de datos

Revisión extensa de la literatura con carácter previo (diseños, resultados positivos y negativos…)
Colaboración con el analista desde el comienzo del estudio
Muestra; diseño, y control de sesgos; diseño de los análisis
Problemas comunes: triple ciego; todas con todas
Costo desglosado por capítulos (proceso de datos, elaboración informe estadístico)
Análisis de los sesgos
Los estudios con resultados negativos también son importantes
No torturar nunca los datos

Anexo: Relación entre pruebas paramétricas y no paramétricas

t-test para muestras independientes –> test de Wilcoxon-Mann-Whitney
t-test para muestras emparejadas –> test de Wilcoxon para muestras emparejadas
ANOVA de una vía –> test de Kruskal-Wallis
Correlación de Pearson –> correlación de Spearman

Anexo: Múltiples pruebas de hipótesis

Se realizan múltiples pruebas de hipótesis cuando:

se comparan varios tratamientos entre si
se utilizan varias resultados para evaluar cambios
Se evalúa el mismo resultado en distintos tiempos con tests estadísticos en cada uno
análisis de subgrupos

Anexo: Múltiples pruebas de hipótesis

Correcciones para tests múltiples
Bonferroni: \(P-corrected = \frac{0.05}{n}\)
Análisis
Objetivo primario –> resultado primario –> análisis primario
Análisis secundarios –> análisis exploratorios

Anexo: Análisis secundarios no protocolizados

Enjoy the result you have found by exploratory data analysis, for you will not find it again

Stephen Senn

Anexo: Análisis de resultados para desenlaces cuantitativos

Tratamiento de 2 ramas
t-test
Wilcoxon
Regresión lineal
Tratamiento con más de 2 ramas
Análisis de la varianza (1 vía)
Regresión lineal

Anexo: Análisis de resultados para desenlaces dicotómicos (0/1)

Tratamiento de 2 ramas
\(\chi^2\)
Regresión logística
Tratamiento con más de 2 ramas
\(\chi^2\)
Regresión logística

Anexo: Análisis de resultados para desenlaces de tiempo hasta evento

Tratamiento de 2 ramas
Kaplan-Meier
Regresión de Cox
Tratamiento con más de 2 ramas
Kaplan-Meier
Regresión de Cox