Publicados los Tutoriales del 7 al 10 del curso de Introducción a la Estadística.

Ya están disponibles para descarga en la página del curso. El tutorial 10, por las urgencias de las clases en la universidad, nos parece un poco más provisional que los anteriores, pero en cualquier caso contiene la información básica necesaria.

Anuncios

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.

Más Coursera

Una entrada mínima para, como ya anuncié que haría, mantener actualizada la información sobre los cursos que me parecen más tentadores.  Aparte de los que aparecían en la anterior entrada sobre este tema, hoy añado:

  • Case-Based Introduction to Biostatistics, de Scott L. Zeger, en la  Johns Hopkins University. La presentación del curso dice, entre otras cosas, ”The course objective is to enable each student to enhance his or her quantitative scientific reasoning about problems related to human health”. 
  • Computing for Data Analysis, de Roger Peng, también en la Johns Hopkins. Para quienes ya saben algo de R (es decir, todos los lectores del blog), este puede ser un curso interesante, con vistas a ir un paso más allá.

 

RStudio

En las entradas precedentes del blog  he usado, hasta ahora, R-Commander como interfaz que simplifica nuestra relación con R, especialmente al principio. Y sigo pensando que es un buen programa para principiantes, y continuaré explicando las opciones que ofrece.  Tengo además una deuda de gratitud con ese programa, porque mis primeros pasos con R habrían sido mucho más difíciles si no hubiera descubierto R-Commander (y el paquete R-UCA de la Universidad de Cádiz).

Sin embargo, en los últimos tiempos me he dejado convencer por otra interfaz gráfica de R, que se llama R-Studio (el enlace es la página principal del proyecto). Para presentarlo, aquí tenemos una captura de pantalla de RStudio, en una sesión de trabajo con el código de la entrada sobre ANOVA para diseños factoriales (con tres factores en este caso).

RStudio001

RStudio ofrece muchas herramientas para el usuario avezado de R, convirtiéndose en un entorno de trabajo muy productivo. Esa es la mayor diferencia con R-Commander, que tal vez está más orientado a proporcionar al usuario una experiencia  similar a la de los programas comerciales clásicos, donde los análisis están disponibles mediante un sistema de menús y cuadros de diálogo, como este:

Rcmdr-UnMenu

En RStudio no encontraremos nada de eso. No hay un menú que nos diga cómo hacer un contraste de hipótesis, un intervalo de confianza. Para eso está R. En ese sentido, RStudio da por sentado que nosotros sabemos pedirle a R que haga esos análisis. Y por eso, para un principiante, puede resultar algo más intimidatorio.

Para el usuario con algo más de experiencia, como decíamos, la situación se invierte, y RStudio se muestra como una de las mejores opciones para el trabajo diario con R. Entre otras, quiero destacar estas características:

  • la posibilidad de trabajar con varios ficheros de comandos R simultáneamente, agrupados por pestañas (en la ventana superior izquierda, en la anterior captura de pantalla).
  • el editor de comandos es, en si mismo, una de las grandes virtudes de RStudio. Es realmente un editor pensado por y para programar en R. No echaremos en falta un editor externo, al que en general teníamos que recurrir con R-Commander.
  • autocompletado de comandos. Si no recordamos la sintaxis exacta de un comando de R, basta con escribir las primeras letras, pulsar la tecla Tab, y RStudio nos ofrece información sobre ese comando. Pero también sobre las variables, ficheros y demás objetos que hayamos creado.
  • gestión de los gráficos de R mucho más eficaz que, por ejemplo, en R-Commander.
  • facilidad de gestión para, por ejemplo, instalar paquetes, consultar ayuda, etcétera.
  • si se dispone de varias versiones de R instaladas en la máquina, es fácil seleccionar la que utilizará en cada sesión.

Pero además, RStudio ofrece la posibilidad de instalación en un servidor (bajo Linux), para posteriormente acceder al programa a través de una interfaz web, desde cualquier máquina, aunque no esté instalado R. Esa interfaz web es esencialmente idéntica a la que hemos mostrado más arriba, como puede verse en esta figura:

RStudio002-InterfazWeb

Esto hace posible (aunque desde luego, no muy cómodo) acceder a RStudio desde un teléfono móvil. Y desde luego es especialmente interesante para usuarios de tablets, (tanto Ipads como Android).

Para finalizar esta primera visita, en la página del proyecto RStudio se puede ver otro producto, Shiny, del que hablaremos en una futura entrada y que permite convertir los análisis que hacemos con R en páginas web interactivas, para que cualquier pueda experimentar con esos análisis, y modificar los parámetros utilizando deslizadores, cajas de diálogo, etc.  Una muestra, en este enlace.

