Acerca del Tema 1 (normal multivariante)


Si considero un vector aleatorio \(\bf X\) con distribución uniforme, ¿la distancia de Mahalanobis de un punto a la esperanza de la distribución será cero porque la densidad uniforme es constante?.

No. En realidad, la distancia de Mahalanobis sólo utiliza la información de varianzas (dispersión) y correlaciones (relación lineal) del vector aleatorio, no toda la información contenida en la densidad. Así que en realidad \(d_M\) sólo toma en cuenta una geometría muy básica de los conjuntos de nivel de la densidad. El hecho de que \(d_M({\bf x},{\bf 0})=d_E({\bf y},{\bf 0})\), donde \(d_E\) es la distancia euclídea e \(\bf y\) es la estandarización multivariante de \(\bf x\), quiere decir que primero centramos \(\bf x\), luego lo medimos en unidades de desviación típica y, luego le quitamos su estructura de correlaciones y finalmente medimos la distancia euclídea del vector resultante, \({\bf y}\), a 0.

Cuando la densidad de \(\bf X\) es normal, entonces eso sí quiere decir que \(d_M({\bf x},{\boldsymbol\mu})\) es más pequeña para vectores \({\bf x}\) en los que la densidad \(f({\bf x})\) es más parecida a \(f({\boldsymbol\mu})\). Pero esto es sólo porque la densidad normal tiene conjuntos de nivel que son elipses cuya forma depende directamente de las varianzas y correlaciones de \({\bf X}\).

Si \({\bf X}=(X_{1},X_{2})'\) fuera un vector con densidad uniforme en el cuadrado \([-1/2,1/2]^2\), entonces las componentes \(X_{1}\) y \(X_{2}\) de \({\bf X}\) serían independientes entre sí y con distribución uniforme en el intervalo \([-1/2,1/2]\). En este caso es fácil comprobar que \({\boldsymbol\mu}= {\bf 0}\) y que la matriz de varianzas-covarianzas es \[ {\boldsymbol\Sigma} = 12 \left( \begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array} \right). \] Por tanto, dado un punto \({\bf x}=(x_{1},x_{2})'\), \(d_M({\bf x},{\boldsymbol\mu})=12\|{\bf x}\|^2\). Es decir, la distancia de Mahalanobis en este caso es proporcional a la distancia euclídea.


