Publicado el Tutorial06 de nuestro curso de Introducción a la Estadística.

Y hoy toca actualizar la información con la publicación de un nuevo tutorial, el sexto de la serie. Además, hemos añadido a nuestra página una lista de las erratas que hemos ido detectando.

Anuncios

Publicado el Tutorial05 de nuestro curso de Introducción a la Estadística.

Aunque hace unos días que nuestros alumnos de Alcalá lo tienen, hasta ahora no había tenido ocasión de anunciarlo como es debido aquí en el blog: ya está disponible el Tutorial05 del curso.

Biomatemáticos en EL PAÍS

Una brevísima nota para hacerme eco de que hoy la tribu de los biomatemáticos se ha asomado a las páginas de EL PÁIS, aunque sea con la triste excusa del caso de contagio de Ébola en Madrid. Supongo que para muchos lectores habrá sido la primera noticia de que exista semejante ocupación. Por si te ha picado la curiosidad, un posible punto de partida (dentro del mundo de la divulgación) podría ser el libro Las matemáticas de la vida, de Ian Stewart.

Una colección de entrevistas interesantes.

Esta entrada va a ser un ejemplo del cambio de ritmo que me he propuesto introducir en este blog. Sin renunciar a entradas mas enjundiosas, vamos a aumentar la frecuencia de publicación de otro tipo de entradas más breves. Especialmente cuando, como es el caso, alguien ha sido tan amable como para darme el trabajo ya hecho. El blog Nada es Gratis es una de mis lecturas habituales, supongo que porque me acerca a algunos temas a los que normalmente tengo pocas oportunidades de asomarme. Hace unos días Jesús Fernández-Villaverde ha publicado allí esta entrada en la que recopila varias entrevistas con gente tan interesante como Hadley Wickham, John Chambers o Yihui Xie. Muy recomendables todos los vídeos que se incluyen en esa entrada.

Estamos de vuelta.

Ha pasado más de un año desde la última entrada del blog. Entonces decía que el blog iba a bajar su ritmo para  dejar crecer al libro de Introducción a la Estadística que estábamos escribiendo (Marcos Marvá y yo mismo).  No pensé, en aquel momento, que eso fuera a suponer una parada tan larga. Pero el libro está prácticamente concluido (la última versión se puede ver en este enlace). Así que hemos vuelto al blog.

En este año he aprendido muchas cosas y he tenido tiempo para pensar sobre el uso futuro que le quiero dar al blog. En particular, a cuál será su relación con el libro.  Las primeras entradas del blog estaban planteadas como una especie de minicurso de introducción a la Estadística usando R, y en particular R-Commander (ver esta entrada y la serie que conduce a ella).   Esa tarea se ha desplazado, durante este año de parada, a los tutoriales del libro. Estos tutoriales acompañan a cada uno de los capítulos del libro, y representan la contraparte computacional y práctica del trabajo más teórico y conceptual que hacemos en el libro. En una próxima entrada me extenderé más sobre la estructura  y características del curso que hemos estado diseñando.

Para el blog, eso se traducirá en una mayor libertad temática. El trabajo formativo se deja para esos tutoriales y, en general, para la web del libro. Aquí vamos a abrir el foco para ir tratando, sin un plan preestablecido, todas las cuestiones que nos parezcan interesantes, dentro del contexto del análisis de datos y la computación científica en general. Seguiremos hablando de R, desde luego. Pero no quiero cerrar la puerta a otras herramientas (en particular, quiero hablar sobre Python) y a temas más generales, incluso con un ojo puesto en la actualidad, porque el análisis de datos nos asalta a diario desde todos los frentes.

Para no perder las buenas costumbres, termino esta entrada con la correspondiente recomendación musical: últimamente no consigo sacarme de la cabeza el tema Tennis Court de Lorde.

 

Compendio de comandos R (y otras cosas)