Shiny  es una forma de avanzar en la idea de mostrar a los demás nuestro trabajo, cada vez de forma más trasparente. Una idea que está en la base del movimiento hacia lo que se ha llamado Investigación reproducible (Reproducible Research, página en la Wikipedia). En esencia, la idea es que los trabajos de investigación, a la hora de publicar sus resultados, debería acompañar la publicación tradicional con el conjunto completo de datos experimentales y los ficheros de código del software que se haya utilizado para el análisis de esos datos.

R tiene una herramienta muy poderosa para contribuir a esta idea, llamada SWeave. Y en esta entrada no puedo dejar de mencionar que RStudio (a diferencia de R-Commander) ofrece mucha ayuda para el trabajo con SWeave. Pero ese un tema que merece ser tratado con mucho más detenimiento. En una próxima entrada volveremos sobre la investigación reproducible (y sus conexiones con la política y la economía, nada menos), sobre SWeave, y sobre la noción de Programación Literaria, que hizo posible todas estas ideas. Y (spoiler alert para profesores) sobre la forma de conectar todo esto con Moodle y usarlo en nuestras clases.

Gracias por la atención.

Tablas de contingencia y contraste de independencia (test chi cuadrado)

Contraste de independencia

El Barómetro del CIS  (Centro de Investigaciones Sociológicas) permite obtener datos sobre las creencias religiosas de la población en España. Una pregunta que puede interesarnos es ¿hay alguna diferencia al respecto entre hombres y mujeres? Vamos a utilizar los datos del Barómetro para intentar contestar esa pregunta.

Por ejemplo, en el mes de enero de 2013 el Barómetro} recoge las respuestas de n=2452 personas sobre sus creencias religiosas. (En realidad son 2483, pero para simplificar vamos a eliminar de nuestra consideración a las 19 mujeres y a los 12 hombres que decidieron no contestar). Observa que, como de costumbre, vamos a usar n para el número total de personas encuestadas. Agrupamos a todos los creyentes de distintas religiones por un lado y a los que se declaran no creyentes o ateos por otro. Y así tenemos una tabla de doble entrada:

Hombres Mujeres Total
Creyentes ?? ?? 1864
No creyentes ?? ?? 588
Total 1205 1247 2452

En la que hemos dejado pendientes los valores interiores de la tabla, y sólo se muestran los valores marginales. Si las creencias religiosas no dependieran del género, ¿cuáles son los valores que esperaríamos ver en esta tabla? Los que se deducen de los marginales, claro. Por ejemplo, el número de hombres creyentes se puede calcular multiplicando el número de hombres por la proporción de creyentes en la población. Esto es:

1205 \cdot \frac{1864}{2452}\approx 916

Con R podemos obtener estos valores esperados fácilmente. Como es la primera vez que hacemos algo así en el blog, vamos a dar los detalles. Empezamos con dos vectores, los de los valores marginales:

creencia=c(1864,588)
genero=c(1205,1247)

Y ahora los multiplicamos como matrices, trasponiendo el vector creencia, que va en vertical en la tabla, y que por tanto hace el papel de vector columna. Para trasponer un vector basta con usar la función t(). El producto matricial se representa mediante el símbolo %*%, así que la tabla (matriz producto) es (mostramos la salida):

> creencia %*% t(genero)
        [,1]    [,2]
[1,] 2246120 2324408
[2,]  708540  733236

Naturalmente, aún tenemos que dividir por el total de personas (el número n). El resultado completo es esta tabla de valores esperados:

> (tablaEsperada=creencia %*% t(genero) / sum(genero))
         [,1]     [,2]
[1,] 916.0359 947.9641
[2,] 288.9641 299.0359

Está bien, pero sería más fácil de interpretar si añadimos nombres a las filas y columnas, y además incluimos los valores marginales con los que habíamos empezado. Hacemos esto con los comandos:

colnames(tablaEsperada)=c("H","M")
rownames(tablaEsperada)=c("CREY","NO_CREY")
(addmargins(tablaEsperada))

y el resultado es mucho más agradable de ver:

                H         M  Sum
CREY     916.0359  947.9641 1864
NO_CREY  288.9641  299.0359  588
Sum     1205.0000 1247.0000 2452

De acuerdo, eso es lo que esperaríamos, si las creencias y el género fueran independientes. ¿Cuáles son los datos reales que se obtuvieron en el Barómetro? Aquí está la tabla de valores observados:

Hombres Mujeres Total
Creyentes 849 1015 1864
No creyentes 356 232 588
Total 1205 1247 2452

Introducimos esta tabla en R directamente usando cbind, y le añadimos la información necesaria para hacerla más útil:

(tablaObservada=rbind(c(849,1015),c(356,232)))
colnames(tablaObservada)=c("H","M")
rownames(tablaObservada)=c("CREY","NO_CREY")
(addmargins(tablaObservada))

El resultado es:

           H    M  Sum
CREY     849 1015 1864
NO_CREY  356  232  588
Sum     1205 1247 2452

Comparando las dos tablas, resulta evidente que los valores observados no coinciden con los  esperados. De hecho, el número de hombres no creyentes es más alto de lo que habíamos estimado a partir de la población en conjunto (y,  lógicamente, el número de mujeres no creyentes es más bajo que la estimación). Pero ese número de hombres no creyentes, ¿es significativamente más alto?

Para contrastar la hipótesis de independencia entre género y creencias, usaremos la función chisq.test de R (se muestra la salida):

> (chisq.Observada=chisq.test(tablaObservada,correct=FALSE))

	Pearson's Chi-squared test

data:  tablaObservada
X-squared = 40.2253, df = 1, p-value = 2.263e-10

Como en otros casos, hemos incluido la opción correct=FALSE para que R (lo haga un poco peor de lo que sabe y) nos devuelva el mismo valor que aparecería en un curso de Introducción a la Estadística.  Se recomienda leer la ayuda de chisq.test para ver lo que implica esa opción.

En este caso el p-valor es muy pequeño, así que con esos datos deberíamos rechazar la hipótesis nula, y concluir que el género y las creencias no son independientes.

El test \chi^2 que estamos usando consiste, para este ejemplo, en calcular el siguiente estadístico:

X^2=\dfrac{(o_{11}-e_{11})^2}{e_{11}}+\dfrac{(o_{12}-e_{12})^2}{e_{12}}+\dfrac{(o_{21}-e_{21})^2}{e_{21}}+\dfrac{(o_{22}-e_{22})^2}{e_{22}}.

donde los valores o_{ij} son los de la tabla observada, y los e_{ij} son los valores esperados. Cada término de la forma

\frac{o_{ij}-e_{ij}}{e_{ij}}

es un residuo, y el estadístico es la suma de los residuos al cuadrado. En R,  al aplicar el test, el programa devuelve,  dentro de la variable chisq.Observada, una lista con todo lo  necesario. Accedemos a las componentes de esa lista con el operador $. Por
ejemplo, las tablas de valores esperados y observados se obtienen así (se muestra la salida):

> (E=chisq.Observada$expected)
               H        M
CREY    916.0359 947.9641
NO_CREY 288.9641 299.0359

> (O=chisq.Observada$observed)
          H    M
CREY    849 1015
NO_CREY 356  232

Y a partir de aquí el cálculo de los residuos es inmediato:

> (residuos=(O-E)/sqrt(E))
                H         M
CREY    -2.214885  2.177266
NO_CREY  3.943532 -3.876553

Lo hemos hecho de esta manera para que quede claro cómo se definen, pero los residuos también se pueden obtener directamente con:

chisq.Observada$res

De cualquier manera, para calcular a mano el estadístico a partir de los residuos, podemos ahora elevar los residuos al cuadrado y sumarlos usando addmargins, de esta manera:

> (tablaChi2=residuos^2)
                H         M
CREY     4.905714  4.740486
NO_CREY 15.551448 15.027663

> addmargins(tablaChi2)
                H         M      Sum
CREY     4.905714  4.740486  9.64620
NO_CREY 15.551448 15.027663 30.57911
Sum     20.457163 19.768148 40.22531

El estadístico es el valor de la esquina inferior derecha de esta tabla:

> (estadistico=addmargins(tablaChi2)[3,3])
[1] 40.22531

Y podemos ver que coincide con lo que obtuvimos antes usando chisq.test. El p-valor se obtiene a partir del estadístico usando pchisq como en cualquier contraste (con un grado de libertad, porque es una tabla 2×2). Usamos la cola derecha, porque los valores del estadístico que contradicen a la hipótesis nula se sitúan ahí:

> (pValor=1-pchisq(estadistico,df=1))
[1] 2.262972e-10

Y como antes, el p-valor coincide con lo que obtuvimos en chisq.test.

