I Conferencia Hispana R-project

30 de Diciembre de 2009

Los días 26 y 27 del mes de noviembre se celebró en la Universidad de Murcia la I Conferencia Hispana R-project. Lamentablemente no pude acercarme a Murcia para este evento pero los organizadores han puesto a nuestra disposición los vídeos de las jornadas. ¡Disfrutadlos!

Más información de cara a la organización de próximas jornadas aquí.

El modelo de ajuste con restricciones “pasa por la media”

23 de Diciembre de 2009

Consideremos el modelo teórico lineal

Y=X \beta + \epsilon

El modelo de ajuste correspondiente a este modelo teórico es

Y=X B + e

Este modelo de ajuste cumple

\bar{Y}=\bar{X} B

Este hecho se suele expresar diciendo que “el modelo de ajuste pasa por la media”. Este hecho resulta de gran utilidad práctica ya que permite el cálculo del término independiente en el caso de que el ajuste del modelo se esté efectuando mediante el método de datos centrados mediante el empleo de la siguiente fórmula:

a= \bar{y} - b_{1} \bar{X_1} - \ldots - b_{k} \bar{X_k}

Pero, ¿qué ocurre cuando al modelo teórico se le añaden restricciones lineales? ¿sigue pasando por la media? Dicho de otra forma: ¿se cumple la siguiente expresión?

\bar{Y}=\bar{X} B_C

Bajo ciertas condiciones, la respuesta es afirmativa y la demostración sumamente sencilla

Sabemos que

B_C = B + (X^t X)^{-1} C^t [C (X^t X)^{-1} C^t]^{-1} (\gamma - C B)

y por tanto

\bar{X} B_C = \bar{X} B +\bar{X} (X^t X)^{-1} C^t [C (X^t X)^{-1} C^t]^{-1} (\gamma - C B)

Para que \bar{Y}=\bar{X} B_C debe ocurrir que

\bar{X} (X^t X)^{-1} C^t [C (X^t X)^{-1} C^t]^{-1} (\gamma - C B) = 0

ya que \bar{Y}=\bar{X} B

Consideremos la matriz

\bar{X} (X^t X)^{-1} C^t

Demostraremos que, bajo ciertas condiciones muy generales, dicha matriz es nula y, en consecuencia, que

\bar{X} (X^t X)^{-1} C^t [C (X^t X)^{-1} C^t]^{-1} (\gamma - C B) = 0

por lo que

\bar{Y}=\bar{X} B_C

En efecto,

\bar{X} (X^t X)^{-1} C^t = H_1 X (X^t X)^{-1} C^t

donde H_1 es la matriz que proyecta los vectores de R^n en la dirección del vector formado por unos.

H_1 = 1 (1^t 1)^{-1} 1^t = \frac {1}{n} 1 1^t

El resultado de premultiplicar cualquier vector de R^n por la matriz H_1 es un vector de R^n cuyos componentes son todos iguales e iguales a la media de los componentes del vector original.

Consideremos ahora la matriz

X (X^t X)^{-1}

Si premultiplicamos esta matriz por X^t, el resultado es la matriz identidad, por lo que todas las columnas de X (X^t X)^{-1} excepto la primera son ortogonales al vector de R^n formado por unos -la primera columna de la matriz X-. Si no fuera así, (X^t X) (X^t X)^{-1} no sería la matriz identidad, lo cual es absurdo. Dicho de otro modo, las sumas de los elementos de todas y cada una de las columnas de la matriz  X (X^t X)^{-1} -excepto la primera- son nulas y, en consecuencia, las medias son nulas.

Como consecuencia de este hecho,

\bar{X} (X^t X)^{-1} = H_1 X (X^t X)^{-1}

es una matriz cuyos componentes son nulos salvo los de la primera columna.

¿Qué debe ocurrir para que la matriz \bar{X} (X^t X)^{-1} C^t sea nula?

Es suficiente con que la primera columna de la matriz C sea nula. En ese caso, al postmultiplicar \bar{X} (X^t X)^{-1} por C^t se obtendrán combinaciones lineales de columnas de ceros, ya que la primera columna de \bar{X} (X^t X)^{-1} -la única que puede no tener ceros- no forma parte de las combinaciones lineales.

¿Cómo podemos expresar la idea de que la primera columna de la matriz C sea nula de una forma más intuitiva? La primera columna de la matriz C recoge los coeficientes del término independiente en las restricciones lineales que se añaden al modelo teórico. Decir que la primera columna de C es nula es lo mismo que decir que las restricciones añadidas al modelo teórico no involucran al término independiente.

En resumen, hemos demostrado que si las restricciones lineales que se añaden a un modelo teórico no incluyen al término independiente -lo que, en la práctica es la situación más habitual- entonces el modelo de ajuste con restricciones “pasa por la media”, es decir, satisface la ecuación

\bar{Y}=\bar{X} B_C

Finalmente, este hecho resulta de gran utilidad práctica para calcular el término independiente en el modelo de ajuste con restricciones si es que se está trabajando con datos centrados ya que se cumple que

a_C= \bar{y} - b_{1C} \bar{X_1} - \ldots - b_{kC} \bar{X_k}

Hierarchical Clustering on Principal Components

30 de Julio de 2009

Al final no me pude acercar a Rennes para asistir a la conferencia uSER2009!. Pero los organizadores han hecho un gran trabajo y han publicado con gran diligencia los resúmenes de las ponencias que allí se presentaron.