Esta entrada será breve. Vamos primero  con la parte más útil. A medida que se va aprendiendo a usar R, resulta conveniente disponer de un compendio o resumen de comandos, funciones, etc., como ayuda a la memoria. Hay varios de estos resúmenes (lo que en inglés de denomina a veces cheat sheet) publicados en la red, y yo suelo tener siempre a mano este en partiicular, elaborado por Tom Short, que me ahorra trabajo a la hora de localizar la sintaxis de un comando que se me resiste. Pero hay otros, que se pueden localizar fácilmente buscando en Google con R cheat sheet, o R reference card, etc.  Conviene mirar unos cuantos, localizar el que se adapte mejor a nuestras necesidades y conocimientos, y tenerlo a mano. No conozco ninguno en español que consiga empaquetar, en un espacio reducido, la cantidad de información de los que existen en inglés. Aprovecho para pedir que, si alguno de los lectores lo conoce, me lo haga saber, y lo reflejaremos aquí.

Las otras cosas: algo sobre planes futuros

La baja actividad del blog de estas semanas no se debe al verano. Los apuntes de las asignaturas de Estadística que imparto van camino de convertirse en un libro, escrito junto con el profesor M.Marvá. Este verano el trabajo se concentra en esa dirección. Queremos, además, escribir un libro que de alguna manera pueda seguir vivo una vez publicado, lo cual significa, entre otras cosas, que este blog se convierta en una extensión natural del libro, que irá acompañado además de un sitio Moodle, y de otros recursos en la red. Por eso, y para definir mejor la relación entre todos estos elementos, mientras el libro avanza el blog frena un poco el paso, pero no se para. En concreto, esta misma semana espero poder publicar la siguiente entrada sobre programación en R, que ya se ha hecho esperar bastante.

QQ-plots y funciones de distribución empírica en R.

Una advertencia antes de empezar: esta entrada es, probablemente, bastante más técnica que la media del blog. Si el lector entiende bien cómo se usa un QQ-plot, seguramente pueda ahorrarse la última parte, donde se intenta explicar cómo se construyen ese tipo de gráficos. Pero en cualquier caso, recomiendo leer la parte sobre funciones de distribución empírica de las muestras, porque esa es una herramienta interesante en si misma.

Y por si luego, metidos en faena, se nos olvida, una buena banda sonora para un discusión elevada como esta puede ser Skytalk, de Beats Antique.

Introducción

Ya hemos hablado de QQ-plots (quantile vs quantile plots) en el blog,en esta entrada sobre ANOVA unifactorial. Un QQ-plot es (aparte de algo que garantiza unas risas el primer día que lo mencionas en clase) una representación gráfica, que sirve para comparar dos distribuciones y ver si coinciden.  Allí vimos un ejemplo, que repito aquí:ANOVA-QQplot-ResiduosLa idea de este tipo de gráficos es que, si las dos distribuciones coinciden, veremos un gráfico muy parecido a una recta. Tanto más parecido cuanto mayor sea la coincidencia.

Cuando hacemos Inferencia, a menudo suponemos que los datos de la muestra proceden de una distribución de cierto tipo, por ejemplo de una distribución normal.  Y en otras ocasiones, al contrario, hemos recogido una muestra y sospechamos que no responde a una distribución normal. ¿Cómo podemos verificar esto? Hay, por supuesto, contrastes diseñados para esta tarea. Pero  la exploración gráfica, como siempre, es una gran ayuda. Y un QQ-plot es una de las formas más sencillas de obtener información en estas situaciones.

En R, obtener un QQ-plot es muy sencillo, y lo veremos en detalle más abajo. Pero en lugar de conformarnos con esto, quiero aprovecharlo como excusa para introducir  la idea de función de distribución empírica de una muestra, y su tratamiento en R. La excusa, desarrollada por completo, va así: “¿si R no nos ofreciera el QQ-plot, cómo podríamos hacerlo a mano?”. Para entender esto necesitamos esa idea de la función de distribución, que por un lado es una idea sencilla, y por otra parte es muy interesante por si misma. Veamos de qué se trata.

