Clase 3_analisis exporatorio

Vista previa del material en texto

Andrés Del Río 31/08/2020 
Clase 3: Análisis de datos exploratorios 
El Análisis de Datos Exploratorios (EDA) es una parte clave de lo que hacemos cuando 
analizamos datos. Comenzamos cada análisis con EDA para familiarizarnos con los datos. 
En primer lugar, queremos comprobar si algunas de las muestras o experimentos 
producen datos inutilizables y sacarlos del análisis. También realizamos EDA al final para 
comprobar si hay resultados sin sentido. De hecho, uno debe realizar EDA durante todo el 
análisis. Por ejemplo, antes de aplicar cualquiera de las técnicas que hemos aprendido, 
queremos asegurarnos de que los datos están de acuerdo con los supuestos necesarios. 
Presentaremos algunas herramientas EDA básicas como histograma, la gráfica Q-Q, 
gráficas de dispersión, diagramas de caja, estratificación, entre otras. 
Histograma 
Introducción 
"El mayor valor de una imagen es cuando nos obliga a darnos cuenta de lo que nunca 
esperábamos ver." -John W. Tukey 
Los sesgos, los errores sistemáticos y la variabilidad inesperada son comunes en los datos 
genómicos. El no descubrir estos problemas a menudo conduce a análisis defectuosos y 
descubrimientos falsos. Por ejemplo, tenga en cuenta que los experimentos a veces fallan 
y no todas las canalizaciones de procesamiento de datos están diseñadas para detectarlos. 
Sin embargo, estas canalizaciones todavía le dan una respuesta y el de los resultados 
finales puede ser difícil o imposible notar un error se hizo. 
Histogramas 
Podemos pensar en cualquier conjunto de datos dado como una lista de números. 
Supongamos que ha medido las alturas de todos los hombres en una población. Imagina 
que tienes que describir estos números a alguien que no tiene idea de cuáles son estas 
alturas, digamos un extraterrestre que nunca ha visitado la tierra. 
library(UsingR) 
x=father.son$fheight 
head(father.son) 
## fheight sheight 
## 1 65.04851 59.77827 
## 2 63.25094 63.21404 
## 3 64.95532 63.34242 
## 4 65.75250 62.79238 
## 5 61.13723 64.28113 
## 6 63.02254 64.24221 
Un enfoque es simplemente enumerar todos los números para que el extraterrestre los 
vea. Aquí hay 20 alturas seleccionadas al azar de 1.078. 
Andrés Del Río 31/08/2020 
round(sample(x,20),1) 
## [1] 70.5 70.9 69.8 67.2 66.0 69.8 65.9 68.5 68.3 69.2 65.8 68.1 66
.6 66.3 
## [15] 69.1 64.5 69.2 62.9 68.3 68.4 
A partir del escaneo a través de estos números comenzamos a tener una idea aproximada 
de cómo se ve toda la lista, pero sin duda es ineficiente. Podemos mejorar rápidamente 
este enfoque creando bins, por ejemplo, redondeando cada valor a su entero más cercano 
e informando del número de individuos en cada bin. Una gráfica de estas alturas se llama 
histograma 
hist(x,breaks=seq(floor(min(x)),ceiling(max(x))),main="",xlab="Height"
) 
 
Mostrar esta trama al alienígena es mucho más informativo que mostrar los números. 
Tenga en cuenta que con esta gráfica simple podemos aproximar el número de individuos 
en cualquier intervalo dado. Por ejemplo, hay alrededor de 70 individuos de más de seis 
pies (72 pulgadas) de altura. 
xs<-seq(floor(min(x)),ceiling(max(x)),0.1) 
plot(xs,ecdf(x)(xs),type="l",xlab="x=Height",ylab="F(x)") 
Andrés Del Río 31/08/2020 
 
Ejercicios de histograma 
 