Esto me reafirma en mi opinión de que el equipo de agrocampus ouest sigue siendo el mejor en lo que a análisis multivariante se refiere. Ahora “amenazan” con mejorar su paquete FactoMineR incluyendo en él funciones que permiten realizar análisis cluster jerárquicos a partir de los resultados obtenidos de un análisis de componentes principales.

Estaré atento a la publicación de la nueva versión.

Distribución log-normal

23 de Julio de 2009

Llevo un par de días leyendo el interesante libro de Rein Taagepera de título Making Social Sciences More Scientific. La aportación fundamental del autor consiste en una comparación entre las formas de trabajo predominates en los ámbitos de las ciencias naturales y las ciencias sociales. Afirma Taagepera que las ciencias sociales son “menos ciencias” que las naturales y que, quizás, la aplicación de los métodos empleados en estas últimas podría mejorar los resultados y la imagen social de aquéllas.

Entre otras cosas, el autor critica la omnipresencia en las ciencias sociales de la distribución normal -y de la media aritmética-. En su opinión, en muchas ocasiones es una mejor opción la distribución log-normal y la media geométrica. El autor nos proporciona unos consejos prácticos para seleccionar la media que mejor se adapta a los datos de que disponemos. Cito:

Geometric means often express the central tendency better than arithmetic means. For the same reason, lognormal data fits often are called for, instead of desperate attempts to fit data into a Procrustean normal distribution. The following advice applies, with some reservations.

  • In the absence of any other information, if a variable can range from minus to plus infinity, a normal distribution is our best  bet, implying that the arithmetic mean is close to the median. (In the presence of further information, the bet may be off.)
  • In the absence of any other information, if a variable can have only positive values, a lognormal distribution is among our best bets, implying that the geometric mean is close to the median. (In the presence of further information, the bet may be off-we may have a gamma distribution or something else.)
  • However, if one tries a normal fit and standard deviation turns out less than one-half of the mean, then one might use this normal distribution. If standard deviation exceeds one-half of the mean, the normal fit should be abandoned in favor of lognormal.
  • If negative values are conceptually excluded but zero values do ocurr, then neither distribution can fit. Neither mean adequately reflects the median, but a pseudo-geometric mean might approximate it.
  • When there are grounds to hesitate between the arithmetic and geometric means, using the median might be the safest way, although it is awkward to calculte.

A este respecto resulta muy interesante el artículo de Limpert et al. (2001) de título Lognormal Distributions across the Sciences: Keys and Clues.

Classification and regression trees

23 de Junio de 2009

Una alternativa muy razonable a los métodos paramétricos de regresión y análisis discriminante es el llamado método CART, siglas de classification and regression trees. Este método, popularizado por Breiman, Friedman, Olshen y Stone en 1984, es de aplicación cuando se trata de predecir el valor de una variable dependiente cuantitativa -caso de la regresión- o cualitativa -caso de la clasificación- a partir del conocimiento de los valores que toman una serie de variables explicativas.

Trataremos en este post el caso de la clasificación, es decir, la situación en la que la variable dependiente es cualitativa, con dos o más modalidades.

En esencia, el método CART parte de un nodo raíz o nodo inicial, que contiene al conjunto de todos los individuos, y establece una partición de dicho conjunto en dos subconjuntos -llamados nodos hijos- siguiendo el criterio de minimizar la heterogeneidad interna de dichos nodos -lo que supone maximizar la heterogeneidad entre ellos-. Este paso se repite con cada uno de los nodos hijos y así sucesivamente hasta satisfacer el criterio de parada -por ejemplo, porque se ha obtenido un nodo completamente homogéneo o de un tamaño menor que un umbral previamente establecido-. Así, el resultado final de esta primera fase del método CART es un árbol binario -cada nodo padre se divide en dos nodos hijos- de nombre árbol máximo.

El conjunto de nodos terminales del árbol máximo constituye una partición del conjunto inicial de individuos. A cada uno de dichos nodos terminales se le asigna una etiqueta igual a la moda en dicho nodo de la variable cualitativa dependiente. Supongamos, por ejemplo, que la variable dependiente analizada es el sexo, con dos modalidades, H y M. Se asignaría la etiqueta H a aquellos nodos terminales en los que haya más H que M y la etiqueta M a los demás. La asignación de la etiqueta a los nodos terminales supone la existencia de un error de clasificación debido a la presencia en nodos terminales de individuos cuya modalidad en la variable dependiente no coincide con la etiqueta asignada a dicho nodo terminal. Naturalmente, el árbol máximo es el árbol con el mínimo error de clasificación.

A partir del árbol máximo comienza el llamado proceso de poda del árbol. Hemos comentado que el árbol máximo tiene un error de clasificación muy pequeño pero sólo a costa de una gran complejidad -un elevado número de nodos-. En la mayoría de las situaciones merece la pena incrementar en pequeña medida el error de clasificación para simplificar la estructura del árbol binario obtenido. Se obtiene de este modo un árbol llamado árbol óptimo -al menos desde el punto de vista del coste-complejidad-.

Es muy importante recordar que el error de clasificación calculado para el árbol óptimo tiende a sobrevalorar el verdadero poder explicativo de las variables independientes para explicar la variable dependiente. La práctica habitual es la de destinar una parte de la muestra total de individuos a “aprender el árbol” y otra parte de la muestra a contrastar la calidad del árbol.

Naturalmente, R dispone de paquetes que permiten la aplicación del método CART. Por ejemplo, el package rpart.