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.

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s