Una posible representación gráfica de este contraste es mediante un gráfico de columnas apiladas (la altura indica el porcentaje), con una columna por  género. Para obtener esa representación empezamos por fabricar una tabla de proporciones (porcentajes), en la que los porcentajes son sobre el total de la muestra. Añadimos las sumas marginales para que sea más fácil entender esta tabla (recordamos que se muestra la salida):

> tablaProporciones1=prop.table(tablaObservada)
> addmargins(tablaProporciones1)
                H          M       Sum
CREY    0.3462480 0.41394780 0.7601958
NO_CREY 0.1451876 0.09461664 0.2398042
Sum     0.4914356 0.50856444 1.0000000

Para que cada columna de la tabla nos de las alturas de las columnas del gráfico, tenemos que dividir la tabla por las sumas por columnas, calculadas con colSums. Además, vamos a redondearla para mejorar la presentación.

> (tablaProporciones2=round(tablaObservada/colSums(tablaObservada),digits=1))
          H   M
CREY    0.7 0.8
NO_CREY 0.3 0.2

Ya podemos usar barplot para dibujar el grafico, con colores y etiquetas en los ejes. Aunque con esto ya se obtiene un gráfico adecuado, la interpretación de un gráfico como este es mucho más fácil, si se rotulan las columnas con el porcentaje de una de las categoría en cada columna. Primero fabricamos los rótulos, y luego los colocamos en sus posiciones con text. Estos son los comandos:

(bp_t2=barplot(tablaProporciones2,col=c("lightseagreen","darkgoldenrod1"),xlab="Género",ylab="Creyente/No creyente",main="Creencias religiosas por género",font=2))

rotulos=paste(format(100*signif(tablaProporciones2[1,],digits=2)),"% creyentes",sep="")

text(bp_t2, tablaProporciones2[1,]+0.05,rotulos , xpd = TRUE, font=2)

Y este es el resultado:

20130514-CreenciasReligiosasPorGenero

Un problema como este, con una tabla de contingencia 2×2, también se puede abordar con un contraste de igualdad de proporciones para dos muestras, como los que vimos en esta entrada. Vamos a hacerlo aquí (de nuevo, con correct=FALSE), para que se vea que R devuelve el mismo valor del estadístico, y el mismo p-valor.

> (creyentesObs=tablaObservada[1,])
   H    M
 849 1015
> prop.test(creyentesObs,genero,correct=F)

	2-sample test for equality of proportions without continuity correction

data:  creyentesObs out of genero
X-squared = 40.2253, df = 1, p-value = 2.263e-10
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.14300581 -0.07577254
sample estimates:
   prop 1    prop 2
0.7045643 0.8139535

Por último, si tenemos dudas sobre las hipótesis de normalidad subyacentes al test \chi^2, podemos hacer un test no paramétrico para confirmar lo anterior.

> fisher.test(tablaObservada)

	Fisher's Exact Test for Count Data

data:  tablaObservada
p-value = 2.745e-10
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.4490685 0.6611693
sample estimates:
odds ratio
 0.5452314

que, en este ejemplo, confirma que el resultado del test es significativo.

Código en R y observaciones finales

En este enlace hay un fichero con el código R que hemos usado en este ejemplo. Naturalmente, el caso 2×2 no agota las posibilidades del test \chi^2 de independencia. Estoy terminando de escribir el correspondiente capítulo de los apuntes de teoría (avisaré en un post cuando acabe), que incluye ejemplos y código para el caso general, con tablas de contingencia con un número mayor de filas y columnas.

Aparte de esto, y como se explica en ese capítulo, existe otra variante del test, que se aplica a problemas de bondad de ajuste, es decir a la pregunta de si los datos observados se corresponden con una distribución teórica. Concretamente, ejemplos de este tipo de preguntas son: ¿cómo podemos saber si un dado está cargado? O ¿corresponden estos datos a una población normal?.  Próximamente le dedicaremos una entrada del blog a este tema. Y, puesto que aquí ha surgido, siquiera sea de pasada, el tema de las operaciones con matrices, también le dedicaremos una entrada a ese tema.

Gracias por la atención.

Y una última novedad (si P. Krugman se atreve, no voy a  ser yo quien se quede atrás): Recomendación musical de esta entrada Headlock de Imogen Heap,

GeoGebra y la enseñanza de la Estadística

Con esta entrada inauguro una nueva sección del blog. Hasta ahora la mayoría de  las entradas tenían a mis alumnos como público objetivo principal (y por extensión a los alumnos de cursos de Introducción a la Estadística basados en R). Pero vamos a ir tratando, en esta y en futuras entradas algunos temas que van dirigidos a los profesores de ese tipo de ese cursos. Naturalmente, cualquiera puede leerlas, pero creo que es de rigor avisar de cuál es el lector que tengo en mente al escribir.