Ejercicios de histograma #1 
Dado el histograma anterior, ¿cuántas personas tienen entre 35 y 45 años? 
Andrés Del Río 31/08/2020 
#sum(age>=35 & age<45) 
A:6 
Aproximación normal 
Si en lugar de los números totales reportamos las proporciones, entonces el histograma es 
una distribución de probabilidad. La distribución de probabilidad que vemos arriba se 
aproxima a una que es muy común en una naturaleza: la curva de campana o distribución 
normal o distribución gaussiana. Cuando el histograma de una lista de números se 
aproxima a la distribución normal podemos utilizar una fórmula matemática conveniente 
para aproximar la proporción de individuos en cualquier intervalo dado 
Pr(a<x<b)=∫ba12πσ2−−−−√exp(−(x−μ)22σ2)dxPr(a<x<b)=∫ab12πσ2exp⁡(−(x−μ
)22σ2)dx 
En este caso, σ y μ se denominan como media y desviación estándar. Si esta aproximación 
se mantiene para nuestra lista, entonces la media de la población y la varianza de nuestra 
lista se pueden utilizar en la fórmula anterior. Para ver esto con un ejemplo recuerde que 
arriba señalamos que 70 individuos o 6% de nuestra población eran más altos que 6 pies. 
La aproximación normal funciona bien: 
1-pnorm(72,mean(x),sd(x)) 
## [1] 0.05806108 
Una característica muy útil de esta aproximación es que sólo se necesita saber μ y σ para 
describir toda la distribución. Todo lo que realmente tenemos que decirle a nuestro amigo 
alienígena es que las alturas siguen una distribución normal con una altura media de 68''y 
una desviación estándar de 3''. A partir de esto podemos calcular la proporción de 
individuos en cualquier intervalo. 
Gráficas cuantitativas cuantitativas (gráficos Q-Q) 
Para corroborar que una distribución teórica, por ejemplo la distribución normal, es de 
hecho una buena aproximación, podemos utilizar parcelas cuantitativas-cuantitativas (qq-
plots). Los cuantiles se entienden mejor considerando el caso especial de los percentiles. El 
percentil p-ésimo de una lista de una distribución se define como el número q que es 
mayor que p% de números (por lo que la inversa de la función de distribución acumulativa 
que definimos anteriormente). Por ejemplo, la mediana del percentil 50 es la mediana. 
Podemos calcular los percentiles para nuestra lista de alturas: 
library(rafalib) 
data(father.son,package="UsingR") 
x <- father.son$fheight 
 
Primer ejemplo de qqplot. 
Andrés Del Río 31/08/2020 
qqnorm(x) 
qqline(x) 
 
Segundo ejemplo de qqplot. Aquí utilizamos la función qqnorm que calcula los cuantiles 
normales teóricos automáticamente. 
Sin embargo, la función qqnorm se traza con una distribución normal estándar. Esta es la 
razón por la que la línea tiene la pendiente popsd(x) e interceptar media(x). 
En el ejemplo anterior, los puntos coinciden muy bien con la línea. De hecho, podemos 
ejecutar simulaciones de Monte Carlo para ver parcelas como esta para los datos que se 
sabe que se distribuyen normalmente. 
n <-1000 
x <- rnorm(n) 
qqnorm(x) 
qqline(x) 
Andrés Del Río 31/08/2020 
 
Ejemplo de la función qqnorm. Aquí lo aplicamos a los números generados para seguir una 
distribución normal. 
También podemos tener una idea de cómo se verán los datos distribuidos no 
normalmente en una gráfica qq. Aquí generamos datos de la distribución t con diferentes 
grados de libertad. Observe que cuanto más pequeños son los grados de libertad, más 
gordos son las colas. Llamamos a estas "colas gordas" porque si trazamos una densidad 
empírica o histograma, la densidad en los extremos sería mayor que la curva teórica. En la 
gráfica qq, esto se puede ver en que la curva es más baja que la línea de identidad en el 
lado izquierdo y más alto en el lado derecho. Esto significa que hay más valores extremos 
de lo predicho por la densidad teórica trazada en el eje X. 
dfs <- c(3,6,12,30) 
mypar(2,2) 
for(df in dfs){ 
 x <- rt(1000,df) 
 qqnorm(x,xlab="t quantiles",main=paste0("d.f=",df),ylim=c(-6,6)) 
 qqline(x) 
} 
Andrés Del Río 31/08/2020 
 
Generamos datos distribuidos en t durante cuatro grados de libertad y hacemos qqplots 
contra cuantiles teóricos normales. 
Ejercicios de QQ-plot 
Descargue este archivo RData en su directorio de trabajo: link. A continuación, cargue los 
datos en R con el siguiente comando: 
load("skew.RData") #Debe tener una matriz dimensional de 1000 x 9 'dat' 
dim(dat) 
## [1] 1000 9 
str(dat) 
## num [1:1000, 1:9] -0.626 0.184 -0.836 1.595 0.33 ... 
names(dat) 
## NULL 
Andrés Del Río 31/08/2020 
Usando gráficas QQ, compare la distribución de cada columna de la matriz con una 
normal. Es decir, utilice qqnorm() en cada columna. Para lograr esto rápidamente, puede 
utilizar la siguiente línea de código para configurar unacuadrícula para trazados de 3x3 a 
9. ("mfrow" significa que queremos una cuadrícula multifigura llenada de fila por fila. Otra 
opción es mfcol.) 
par(mfrow = c(3,3)) 
#A continuación, puede utilizar un for loop, para recorrer las columnas y mostrar una gráfica qqnor
m() a la vez. 
for (i in 1:9) { 
 x <- dat[,i] 
 qqnorm(x,xlab="quantiles",main=paste0("Col No=",i)) 
 qqline(x) 
} 
 