Función de distribución empírica de una muestra

En el fichero empiricalCDF-01.csv hay una muestra de datos, de tamaño n=10Es una muestra pequeña, pero la elegimos así para que las gráficas que vamos a dibujar resulten más claras. Más abajo trabajaremos con muestras de mayor tamaño. 

Uno de los valores de esa muestra es 4.6550. Supongamos que queremos responder a esta pregunta: si elegimos un valor de la muestra al azar (de manera que todos los valores tienen la misma probabilidad de ser elegidos) ¿cual es probabilidad de que el valor elegido sea menor o igual que 4.6550?  El procedimiento – bastante sencillo – para responder, es este:

  • Ordenamos los valores de la muestra. Se obtiene:

    x1 = 0.2849
    x2 = 0.5262
    x3 = 1.5600
    x4 = 2.7790
    x5 = 3.8230
    x6 = 4.6550
    x7 = 5.0150
    x8 = 5.3990
    x9 = 5.6570
    x10 = 9.9890
    

    Obsérvese que el valor 4.6550 ocupa la posición número seis en esta lista ordenada.

  • Puesto que son equiprobables, la probabilidad de elegir cada uno de ellos es 1/n, es decir, 1/10. La probabilidad de elegir un valor menor o igual que x6=4.6550 es, por lo tanto, 6/10. En general, si el valor xk ocupa la posición número k en la lista, la probabilidad de elegir un valor menor o igual que xk es k/n. Diremos que

F(xk)=\frac{k}{n}

Una vez entendido esto, la función de distribución empírica de la muestra F(x) se define de forma muy parecida, extendiendo esa idea a un número cualquiera, que no tiene por que formar parte de la muestra. Concretamente,

F(x)=\displaystyle\sum_{xk \leq x}\frac{1}{n}

Y, con menos notación pero tal vez más claridad, si x está, por ejemplo, entre el quinto y el sexto valor de una muestra (ordenada) de 17 valores, entonces F(x)=5/17.

En R disponemos de la función ecdf (empirical cumulative distribution function, ver aquí el fichero de ayuda de esa función), para obtener la función de distribución empírica de una muestra cualquiera. En nuestro caso, si el fichero  empiricalCDF-01.csv  está en el Escritorio de la máquina, para leer los datos, definir y la función y dibujarla (con algunos adornos), ejecutaríamos este código:

# En Windows, las siguientes lineas
# establecen el Escritorio del usuario actual como directorio de trabajo

(dirTrabajo=paste("C:/Users/",Sys.info()["user"],"/Desktop",sep=""))
setwd(dirTrabajo)

# En sistemas Mac o Linux alguna de estas lineas (descomentar) debería tener el mismo efecto.

# setwd("~/Desktop")
# setwd("~/Escritorio")

(datos=scan(file="empiricalCDF-01.csv",sep=" "))

# Aunque en general no es necesario, en este ejemplo vamos a reordenarlos

(datos=sort(datos))

# Ahora definimos y dibujamos la gráfica de la función de distribución empírica:

ECDF=ecdf(datos)
plot(ECDF,col="red",lwd=3,xlab="datos",ylab="")
points(datos,rep(0,n),col="blue",pch=20,cex=2)
segments(datos,rep(0,n),datos,sapply(datos,ECDF),col="blue",lwd=3)

Como puede verse, hemos llamado ECDF a la función. No es una buena idea llamarla F, porque R interpreta ese símbolo como una abreviatura del valor lógico FALSE, y se pueden generar errores inesperados. La gráfica que se obtiene es esta: empiricalCDF01Como se ve, es una gráfica en forma de escalera, ascendente de izquierda a derecha, con un escalón (todos de altura 1/n) sobre cada uno de los valores de la muestra, que sube desde 0 hasta 1 a medida que recorremos los valores de la muestra. Estas son características generales de cualquier función de distribución empírica (entrada en la Wikipedia, en inglés.).