Por razones que explicaré más abajo, hasta ahora me he resistido todo lo que he podido a hablar de GeoGebra (enlace a la página oficial; es interesante además ver la página de la Wikipedia) en este blog. ¡Y sólo aquí! En realidad hace tiempo que hablo y no callo sobre el tema. Es, sin la menor duda, uno de los mejores programas para la enseñanza de las Matemáticas, desde la enseñanza primaria hasta la universidad. Creado por Markus Hohenwarter en 2001, ha ido ganando usuarios a un ritmo vertiginoso, sólo superado por el aún más vertiginoso ritmo al que el programa ha ido ganando en capacidades. Inicialmente, GeoGebra parecía ser uno más de los entornos de Geometría Dinámica (como Cabri, Cinderella o Sketchpad, por citar algunos). Pero la genialidad de su creador consistió en incorporar una doble representación de los objetos matemáticos, geométrica, pero también algebraica (de ahí el nombre). En la siguiente figura puede verse la interfaz habitual de GeoGebra (en la versión 4.2, la última oficial a fecha de hoy), con algunos objetos: una función, una recta, una circunferencia, puntos, etc. Y a la vez, en la vista algebraica, a la izquierda, aparece la descripción de esos objetos en una notación que, se acerca mucho (más con cada nueva versión de GeoGebra) a la notación algebraica habitual.

GeoGebraInterfaz

La vista gráfica de GeoGebra es una “pizarra con superpoderes”, puesto que los objetos tienen un carácter dinámico. Pero ante todo, el fundamento de su utilidad para la enseñanza es que mantiene la unidad de la representación habitual de los objetos matemáticos, a la que todos nos hemos acostumbrado en los libros de texto y en las clases. El usuario, con conocimientos de Matemáticas, que mira la ventana gráfica de GeoGebra, no necesita “traducción” de ningún tipo para reconocer la situación que tiene delante, porque el lenguaje es el habitual.

Conviene además recordar que GeoGebra es gratuito, de código abierto, y que cuenta con una auténtica legión de seguidores entusiastas.

¿Y GeoGebra tiene algo que ver con la Estadística? Hace unos años la respuesta hubiera sido: no mucho, salvo a costa de que el usuario hiciera bastante parte del trabajo por su cuenta. Pero ya decíamos que GeoGebra no deja de crecer. Y en las últimas versiones se han ido incorporando cada vez más características, que convierten al programa en una herramienta a tener muy en cuenta en la enseñanza de la Estadística. Vamos con algunas de ellas.

Los creadores de GeoGebra tuvieron, en algún momento de inspiración, la original idea de añadir, junto a las vistas geométrica y algebraica, una “vista de hoja de cálculo”, que se aprecia en la parte derecha de esta figura:

GeoGebraHojaCalculo

La hoja de cálculo es un primer paso hacia la Estadística, como bien sabemos. En el caso de GeoGebra, no obstante, lo que hace especial la hoja de cálculo es que el contenido de las celdas pueden ser objetos geométricos (puntos, rectas, curvas, etcétera), a los que se pueden aplicar operaciones en la forma habitual en una hoja de cálculo. Por ejemplo, si tenemos una columna de puntos, es posible marcar esa columna y trasladar todos los puntos simultáneamente mediante un cierto vector, obteniendo los puntos resultantes en otra columna. Es una forma muy interesante de trabajar, que da mucho juego en la enseñanza. Pero que nos desvía un poco de la temática del blog. La hoja de cálculo es una condición necesaria, pero no suficiente, para que GeoGebra juegue un papel en la enseñanza de la Estadística.

El paso definitivo (al menos para mi), sin embargo, se ha dado en las últimas versiones (esencialmente, a partir de la versión 4).  Ahora en uno de los menús de GeoGebra disponemos de la herramienta Cálculo de Probabilidades:GeoGebraMenuProbabilidad

Al pulsar sobre esa herramienta se abre una nueva ventana, que ofrece una gran cantidad de posibilidades para el trabajo con las distribuciones de probabilidad (discretas y continuas) más habituales en los cursos de introducción a la Estadística:

GeoGebraVentanaCalculosProbabilidad

Al abrir la ventana nos recibe la distribución Z (¿Cuál, si no?) Podemos usar los formularios y botones de esta ventana para estudiar problemas de probabilidad directa o inversa (cálculo de cuantiles o percentiles), tanto con colas como con intervalos. Además, desde luego, podemos modificar la media \mu y desviación típica \sigma de la normal. También podemos, en lugar de la función de densidad, pulsar el pequeño botón que aparece al lado del nombre de la distribución, y cambiar la función de densidad por la de distribución, que se muestra en esta figura:

GeoGebraNormalFuncionDistribucion

Esto mismo es posible con muchas otras distribuciones:

  • continuas: t de Student, \chi^2, F de Fischer, exponencial, Cauchy, Weibull, Gamma, LogNormal y Logística.
  • Discretas: Binomial, Pascal, Poisson, Hipergeométrica.

Además, y de acuerdo con el espíritu dinámico de GeoGebra, podemos deslizar los puntos del eje x en la figura, y ver esos cambios reflejados automáticamente en los correspondientes valores de probabilidad. Dejo aquí enlazado un pequeño vídeo con una captura de ese tipo de manipulaciones dinámicas, que creo que ilustra los útil que puede llegar a ser GeoGebra en nuestras clases (se recomienda ver a pantalla completa).

Distribuciones de probabilidad con GeoGebra

Finalmente, en este breve recorrido de las posibilidades de GeoGebra, no quiero olvidarme de mencionar que (al margen de esta herramienta específica para las distribuciones más conocidas) GeoGebra ha incorporado una gran cantidad de comandos relacionados con la Probabilidad y la Estadística. En estas dos capturas de pantalla pueden verse las listas de comandos de GeoGebra para esos dos temas. De hecho la lista de comandos estadísticos es ya tan extensa que no cabe en esta captura…

GeoGebraComandosProbabilidad                                        GeoGebraComandosEstadistica

En cualquier caso, el lector interesado podrá encontrar mucha más información aquí.

Gracias por la atención.

Anova para diseños factoriales con 2 factores (efectos fijos).

Editado el 10/5/2013 para añadir, al final,  código del ANOVA con 3 factores

Otra larga pausa sin publicar. Como en la canción de Lennon, “Life is what happens to you while you’re busy making other plans“.  Tengo varias entradas del blog y dos capítulos más de los apuntes, preparados pero que no acaban de ver la luz, porque varias cosas se cruzan por el camino. Y la entrada de hoy es, de hecho, una de esas cosas que se cruzan. Varios alumnos me han pedido ayuda en relación con un análisis de datos que involucra una tabla Anova con dos factores (y efectos fijos). En mi – ingente – lista de tareas para el blog está previsto tratar ese tema con más profundidad, a lo largo de varias entradas, y con algún artista invitado, si puede ser.  Porque las generalizaciones del ANOVA que vemos en el curso son, de hecho, una de las técnicas estadísticas más utilizadas en la práctica. Pero la realidad impone sus ritmos, así que vamos a hacer una introducción, “quick and dirty” al tema; el lector queda avisado: esta no es, ni mucho menos, la historia completa.

Lo que sigue se inspira en el Ejemplo 12A del Capítulo 12 (pág. 334) de “Biostatistical Design and Analysis Using R“, de M. Logan. A su vez, ese ejemplo se basa en los datos del Ejemplo 9.4, del Capítulo 9 (pág. 224) de “Experimental Design and Data Analysis for Biologists“, de G.P. Quinn y M.J. Keough. En este ejemplo, se usan los datos de un estudio (Quinn, G. P. (1988). “Ecology of the intertidal pulmonate limpet Siphonaria Diemenensis”, Journal of experimental marine biology and ecology, vol. 117, no2, pp. 115-136.) sobre los efectos que tienen

  1. la estación del año y
  2. la densidad de ejemplares adultos

sobre la producción de huevos en la especie Siphonaria Diemenensis.

El fichero de datos que vamos a usar es quinn.csv. Lo leemos con R (esta vez, como novedad, en lugar de descargarlo primero, lo vamos a descargar directamente desde R) y le echamos un vistazo, no es excesivamente largo. 

> fichero="http://www.fernandosansegundo.es/Estadistica/quinn.csv"
> (quinn=read.table(fichero,header=T,sep=","))
 DENSITY SEASON EGGS
1 8 spring 2.875
2 8 spring 2.625
3 8 spring 1.750
4 8 summer 2.125
5 8 summer 1.500
6 8 summer 1.875
7 15 spring 2.600
8 15 spring 1.866
9 15 spring 2.066
10 15 summer 0.867
11 15 summer 0.933
12 15 summer 1.733
13 30 spring 2.230
14 30 spring 1.466
15 30 spring 1.000
16 30 summer 1.267
17 30 summer 0.467
18 30 summer 0.700
19 45 spring 1.400
20 45 spring 1.022
21 45 spring 1.177
22 45 summer 0.711
23 45 summer 0.356
24 45 summer 0.711