#plotting múltiples histogramas 
par(mfrow = c(3,3)) 
#A continuación, puede utilizar un for loop, para recorrer las columnas y mostrar una gráfica qqn
orm() a la vez. 
for (i in 1:9) { 
 x <- dat[,i] 
Andrés Del Río 31/08/2020 
hist(x,xlab="X",main=paste0("Col.No=",i)) 
} 
 Identifique las dos columnas que están sesgadas. 
Examine cada una de estas dos columnas utilizando un histograma. Tenga en cuenta qué 
columna tiene "sesgo positivo", es decir, el histograma muestra una cola larga a la derecha 
(hacia valores más grandes). Tenga en cuenta qué columna tiene "sesgo negativo", es 
decir, una cola larga a la izquierda (hacia valores más pequeños). Tenga en cuenta que el 
sesgo positivo parece una curva de modelado hacia arriba en una gráfica qqnorm(), 
mientras que el sesgo negativo parece una curva de forma hacia abajo. 
Puede utilizar la siguiente línea para restablecer el gráfico para mostrarlo uno a la vez: 
par(mfrow-c(1,1)) 
Ejercicios de trazado QQ #1 
¿Qué columna tiene sesgo positivo (una cola larga a la derecha)? 
par(mfrow=c(1,1)) 
hist(dat[,4]) 
Andrés Del Río 31/08/2020 
 A:
4 
Ejercicios de trazado QQ #2 
¿Qué columna tiene sesgo negativo (una larga cola a la izquierda)? 
hist(dat[,9]) 
 A:9 
Andrés Del Río 31/08/2020 
Boxplots 
Los datos no siempre se distribuyen normalmente. Los ingresos son un ejemplo 
ampliamente citado. En estos casos, la desviación media y estándar no son 
necesariamente informativas, ya que no se puede inferir la distribución de solo estos dos 
números. Las propiedades descritas anteriormente son específicas de la normal. Por 
ejemplo, la distribución normal no parece ser una buena aproximación para la 
compensación directa para 199 CEOs de los Estados Unidos en el año 2000. 
data(exec.pay,package="UsingR") 
mypar(1,2) 
hist(exec.pay) 
qqnorm(exec.pay) 
qqline(exec.pay) 
 
Histograma y QQ-gráfico de pago ejecutivo. 
Además de las gráficas qq, un resumen práctico de los datos es calcular 3 percentiles: 25-o, 
50-th (la mediana) y el 75-o. Una gráfica de caja muestra estos 3 valores junto ±con un 
rango de los puntos dentro de la medianade 1,5 (percentil 75 -percentil 25). Los valores 
fuera de este rango se muestran como puntos y a veces se denominan valores atípicos.. 
boxplot(exec.pay, ylab="10,000s of dollars", ylim=c(0,400)) 
Andrés Del Río 31/08/2020 
 
Una simple gráfica de caja de la paga ejecutiva. 
Aquí mostramos sólo una gráfica de caja. Sin embargo, uno de los grandes beneficios de 
las gráficas de caja es que podríamos mostrar fácilmente muchas distribuciones en una 
parcela, a reenciéndolas, una al lado de la otra. Veremos varios ejemplos de esto a lo largo 
del libro. 
Ejercicios de Boxplot 
El conjunto de datos InsectSprays mide los recuentos de insectos en unidades 
experimentales agrícolas tratadas con diferentes insecticidas. Este conjunto de datos se 
incluye en R y puede examinarlo escribiendo: 
head(InsectSprays) 
## count spray 
## 1 10 A 
## 2 7 A 
## 3 20 A 
Andrés Del Río 31/08/2020 
## 4 14 A 
## 5 14 A 
## 6 12 A 
Pruebe dos formas equivalentes de dibujar diagramas de caja en R, utilizando el conjunto 
de datos InsectSprays. A continuación se muestra el pseudocódigo, que debe modificar 
para trabajar con el conjunto de datos InsectSprays. 
1. Usando split: 
boxplot(split(values, factor)) 
1. utilizando una fórmula: 
boxplot(values ~ factor) 
Ejercicios de diagrama de caja #1 
¿Qué spray parece el más eficaz (tiene la mediana más baja)? 
#Usando split 
boxplot(split(InsectSprays$count, InsectSprays$spray)) 
 