La función ECDF que hemos definido en R es, por supuesto, un objeto de tipo function (usar class(ECDF) para comprobarlo, y ver esta entrada si se necesita más información sobre funciones). Así que podemos calcular ECDF(x) para cualquier valor x. ¿Cómo interpretar estos valores ECDF(x)? Sabemos como hacerlo para los valores de la muestra, pero ¿y para un valor x cualquiera? Si calculamos ECDF(x) y obtenemos 0.36, podemos interpretar que el 36% de los valores de la muestra son menores o iguales que x. Visto así, esta noción está estrechamente emparentada con la idea de percentil (o cuantil).

¿Qué tienen que ver estas funciones de distribución empíricas con las funciones de distribución teóricas de una variable aleatoria (de las que ya hemos hablado en la duodécima sesión)? Recordemos que si X es una variable aleatoria (discreta o continua), su función de distribución se define como:

F(x)=P(X\leq x)

Es decir, la probabilidad (según la distribuye la variable X) de la cola izquierda definida por el valor x. Por ejemplo, la siguiente gráfica muestra la función de distribución de una variable Z de tipo normal estándar N(0,1) que, en R, no es otra que nuestra vieja amiga pnormgrafica_pnorm

Como se ve, una función que sube, de izquierda a derecha, desde el 0 hasta el 1. ¿Qué sucede si tomamos una muestra de, pongamos 200 valores de una distribución normal N(\mu,\sigma), los tipificamos, después los usamos para definir una función de distribución empírica, y finalmente la superponemos sobre esta teórica? Para hacerlo, basta con ejecutar este código:

set.seed(2013)

#Generamos unos datos para que sirvan de ejemplo de cdf
n=200
datos=rnorm(n)
(xbar=mean(datos))
(desvt=sd(datos))
datosNormalizados=(datos-xbar)/desvt

ft=function(x){pt(x,df=2)}

plot(pnorm,from=-5,to=5,col="red",lwd=3)

plot(ft,from=-5,to=5,col="magenta",lwd=3,add=TRUE,lty=2)
ECDF=ecdf(datosNormalizados)
plot(ECDF,col="blue",lwd=1,xlab="datos",ylab="",add=TRUE)

Hemos incluido la orden set.seed(2013) para que el lector obtenga exactamente la misma muestra, y por tanto los mismos resultados, que nosotros. Si se desea obtener muestras normales distintas, comentar esa línea. Además, hemos incluido (a trazos) la gráfica de la función de distribución teórica de una t de Student con dos grados de libertad, para que quede claro que, en este caso, hay una clara diferencia con la muestra. El resultado es esta figura: ECDFvsDF-normal

Y como se ve, la distribución empírica, de una muestra de este tamaño, sigue muy fielmente a la distribución teórica. El mensaje que hay que extraer de este ejemplo es que, en general, para muestras grandes, la función de distribución empírica caracteriza bastante bien a la distribución teórica subyacente. Y como se trataba de tener una forma de comparar distribuciones, esto nos deja a las puertas del QQ-plot.

Cómo construir un QQ-plots en R “a mano”

Empecemos con una muestra (ordenada)

x_1 \leq x_2 \leq \cdots \leq x_k

suficientemente grande, que procede de una distribución normal N(\mu,\sigma). La tipificación consiste en calcular:

x\rightarrow \displaystyle\frac{x-\mu}{\sigma}= \displaystyle\frac{1}{\sigma}\cdot x - \displaystyle\frac{\mu}{\sigma}

Es decir, que desde el punto de vista matemático, es una transformación de la forma:

x\rightarrow y=a+b\cdot x