Este conjunto de datos es típico de un ANOVA con dos factores (y efectos  fijos). Tenemos una variable (respuesta) cuantitativa, que es la producción de huevos (la tercera columna de los datos). Y queremos estudiar la relación de esa variable con dos variables (explicativas), que son la la densidad de adultos y la estación (primera y segunda columnas, respectivamente).  

Se trata de dos variables cualitativas o factores. Esto es especialmente evidente en el caso de la variable SEASON, que tiene dos niveles (spring y summer). Podemos confirmarlo preguntándole a R de qué tipo de objeto se trata, así:

> class(quinn$SEASON)
[1] "factor"

Por otro lado, si hacemos lo mismo con DENSITY:

> class(quinn$DENSITY)
[1] "integer"

veremos que R ha entendido esta variable como un vector de enteros. Pero nosotros queremos tratar a DENSITY como un factor con cuatro niveles, que son 8, 15, 30 y 45. Esta es una situación típica, cuando tenemos un factor cuyos niveles se representan numéricamente. Para asegurarnos de que las cosas son como queremos que sean, hacemos:

> quinn$DENSITY = factor(quinn$DENSITY)
> class(quinn$DENSITY)
[1] "factor"

y, como puede verse, DENSITY pasa a ser considerado como un factor. Ya he avisado de que no vamos a tratar el tema en profundidad, pero merece la pena subrayar que, como se explica en el texto de Quinn y Keough, estamos considerando ambos como efectos fijos (y no como efectos aleatorios). En particular, en el diseño de este experimento, esos valores se han fijado por el experimentador  (si alguien quiere repetir este experimento, volverá a utilizar los valores 8, 15, 30 y 45 de DENSITY, y no otras densidades al azar; con densidades al azar sería otro diseño experimental). Para entender con más detalle esta diferencia entre efectos fijos y aleatorios recomiendo:

(a) leer los libros citados (o, para mis alumnos de Alcalá, el capítulo 10 de “Estadística para Biología y Ciencias de la Salud”, de S.Milton, que está disponible en nuestra biblioteca)

(b) esperar a que terminemos de escribir el capítulo correspondiente de los apuntes :).

¡No es una cuestión menor! El análisis de los datos es distinto dependiendo de esta diferencia entre efectos fijos y aleatorios. Vamos a ocuparnos del caso de efectos fijos porque es el más sencillo para empezar.

El modelo que vamos a estudiar, la relación de la variable (respuesta) EGGS con las variables (explicativas) DENSITY  y SEASON, se describe así en R:

modelQuinn=lm(EGGS ~ DENSITY * SEASON,data=quinn)

Esta notación es una generalización de la notación ~ que ya vimos al tratar el análisis unifactorial.  El asterisco * que aparece entre DENSITY y SEASON le indica a R que estamos interesados en la posible interacción entre esas dos variables. Volveremos sobre esto más abajo.

Otra pregunta importante que tenemos que hacernos, al principio del análisis, es si el diseño es equilibrado (balanced). Es decir, si todas las muestras son del mismo tamaño. Debemos tener en cuenta, entre otras consideraciones, que las hipótesis que justifican la validez del ANOVA son especialmente  sensibles a la existencia de muestras de tamaños distintos (como ya mencionamos al hablar  del análisis unifactorial). En nuestro ejemplo podemos comprobar si el diseño es equilibrado así:

> replications(modelQuinn,data=quinn)
DENSITY SEASON DENSITY:SEASON
6 12 3
> (isBalanced=!is.list(replications(modelQuinn, data=quinn)))
[1] TRUE

La primera instrucción usa la función replications para obtener un vector, o una lista, de las veces que cada nivel de cada factor aparece repetido en el experimento.  De hecho,  si el experimento es no-equilibrado, la función replications devuelve como resultado una lista, no un vector. La segunda instrucción utiliza esto para preguntarle a R si el diseño es equilibrado o no. Para ver esto con algo más de detalle, en este enlace hay un fichero de instrucciones que, al ejecutarlo, descargará una versión modificada de los datos, no equilibrada. Es una buena idea ejecutar este código y ver la diferencia en los resultados de replications.

Ahora que ya sabemos que el diseño es equilibrado, construimos la tabla ANOVA:

> quinn.aov = aov(EGGS ~ SEASON * DENSITY, data = quinn)
> summary(quinn.aov)
Df Sum Sq Mean Sq F value Pr(>F)
SEASON 1 3.250 3.250 17.842 0.000645 ***
DENSITY 3 5.284 1.761 9.669 0.000704 ***
SEASON:DENSITY 3 0.165 0.055 0.301 0.823955
Residuals 16 2.915 0.182
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

El resultado (y los propios comandos que hemos utilizado) son similares a los del caso unifactorial. Hay una línea para cada factor, y una para el término de interacción (SEASON:DENSITY). Los p-valores, y los códigos de asteriscos que los acompañan, nos indican que:

  • se obtienen resultados significativos para las dos variables SEASON y DENSITY, cada una por separado. Es decir, que la variable EGGS depende, en efecto, de estas dos variables.
  • el término de interacción no es significativo: no hay interacción entre las dos variables.

Podemos ilustrar esto último gráficamente mediante la función interaction.plot de R:

interaction.plot(quinn$DENSITY, quinn$SEASON, quinn$EGGS)

que produce este resultado:ANOVA2-InteractionPlot

El factor DENSITY se muestra en el eje horizontal, la respuesta EGGS en el eje vertical, y cada uno de los dos factores de SEASON se muestra como una línea que conecta los correspondientes valores. En un caso ideal, la ausencia completa de interacción correspondería a dos líneas perfectamente paralelas en este gráfico (y la existencia de interacción es evidente cuando las líneas se cruzan). En la práctica, el paralelismo en general está lejos de ser perfecto. No obstante, la gráfica de este ejemplo sí parece indicar que no existe interacción entre ambos factores. Y en cualquier caso, tenemos los resultados de la tabla ANOVA para corroborarlo.

El ejemplo 12B del libro de Logan ilustra otro experimento de Quinn en el que, para las mismas variables, de hecho, existe interacción. Para no repetir toda la discusión, vamos a enlazar un fichero de código R que corresponde al mismo análisis, pero para ese segundo conjunto de datos.  En este caso la tabla ANOVA es:

Df Sum Sq Mean Sq F value Pr(>F)
SEASON 1 17.148 17.148 119.845 1.34e-07 ***
DENSITY 2 4.002 2.001 13.984 0.000732 ***
SEASON:DENSITY 2 1.691 0.845 5.908 0.016363 *
Residuals 12 1.717 0.143

que muestra que el término de interacción es, en este caso, significativo. La gráfica obtenida con interaction.plot es esta: ANOVA2-InteractionPlot-1

en la que la falta de paralelismo entre ambas líneas resulta bastante más evidente.

Volviendo al ejemplo 12A, otra representación gráfica interesante se consigue utilizando el paquete effects. Sin entrar en detalles, el código para este ejemplo sería:

library(effects)
Quinn.eff=effect("DENSITY*SEASON",modelQuinn,confidence.level=0.90)
plot(Quinn.eff)

y el resultado es:

ANOVA2-effectPlotSe trata esencialmente de otra manera de presentar los resultados que ya obtuvimos en el interaction.plot, junto con intervalos de confianza para la media en cada uno de los niveles del factor DENSITY.

Añadido el 10/5/2013: varios de mis alumnos me han explicado que en sus problemas aparecen tres factores. No puedo entrar aquí en muchos detalles, que quedan para una  futura entrada del blog, pero voy a dejar enlazado un fichero de código R (que carga los datos necesarios) correspondiente a un ejemple de este tipo de análisis (ANOVA, 3 factores, efectos fijos, diseño equilibrado). El problema consiste en estudiar el rendimiento de tres variedades de centeno (I,II,III), con dos fertilizantes (A y B), y dos niveles de agua

La tabla de datos original del problema era esta

I II III
A B A B A B
Bajo 18.6 10.2 19.4 14.7 15.9 20.9
18.8 10.5 18.9 17.8 16.5 21.0
15.8 8.5 18.4 15.6 17.2 21.3
17.9 10.6 18.9 16.5 16.5 20.4
Alto 12.3 12.3 13.2 10.2 19.4 13.2
10.9 10.7 15.5 8.5 21.2 13.5
10.7 12.0 11.0 7.1 20.4 12.8
11.8 11.8 11.8 6.5 20.3 15.2

Y una de las primeras tareas del código R consiste en convertir estos datos en un data.frame manejable desde R. El problema procede de Statgraphics Plus 4, de J.Llovet, D. Delgado, J. Martínez. Pág. 237. Ed. Anaya Multimedia, ISBN 9788441509504. Es muy posible que haya que instalar en R algunas librerías adicionales para poder ejecutar el código completo. Hay instrucciones de cómo hacer eso en esta entrada.

Gracias por la atención.