#Usando formula 
boxplot(InsectSprays$count ~ InsectSprays$spray) 
 A
:C 
Andrés Del Río 31/08/2020 
Boxplot Exercises #2 
Consideremos una muestra aleatoria de remates del Maratón de Nueva York en 2002. Este 
conjunto de datos se puede encontrar en el paquete UsingR. Cargue la biblioteca y, a 
continuación, cargue el conjunto de datos nym.2002. 
library(dplyr) 
data(nym.2002, package="UsingR") 
Utilice diagramas de caja e histogramas para comparar los tiempos de acabado de 
hombres y mujeres. ¿Cuál de las siguientes opciones describe mejor la diferencia? 
data(nym.2002, package="UsingR") 
str(nym.2002) 
## 'data.frame': 1000 obs. of 5 variables: 
## $ place : num 3592 13853 12256 10457 9686 ... 
## $ gender: Factor w/ 2 levels "Female","Male": 2 1 2 1 2 2 1 2 2 2 
... 
## $ age : num 52 40 31 33 33 40 30 27 42 48 ... 
## $ home : Factor w/ 165 levels "","--","AB","ABE",..: 50 107 47 91 
107 103 23 50 100 15 ... 
## $ time : num 217 273 265 256 252 ... 
head(nym.2002) 
## place gender age home time 
## 3475 3592 Male 52 GBR 217.4833 
## 13594 13853 Female 40 NY 272.5500 
## 12012 12256 Male 31 FRA 265.2833 
## 10236 10457 Female 33 MI 256.1500 
## 9476 9686 Male 33 NY 252.2500 
## 1720 1784 Male 40 NJ 201.9667 
library(dplyr) 
males<-filter(nym.2002, gender=="Male") 
females<-filter(nym.2002, gender=="Female") 
 
library(rafalib) 
mypar(1,3) 
 
boxplot(females$time, males$time) 
hist(females$time,xlim=c(range( nym.2002$time))) 
hist(males$time,xlim=c(range( nym.2002$time))) 
Andrés Del Río 31/08/2020 
 A:Male and females have similar right skewed distributions with the former, 20 