cuya gráfica de y frente a x es una recta. Es importante entender esto, porque esa es la razón última por la que en el QQ-plot esperamos ver algo parecido a una recta, si hemos partido de una muestra normal. Hemos empezado con los valores x, pero ahora, después de tipificar, tenemos unos valores

y_1 \leq y_2 \leq \cdots \leq y_k

que podemos considerar como una muestra de Z (la normal N(0,1)). Eso significa que su función de distribución empírica es la escalera determinada por los puntos:

(y_1,\frac{1}{n}), (y_2,\frac{2}{n}), \ldots , (y_n,\frac{n}{n})

Y hemos visto (esto es muy importante, revisar la última gráfica) que la escalera que forman estos puntos se parece mucho a la función de distribución teórica de Z, y que esa distribución teórica es lo que en R llamamos pnorm.

La función qnorm es la inversa de pnorm. Ahora, si empezamos con las probabilidades:

\frac{1}{n}, \frac{2}{n}, \ldots , \frac{n}{n}

y calculamos los valores:

y'_1 = qnorm(\frac{1}{n}), y'_2 = qnorm(\frac{2}{n}),…,y'_n = qnorm(\frac{n}{n}),

entonces es de esperar que los valores

y'_1, y'_2,\ldots , y'_n

sean muy parecidos a los y_i.

Esperamos que el lector no haya sufrido demasiado, y no se haya perdido en la discusión. Es un poco complicada, pero ya hemos llegado. Si los y'_i se parecen mucho a los y'_i, y la gráfica de los x_i frente a los y_i  era una recta, entonces la gráfica de los x_i frente a los y'_i se debe parecer mucho a una recta. Esta última gráfica es nuestro QQ-plot “artesanal”.

A continuación incluyo el código en R para este proceso que hemos descrito. Pero antes, es necesario hacer algunos comentarios de advertencia: el código usa, además,  la función qqnorm de R, para comparar el resultado con el proceso artesanal que hemos descrito. Y hay buenas razones para que ambos no coincidan exactamente. En particular, nosotros hemos usado 1/n,…,n/n como probabilidades, mientras que R usa un conjunto distinto (ver el Capítulo 2 del libro Testing for Normality, de H.C.Thode. ) Además, R coloca la muestra original en el eje vertical (es la práctica habitual en un QQ-plot), y nosotros lo hemos colocado en el horizontal porque, de lo contrario, la discusión anterior habría sido aún más complicada de entender… eso nos obliga a hacer algunas manipulaciones a partir del resultado de qqnorm, para colocar las cosas como queremos.

set.seed(2013)
#Generamos una muestra con datos de una distribucion normal,
n=200
x=rnorm(n,mean=7,sd=1.5)

# la ordenamos,
x=sort(x)

# calculamos su media y cuasidesviacion tipicas,
(xbar=mean(x))
(desvt=sd(x))

# y hacemos una especie de tipificacion.
y=(x-xbar)/desvt

# Ahora empezamos con las probabilidades,
probs=(1:n)/n

# y las invertimos usando qnorm.

y2=qnorm(probs)

# El QQ-plot artesanal es la grafica x frente a y2.
# Para comparar con lo que hace R, usamos qqnorm

qqplot=qqnorm(x,plot.it=FALSE)

# Ademas R invierte los ejes (cambia x por y)
# asi que tenemos que colocarlos al reves para poder comparar
qqx=qqplot$x
qqy=qqplot$y

# En azul el QQ-plot artesanal, y en rojo el de R.
plot(x,y2,col="blue",pch=3)
points(qqy,qqx,col="red") # Observar el cambio de orden en los ejes

Y el resultado es esta gráfica. qqplot-artesanalen la que el resultado de R se muestra en rojo, y el nuestro en azul. Como puede verse, ambos son muy similares, y ambos además se parecen mucho a una recta, como esperábamos al haber partido de una muestra normal.

Muchas gracias por la atención, especialmente en una entrada larga y seguramente más técnica que la media.