¿Cómo es posible que el vector \(\mathbf X=(X_1,X_2)'\) tenga marginales normales y no sea normal bivariante?

Pongamos un ejemplo de distribución bivariante no gaussiana con marginales normales. Sea \(X_1\) una v.a. N(0,1) y sea \[ X_2 = \left\{ \begin{array}{rl} -X_1 & \mbox{si $-1\leq X_1\leq 1$;} \\ X_1 & \mbox{en otro caso.} \end{array} \right. \]

  1. Calculando la función de distribución de \(X_2\), \(F_{X_{2}}(x_{2})=\mathbb P\{X_{2}\leq x_{2}\}\), es fácil ver que \(X_2\) sigue una distribución \(N(0,1)\): \[ {\mathbb P}\{X_{2}\leq x_{2}\} = \left\{ \begin{array}{ll} {\mathbb P} \{ X_{1} \leq x_{2} \} = \Phi(x_{2}) & \mbox{si $x_2< -1$} \\ {\mathbb P} \{ X_{1} < -1 \} + \mathbb P \{ x_2\leq X_{1}\leq 1\} = 1-\Phi(1)+\Phi(1)-\Phi(x_2) = \Phi(x_2) & \mbox{si $-1 \leq x_2\leq 1$} \\ {\mathbb P} 1-\{ X_{1} \geq x_{2} \}=1-(1-\Phi(x_2)) =\Phi(x_2) & \mbox{si $x_2> 1$} \end{array} \right. \] donde \(\Phi\) denota la función de distribución normal estándar.

  2. Calculando la probabilidad de que \(X_1-X_2\) sea igual a 0 y viendo que no es 0, deducimos que la distribución conjunta de \(X_1\) y \(X_2\) no es normal bivariante, porque, si lo fuera, entonces \(X_1-X_2\) también debería ser normal univariante y en este caso \(\mathbb P\{X_1-X_2=0\}=0\): \[ \mathbb P\{X_1-X_2=0\} = \mathbb P \{X_1< -1\} + \mathbb P \{X_1> 1\} = 2(1-\Phi(1)) \neq 0. \]

Lo que no he calculado es \(\mbox{Cov}(X_1,X_2)\), pero creo que no será 0. Me faltaría pensar un ejemplo de dos marginales univariantes \(X_1\) y \(X_2\) que sean \(N(0,1)\) y con \(\mbox{Cov}(X_1,X_2)=0\), pero tal que \(\mathbf X=(X_1,X_2)'\) no sea normal bivariante. Si a alguien se le ocurre, que me lo diga.


En la descomposición espectral de una matriz de covarianzas \(\boldsymbol\Sigma = \bf C \bf D \bf C'\), siendo \(\bf C\) la matriz ortogonal cuyas columnas son los autovectores ortonormalizados de \(\boldsymbol\Sigma\) y \(\bf D\) la matriz diagonal con los correspondientes autovalores de \(\boldsymbol\Sigma\), ¿es \(\bf D\) la matriz de Jordan?.

Agradezco a Ernesto Girondo que me resolviera la pregunta. La respuesta es sí.


Acerca del Tema 2 (Contrastes no paramétricos)


Queremos comprobar si la distribución de unas localizaciones en una ventana de observación es uniforme. Para ello dividimos la ventana superponiendo una malla que da lugar a una partición de rectángulos \(A_1,\ldots,A_k\), todos de la misma área. Definimos la v.a. cualitativa \(X\) = ‘’rectángulo en el que se observa la localización’‘, cuyo espacio muestral es \(A_1,\ldots,A_k\). También podemos definir la v.a. cuantitativa \(Y\) = “Número de localizaciones observadas en un rectángulo de la malla’’, cuyo espacio muestral es 0,1,2,… ¿Se usa un contraste de bondad de ajuste de la v.a. \(X\) a una uniforme o se contrasta si todas las \(Y\) sigue una distribución de Poisson?

Las dos opciones son correctas, pero depende de si tenemos los datos originales y, si no, de cómo están agrupados los datos. Si dividimos la ventana de observación de las localizaciones observadas mediante una retícula muy fina entonces seguramente funcione mejor el contraste de bondad de ajuste a una Poisson. Si la retícula tiene rectángulos muy amplios en cada uno de los cuales pueden caer muchas observaciones, entonces es mejor usar el contraste a una uniforme.

Quizá se entienda mejor con un ejemplo.

Genero \(n\) observaciones de una distribución uniforme en el cuadrado unidad.

n <- 100
Datos <- cbind(runif(n),runif(n))
par(pty="s")
plot(Datos,xlab="",ylab="",pch=20,xlim=c(0,1),ylim=c(0,1))

Separo los datos según una malla “grosera” de 3x3. La función quadratcountnos devuelve la frecuencia observada de cada rectángulo.

library(spatstat)
Datosppp <- ppp(x=Datos[,1],y=Datos[,2],xrange=c(0,1),yrange=c(0,1))
Q <- quadratcount(Datosppp,nx=3,ny=3)
Q
##                x
## y               [0,0.333) [0.333,0.667) [0.667,1]
##   [0.667,1]             9             8        13
##   [0.333,0.667)        11            11        15
##   [0,0.333)            10            10        13
plot(Q)

Contraste de bondad de ajuste de \(X\) a una uniforme: \[ H_{0}: \; {\mathbb P}\{X=A_j\} = \frac{1}{k}, \;\; j=1,\ldots,k. \]

vectorQ <- matrix(Q,nrow=1)
chisq.test(vectorQ)
## 
##  Chi-squared test for given probabilities
## 
## data:  vectorQ
## X-squared = 3.5, df = 8, p-value = 0.8992

Vemos que no podemos rechazar la hipótesis nula, que es la decisión acertada.


Contraste de bondad de ajuste de \(Y\) a una Poisson(\(\lambda\)): \[ H_{0}: \; Y \sim \mbox{Poisson}(\lambda) \;\; \mbox{para algún } \lambda>0. \] Estimamos el parámetro \(\lambda\) de la Poisson:

N <- length(vectorQ) # Número de rectángulos en la malla
lambdaEst = sum(vectorQ)/N # Número promedio de localizaciones por rectángulo de la malla
lambdaEst
## [1] 11.11111

Valores que ha tomado la variable \(Y\) en la muestra (es un valor de \(Y\) por cada rectángulo de la malla).

sort(vectorQ)
## [1]  8  9 10 10 11 11 13 13 15

Vemos que, aunque el espacio muestral de la Poisson es \({\mathbb N}\), la cantidad de valores observados de \(Y\) es pequeña, porque la partición es grosera. Sin embargo, nosotros tenemos que calcular las frecuencias observadas y esperadas (bajo \(H_0\)) de muchos valores del espacio muestral, tanto si los hemos observado como si no. Vamos a calcular las frecuencias de los valores 0,1,2,…,\(M\) = máximo valor observado de \(Y\).

M <- max(vectorQ)
# Frecuencias observadas y esperadas para los valores 0,1,2,...,M
obs <- rep(NA,M+1)
esp <- rep(NA,M+1)
for (i in 1:(M+1)){
  obs[i] <- sum(vectorQ == (i-1))
  esp[i] <- n * dpois(i-1,lambdaEst)
}
esp[M+1] <- n-sum(esp[1:M])
Estchi2 <- sum(((obs-esp)^2)/esp)
Estchi2
## [1] 83.3726
pvalor <- pchisq(Estchi2,df=(M+1)-1-1,lower.tail=F)
pvalor
## [1] 6.669131e-12

Vemos que el p-valor es muy pequeño, por lo que rechazamos la hipótesis nula, que es una decisión incorrecta.


Ahora repito el proceso anterior con una malla “fina” de 10x10. El contraste de bondad de ajuste de \(X\) a la uniforme ni siquiera es válido, porque las frecuencias observadas son muchas menores que 5. Por el contrario, el contraste de bondad de ajuste de \(Y\) a la Poisson funciona estupendamente (en general): el p-valor no es significativo.

Q <- quadratcount(Datosppp,nx=10,ny=10)
Q
##            x
## y           [0,0.1) [0.1,0.2) [0.2,0.3) [0.3,0.4) [0.4,0.5) [0.5,0.6) [0.6,0.7)
##   [0.9,1]         1         0         1         3         1         1         2
##   [0.8,0.9)       0         1         0         0         1         0         0
##   [0.7,0.8)       2         1         1         1         0         0         1
##   [0.6,0.7)       1         0         1         0         1         0         4
##   [0.5,0.6)       1         2         2         0         3         0         1
##   [0.4,0.5)       1         0         1         1         1         1         1
##   [0.3,0.4)       0         0         1         1         1         0         1
##   [0.2,0.3)       2         0         2         2         1         1         0
##   [0.1,0.2)       0         0         0         1         0         0         1
##   [0,0.1)         1         2         3         1         1         0         2
##            x
## y           [0.7,0.8) [0.8,0.9) [0.9,1]
##   [0.9,1]           3         1       0
##   [0.8,0.9)         1         4       0
##   [0.7,0.8)         1         0       0
##   [0.6,0.7)         0         0       2
##   [0.5,0.6)         1         3       0
##   [0.4,0.5)         2         1       3
##   [0.3,0.4)         0         1       3
##   [0.2,0.3)         4         1       0
##   [0.1,0.2)         0         3       1
##   [0,0.1)           0         0       2
vectorQ <- matrix(Q,nrow=1)
chisq.test(vectorQ)
## Warning in chisq.test(vectorQ): Chi-squared approximation may be incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  vectorQ
## X-squared = 108, df = 99, p-value = 0.252
N <- length(vectorQ) # Número de rectángulos en la malla
lambdaEst = sum(vectorQ)/N # Número promedio de localizaciones por rectángulo de la malla
lambdaEst
## [1] 1
sort(vectorQ)
##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4
M <- max(vectorQ)
# Frecuencias observadas y esperadas para los valores 0,1,2,...,M
obs <- rep(NA,M+1)
esp <- rep(NA,M+1)
for (i in 1:(M+1)){
  obs[i] <- sum(vectorQ == (i-1))
  esp[i] <- n * dpois(i-1,lambdaEst)
}
esp[M+1] <- n-sum(esp[1:M])
obs
## [1] 37 40 12  8  3
esp
## [1] 36.787944 36.787944 18.393972  6.131324  1.898816
Estchi2 <- sum(((obs-esp)^2)/esp)
Estchi2
## [1] 3.712438
pvalor <- pchisq(Estchi2,df=(M+1)-1-1,lower.tail=F)
pvalor
## [1] 0.2942367

Acerca del Tema 4 (Clasificación)


La cuestión es que no entiendo qué muestra R en Deviance Residuals. Por ejemplo, en el ejemplo de la esclerosis múltiple nos muestra esto:

Deviance Residuals:
Min 1Q Median 3Q Max
-1.31250 -0.45006 -0.21017 0.01804 2.63976

En principio creía que mostraba información sobre la desviación en el individuo \(i\) (definida como \(D^2_i = -2[y_i \, \log \eta(x_i) + (1 - y_i) \, \log(1 - \eta(x_i))]\): su valor mínimo, su mediana, etc. Sin embargo, creo que estos valores son siempre positivos, por lo que no tiene sentido que su mediana sea un valor negativo.

Cada deviance residual del output de Res la raíz cuadrada de la desviación \(D^2_i\), con un signo asociado: +1 si \(y_i=1\), -1 si \(y_i=0\). Los valores de estos residuos los puedes extraer con la función resid:

reglog = glm(V6~V1+V2+V3+V4+V5,data=Datos,family="binomial")

summary(resid(reglog))

Min. 1st Qu. Median Mean 3rd Qu. Max. -1.31250 -0.45006 -0.21017 -0.08950 0.01804 2.63976


El comando ClasFisher$scaling ¿da el vector \({\mathbf w}_F\)? Por otro lado, ¿qué significa el output que da summary(ClasFisher)? Y por último, no entiendo muy bien el procedimiento leave-one-out, ni para qué sirve ni la forma en que lo implementas en R.

ClasFisher$scaling no da el vector \({\mathbf w}_F\) exactamente, sino un vector proporcional. Pero no sé decir por qué. Ni siquiera tiene necesariamente longitud 1. De todas formas, para el procedimiento de clasificación no importa la longitud por la expresión de la regla de decisión (que es una desigualdad en la que en ambos lados se multiplica escalarmente por ese vector de \({\mathbf w}_{F}\)).

¡Yo tampoco sé qué proporciona el output de summary(ClasFisher)! No lo hemos visto en clase y lo he mirado ahora pero no sé qué da. En cualquier caso, no es nada interesante. No merece la pena que perdamos el tiempo en averiguar en qué consiste.

El procedimiento leave-one-out (LOO) es un caso particular de una importante metodología práctica (computacional) que se utiliza en Estadística para comprobar si un procedimiento (típicamente de clasificación o regresión) funciona bien o no. LOO es un caso particular de validación cruzada (o CV, por sus siglas en inglés, cross-validation). CV consiste en un procedimiento iterativo que, en cada iteración del bucle, separa la muestra en dos submuestras (una grande y una pequeña). Con la grande construye un clasificador (si nuestro problema estadístico es la clasificación) y lo utiliza para clasificar la submuestra pequeña. Entonces, para la submuestra pequeña se comparan las clases que el clasificador ha asignado a los individuos con las clases a las que pertenecen realmente. Se guarda la información del número de clasificaciones erróneas y el tamaño de la muestra pequeña en cada iteración. Al terminar el procedimiento iterativo (un bucle for en R) se utiliza la información guardada para proporcionar una estimación de la proporción o tasa de clasificaciones erróneas de ese clasificador (en nuestro caso, el de Fisher o la regresión logística).

Seguramente se entiende mejor en el caso del algoritmo LOO. En este caso, en la iteración \(i\) del bucle for se toma la observación \(i\)-ésima de la muestra. Esta observación es la submuestra pequeña y las restantes \(n-1\) observaciones son la submuestra grande. Con la submuestra grande construyo el clasificador y con éste clasifico la observación \(i\). Acumulo la clase en la \(i\)-ésima componente de un vector de tamaño \(n\) que he inicializado antes del bucle. Al finalizar el bucle tengo un vector de las poblaciones que el clasificador ha asignado a cada individuo de la muestra. Comparo la clase decidida por el clasificador con la clase a la que realmente pertenece cada individuo. Con esto determino la proporción de fallos (tasa de error) del clasificador. Si en lugar del usar LOO uso otro procedimiento de CV, esta tasa de error puede cambiar (esperemos que poco).

Respecto a cómo llevar a cabo LOO en la práctica para el clasificador de Fisher con R, conviene aprovechar el output de la función lda del paquete MASS. El resultado final es el mismo que en el algoritmo que expliqué arriba.

  1. Primero se hace el bucle en el que clasifico cada individuo usando como muestra de entrenamiento las \(n-1\) restantes observaciones. Esto se consigue con la opción CV=TRUE en la función lda.

ClasFisher = lda(X,Y,CV=TRUE)

  1. Calculo el tamaño muestral

n = nrow(X) # Tamaño muestral

  1. Comparo la clase a la que realmente pertenece cada individuo de la muestra (que están guardadas en el vector Y) con la clase asignada por el clasificador, que se obtiene con ClasFisher$class. Obtengo la proporción de clasificaciones erróneas.

sum(Y != ClasFisher$class)/n


Acerca de R


¿Cómo divido la pantalla de gráficos de R en una cuadrícula?

Hay dos maneras. Supongamos que queremos dividirla en una cuadrícula de 2 filas y tres columnas. Con layout le tenemos que dar una matriz con componentes de 1 a 6 en 2 filas y 3 columnas:

X = matrix(rnorm(20),ncol=2) # Datos
layout(matrix(c(1:6),2,3)) # Divido la pantalla de gráficos
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
layout(matrix(1)) # Deshago la división de la pantalla

Segunda manera: con mfrow directamente le indico el número de filas y columnas.

X = matrix(rnorm(20),ncol=2) # Datos
par(mfrow=c(2,3)) # Divido la pantalla de gráficos
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
plot(X,pch=21)
par(mfrow=c(1,1)) # Deshago la división de la pantalla


Cuando escribo set.seed luego siempre se genera el mismo número aleatorio, aunque no vuelva a ejecutar set.seed.

No. Ha debido de ser que seguías ejecutando el set.seed aunque no te dieras cuenta. Por ejemplo,

set.seed(20)
rnorm(1)
## [1] 1.162685
rnorm(1)
## [1] -0.5859245
rnorm(1)
## [1] 1.785465

Lo que sí es cierto es que la semilla fijada determina (secuencialmente) todos los números aleatorios que se generen a partir de ese momento. Es decir, si vuelvo a escribir las cuatro líneas de código de arriba volveré a obtener exactamente los mismo tres números aleatorios.

Normalmente al iniciar una nueva sesión de R, el programa internamente toma una semilla nueva basada, entre otras cosas, en la hora. Para deshacer la orden de fijar la semilla se pueden usar cualquiera de los dos comandos siguientes

rm(.Random.seed)
set.seed(NULL)


¿Cómo convierto una matriz de números en una tabla de datos?

Con la función data.frame. Por ejemplo,

v = 1:6
M = matrix(v,nrow=3) # Matriz numérica 3x2
Datos = data.frame(M) # Convertimos la matriz M en una tabla de datos (objeto data.frame)
names(Datos) = c("a","b") # Le damos nombre a las columnas de la tabla