minutes shifted to the left. 
Scatterplots y Correlación 
Los métodos descritos anteriormente se relacionan con variables univariadas. En las 
ciencias biomédicas, es común estar interesado en la relación entre dos o más variables. 
Un ejemplo clásico son los datos de altura padre/hijo utilizados por Francis Galton para 
entender la herencia. Si tuviéramos que resumir estos datos, podríamos usar los dos 
promedios y dos desviaciones estándar ya que ambas distribuciones están bien 
aproximadadas por la distribución normal. Este resumen, sin embargo, no describe una 
característica importante de los datos. 
data(father.son,package="UsingR") 
x=father.son$fheight 
y=father.son$sheight 
plot(x,y,xlab="Father's height in inches",ylab="Son's height in inches
",main=paste("correlation =",signif(cor(x,y),2))) 
https://en.wikipedia.org/wiki/Francis_Galton
Andrés Del Río 31/08/2020 
 
Alturas de las parejas de padre e hijo conspiraron entre sí. 
La parcela de dispersión muestra una tendencia general: cuanto más alto es el padre, más 
alto es el hijo. Un resumen de esta tendencia es el coeficiente de correlación, que en estos 
casos es 0.5. Motivaremos esta estadística tratando de predecir la altura del hijo usando la 
altura del padre. 
Estratificación 
Supongamos que se nos pide que adivinemos la altura de los hijos seleccionados al azar. La 
altura promedio, 68,7 pulgadas, es el valor con la proporción más alta (ver histograma) y 
sería nuestra predicción. Pero, ¿y si nos dicen que el padre mide 72 pulgadas de alto, sí 
que adivinamos 68.7? 
El padre es más alto que el promedio. Específicamente, es 1.75 desviaciones estándar más 
altas que el padre promedio. Entonces, ¿debemos predecir que el hijo es también 1.75 
desviaciones estándar más altas? Resulta que esto sería una sobreestimación. Para ver 
esto, miramos a todos los hijos con padres que son de unos 72 pulgadas. Hacemos esto 
estratificando las alturas del padre. 
groups <- split(y,round(x)) 
boxplot(groups) 
Andrés Del Río 31/08/2020 
 
Caja de alturas de hijo estratificadas por las alturas del padre. 
print(mean(y[ round(x) == 72])) 
## [1] 70.67719 
La estratificación seguida de diagramas de caja nos permite ver la distribución de cada 
grupo. La altura promediode los hijos con padres de 72 pulgadas de alto es de 70,7 
pulgadas. También vemos que las medianas de los estratos parecen seguir una línea recta 
(recuerde que la línea media en la gráfica de caja muestra la mediana, no la media). Esta 
línea es similar a la línea de regresión,con una pendiente que está relacionada con la 
correlación, como aprenderemos a continuación. 
Correlación de Spearman 
Al igual que la desviación media y estándar no son buenos resúmenes cuando los datos no 
son bien aproximados por la distribución normal, la correlación no es un buen resumen 
cuando los pares de listas no son aproximados por la distribución normal bivariada. Los 
ejemplos incluyen casos en los que la variable está relacionada con otra mediante una 
función parabólica. Otro ejemplo más común es causado por valores atípicos o valores 
extremos. 
a=rnorm(100);a[1]=10 
b=rnorm(100);b[1]=11 
plot(a,b,main=paste("correlation =",signif(cor(a,b),2))) 
cor(x, method= #c(‘pearon’,’spearman’), hay varias mas… 
Andrés Del Río 31/08/2020 
 
En el ejemplo anterior los datos no están asociados, pero para un par ambos valores son 
muy grandes. La correlación aquí es de aproximadamente 0.5. Esto es impulsado por sólo 
ese punto como sacarlo reduce a la correlación a aproximadamente 0. Un resumen 
alternativo para los casos con valores atípicos o valores extremos es la correlación de 
Spearman que se basa en rangos en lugar de los propios valores. 
Ejercicios de diagrama de dispersión 
Consideremos una muestra aleatoria de remates del Maratón de Nueva York en 2002. Este 
conjunto de datos se puede encontrar en el paquete UsingR. Cargue la biblioteca y, a 
continuación, cargue el conjunto de datos nym.2002. 
Aquí usaremos las gráficas que hemos aprendido para explorar un conjunto de datos: 
algunas estadísticas sobre una muestra aleatoria de corredores del Maraton de la ciudad 
de Nueva York en 2002. Este conjunto de datos se puede encontrar en el paquete UsingR 
(utilizado en la evaluación anterior). Cargue la biblioteca y, a continuación, cargue el 
conjunto de datos nym.2002 con la siguiente línea: 
data(nym.2002, package="UsingR") 
head(nym.2002) 
## place gender age home time 
## 3475 3592 Male 52 GBR 217.4833 
## 13594 13853 Female 40 NY 272.5500 
## 12012 12256 Male 31 FRA 265.2833 
## 10236 10457 Female 33 MI 256.1500 
## 9476 9686 Male 33 NY 252.2500 
## 1720 1784 Male 40 NJ 201.9667 
Ejercicios de trazado de dispersión #1 
Andrés Del Río 31/08/2020 
Utilice dplyr para crear dos nuevos marcos de datos: machos y mujeres, con los datos de 
cada género. Para los hombres, ¿cuál es la correlación de Pearson entre la edad y el 
tiempo para terminar? 
library(dplyr) 
males <- filter(nym.2002, gender=="Male") 
females <- filter(nym.2002, gender=="Female") 
 
cor(males$age,males$time) 
## [1] 0.2432273 
Un: 
Ejercicios de trazado de dispersión #2 
Para las mujeres, ¿cuál es la correlación de Pearson entre la edad y el tiempo para 
terminar? 
cor(females$age,females$time) 
## [1] 0.2443156 
A:0.2443156 
Ejercicios de trazado de dispersión #3 
Si interpretamos estas correlaciones sin visualizar los datos, llegaríamos a la conclusión de 
que cuanto más viejos seamos, más lentos corremos maratones, independientemente del 
género. Mira las gráficas de dispersión y las gráficas de cajas de los tiempos estratificados 
por grupos de edad (20-25, 25-30, etc.). Después de examinar los datos, ¿qué es una 
conclusión más razonable? 
library(rafalib) 
mypar(2,2) 
plot(females$age, females$time) 
plot(males$age, males$time) 
group <- floor(females$age/5) * 5 
boxplot(females$time~group) 
group <- floor(males$age/5) * 5 
boxplot(males$time~group) 
Andrés Del Río 31/08/2020 
 
R: Los tiempos de finalización son constantes hasta algún momento entre 50-60, luego nos 
ponemos más lentos.