Manual de Estadística Multivariante

MANUAL ABREVIADO DE
ESTADÍSTICA MULTIVARIANTE
Jesús Montanero Fernández

Introducción
El presente manual autoeditado y revisado periódicamente pretende constituir una introduc-
ción a las técnicas clásicas de la Estadística Multivariante, con breves incursiones en métodos
más novedosos. Está concebido en principio como apoyo a la docencia en una asignatura de
unas sesenta horas lectivas, correspondiente al cuarto curso del Grado en Estadística de la Uni-
versidad de Extremadura, e impartida por tanto a alumnos con cierta formación matemática
en general y estadística en particular. Consta de un capítulo inicial enfocado mayormente a
la comprensión de los modelos lineales univariante y multivariante, que se tratan en los dos
capítulos siguientes. El estudio del modelo lineal es interesante de por sí para cualquier esta-
dístico, pero en nuestro contexto debe entenderse mayormente como una herramienta teórica
que nos permite comprender mejor el problema de clasificación, que se afronta en el capítulo
4. Se trata del tema más interesante desde el punto de vista práctico, junto con el capítulo 5,
que aborda un problema genuinamente multivariante como es el de simplificar observaciones
multidimensionales con el fin último de entender datos complejos. El último capítulo está de-
dicado el análisis de conglomerados, técnica que puede estudiarse también desde el punto de
vista de la Minería de Datos.
La intención a la hora de de elaborar este manual ha sido exponer los contenidos de ma-
nera breve, pero indicando al lector referencias bibliográficas oportunas para profundizar en
el estudio de la Estadística Multivariante. En particular, no se incluyen las demostraciones
de los resultados teóricos. Algunas son asequibles y se proponen como ejercicio; otras pueden
encontrarse en la bibliografía recomendada, por ejemplo en los manuales 56 y 59 del Servicio
de Publicaciones de la UEx, que podrían considerarse como versiones extendidas del presente
volumen. Tampoco pretende ser exhaustivo. De hecho, ciertas técnicas que podemos catalo-
gar de multivariantes, como el análisis de componentes principales no lineal, el escalamiento
multidimensional, el análisis de correspondencias y las redes neuronales se introducen a título
meramente informativo. El lector interesado puede encontrar información más amplia sobre
todas ellas en Hastie et al. (2008), Gifi (1990), Hair et al. (1999) y Uriel y Aldás (2005). Debe
mencionarse también que, de todas las referencias indicadas en la bibliografía, los que han
tenido una influencia más patente en la redacción de este volumen ha sido Arnold (1981) y
Hatie et al. (2008).
Por último, hacemos constar que los diferentes gráficos y tablas de resultados que aparecen
a lo largo del volumen han sido obtenidos mediante los programas SPSS y R. De hecho, se
incluyen algunas capturas de pantallas para ilustrar la ejecución de las técnicas más relevantes
con SPSS.
Badajoz, diciembre de 2018

Índice general
1. Preliminares 7
1.1. Notación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2. Principales parámetros probabilísticos . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3. Principales parámetros muestrales . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4. Regresión lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5. Nociones básicas de Álgebra Lineal . . . . . . . . . . . . . . . . . . . . . . . . . 15
2. Modelo lineal multivariante 19

2.1. Distribución normal multivariante . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2. Modelo lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.1. Distribuciones asociadas al modelo . . . . . . . . . . . . . . . . . . . . . 22
2.2.2. El modelo y ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2.3. Estimación y contraste de hipótesis . . . . . . . . . . . . . . . . . . . . . 26
2.3. Modelo multivariante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1. Distribuciones asociadas al modelo . . . . . . . . . . . . . . . . . . . . . 31
2.3.2. El modelo y ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3.3. Estimación y contraste de hipótesis . . . . . . . . . . . . . . . . . . . . . 34
3. Aplicaciones del modelo 37

3.1. Inferencia para una media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2. Inferencia para dos medias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3. Manova de una vía . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.1. Ejes discriminantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.4. Regresión multivariante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4.1. Contraste total: análisis de correlación canónica . . . . . . . . . . . . . . 45
3.4.2. Contrastes parciales: método Lambda de Wilks . . . . . . . . . . . . . . 47
3.5. Análisis de perfiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4. Problema de clasificación 53
4.1. Planteamiento general . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2. Análisis Discriminate Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2.1. LDA y ejes discriminantes . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2.2. Estrategia cuadrática de Fisher . . . . . . . . . . . . . . . . . . . . . . . 60
4.3. Métodos alternativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.1. Regresión logística . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.2. Vecino más próximo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3.3. Árbol de decisión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5
Introducción
4.3.4. Validación de estrategias: . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.3.5. Indicaciones sobre Redes Neuronales . . . . . . . . . . . . . . . . . . . . 69
5. Reducción dimensional 73
5.1. Componentes principales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.2. Justificación geométrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.3. Representación de observaciones y variables . . . . . . . . . . . . . . . . . . . . 77
5.3.1. Representación de observaciones . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.2. Representación de variables . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.3. Representación conjunta de observaciones y variables . . . . . . . . . . . 80
5.3.4. Rotación de ejes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.4. Análisis factorial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4.1. Modelos con factores latentes . . . . . . . . . . . . . . . . . . . . . . . . 84
5.5. Indicaciones sobre Análisis de Correspondencias . . . . . . . . . . . . . . . . . . 85
5.6. Multicolinealidad y PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6. Análisis de conglomerados 91
6.1. Método de k-medias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.2. Método jerárquico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.3. Algoritmo EM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.4. Análisis de conglomerados bietápico . . . . . . . . . . . . . . . . . . . . . . . . . 97
Índice alfabético 101
6
Capítulo 1
Preliminares
En este capítulo intentaremos fijar la notación, así como definir e interpretar conceptos
fundamentales en el contexto de la Estadística Multivariante, muchos de los cuales deben ser
conocidos. También llevaremos a cabo un breve repaso de Álgebra Lineal.
1.1. Notación
En general, solemos manejar en estadística dos tipos de lenguajes: probabilístico y muestral.
El primero sirve para expresar las propiedades de la población objeto del estudio, entendiendo
población en un sentido amplio; el segundo se utiliza para expresar las propiedades de una
muestra de n datos extraídos, se supone que aleatoriamente, de dicha población.
El marco formal en el que se desarrolla el estudio poblacional es el espacio L2 de funciones
reales de cuadrado integrable, definidas sobre cierto espacio de probabilidad. Queremos decir
que las variables aleatorias que estudiemos se identificarán con elementos de L2 . El estudio
muestral tiene lugar en el espacio Euclídeo Rn , es decir que, dada una variable aleatoria X ∈ L2 ,
una muestra aleatoria de tamaño n de dicha variable se identificará con un vector X de Rn ,
cuyas componentes Xi serán las distintas mediciones de la misma. Obsérvese que hemos utilizado
distintas fuentes de letra para denotar ambos conceptos, norma que intentaremos seguir en la
medida de lo posible.
En el contexto del análisis multivariante, X puede denotar con frecuencia un vector aleatorio
p-dimensional de componentes X[1], . . . , X[p]. En tal caso, una muestra aleatoria de tamaño n
para dicho vector aleatorio se expresará mediante la matriz X ∈ Mn×p que descompone así:
X01
   
X1 [1] . . . X1 [p]
X = (X[1], . . . , X[p]) =  .. ..  =  ..  (1.1)
 . .   . 
Xn [1] . . . Xn [p] X0n
A título de ejemplo, en el cuadro 1.1 se expone una muestra de tamaño n = 38 de un vector
aleatorio de dimensión p = 8. Los datos corresponden a medidas de la motilidad de esperma-
tozoides en moruecos y fueron recogidos por J.A. Bravo en el CENSYRA de Badajoz.
L2 forma parte de una categoría de espacios que generalizan el concepto de espacio Euclídeo
por estar también dotados de un producto interior. Concretamente, dados f, g ∈ L2 , se define
hf, gi = EP [f · g] (1.2)
7
Capítulo 1 Notación
EP se entiende como el funcional que asigna a cada variable aleatoria su integral respecto a la
probabilidad P definida en el espacio de origen. El subíndice P suele omitirse. En el contexto
estadístico en podemos considerar en Rn el siguiente producto interior proporcional al conocido
como producto escalar:
n
1X
ha, bi = ai bi (1.3)
n i=1
En ambos espacios, los respectivos productos inducen sendas normas (al cuadrado), definidas
en general mediante kak2 = ha, ai y, en consecuencia, sendas métricas basadas en la norma al
cuadrado de las diferencias:
d2 (X, Y ) = E[(X − Y )2 ], X, Y ∈ L2 (1.4)
1X
d2n (X, Y) = (Xi − Yi )2 , X, Y ∈ R n (1.5)
n i
La segunda es, salvo una homotecia, la distancia Euclídea al cuadrado en Rn . El uso de estas
distancias para cuantificar errores se asocia al denominado método de Mínimos Cuadrados. Por
otra parte, del producto interior se deriva a su vez una noción de ortogonalidad o perpendicu-
laridad. En Rn decimos que a y b son ortogonales entre sí cuando ha, bi = 0, en cuyo caso se
denota a ⊥ b. En L2 se define de manera análoga.
Proyección ortogonal: La noción de perpendicularidad se relacionará bajo ciertas condi-

ciones con los conceptos estadísticos de incorrelación y de independencia. Además, da pie a
considerar un tipo de función lineal denominada proyección ortogonal. Concretamente, si V es
un subespacio lineal cerrado del espacio E (E en nuestro caso puede tratarse de L2 o de Rn ),
se define PV como la aplicación que asigna a cada elemento e del espacio el único elemento de
V tal que e − PV e ⊥ V , en cuyo caso la distancia entre e y PV e es la mínima posible entre
e y un elemento de V . Si V1 y V2 son dos subespacios ortogonales de E, se verifica que que
PV1 ⊕V2 = PV1 + PV2 . Además, kPV1 ⊕V2 ek2 = kPV1 ek2 + kPV2 ek2 . Para V ⊂ Rn , dado que PV es
una aplicación lineal se identificará con una matriz n × n que se denotará de la misma forma.
Ejercicio 1. Dado V ⊂ Rm , probar que tr(PV ) = dim V , y que todos lo elementos de la
diagonal de PV pertenecen al intervalo [0, 1].
6

E

e − PV e

e

1

t PV e
0

V
Figura 1.1: Proyección ortogonal
8
Capítulo 1 Principales parámetros probabilísticos
La colección de resultados teóricos conocida como Leyes de los Grandes Números y Teore-
ma Central del Límite establecen una clara conexión entre los espacios Rn y L2 , si entendemos
X ∈ Rn como una muestra aleatoria simple de una variable aleatoria X ∈ L2 . La idea germi-
nal de estos resultados y, por lo tanto, la conexión entre Estadística y Probabilidad, podría
ilustrarse mediante el triángulo de Pascal y enunciarse en términos intuitivos así: si asumimos
una completa simetría en la elección de n individuos (equiprobabilidad e independencia), la
media de sus respectivos datos se estabiliza conforme n aumenta, debido a la conmutatividad
(invarianza ante permutaciones) de la operación suma.
En todo caso, lo más importante en esta sección es resaltar que todos las definiciones
en L2 expresadas en términos del producto interior pueden traducirse automáticamente al
lenguaje muestral e interpretarse de manera completamente análoga. Por ello, en este capítulo
nos centraremos principalmente en el estudio de los parámetros probabilísticos o poblacionales
(salvo en la sección 1.3, donde se fijará la notación de los parámetros muestrales), dejando como
ejercicio para el lector el estudio paralelo en términos muestrales. Por lo general seguiremos
la costumbre habitual de expresar los parámetros probabilísticos mediante letras griegas y sus
homólogos muestrales con notación latina.
Si X es una familia de k elementos, bien sean de L2 o de Rn (en el segundo caso puede
identificarse con una matriz n × k), se denota por hX i su expansión lineal. En el espacio L2
se denotará por 1 la variable aleatoria con valor constante 1, siendo entonces h1i el subespacio
unidimensional de las funciones constantes en L2 ; se denotará por h1i⊥ su ortogonal, que es
un hiperplano de L2 . Análogamente, se denotará por 1n al vector de Rn cuyas componentes
son todas 1, siendo por tanto h1n i la recta de los vectores constantes y h1n i⊥ su ortogonal, de
dimensión (n − 1).
1.2. Principales parámetros probabilísticos

En esta sección definiremos los parámetros relacionados con los momentos de orden uno y
dos. Con ello estamos centrando indirectamente nuestro estudio en el ámbito de la distribución
normal y de las relaciones de tipo lineal.
Media: Primeramente definimos la media de una variable aleatoria X como su esperanza, es

decir, su integral respecto a la probabilidad considerada. Se denota por E[X] o por la letra µ,
acompañada si es necesario de un subíndice aclaratorio. Si X es un vector p-dimensional, su
media es el vector p-dimensional compuesto por las respectivas medias, y se denotará de forma
idéntica.
Varianza: Dada una variable aleatoria X ∈ L2 de media µ, se define su varianza mediante
var[X] = E[(X − µ)2 ] (1.6)
denotándose también por la letra σ 2 . La raíz cuadrada positiva de la varianza se denomina

desviación típica. Nótese que la varianza está bien definida al ser X de cuadrado integrable.
De hecho, puede expresarse mediate
var[X] = E[X 2 ] − E[X]2 (1.7)
Ejercicio 2. Probar (1.7).
9
Ejercicio 3. Probar que µ = Ph1i X.

Del ejercicio 3 se deduce que la media µ de una variable aleatoria X es la función constante
más próxima en términos de la distancia (1.4). Así pues µ se interpreta como la constante más
próxima a X. La diferencia X − µ ∈ h1i⊥ se denomina variabilidad de X. La distancia respecto
a esa constante más próxima es precisamente la varianza:
var[X] = d2 (X, E[X]) (1.8)
Varianza total: Si X es un vector aleatorio p-dimensional de componentes X[1], . . . , X[p], se

define la varianza total de X mediante
p
X
varT [X] = var X[j] (1.9)
j=1
Este parámetro puede interpretarse en términos de la distancia d2[p] definida en el espacio de

los p-vectores aleatorios con componentes en L2 mediante
d2[p] (X, Y ) = EP kX − Y k2Rp

(1.10)
Ejercicio 4. Probar que E[X] es el vector aleatorio constante que más se aproxima a X en
términos de la distancia (1.10) y que, además,
varT [X] = d2[p] (X, E[X]) (1.11)
Covarianza: Dado un vector aleatorio p-dimensional X, se define la covarianza entre dos

componentes cualesquiera X[i] y X[j] del mismo como el producto interior de sus respectivas
variabilidades, es decir,

cov X[i], X[j] = hX[i] − µi , X[j] − µj i (1.12)
denotándose también por σij . Se trata de una generalización de la varianza, pues σii = σi2 ,
que describe, según veremos en la próxima sección, el grado de relación lineal existente entre
las variabilidades, es decir, el grado de relación afín existente entre las variables originales.
Se dice que dos variables son incorreladas cuando su covarianza es nula, es decir, cuando sus
variabilidades son ortogonales.
Ejercicio 5. Probar que −σi σj ≤ σij ≤ σi σj
Coeficiente de correlación: La desigualdad anterior invita a definir el denominado coefi-

ciente de correlación lineal
σij
ρij = (1.13)
σi σj
que tiene la virtud de ser adimensional y estar comprendido en todo caso entre -1 y 1. La
incorrelación se identifica con ρij = 0. Procuraremos utilizar los subíndices sólo cuando sea
estrictamente necesario.
Dado un vector aleatorio p-dimensional X, las posibles covarianzas componen una matriz
simétrica que puede definirse mediante
Cov[X] = E[(X − µ)(X − µ)0 ] (1.14)
10
cuya diagonal está compuesta por las diferentes varianzas. Suele denotarse por la letra Σ. Lo
mismo ocurre con los coeficientes de correlación, que componen una matriz de correlaciones
p × p simétrica cuya diagonal está compuesta por unos.
Ejercicio 6. ¿Por qué es simétrica Σ? ¿Por qué la diagonal de la matriz de correlaciones está
compuesta por unos?
Es muy frecuente contemplar transformaciones de un vector aleatorio del tipo X̃ = A0 X +b,
con A ∈ Mp×m y b ∈ Rm .
Ejercicio 7. Probar que, en ese caso, el vector m-dimensional X̃ verifica
E[X̃] = A0 E[X] + b, Cov[X̃] = A0 Cov[X]A (1.15)
También es frecuente considerar una partición del vector aleatorio p-dimensional X en dos
vectores X1 y X2 de dimensiones p1 y p2 , respectivamente, lo cual da lugar a su vez a particiones
obvias de la media y la matriz de covarianzas:

X1 µ1 Σ11 Σ12
X= , µ= , Σ= (1.16)
X2 µ2 Σ21 Σ22
En el caso particular p1 = 1, es decir, cuando X1 es una variable aleatoria real y X2 un vector
aleatorio (p − 1)-dimesional, la descomposición de Σ será de la forma
2
σ1 Σ12
Σ= (1.17)
Σ21 Σ22
En tal caso cabe definir el coeficiente de correlación lineal múltiple (al cuadrado) entre X1 y
X2 mediante
Σ12 Σ−1
22 Σ21
ρ212 = (1.18)
σ12
Se trata de una generalización del coeficiente de correlación simple (al cuadrado) que interpre-
taremos en la siguiente sección.
1.3. Principales parámetros muestrales

Aunque los conceptos que introduciremos a continuación son completamente análogos a
los propuestos en la sección anterior (con la salvedad de que se definen en esta ocasión en el
espacio Euclídeo Rn , correspondiente a las muestras aleatoria de tamaño n), conviene dedicarle
al menos un pequeño espacio aunque sea sólo por aclarar la notación.
En el caso univariante (p = 1), dado un vector X ∈ Rn definimos su media y varianza
mediante n n
1X 2 1X
x= xi , s = ( xi − x) 2 (1.19)
n i=1 n i=1
Estos parámetros gozan de interpretación análoga a la de sus homólogos probabilísticos. En
general, si X = (X[1], . . . , X[p]) denota una matriz de dimensiones n × p, x denotará al vector de
Rp (matriz de dimensiones p × 1) cuyas componentes son las medias de las diferentes columnas
de X ordenadas, es decir, x = (x[1], . . . , x[p])0 . En ese caso, se denotará por X la matriz n × p
definida mediante 1n · x0 , es decir, la que consiste en reemplazar cada dato de X por la media
aritmética de su repectiva columna.
11
Capítulo 1 Regresión lineal
Covarianza muestral: Dados X, Y ∈ Rn , se denota
sxy = n−1 hX − X, Y − Yi
= n−1 (X − X)0 (Y − Y)
n
X
−1
= n (Xi − x)(Yi − y)
i=1
En general, si X denota una matriz dimensiones n×p, se define la matriz de covarianzas muestral
S = n−1 (X − X)0 (Y − Y) ∈ Mp×p (1.20)

Si se contempla una partición concreta de las columnas de X, la matriz S se descompondría de
manera análoga a (1.17), lo cual da pie a la definición de nuevos parámetros como el coeficiente
de correlación múltiple.
Coeficientes de correlación simple y múltiple: Dados X1 ∈ Rn y X2 ∈ Mn×q se define

2
R12 o R2 a secas mediante
−1
S12 S22 S21
R2 = 2
(1.21)
s1
En el caso particular q = 1 estaríamos hablando del cuadrado del coeficiente de correlación
simple, que se define mediante r = sxy /sx sy . Por último, si X es una matriz p × p se denota
por R la matriz compuesta por los coeficientes de correlación lineal simple entre sus diferentes
columnas, que puede expresarse matricialmente mediante R = S −1/2 · S · S −1/2 .
Análogamente a (1.15), las transformaciones lineales de un vector aleatorio p-dimensional
X, es decir, del tipo A0 X, con A ∈ Mp×m , se traducen en términos de una muestra de tamaño
X ∈ Mn×p de tamaño n en una transformación del tipo XA.
Ejercicio 8. Probar que, en ese caso, se verifica
X · A = X · A, Sx·A = A0 Sx A (1.22)
1.4. Regresión lineal

Consideremos un vector aleatorio X descompuesto en X1 y X2 según (1.16) con p1 = 1, es
decir, tenemos una variable aleatoria real X1 y un vector X2 p2 -dimensional. Nuestra intención
es explicar la variabilidad de X1 como función lineal de la variabilidad de X2 , en la medida de
lo posible. Por lo tanto, buscamos el vector β ∈ Rp2 que alcance el siguiente mínimo, en cuyo
caso se denominará solución mínimo-cuadrática:
mı́n kX1 − E[X1 ] − b0 (X2 − E[X2 ])k2 : b ∈ Rp2

(1.23)
Ecuación de regresión lineal: La solución se obtiene proyectando ortogonalmente el vector

aleatorio X1 − E[X1 ] sobre el subespacio hX2 − E[X2 ]i, como indica la figura 1.2. Se trata pues
de buscar el vector β tal que
X1 − E[X1 ] − β 0 (X2 − E[X2 ]) ⊥ hX2 − E[X2 ]i (1.24)
12
E = h1i⊥ 6

X1 − E[X1 ] X1 − E[X1 ] − β 0 (X2 − E[X2 ])

1

v
β 0 (X2 − E[X2 ])
0

hX2 − E[X2 ]i
Figura 1.2: Ecuación de regresión lineal
En el caso p2 = 1 (regresión simple), el problema se reduce a encontrar el escalar β ∈ R tal

que
hX1 − E[X1 ], X2 − E[X2 ]i = β · hX2 − E[X2 ], X2 − E[X2 ]i (1.25)
−1
que, según (1.12), es β = σ22 σ21 .
Ejercicio 9. Probar que, en general, la ortogonalidad (1.24) se alcanza con el vector
β = Σ−1
22 Σ21 (1.26)
Por otra parte, si se define

α = E[X1 ] − β 0 E[X2 ], (1.27)
se verifica que X1 = α + β 0 X2 + E, donde E, denominado variabilidad parcial de X1 dado X2 ,

tiene media nula. Respecto a var[E], denominada varianza parcial, dado que kX1 −E[X1 ]k2 = σ12
y kβ 0 (X2 − E[X2 ])k2 = Σ12 Σ−1
22 Σ21 , se deduce que la proporción de variabilidad de X1 explicada
linealmente por la variabilidad de X2 es ρ212 , definido en (1.18). En definitiva, tal y como se
ilustra en la figura 1.3,
var[E] = σ12 − Σ12 Σ−1

22 Σ21 (1.28)
2 2
= σ1 (1 − ρ12 ) (1.29)
13
X1 − E[X1 ]
E = h1i⊥ 6

σ12 σ12 (1 − ρ212 )

1

v
σ12 ρ212
0

hX2 − E[X2 ]i
Figura 1.3: Descomposición de la varianza
Coeficientes de correlación parcial Razonando de manera análoga en el caso general

p1 ≥ 1, obtenemos que la matriz de covarianzas de E, denominada matriz de covarianzas
parciales, es la siguiente
Cov[E] = Σ11 − Σ12 Σ−1
22 Σ21 (1.30)
que se denota Σ11·2 . Si p1 > 1, los coeficientes de correlación asociados a esta matriz se de-
nominan correlaciones parciales entre las componentes de X1 dado el vector aleatorio X2 , que
equivalen pues a las correlaciones simples entre las componentes del vector aleatorio E. Así
pues, si consideramos dos componentes de X1 [i] y X1 [j] de X1 , su correlación parcial al cua-
drado dado X2 , que se denota por ρ2ij·2 , se interpreta como la proporción de la variabilidad
parcial de X1 [i] dado X2 explicada linealmente por la variabilidad parcial de X1 [j] dado X2 , y
viceversa. En un lenguaje más intuitivo podríamos definirlo como la proporción de variabilidad
que comparten “en exclusiva”, al margen de lo ya explicado por X2 .
Asimismo, podríamos generalizar el coeficiente de correlación múltiple como la matriz
Σ−1 −1
11 Σ12 Σ22 Σ21 , pero de dicha matriz sólo nos interesarán sus autovalores, que denominare-
mos en el capítulo 3 coeficientes de correlación canónica.
Ejercicio 10. Probar que α + β 0 X2 es la función afín de X2 que minimiza la distancia respecto
a X1 .
Incorrelación e independencia: Se dice que X1 y X2 son incorreladas cuando Σ12 = 0,

lo cual equivale a β = 0 o, si p1 = 1, ρ12 = 0. Se traduce por tanto en la imposibilidad de
14
Capítulo 1 Nociones básicas de Álgebra Lineal
explicar parte alguna de la variabilidad de X1 como función lineal de la variabilidad de X2 .

Geométricamente puede definirse así:
X1 y X2 incorreladas ⇔ hX1 − E[X1 ]i ⊥ hX2 − E[X2 ]i (1.31)
La independencia supone sin embargo una propiedad estrictamente más fuerte que la incorre-
lación. Efectivamente, puede ocurrir que entre X1 y X2 no se dé relación afín alguna pero que,
sin embargo, exista entre ambas una relación de otro tipo, que podría ser incluso funcional.
E[X1 |X2 ] ◦ X2 es la función medible de X2 que mejor se aproxima a X1 según la métrica (1.4) y
no podemos en general afirmar que se trate de una función afín. Eso sí ocurre bajo el supuesto
de (p1 + p2 )-normalidad, como veremos en el próximo capítulo. En ese caso, debe verificarse
entonces E[X1 |X2 ] ◦ X2 = α + βX2 , con α y β definidas como antes.
El concepto probabilístico de independencia lo suponemos conocido. Desde un punto de
vista geométrico, podría definirse como sigue: primeramente, dado un vector k-dimensional
Y con componentes en L2 , denótese por M(Y ) el espacio de las variables en h1i⊥ que son
funciones medibles de Y . En tal caso, se verifica
X1 y X2 independientes ⇔ M(X1 ) ⊥ M(X2 ) (1.32)
Al igual que podemos establecer una conexión entre el coeficiente de correlación y el concep-
to de independencia, podríamos asociar el coeficiente de correlación parcial al de independencia
condicional.
Ejercicio 11. Probar (1.31) y (1.32). Deducir entonces que la independencia implica incorre-
lación.
Ejercicio 12. Interpretar los parámetros muestrales definidos en la sección 3 en los términos
de la sección 4.
1.5. Nociones básicas de Álgebra Lineal

Aparte de los conceptos introducidos en la primera sección debemos destacar algunas no-
ciones y resultados propios del Álgebra Lineal que se manejan con frecuencia en nuestra teoría.
Hemos de tener presente en todo momento tres observaciones: primero que, fijada una base
vectorial en Rm , las aplicaciones lineales de Rm en Rm se identifican con las matrices cuadradas
de orden m; segundo, que una vez fijado un orden de lectura, el conjunto Mn×p de matrices
de dimensión n × p se identifica con Rnp ; tercero, que dicha identificación permite definir un
producto interior en Mn×p mediante
hA, Bi = tr(A0 B) (1.33)

X
= aij bij (1.34)
i,j
Este producto interior permite generalizar la distancia (1.5) al conjunto Mn×p mediante:
d2n,p (A, B) = n−1 tr[(A − B)0 (A − B)] (1.35)

n
X
−1
= n kai − bi k2Rp (1.36)
i=1
15
donde a0i y b0i denotan las filas de A y B, respectivamente. Esta distancia generalizada puede
entenderse a su vez como una versión muestral de la distancia (1.10). Entre otras propiedades,
podemos destacar que tr(A0 B) = tr(B 0 A) y que, si A, B, C son matrices cuadradas de orden
m, se verifica que tr(ABC) = tr(CAB) = tr(BCA).
Ejercicio 13. Probar (1.34) y (1.36).
Ejercicio 14. Dada una matriz de datos X ∈ Mn×p , probar que la varianza total muestral
de X, definida de manera análoga a (1.9) como la suma de las varianzas muestrales de sus
p-componentes, verifica
s2T = d2n,p (X, X) (1.37)
Matriz positiva: En el conjunto de matrices cuadradas m × m, podemos definir el siguiente

preorden que generaliza el orden natural en R: decimos que A ≥ B cuando x0 Ax ≥ x0 Bx para
todo x ∈ Rm . Así mismo, decimos que A > B cuando la desigualdad anterior es estricta si
x 6= 0. En consecuencia, A ≥ 0 cuando x0 Ax ≥ 0 para todo x ∈ Rm , en cuyo caso se dice que
A es semidefinida positiva. Si A > 0 se dice definida positiva.
Distancia de Mahalanobis: Dada una matriz A ∈ Mm×m simétrica y positiva podemos

definir en Rm la distancia de Mahalanobis DA
2
mediante
2
DA (x, y) = (x − y)0 A−1 (x − y), x, y ∈ Rm (1.38)
Se trata de una generalización de la métrica Euclídea, que se obtendría en el caso A = Id.
Matriz ortogonal: Se dice que una matriz Γ ∈ Mm×m es ortogonal cuando sus columnas
constituyen una base ortonormal de Rm , es decir, cuando Γ0 = Γ−1 . El conjunto de matrices
ortogonales de orden m se denotará por Om .
Matriz idempotente: Se dice que una matriz A ∈ Mm×m es idempotente cuando A2 = A.

Puede probarse que, si V es un subespacio lineal de Rm y B ∈ Mm×dim V es una base de V
(entendemos con esto que las columnas de B constituyen una base de V , es decir, que V = hBi),
entonces la matriz PV que se identifica con la proyección ortogonal sobre V puede calcularse
mediante
PV = B(B 0 B)−1 B 0 (1.39)
Se trata pues de una matriz simétrica e idempotente.
Ejercicio 15. Probar (1.39). Es más, probar que una matriz A ∈ Mm×m simétrica e idempo-
tente se identifica con la proyección ortogonal sobre V = hAi.
Autovalores y autovectores: Dada una matriz A ∈ Mm×m , se dice que δ ∈ R es un

autovalor real de A y γ ∈ Rm un autovector asociado cuando se verifica que Ae = δ · γ. En
tal caso, δ debe ser necesariamente una raíz del polinomio p(x) = |A − x · Id| y hγi debe estar
incluido en ker(A − δ · Id). Puede probarse que, si A es simétrica, las m raíces de p(x) son
reales, lo cual equivale a la existencia de m autovalores reales contados con su multiplicidad. El
siguiente resultado, conocido como teorema de diagonalización de una matriz simétrica, aclara
la estructura de la familia de autovectores asociados.
16
Teorema 1.5.1. Dada una matriz A ∈ Mm×m simétrica, si ∆ denota la matriz diagonal
compuesta por los autovalores δ1 , . . . , δm de A ordenados de mayor a menor y contados con su
multiplicidad, existe una matriz Γ ∈ Om , cuyos vectores columnas se denotan por γ1 , . . . , γm ,
tal que
A = Γ∆Γ0 (1.40)
Se verifica además que δ1 = máx{γ 0 Aγ : kγk = 1}, que se alcanza con γ = γ1 , y que, para
todo j = 2, . . . , m, δj = máx{γ 0 Aγ : kγk = 1, γ ⊥ hγ1 , . . . , γj−1 i}, alcanzándose con γ = γj .
Del teorema se sigue directamente que las columnas de Γ constituyen una base ortonormal
de autovectores asociados a los correspondientes autovalores. También podemos de deducir de
(1.40) que ∆ = Γ−1 AΓ. Por lo tanto, la aplicación lineal identificada con la matriz A para la
base vectorial original admite una expresión diagonal respecto a una base ortonormal canónica
de autovectores. Es decir, el cambio a la base de autovectores permite expresar la matriz de
forma sencilla. A modo de ejemplo, podemos utilizar ese procedimiento para demostrar las
siguientes propiedades;
Ejercicio 16. Dada una matriz simétrica A, probar:
(i) Si A es simétrica, su rango coincide con el número de autovalores no nulos.
(ii) Si A ≥ 0, sus autovalores son todos no negativos. Si A > 0, son todos estrictamente
positivos.
(iii) Si A ≥ 0, existe una matriz simétrica A1/2 tal que A = A1/2 A1/2 . Si A > 0, existe
también una matriz simétrica A−1/2 tal que A−1 = A−1/2 A−1/2 .
(iv) Si A ≥ 0, existe una matriz X con las mismas dimensiones tal que A = X 0 X.
(v) La traza de una matriz simétrica es la suma de sus autovalores y el determinante, el

producto de los mismos.
(vi) La inversa de una matriz simétrica positiva también es positiva.
(vii) La matriz PV de la proyección ortogonal sobre el subespacio lineal V ⊂ Rn consta de

tantos autovalores 1 como dim V , siendo nulos los n − dim V restantes.
A partir del teorema 1.5.1 y del ejercicio 1 podemos probar el siguiente resultado en el cual
se fundamenta el capítulo 5:
Lema 1.5.2. En las condiciones del teorema 1.5.1 y dado k ≤ m, si Γ1 es la matriz con los
autovectores asociados a los k primeros autovalores de A, se verifica que
k
X
máx{tr(B 0 AB) : B ∈ Mm×k , B 0 B = Id} = δi (1.41)
i=1
y se alcanza en B = Γ1 .
17
Cuadro 1.1: Matriz de datos muestra tamaño n = 38 y dimensión p = 8
18
Capítulo 2
Modelo lineal multivariante
En este capítulo expondremos los aspectos más generales del modelo lineal normal multiva-
riante. Previamente, estudiaremos con brevedad las distribuciones de probabilidad relacionadas
con este modelo así como el modelo lineal normal (univariante) que pretende generalizar.
2.1. Distribución normal multivariante

La distribución normal multivariante p-dimensional o p-normal se trata de una generali-
zación natural de la distribución normal que servirá como hipótesis de partida en el modelo
estadístico objeto de estudio.
Dados µ ∈ Rp y Σ ∈ Mp×p simétrica y semidefinida positiva, se dice que un vector aleatorio
X p-dimensional sigue un modelo de distribución Np (µ, Σ) cuando admite la siguiente función
característica:
0 1 0
ϕX (t) = exp it µ − t Σt , t ∈ Rp (2.1)
2
En ese caso se denota X ∼ Np (µ, Σ) y puede comprobarse trivialmente que generaliza la distri-
bución normal unidimensional. Vamos a enunciar a continuación las propiedades fundamentales
de esta distribución. Las dos siguientes se siguen de las propiedades de la función característica.
Proposición 2.1.1. Si X ∼ Np2 (µ, Σ), A ∈ Mp1 ×p2 y b ∈ Rp1 , entonces
AX + b ∼ Np1 (Aµ + b, AΣA0 ) (2.2)
Proposición 2.1.2. Si Z[1], . . . , Z[p] iid N(0,1), entonces Z = (Z[1], . . . , Z[p])0 ∼ Np (0, Id)
A partir de las dos propiedades anteriores podemos construir cualquier vector normal:
Proposición 2.1.3. Dados µ y Σ como en la definición, si consideramos el vector aleatorio Z
anterior, la descomposición Σ = Γ∆Γ0 y se denota A = Γ∆1/2 , se sigue que AZ + µ ∼ Np (µ, Σ).
En consecuencia, se sigue de (1.15) el siguiente resultado:
Proposición 2.1.4. Si X ∼ Np (µ, Σ), E[X] = µ y Cov[X] = Σ.
También es consecuencia de la proposición 2.1.1 que, si X ∼ N (µ, Σ), cada componente
X[i] de X sigue un modelo de distribución N (µi , σi2 ). Sin embargo, el recíproco no es cierto.
Hemos de tener en cuenta que la componente X[i] puede obtenerse mediante e0i X, siendo ei el
vector unidad en el eje de coordenadas i-ésimo, y que la siguiente afirmación puede probarse
con relativa facilidad:
19
Capítulo 2 Distribución normal multivariante
Proposición 2.1.5. Dado un vector aleatorio p-dimensional X, cualquiera de las condiciones

siguientes garantizan la p-normalidad del mismo:
(i) a0 X es 1-normal, para todo a ∈ Rp .
(ii) Sus componentes son todas normales e independientes entre sí.

El siguiente resultado puede probarse también a través de la función característica y esta-
blece la equivalencia entre incorrelación e independencia bajo la hipótesis de normalidad.
Proposición 2.1.6. Si descomponemos un vector (p1 +p2 )-normal X con matriz de covarianzas
Σ en X1 de dimensión p1 y X2 de dimensión p2 , entonces X1 y X2 son independientes s, y sólo
si, Σ12 = 0.
Si la matriz de covarianzas es estrictamente positiva, la distribución p-normal es dominada
por la medida de Lebesgue en Rp . Teniendo en cuenta las proposiciones 2.1.1, 2.1.2 y el teorema
del cambio de variables, podemos obtener la densidad de dicha distribución:
Proposición 2.1.7. Si X ∼ Np (µ, Σ) con Σ > 0 admite la siguiente función de densidad:

1 1 0 −1
f ( x) = p exp − (x − µ) Σ (x − µ) , x ∈ Rn . (2.3)
p
(2π) |Σ| 2
z
Figura 2.1: Función de densidad distribución 2-normal
Ejercicio 17. Probar las siete proposiciones anteriores.

Nótese que en la función de verosimilitud determinada por (2.3) la observación x y los
parámetros (µ, Σ) que caracterizan la distribución de probabilidad se relacionan a través de la
distancia de Mahalanobis DΣ2 (x, µ). Concretamente, para cada k ∈ [0, [(2π)p |Σ|]−1/2 ], la región
de los puntos {x ∈ Rp : f (x) = k}, es decir, aquéllos cuya densidad es igual a k, es el elipsoide
siguiente:
{x ∈ Rp : DΣ2 (x, µ) = k̃} (2.4)
p
para k̃ = −2 log k (2π)p |Σ| .
20
Capítulo 2 Distribución normal multivariante
2
1
0
y
-1
-2
-3
-3 -2 -1 0 1 2
Figura 2.2: Contornos distribución 2-normal
En la figura 2.1 se aprecia la función de densidad de la distribución bidimensional

0 1 0.8
N2 , (2.5)
0 0.8 1
mientras que en la figura 2.2 se podemos ver un diagrama de dispersión con una muestra
aleatoria simple de tamaño n = 150 de dicha distribución en la que aparecen marcados dos
contornos elípticos de la misma.
Consideremos un vector aleatorio X (p1 + p2 )-normal que descompone de la forma

X1 µ1 Σ11 Σ12
X= ∼ Np1 +p2 , (2.6)
X2 µ2 Σ21 Σ22
El siguiente resultado puede probarse teniendo en cuenta el hecho conocido de que la densidad
de la distribución condicional P X1 |X2 puede calcularse mediante
fX1 ,X2 (x1 , x2 )

fX1 |X2 =x2 (x1 ) = (2.7)
fX2 (x2 )
Proposición 2.1.8. Si Σ22 > 0, se verifica
X1 |X2 = x2 ∼ Np1 (α + β 0 x2 , Σ11·2 ) (2.8)
con α, β y Σ11·2 definidas según (1.27), (1.26) y (1.30), respectivamente.
Como consecuencia se deduce que, bajo el supuesto de normalidad, E[X1 |X2 ]◦X2 = α+β 0 X2 .
Es más, podemos garantizar que
X1 = α + β 0 X2 + E, E y X2 independientes, E[E] = 0, Cov[E] = Σ11·2 (2.9)
Esta afirmación puede probarse también teniendo en cuenta la proposición 2.1.6, (1.31) y (1.32).
En definitiva, establece una clara conexión entre los conceptos de normalidad y linealidad.
21
Capítulo 2 El modelo lineal
Ejercicio 18. Si X denota un vector 2-normal siguiendo un modelo de distribución (2.5),

razonar qué modelo de distribución sigue en cada caso el vector Y indicando, si procede, su
función de densidad:
(a) Y [1] = 1 + 2X[1] + 3X[2]; Y [2] = 4 − X[1] + X[2]
(b) Y [1] = 2 + 5X[1] − 4X[2]
(c) Y [1] = 1 + 2X[1] + 3X[2]; Y [2] = 4 − 4X[1] − 6X[2]
Ejercicio 19. Simular de manera aproximada una muestra de tamaño n = 200 de la distribu-
ción (2.5).
Desde el punto de vista estadístico, podemos proponer tests para contrastar la hipótesis
inicial de normalidad multivariante. En Bilodeau y Brenner (1999) se recoge un test que se
basa en el hecho de que, para una muestra aleatoria simple de tamaño n de una distribución
p-normal, las distancias de Mahalanobis entre las observaciones y la media aritmética de la
misma dada la matriz de covarianzas muestral siguen una distribución de la familia Beta y
tienden a la incorrelación conforme aumenta el tamaño de muestra. Desde una perspectiva
eminentemente práctica, si realmente tenemos la intención de utilizar alguno de los procedi-
mientos de tipo paramétrico que expondremos a continuación, resulta más realista comprobar
que los diagramas de dispersión entre las diferentes componentes revelan al menos relaciones
de tipo lineal, estando muy pendiente de la presencia de sesgos, que pueden conducirnos a
transformar las variables originales, o fragmentaciones de los datos, que pueden conducirnos a
introducir factores cualitativos en el modelo.
2.2. Modelo lineal

Antes de abordar el estudio del modelo lineal multivariante, repasaremos muy brevemente
el modelo lineal en dimensión 1, empezando por las distribuciones de probabilidad asociadas
al mismo.
2.2.1. Distribuciones asociadas al modelo

No pretendemos aquí describir con detalle las diferentes distribuciones pues se suponen
conocidas (ver por ejemplo Nogales (1998)), sino interpretarlas desde un punto de vista geo-
métrico.
Distribución normal esférica: El punto de partida del modelo es la distribución normal

multivariante esférica, que se obtiene cuando la matriz de covarianzas del vector es de la forma
Σ = σ 2 Id, para algún σ 2 > 0. Efectivamente, puede comprobarse que, en ese caso, la distancia
de Mahalanobis DΣ2 es, salvo una homotecia, la distancia Euclídea, por lo que los elipsoides
(2.4) son en particular esferas cuyo centro es la media. En virtud de la proposición 2.1.5-(ii),
Y ∼ Nn (µ, σ2 Id) si, y sólo si, sus componentes Y1 , . . . , Yn son independientes, normales y con
idéntica varianza (homocedásticos). Es decir, que esta distribución está asociada a una muestra
de tamaño n en sentido amplio, de ahí la notación utilizada.
22
psico1
psico2
psico3
psico4
psico5
psico6
psico1 psico2 psico3 psico4 psico5 psico6
Figura 2.3: Simulación muestra n = 100 distribución 6-normal
Distribución χ2 : Cuando la media es 0, la distribución normal esférica tiene además la

particularidad de ser invariante ante cualquier transformación mediante una matriz ortogonal,
es decir, que la verosimilitud de cada observación Y depende exclusivamente de kYk2 . Eso explica
nuestro interés en estudiar la distribución de kYk2 bajo la hipótesis de normalidad esférica.
Efectivamente, si Y ∼ Nn (0, Id), se dice que kYk2 ∼ χ2n ; la anterior distribución denominada
χ2 central puede generalizarse si consideramos la norma Euclídea al cuadrado de un vector
Y ∼ Nn (µ, Id), que se distribuye según un modelo χ2n (δ), con δ = kµk2 . Si Y ∼ Nn (µ, σ2 Id),
entonces σ −2 kYk2 / ∼ χ2n (δ) con δ = kµk2 /σ 2 , lo cual se denota por kYk2 ∼ σ 2 χ2n (δ).
Ejercicio 20. Probar que, en general, si Y ∼ Nn (µ, σ 2 Id) y E ⊂ Rn ,
1 1
2
kPE Yk2 ∼ χ2dim E (δ), δ= 2
kPE µk2 (2.10)
σ σ
Por otra parte, se verifica además que E [kPE Yk2 ] = (dim E)σ 2 (1 + δ). Como caso particular,
si µ ∈ E ⊥ , entonces kPE Yk2 ∼ σ 2 χ2dim E .
Ejercicio 21. Probar que, dados E1 ⊥ E2 y X ∼ Nn (µ, σ 2 Id), se verifica que kPEi Yk2 ∼
σ 2 χ2dim Ei (kPEi µk2 /σ 2 ),
para i = 1, 2, y son independientes. Página 1
23
Distribución F : Este modelo probabilístico de distribución surge de manera natural de la

aplicación del Principio de Invarianza en el problema estadístico de contraste de una hipóte-
sis lineal para la media que veremos más adelante. Efectivamente, en la resolución de dicho
problema nos vemos avocados a considerar, dado un vector Y ∼ Nn (µ, σ 2 Id), el cociente entre
kPE1 Yk2 y kPE2 Yk2 para ciertos subespacios E1 y E2 ortogonales entre sí y tales que µ ∈ E2⊥ .
No obstante, ambos términos se normalizan dividiéndolos por sus respectivas dimensiones, de
manera que la distribución F se obtiene mediante
kPE1 Yk2 / dim E1 kPE1 µk2

∼ Fdim E1 ,dim E2 (δ), δ= (2.11)
kPE2 Yk2 / dim E2 σ2
Nótese que el cociente entre las medias del numerador y el denominador es (1 + δ) y, por lo
tanto, 1 cuando δ = 0. La distribución m · Fm,n converge a χ2m cuando n tiende a infinito.
E2 ⊂ hµi⊥
PE2 Y 6 PE1 ⊕E2 Y
'$
0
u
E1
-
PE1 Y
&%
E1 ⊕ E2 ⊂ Rn
Figura 2.4: Interpretación geométrica de la distribución F
Distribución t: En esta sección interpretaremos la conocida distribución t-Student como un

caso particular de la distribución F . Concretamente, decimos que una variable real t sigue un
modelo de distribución tm (δ) cuando es simétrico respecto a 0 y tal que t2 ∼ F1,m (δ). De esta
forma, nos encontraremos con dicha distribución cuando operemos como en la figura 2.4 con
dim E1 = 1.
2.2.2. El modelo y ejemplos

El modelo lineal (normal) consiste en una estructura o experimento estadístico en Rn donde
la distribución de probabilidad es normal esférica Nn (µ, σ 2 Id). No se impone ninguna condición
respecto al parámetro σ 2 pero si se impone una restricción de tipo lineal para el parámetro µ,
pues se supondrá por hipótesis que µ ∈ V para un cierto subespacio lineal conocido V ⊂ Rn .
Se denota mediante
Y ∼ Nn (µ, σ2 ), µ ∈ V, σ 2 > 0 (2.12)
24
La restricción lineal µ ∈ V vendrá dada, bien por la presencia de factores cualitativos, bien por
la relación lineal respecto a otras variables numéricas con valores conocidos.
Si una matriz X ∈ Mn×dim V constituye una base de V , podemos parametrizar el mode-
lo (2.12) a través de las coordenadas β de µ respecto a X, es decir, Y ∼ Nn (Xβ, σ 2 Id), o
equivalentemente,
Y = Xβ + E, β ∈ Rdim V , E ∼ Nn (0, σ2 Id), σ2 > 0 (2.13)
Enunciaremos a continuación cuatro ejemplo de problemas estadísticos que se formalizan me-

diante el modelo lineal:
Ejemplo 1. [Muestra aleatoria simple de una distribución normal] Consideremos Y1 , . . . , Yn iid

N (ν, σ 2 ). En ese caso, el vector aleatorio Y = (Y1 , . . . , Yn )0 sigue una distribución Nn (µ, σ 2 Id)
con µ ∈ V = h1n i y σ 2 > 0.
Ejemplo 2. [Muestras independientes de distribuciones normales con idéntica varianza] Consi-

deremos ahora, para i = 1, 2, sendas muestras independientes entre sí Yi1 , . . . , Yini iid N (µi , σ 2 ).
Si se denota n = n1 + n2 e Y = (Y11 , . . . , Y2n2 )0 , se verifica que Y ∼ Nn (µ, σ 2 Id) con σ 2 > 0 y
µ ∈ V = hv1 , v2 i, donde v1 denota el vector de Rn cuyas n1 primeras componentes son 1 y el
resto 0. De manera análoga se define v2 .
Ejemplo 3. [Diseño completamente aleatorizado] Se trata de una generalización del problema

anterior para r ≥ 2 muestras independientes Yi1 , . . . , Yini iid N (µi , σ 2 ). En este caso, si n =
0 2 2
P
i ni e Y = (Y11 , . . . , Yrnr ) , se verifica que Y ∼ Nn (µ, σ Id) con σ > 0 y µ ∈ V = hv1 , . . . , vr i.
Ejemplo 4. [Regresión lineal múltiple] Supongamos que se recogen n observaciones indepen-

dientes que pueden calcularse mediante una relación afín con los valores de otras q variables
numéricas controladas en el experimento, salvo errores independientes, normalmente distribui-
dos y homocedásticos. Es decir,
Y1 = β0 + β1 Z1 [1] + . . . + βq Z1 [q] +E1

.. .. .. .. .. .. (2.14)
. . . . . .
Yn = β0 + β1 Zn [1] + . . . + βq Zn [q] +En
donde E1 , . . . , En iid N (0, σ 2 ). Si se denota

       
Y1 1 Z1 [1] . . . Z1 [q] E1 β0
.  . .. ..  , E =  ...  , β =  ... 
 ..  ,
Y=  ..
X= (2.15)
   
. . 
Yn 1 Zn [1] . . . Zn [q] En βq
el modelo puede expresase de la forma (2.13) con V = hXi. En lo sucesivo se denotará β =

(β1 , . . . , βq )0 ∈ Rq . Este vector expresa la influencia de los vectores explicativos Z[1], . . . , Z[q] en
la predicción de la respuesta Y. Así mismo, se denota por Z la matriz que resulta al eliminar de
X el término independiente 1n . De esta, forma, el modelo puede expresarse también mediante
Y = β0 · 1n + Zβ + E (2.16)
25
Regresión respecto a variables dummys: Cualquiera de los problemas considerados an-

teriormente puede entenderse como un problema de regresión lineal, es decir, pueden para-
metrizarse de la forma (2.13) para una base X con término independiente. Así, en el caso del
ejemplo 3, podemos considerar entre otras posibilidades la matriz X = (1n , v1 , . . . , vr−1 ). Con
esta parametrización particular, la relación entre µ y β es la siguiente:
β0 = µr , βj = µj − µr , j = 1, . . . , r − 1 (2.17)
Ejercicio 22. Probar (2.17). Indicar así mismo cómo se relacionaría µ con β si consideráramos
la base natural X̃ = (1n , v1 − vr , . . . , vr−1 − vr ).
Los vectores Z[1], . . . , Z[r − 1] de X en la parametrización anterior recogen valores concretos
de unas variables denominadas dummys que indican la muestra o categoría a la que pertenece
cada dato. Que las medias µ1 , . . . , µr sean idénticas, es decir, que las muestras procedan de una
única distribución común, equivale a que β sea nulo, independientemente de la parametrización
particular considerada. En otras palabras, la ausencia de relación entre el factor cualitativo que
distingue las muestras con la variable numérica Y equivale a la ausencia de relación de ésta
con las variables numéricas dummys.
Ejercicio 23. Desarrollar con detalle los modelos asociados a los cuatro ejemplos anteriores.
2.2.3. Estimación y contraste de hipótesis

Dado que suponemos µ ∈ V y que, al seguir Y un modelo de distribución n-normal, es
más verosímil que la observación Y sea próxima a la media que lo contrario, parece razonable
estimar µ mediante un vector de V próximo a Y. De hecho, definimos el estimador
µ̂ = PV Y (2.18)
En tal caso, resulta también razonable estimar σ 2 mediante la distancia (1.5) entre Y y µ̂, es
2 −1 2
decir, σ̂M V = n kPV ⊥ Yk . Puede probarse que ambos estimadores son independientes y que
constituyen un estadístico suficiente y completo. Se sigue entonces del teorema de Lehmann-
Scheffe, que µ̂ es el estimador insesgado de mínima varianza (EIMV) de µ. También puede
2
probarse a partir de (2.3) que (µ̂, σ̂M V ) constituyen un estimador de máxima verosimilitud
2 2
(EMV) de (µ, σ ). Sin embargo, σ̂M V no es insesgado, de ahí que se proponga el siguiente
estimador que es, según el teorema de Lehmann-Scheffe, EIMV:
1
σ̂ 2 = kP ⊥ Yk2 (2.19)
n − dim V V
Si el modelo está parametrizado de la forma (2.13), podemos estimar β como las coordenadas
del estimador de µ, es decir:
β̂ = (X0 X)−1 X0 Y (2.20)
En definitiva, los estimadores de µ y σ 2 pueden entenderse geométricamente según la figura 2.4
con E1 = V y E2 = V ⊥ .
Ejercicio 24. Obtener los estimadores µ̂ y σ̂ 2 para los ejemplos 1 y 2.
Ejercicio 25. Obtener µ̂ para el ejemplo 3. Probar que, en dicho ejemplo, el EIMV de σ 2 es
i r n
2 1 XX
σ̂ = (Yij − yi· )2 (2.21)
n − r i=1 j=1
26
Ejercicio 26. Probar que, en el ejemplo 4, podemos estimar β a partir de las medias aritméticas
del vector Y y la matriz Z, así como de la matriz de covarianza muestral conjunta mediante
−1
β̂ = Szz Szy , β̂0 = y − z0 β̂ (2.22)
Para resolver el anterior ejercicio se aconseja utilizar un razonamiento análogo al que se ilustra
en la figura 1.2, pero en términos muestrales, y tener en cuenta la descomposición ortogonal
hXi = h1n i ⊕ hZ − Zi.
2 1 2
Ejercicio 27. Probar que, en el ejemplo 4, σ̂M V puede relacionarse con la varianza sy del
2
vector Y y el coeficiente de correlación múltiple al cuadrado R de Y respecto a Z, definido en
(1.21), mediante
2 −1 2 2
σ̂M V = n (n − 1)sy (1 − R ) (2.23)
El problema de contraste de hipótesis relativas al parámetro σ 2 no será expuesto aquí debido
a que los tests que los resuelven son sensibles ante la violación del supuesto de normalidad.
No ocurre lo mismo con el test F o anova que resuelve el contraste de hipótesis de tipo lineal
sobre el parámetro µ pues, tal y como se prueba en Arnold (1981), es asintóticamente válido
aunque no se verifique el supuesto de normalidad. Además, es relativamente robusto ante la
heretocedasticidad. Lo mismo ocurre en el modelo multivariante.
Anova: Nos ocuparemos pues del contraste de hipótesis tipo H0 : µ ∈ W , para algún subes-
pacio lineal W ⊂ V . Veamos ejemplos de hipótesis de este tipo:
Ejercicio 28. En el ejemplo 1 podemos contrastar si la media ν de la distribución es nula.
Probar que se corresponde con W = 0.
Ejercicio 29. En los ejemplos 2 y 3 podemos contrastar si todas las muestras consideradas pro-
vienen de una misma distribución de probabilidad. Probar que en ambos casos se corresponde
con W = h1n i.
Ejercicio 30. En el ejemplo 4 podemos contrastar si los vectores explicativos Z[1], . . . , Z[q]
no intervienen en la explicación de Y, lo cual equivale a β = 0. Probar que se corresponde con
W = h1n i. Dicho contraste se denomina total.
Ejercicio 31. En las condiciones del ejemplo 4 podemos contrastar también hipótesis del tipo
βj = 0. Probar que, por ejemplo, en el caso j = q, se corresponde con W = h1n , Z[1], . . . , Z[q −
1]i. Dicho contraste se denomina parcial.
Si se denota V |W = W ⊥ ∩ V, la hipótesis inicial H0 : µ ∈ W equivale a PV |W µ = 0. Ello
invita a descomponer Rn en tres subespacios ortogonales: Rn = W ⊕ V |W ⊕ V ⊥ . De dicha
descomposición se deriva la siguiente descomposición ortogonal del vector de observaciones:
Y = PW Y + PV |W Y + PV ⊥ Y (2.24)
Si W = 0, como en el ejercicio 28, la descomposisión (2.24) se reduce a los dos últimos sumandos.
Teniendo en cuenta el Principio de Invarianza y el hecho de que (µ̂, σ̂ 2 ) es suficiente, se
sigue que la decisión correspondiente al contraste de H0 debe depender de la observación Y
a través del cociente entre kPV |W Y k2 y kPV ⊥ Y k2 . Si dividimos ambos números positivos por
Pn
1
Nos referimos al estimador insesgado de la varianza s2y = (n − 1)−1 Yi − y)2 .
i=1 (
27
sus respectivos grados de libertad obtenemos la distribución F . En definitiva, consideramos el

estadístico de contraste
n − dim V kPV |W Yk2
F (Y) = · (2.25)
dim V |W kPV ⊥ Yk2
o, equivalentemente,
kPV |W Yk2 /(dim V |W )
F (Y) = (2.26)
σ̂ 2
que, según (2.11), sigue en general un modelo de distribución Fdim V |W,n−dim V (δ), con δ =
kPV |W µk2 /σ 2 . En particular, bajo la hipótesis inicial sigue una distribución Fdim V |W,n−dim V .
Teniendo en cuenta el Principio de Máxima Verosimilitud, construimos el denominado test F o
anova de manera que se rechace la hipótesis inicial si el estadístico F toma valores extremos, los
cuales serán tanto más verosímiles cuanto más nos alejemos de la hipótesis inicial (δ = 0). Del
lema fundamental de Neyman-Pearson se deduce que el test así construido es UMP-invariante;
además, es el test de la razón de verosimilitudes (TRV). Desde el punto de vista geométrico
puede entenderse según la figura 2.4 con E1 = V |W y E2 = V ⊥ .
El caso de mayor interés práctico es W = h1n i (como en los ejercicios 29 y 30), en el cual
la decomposión (2.24) de kYk2 es la siguiente:
kYk2 = kPh1n i Yk2 + kPh1n i⊥ Yk2 (2.27)

2 2 2
= ny + kPV |h1n i Yk + kPV ⊥ Yk (2.28)
Ejercicio 32. Probar que, si consideramos cualquier matriz Z tal que de h(1n |Z)i = V y R2
denota el coeficiente de correlación múltiple entre Y y Z, como en (2.23), la igualdad (2.28)
puede expresarse también así
kYk2 = ny2 + (n − 1) R2 s2y + (1 − R2 )s2y

(2.29)
Nótese que el Principio de Invarianza nos conduce a desentendernos del primer sumando,
correspondiente a un vector en h1n i, y concentranos en el cociente entre el segundo y el tercero,
tal y como se ilustra mediante la figura 2.5, que es la versión muestral de la figura 1.3.
Este razonamiento nos ayuda a comprender que el valor del coeficiente de correlación li-
neal múltiple R2 depende exclusivamente del subespacio V (conteniendo a h1n i) sobre el que
deseemos proyectar las observaciones, independientemente de la matriz (1n |Z) que escojamos
como base del mismo, lo cual invita a entender de manera más general la definición propuesta
en (1.21). Concretamente, dados Y ∈ Rn y V ⊂ Rn conteniendo a h1n i, se puede definir el
coeficiente de correlación de Y respecto a V mediante
kPV |h1n i Yk2
R2 = (2.30)
kPh1n i⊥ Yk2
En el caso particular dim V |W = 1, es decir, cuando W es un hiperplano de V , el estadístico

de contrate se denota por t2 en lugar de F pues se confronta con la distribución t2n−dim V , dando
lugar a lo que conocemos como test de Student.
Ejercicio 33. Relacionar la descomposición (2.28) con los términso de la tabla 2.1.
28
E = h1n i⊥ Y−Y
6

s2y (1 − R2 )s2y

1

t R2 s2

y
0

V |h1n i
Figura 2.5: Descomposición de la varianza muestral
Ejercicio 34. Resolver el contrate de la hipótesis inicial H0 : ν = 0 en el ejemplo 1; resolver

el contrate de la hipótesis inicial H0 : µ1 = µ2 en el ejemplo 2.
Ejercicio 35. Probar que el test anova que resuelve el contraste H0 : µ1 = . . . = µr en el
ejemplo 3 consiste en confrontar con la distribución Fr−1,n−r el estadístico de contraste
SCH/(r − 1)
F = , (2.31)
SCE/(n − r)
donde
X
SCH = ni (Yi· − y·· )2 , (2.32)
i
XX
SCE = (Yij − yi· )2 (2.33)
i j
Ejercicio 36. Probar que el test anova que resuelve el contrate para W = h1n i consiste en
confrontar con la distribución Fdim V −1,n−dim V el estadístico de contraste
n − dim V R2
F = · . (2.34)
dim V − 1 1 − R2
con R2 definido según (2.30). En particular, para contrastar H0 : β = 0 en el ejemplo 4 se

utiliza el estadístico
n − (q + 1) R2
F = · . (2.35)
q 1 − R2
¿Qué sucede en el caso particular q = 1?
Ejercicio 37. En las condiciones del ejemplo 4, ¿qué distribución sigue bajo la hipótesis inicial
H0 : βq = 0 el estadístico de contraste del anova correspondiente? Probar que el estadístico
de contraste correspondiente puede expresarse de manera análoga a (2.35) pero mediante el
29
coeficiente de correlación parcial muestral entre Y y Z[q], dados Z[1], . . . Z[q − 1], que se denota
por ry,z[q]·(z[1],...,z[q−1]) , o abreviadamente por rq∗ , mediante
2
rq∗
F = [n − (q + 1)] · 2
(2.36)
1 − rq∗
Ejercicio 38. Probar que2 el estadístico de contraste (2.37) puede expresarse también a través
de β̂q y del coeficiente de correlación múltiple de Z[q] respecto al resto de variables explicativas
Z[1], . . . , Z[q − 1], que se denota abreviadamente por Rq∗ , mediante
β̂q2 2
F =n· · (1 − Rq∗ ) (2.37)
σ̂ 2 · s−2
z[q]
2
El término 1−Rq∗ que aparece a la derecha en (2.37) se denomina tolerancia y su inverso, factor
de inflación de la varianza (FIV). (2.37) y (2.36) vienen a expresar pues , en términos intuitivos,
que una redundancia lineal entre las variables explicativas, que a su vez puede asociarse a bajas
correlaciones parciales, resta fiabilidad a las estimaciones de los coeficientes de regresión.
Por otra parte, teniendo en cuenta que, en las condiciones del ejemplo 3, la hipótesis inicial
H0 : µ1 = . . . = µr equivale a β = 0 para cualquier parametrización del modelo mediante
variables dummys, se sigue de (2.35) que la decisión al respecto se puede expresar a través de
la correlación múltiple R2 de Y con dichas variables dummys, independientemente de cuáles se
elijan. Esto hecho, que tendrá lugar igualmente en el modelo multivariante, justifica el estudio
de los coeficientes de correlación canónicos.
En la salida de SPSS recogida en el cuadro 2.1 podemos apreciar muchos de los ingredientes
estudiados en la sección.
Pruebas de los efectos inter-sujetos
Variable dependiente: sepleng
Suma de
cuadrados Media
Fuente tipo III gl cuadrática F Significación
Modelo corregido 63,212a 2 31,606 119,265 ,000
Intersección 5121,682 1 5121,682 19326,505 ,000
species 63,212 2 31,606 119,265 ,000
Error 38,956 147 ,265
Total 5223,850 150
Total corregida 102,168 149
a. R cuadrado = ,619 (R cuadrado corregida = ,614)
Cuadro 2.1: Tabla anova; sepleng vs especies en irisdata.sav
Ejercicio 39. Construye mediante SPSS dos variables dummys para distinguir las tres especies
de flores de irisdata y comprueba que el coeficiente de correlación múltiple R2 entre sepleng y
dichas variables es el que aparece en la tabla 2.1.
2
Ver Manual de Modelos Lineales, Montanero (2008).
30
Capítulo 2 Modelo multivariante
Ejercicio 40. En el cuadro 2.2 aparece el resultado del anova para comparar los valores medios
de glucemia de cuatro categorías de recién nacidos (control, respiratoria, metabólica y mixta)
mediante un diseño equilibrado. Relacionar los valores que aparecen en dicha tabla con las
columnmas de la matriz de datos del cuadro 2.6, donde a los datos originales se les ha añadido
las proyecciones relacionadas con la descomposición (2.24), las sumas de cuadrados correspon-
dientes a las descomposición (2.28) y las variables dummys asociadas a la parametrización del
ejercicio 22.
Pruebas de los efectos inter-sujetos

Variable dependiente:Nivel de glucemia en el cordón umbilical
Suma de
cuadrados Media
Origen tipo III gl cuadrática F Sig.
a
Modelo corregido 643,930 3 214,643 4,557 ,009
Intersección 170156,250 1 170156,250 3612,563 ,000
Enfermedad 643,930 3 214,643 4,557 ,009
Error 1507,240 32 47,101
Total 172307,420 36
Total corregida 2151,170 35
a. R cuadrado = ,299 (R cuadrado corregida = ,234)
Cuadro 2.2: Tabla anova; glucemia vs enfermedad en acidosis-SPSS.sav
2.3. Modelo multivariante

Una vez repasado el modelo lineal univariante estamos en condiciones de generalizarlo al
caso multivariante, en el cual no contamos con una sino con p variables respuesta. Previamente,
debemos introducir con brevedad las distribuciones de probabilidad asociadas al nuevo modelo.
Para un estudio más detallado, consultar Arnold (1981), Anderson (1958) y Mardia et al. (1979).
En lo que sigue supondremos en todo momento Σ > 0 y n ≥ p.
2.3.1. Distribuciones asociadas al modelo

Seguimos el mismo esquema que en el caso unidimensional, con la salvedad de que no existe
una distribución que generalice unívocamente la distribución F . Debemos tener en cuenta que,
en nuestro modelo estadístico, la observación es una matriz Y ∈ Mn×p de datos como la que
aparece en el cuadro 1.1, que se denotará como en (1.1).
Distribución normal matricial: Se trata de la distribución de partida del modelo, al igual

que la normal esférica lo era en el caso univariante. Dados µ ∈ Mn×p y Σ ∈ Mp×p simétrica
y definida positiva, se dice que Y ∼ Nn,p (µ, Id, Σ) cuando Yi ∼ Np (µi , Σ), i = 1, . . . , n, siendo
todas independientes. Esta distribución es un caso particular de otra más general que se trata
31
Figura 2.6: Proyecciones en acidosis-SPSS.sav
con detalle en Arnold (1981). La función de densidad se define, para cada matriz X ∈ Mn×p ,
mediante
1 1 −1 0
f ( x) = exp − tr[(X − µ)Σ (X − µ) ] (2.38)
(2π)np |Σ|n/2 2
Distribución de Wishart: Generaliza la distribución χ2 . Dado Y ∼ Nn,p (µ, Id, Σ), puede
probarse que la distribución de Y0 Y depende de µ a través de µ0 µ. Teniendo en cuenta eso y
dado E ⊂ Rn , se define la distribución de Wishart mediante Y 0 PE Y ∼ Wp (dim E, δ, Σ), con
δ = µ0 PE µ. Si δ = 0 se denota Wp (dim E, Σ). Las propiedades de la distribución de Wishart
son por completo análogas a la de la distribución χ2 y se estudian con detalle en Arnold (1981).
Ejercicio 41. Comprobar que W1 (m, δ, σ 2 ) = σ 2 χ2m (δ/σ 2 )
Distribución T 2 de Hotelling: Dados X ∼ Np (ν, Σ) y W ∼ Wp (m, Σ) independientes, se

define la distribución T 2 -Hotelling mediante
mX 0 W −1 X ∼ Tp,m
2
(δ), δ = ν 0 Σ−1 ν (2.39)
32
2
En el caso δ = 0 se denota Tp,m . En Arnold (1981) se prueba que esta distribución no es en
esencia nueva, sino que se identifica, salvo un factor escala, con un modelo tipo F , lo cual
garantiza que está bien definida. Concretamente
2 mp
Tp,m (δ) = Fp,m−p+1 (δ) (2.40)
m−p+1
2
En particular, se verifica que T1,m = t2m , por lo que debemos entender la distribución T 2 una
generalización en sentido estadístico de la distribución t2 . Es decir, que se utilizará en aquellos
problemas multivariantes cuyos análogos univariantes precisen de la distribución t-Student,
concretamente, en el contraste de hipótesis del tipo H0 : µ ∈ W con dim V |W = 1. Veremos
que en tales casos el estadístico de contraste puede entenderse geométricamente como una
2
distancia de Mahalanobis. Además, puede probarse que Tp,m converge en distribución a χ2p
conforme m tiende a infinito.
Distribuciones de Wilks, Lawley-Hotelling, Roy y Pillay: Pueden entenderse como

cuatro formas diferentes de generalizar la distribución F en el caso multivariante. Se estudian
con detalle en Arnold (1981). Al igual que ocurre con la distribución F , convergen en distribu-
ción a χ2p·dim V |W conforme aumenta el segundo grado de libertad, por lo cual omitiremos aquí
su estudio.
2.3.2. El modelo y ejemplos

Dada una matriz A ∈ Mn×p y E ⊂ Rn , se denota A ∈ E cuando cada columna de A
pertenece al subespacio E. Dicho esto, el modelo lineal normal multivariante viene dado por
una matriz de datos Y ∼ Nn,p (µ, Id, Σ), con Σ > 0 y la restricción µ ∈ V para algún V ⊂ Rn
conocido. Por lo tanto, Y constituye una matriz como la que aparece en el cuadro 1.1 que
recoge una muestra (en sentido amplio) de n observaciones Yi ∼ Np (µi , Σ) independientes. Si
consideramos una base X de V , el modelo puede parametrizarse también de la forma
Y = Xβ + E, E ∼ Nn,p (0, Id, Σ), β ∈ Mdim V ×p , Σ > 0 (2.41)
Los cuatro problemas univariantes (ejemplos 1-4) considerados en el apartado 2.2.2 se genera-
lizan al caso multivariante dando lugar a los siguientes problemas estadísticos multivariantes
que se estudiarán con más detalle en el siguiente capítulo. Basta tener en cuenta que la variable
respuesta Y se convierte en este caso en un vector respuesta p-dimensional de componentes
Y [1], . . . , Y [p].
Ejemplo 5. [Muestra aleatoria simple de una distribución p-normal] Consideremos Y1 , . . . , Yn
iid Np (ν, Σ). En ese caso, la matriz aleatoria Y = (Y1 , . . . , Yn )0 sigue un modelo de distribución
Nn,p (µ, Id, Σ) con µ ∈ V = h1n i y Σ > 0. Efectivamente, se verifica que cada columna µ[j] de
µ, que corresponde a la componente Y [j] del vector Y , pertenece a V .
Ejemplo 6. [Muestras independientes de p-normales con idénticas matrices de covarianzas]
Consideremos, para i = 1, 2, sendas muestras independientes Yi1 , . . . , Yini iid Np (µi , Σ). Si se
denota n = n1 + n2 e Y = (Y11 , . . . , Y2n2 )0 , se verifica que Y ∼ Nn,p (µ, Id, Σ) con Σ > 0 y
µ ∈ V = hv1 , v2 i.
Ejemplo 7. [Diseño completamente aleatorizado multivariante] Se generaliza el caso univarian-
te como en los ejemplos 5 y 6.
33
Ejemplo 8. [Regresión lineal multivariante] A diferencia del ejemplo 4 univariante, se pretende

explicar p variables respuesta, Y [1], . . . , Y [p], a partir de q variables explicativas, lo cual nos
lleva a un modelo tipo (2.41) donde Y es la matriz n × p de observaciones respuesta, expresada
como en (1.1), E la matriz n × p de errores, X es la misma matriz que aparece en (2.15) y β es
la matriz (q + 1) × p siguiente
   0 
β0 [1] . . . β0 [p] β0
 β1 [1] . . . β1 [p]   β 0 
  1 
β =  .. ..  =  ..  (2.42)

 . ... .   . 
βq [1] . . . βq [p] βq0
Para cada coeficiente βi [j], el subíndice i y el índice entre corchetes j indican, respectivamente,
a qué vector explicativo y a qué vector respuesta hace referencia. La primera fila, relativa al
término independiente, se denota por β0 , y el resto de la matriz por β.
Al igual que en el caso univariante, un problema como el del ejemplo 7 puede parametrizarse
de idéntica forma mediante variables dummys para convertirse en un problema de regresión
lineal multivariante, donde el contraste de la igualdad de las r medias equivale al contraste
total de la hipótesis β = 0.
Estos cuatro problemas se abordarán con más detalle en el siguiente capítulo. A continuación
estudiaremos brevemente la solución teórica a los problemas de estimación y contraste de
hipótesis.
2.3.3. Estimación y contraste de hipótesis

Los estimadores de µ y σ 2 en el modelo univariante pueden generalizarse de manera natural
mediante
µ̂ = PV Y, (2.43)
1
Σ̂ = Y0 P V ⊥ Y (2.44)
n − dim V
Puede probarse que, así definidos, µ̂ y Σ̂ son EIMV de µ y Σ y que, si reemplazamos en Σ̂

el denominador n − dim V por n, constituyen el EMV. El estimador µ consiste en estimar la
media de las distintas componentes por separado. Si el modelo está parametrizado de la forma
(2.41), el estimador de β será igualmente
β̂ = (X0 X)−1 X0 Y (2.45)
Ejercicio 42. Probar que (n − dim V )Σ̂ ∼ Wp (n − dim V, Σ)

En lo referente al contraste de hipótesis tipo H0 : µ ∈ W , las afirmaciones de los ejercicios
(28), (29) y (30) pueden extenderse trivialmente al caso multivariante. El test que resuelve el
contraste se denomina manova.
Manova: Consideraremos nuevamente la descomposición ortogonal Rn = W ⊕ V |W ⊕ V ⊥ ,

que da pie a definir las siguientes matrices p × p simétricas y semidefinidas positivas:
S 1 = Y 0 PW Y , S2 = Y0 PV |W Y, S 3 = Y0 P V ⊥ Y (2.46)
34
Ejercicio 43. Probar que,las tres matrices aleatorias siguen las siguientes distribuciones, in-
dependientes entre sí:
S1 ∼ Wp (dim W, µ0 PW µ, Σ), S2 ∼ Wp (dim V |W, µ0 PV |W µ, Σ), S3 ∼ Wp (n − dim V, Σ) (2.47)
Las matrices Si , i = 1, 2, 3 en (2.46), de dimensión p × p y semidefinidas positivas, pueden

expresarse mediante Z0i Zi , siendo Zi = G0i Y y Gi una matriz de n filas cuyas columnas constitu-
yan una base ortonormal de W , V |W y V ⊥ , respectivamente. En lo sucesivo prescindiremos de
S1 por argumentos de invarianza. Es más, la aplicación del Principio de Invarianza de manera
análoga al caso univariante nos conduce a desechar todo test cuyo estadístico de contraste no
pueda expresarse en función de los autovalores de S−1 3 S2 , que se denotan de mayor a menor
mediante t1 , . . . , tp . Si los consideramos como vector aleatorio su distribución exacta se calcula
partiendo de (2.47). En lo sucesivo, se denotará
b = mı́n{p, dim V |W } (2.48)
Ejercicio 44. Probar que p autovalores reales y no negativos de la matriz simétrica y semidefi-
0
nida postiva Z2 S−1 −1
3 Z2 los son a su vez de la matriz S3 S2 . Además, si p > b, los últimos p − b au-
tovalores de ambas matrices son necesariamente nulos. Lo mismo puede decirse de (θ1 , . . . , θp ),
definidos como los autovalores de la matriz Σ−1 µ0 PV |W µ, de los cuales (t1 , . . . , tb ) pueden con-
siderarse estimadores. Nótese que la hipótesis inicial H0 : µ ∈ W equivale a θ1 = . . . = θb = 0.
Así pues, el Principio de Invarianza nos conduce a considerar sólo los tests que se constru-
yan a partir de los b primeros autovalores de S−13 S2 , (t1 , . . . , tb ). En el capítulo 3 se verá un
interpretación precisa de estos autovalores. Sólo en el caso b = 1 estaremos en condiciones de
formular directamente un test basado en la distribución de t1 . Se da tal situación cuando p = 1
o dim V |W = 1:
(i) Si p = 1 las matrices de (2.46) son números positivos y t1 es, salvo una constante, el
estadístico F . Se trata pues del propio anova.
(ii) Si dim V |W = 1 puede probarse que t1 sigue, salvo una constante, una distribución T 2 -
Hotelling, lo cual permite formular un test UMP-invariante y de razón de verosimilitudes.
Si, además, p = 1, estaremos hablando del test de Student.
Dado que en el caso b > 1 el Principio de Invarianza no propicia una simplificación completa
de la información, el problema se ha abordado históricamente acogiéndose a otros diferentes
principios estadísticos que conducen a respectivas soluciones razonables, que pueden expresarse
a partir de los mencionados autovalores, es decir, que son invariantes. De esta manera aparecen
en la literatura estadística cuatro tests diferentes: el test de Wilks, que es el TRV y, por lo tanto,
responde al Principio de Máxima Verosimilitud; el test de Lawley-Hotelling, que responde al
Principio de Sustitución (o método plug-in); el test de Pillay, que sigue la línea del test anterior
pero utilizando los coeficientes de correlación canónica (ver apartado 3.4.1), y por último el
de Roy, que obedece al Principio de Unión-Intersección3 . La formulación precisa de estos tests
como funciones de los autovalores t1 , . . . , tb se comentará brevemente en el apartado 3.4.1. No
obstante, nos centraremos aquí en el test de Wilks por dos razones: por obedecer al Principio de
Máxima Verosimilitud que es, posiblemente, el más intuitivo en Estadística, y porque facilita
3
Los detalles de los principios mencionados los podemos encontrarlos en Arnold (1981) o en el Manual 59
citado en la bibliografía.
35
el algoritmo de selección de variables en regresión lineal, lo cual es especialmente interesante

en el análisis discriminante lineal.
De (2.38) se sigue que el estadístico de contraste del test de Wilks, es decir, la razón de
verosimilitudes, es la siguiente:
|S3 |
λ(Y) = (2.49)
|S2 + S3 |
Ejercicio 45. Probar que λ(Y) puede expresarse a través de t1 , . . . , tb mediante
b
Y
λ(Y) = (1 + ti )−1 (2.50)
i=1
Se demuestra en Arnold (1981) o, también, en el Manual 59, que bajo la hipótesis nula,
−(n− dimV ) log λ converge en distribución a χ2p·dim V |W cuando n tiende a infinito. Este resultado
es incluso cierto aunque no se respete el supuesto de normalidad, siempre y cuando el diseño de
la muestra respete ciertas condiciones razonables. Los otros tres tests recogidos por la literatura
(Lawley-Hotelling, Pillay y Roy) también convergen asintóticamente a la distribución χ2p·dim V |W .
En definitiva, para muestras de gran tamaño utilizaremos la distribución χ2 como referencia,
aunque el programa SPSS puede trabajar con otras aproximaciones a la distribución F . En el
cuadro 2.3 podemos apreciar un esquema explicativo.
Exacta p=1 p≥1 Asintótica p=1 p≥1

n↑
dim V |W = 1 t2n−dim v T 2p,n−dim v −→ dim V |W = 1 χ21 χ2p
dim V |W ≥ 1 F dim v|w,n−dim v Wilksp,dim v|w,n−dim v dim V |W ≥ 1 χ2dim v|w χ2p·dim v|w
Cuadro 2.3: Distribuciones en los contrastes de la media en función de las dimensiones del modelo
lineal V y de la hipótesis inicial W , del número p de variables y del tamaño n de la muestra.
También se recogen en Arnold (1981), Dillon y Goldstein (1984), Flury (1996) y Rencher
(1995), entre otras referencias, diversos tests para contrastes de hipótesis relativos a la matriz
de covarianzas implementados en los programas estadísticos, como el test M de Box, el de
esfericidad de Barlett y algunos otros, que no abordamos aquí por brevedad y dado que son
sensibles ante la violación del supuesto de normalidad.
36
Capítulo 3
Aplicaciones del modelo
En este capítulo desarrollaremos los cuatro problemas estadísticos formulados en los ejem-
plos 5-8 de la página 33 del capítulo anterior, cuyo denominador común es que se formalizan
mediante el modelo lineal multivariante. Añadimos además un apartado dedicado al análisis de
correlación canónica, relacionado directamente con el problema de regresión lineal multivarian-
te, y una sección dedicada al análisis de perfiles, relacionado con los tres problemas restantes.
Por último, ilustraremos con un ejemplo algunas de las técnicas estudiadas. En los distintos
casos se aplicarán los métodos teóricos de estimación y contraste de hipótesis expuestos en el
capítulo anterior. Se da por supuesto que el lector conoce ya las técnicas univariante análogas
(test de Student para muestras independientes y relacionadas, anova de una vía y estudio de
regresión lineal múltiple), que puede consultar, por ejemplo, en Peña (2010). A lo largo del
capítulo se hará uso del siguiente resultado, comúnmente conocido como teorema de los mul-
tiplicadores finitos de Langrange, que se deduce del teorema de la función implícita y permite
obtener valores extremos para una función definida en Rp bajo una serie de restricciones.
Lema 3.0.1. Sean k < p enteros y φ y f funciones derivables de Rp en R y Rk , respectivamente,
tales que existe máx{φ(x) : f (x) = 0} alcanzándose en c ∈ Rp tal que 5f (c) 6= 0. Entonces,
existe η ∈ Rk tal que 5(φ − η 0 f )(c) = 0.
3.1. Inferencia para una media

Desarrollamos aquí el ejemplo 5 de la página 33. Partimos pues de una muestra aleatoria
simple de una distribución p-normal, es decir,
Y1 , . . . , Yn iid Np (ν, Σ) (3.1)
de tal forma que la matriz de datos Y sigue un modelo de distribución Nn,p (µ, Id, Σ) con
µ ∈ V = h1n i y Σ > 0. Denótese por y el vector de medias (y[1], . . . , y[p])0 y por S la matriz
de covarianzas muestral. Podemos probar entonces los siguientes resultados.
Proposición 3.1.1. Los EIMV de µ y Σ son µ̂ = 1n · y0 y Σ̂ = n
n−1
S, respectivamente.
Proposición 3.1.2. n(y − ν)0 Σ̂−1 (y − ν) ∼ Tp,n−1
2
De la proposición 3.1.2 se sigue que el siguiente conjunto de Rp es una región de confianza

a nivel 1 − α para la media ν.
n o
p 0 −1 2,α
Eα (Y) = x ∈ R : n(y − x) Σ̂ (y − x) ≤ Tp,n−1 (3.2)
37
Capítulo 3 Inferencia para una media
Esta región geométrica es un elipsoide cuyo centro es y y cuya forma viene dada por Σ̂. Si
pretendemos contrastar la hipótesis inicial H0 : ν = 0, que equivale a µ ∈ W = 0, la proposición
2
3.1.2 invita confrontar con la distribución Tp,n−1 el estadístico de contraste
T 2 (Y) = n · y0 Σ̂−1 y (3.3)
Tanto el elipsoide (3.2) como el estadístico de contraste (3.3) pueden expresarse en términos
de la distancia de Mahalanobis DΣ̂2 definida en (1.38). Concretamente,
T 2 (Y) = n · DΣ̂2 (y, 0) (3.4)
Éste es precisamente el test UMP-invariante y de razón de verosimilitudes que se propone en el

capítulo anterior para este caso particular, donde se da la circunstancia de que dim V |W = 1.
Concretamente, con la notación introducida en dicho capítulo se verifica que t1 y θ1 , entendidos
como funciones sobre el espacio de observaciones y el de parámetros probabilísticos, son iguales
salvo escalares a DΣ̂2 (y, 0) y DΣ2 (µ, 0), respectivamente.
Ejercicio 46. Probar que Σ̂−1 > 0
Ejercicio 47. Probar que, en , p = 1, el test (3.3) es el de Student para una muestra.
Ejercicio 48. En el cuadro 3.1 se muestra el resultado de aplicar el test 3.3 con tres variables
y 36 individuos (se muestra el correspondiente cuadro de diálogo de SPSS). Interpretar la tabla
en función de los conocimientos teóricos.
Por otra parte, del Teorema Central el Límite y la Ley Débil de los Grandes Números se
sigue:
Proposición 3.1.3. Si Y1 , . . . , Yn iid con media ν y componentes en L2 , entonces se verifica la

siguiente convergencia en distribución:
lı́m nDΣ̂2 (y, ν) = χ2p (3.5)

n→∞
38
Capítulo 3 Inferencia para dos medias
Cuadro 3.1: Tabla Manova una muestra
b
Contrastes multivariados
Gl de la
Efecto Valor F hipótesis Gl del error Sig.
a
Intersección Traza de Pillai ,802 44,582 3,000 33,000 ,000
a
Lambda de Wilks ,198 44,582 3,000 33,000 ,000
a
Traza de Hotelling 4,053 44,582 3,000 33,000 ,000
a
Raíz mayor de Roy 4,053 44,582 3,000 33,000 ,000
a. Estadístico exacto
b. Diseño: Intersección
Cuadro 3.2: Cuadro de diálogos Manova para una muestra
Este resultado otorga validez asintótica al test propuesto aunque no se verifique el supuesto
de normalidad. Nótese también que podemos construir una región de confianza a nivel 1 − α sin
utilizar técnicas multivariantes, calculando para cada componente del vector respuesta Y un
intervalo de confianzas a nivel 1 − α∗ y componiendo entonces un rectángulo en dimensión p.
El valor de α∗ puede determinarse mediante de manera conservadora mediante la desigualdad
de Bonferroni:
\m Xm
P Aci

P Ai ≥ 1 − (3.6)
i=1 i=1
y[2] y
y[1]
Figura 3.1: Rectángulo y elipse de confianza
El elipsoide (3.2) delimita una región del espacio de menor tamaño que el del rectángulo,
siendo mayor su diferencia cuanto mayor sea la correlación entre las variables. Ello es debido a
que el método univariante no hace uso en ningún momento de las covarianzas y, por lo tanto,
emplea menos información que el multivariante. Si las componentes del vector aleatorio Y
fueran incorreladas (independientes bajo el supuesto de p-normalidad) el rectángulo anterior
podría construirse sin recurrir a la desigualdad de Bonferroni (3.6) y tendría un área similar al
de la elipse, cuyos ejes coincidirían con los ejes de coordenadas. En ese caso no procedería el
uso de métodos multivariantes. Página 1
39
3.2. Inferencia para dos medias

En esta sección desarrollamos el ejemplo 6 de la página 33. Se trata pues de estudiar la
posible relación entre un vector respuesta p-dimensional Y y un factor cualitativo que distingue
dos categorías. Partimos de dos muestras independientes de sendas distribuciones p-normales
con matriz de covarianzas común

Y11 , . . . , Y1n1 iid Np (µ1 , Σ)
(3.7)
Y21 , . . . , Y2n2 iid Np (µ2 , Σ)
La matriz de datos Y sigue un modelo de distribución Nn1 +n2 ,p (µ, Id, Σ) con µ ∈ V = hv1 , v2 i
y Σ > 0.
Ejercicio 49. Construir los EIMV de µ y Σ a partir de las medias aritméticas y1 , y2 ∈ Rp de
ambas muestras.
Ejercicio 50. Probar que
n1 n2
· DΣ̂2 (y1 , y2 ) ∼ Tp,n
2
1 +n2 −2
(θ), θ = DΣ2 (µ1 , µ2 ) (3.8)
n1 + n2
Si pretendemos contrastar la hipótesis inicial H0 : µ1 = µ2 , (3.8) invita a confrontar con la
2
distribución Tp,n 1 +n2
el estadístico de contraste
n1 n2
T 2 (Y) = · DΣ̂2 (y1 , y2 ) (3.9)
n1 + n2
En eso consiste precisamente el test UMP-invariante y de razón de verosimilitudes que se
propone en el capítulo anterior para este caso particular, donde se da la circunstancia también
de que dim V |W = 1. Concretamente, se verifica que t1 y θ1 , entendidos como funciones sobre
el espacio de observaciones y el de parámetros probabilísticos, son iguales salvo escalares a
DΣ̂2 (y1 , y2 ) y θ, respectivamente.
Como en la sección anterior, estamos también en condiciones de garantizar la validez asintó-
tica del test aunque no se verifique el supuesto de p-normalidad si n1 , n2 → ∞; también podemos
garantizarla aunque no se verifique el supuesto de homocedasticidad si, además, n1 /n2 → 1. Si
p = 1 el test propuesto es el conocido test de Student para dos muestras independientes.
La hipótesis H0 : µ1 = µ2 podría contrastarse prescindiendo de técnicas multivariantes
aplicando de manera independiente sendos tests de Student para cada una de las p componentes
del vector respuesta Y . En ese caso, los niveles de significación de cada test deberían calcularse
de manera conservadora mediante la desigualdad de Bonferroni. No podemos descartar que
el método multivariante (3.9) aprecie diferencias significativas entre ambas medias mientras
que ninguno de los tests de Student univariantes logre diferenciar las componentes de las
mismas. Hemos de ser conscientes, nuevamente, de que el método multivariante hace uso de la
información que aportan las covarianzas, lo cual no se tiene en cuenta en ninguno de los p tests
de Student.
Ejercicio 51. Interprétese en los términos anteriores la tabla que aparece en el cuadro 3.4,
obtenida según se indica en el cuadro 3.3, que corresponde a un estudio comparativo efectuado
a 25 lobos en los se que relacionan de 9 variables numéricas con el sexo.
40
Cuadro 3.3: Cuadro de diálogos Manova para dos muestras
b
Contrastes multivariados
Gl de la
Efecto Valor F hipótesis Gl del error Sig.
a
Intersección Traza de Pillai 1,000 6886,561 9,000 15,000 ,000
a
Lambda de Wilks ,000 6886,561 9,000 15,000 ,000
a
a
Raíz mayor de Roy 4131,937 6886,561 9,000 15,000 ,000
a
Sex Traza de Pillai ,784 6,038 9,000 15,000 ,001
a
Lambda de Wilks ,216 6,038 9,000 15,000 ,001
a
a
Raíz mayor de Roy 3,623 6,038 9,000 15,000 ,001
a. Estadístico exacto
b. Diseño: Intersección + Sex
Cuadro 3.4: Tabla Manova dos muestras
Nótese, por otra parte, que la j-ésima componente del vector respuesta, Y [j] es la proyec-
ción del vector Y sobre el j-ésimo eje de coordenadas. Si ej denota un vector unitario que lo
determina, podemos expresar Y [j] = e0j Y . En general, para cada eje hai con kak = 1, podemos
considerar la proyección a0 Y sobre hai que da lugar, teniendo en cuenta (1.15), a dos muestras
independientes
0
a Y11 , . . . , a0 Y1n1 iid N1 (a0 µ1 , a0 Σa)
(3.10)
a0 Y21 , . . . , a0 Y2n2 iid N1 (a0 µ2 , a0 Σa)
y a una hipótesis inicial H0a : a0 µ1 = a0 µ2 , que puede contrastarse a partir de los datos pro-
yectados mediante el test de Student. Concretamente, se confronta con la distribución tn1 +n2 −2
el estadístico de contrate thai (Y) definido como t(Ya). Conocido Y, debe existir necesariamente
un eje ha1 i que aporte un valor máximo para thai (Y). Teniendo en cuenta (1.15) y el lema 3.0.1
obtenemos la solución concreta1
(n1 − 1)S1 + (n2 − 1)S2

ha1 i = Sc−1 (y1 − y2 ), Sc = (3.11)
n1 + n2 − 2
1
Sc es el EIMV de Σ que se pedía calcular en el ejercicio 49.
41
Capítulo 3 Manova de una vía
Es más, si se denota
Wij [1] = a01 Yij , i = 1, 2, j = 1, . . . , ni (3.12)
se verifica entonces que t2 (W[1]) = T 2 (Y). En ese sentido podemos afirmar que distinguir las
dos muestras en dimensión p es equivalente a distinguirlas en dimensión 1 sobre el eje ha1 i,
que se denomina (primer) eje discriminante. El vector de proyecciones W[1] = Ya1 se denomina
vector de las (primeras) puntuaciones discriminantes. En la figura 3.2 el eje discriminante se
representa con líneas discontinuas:
Figura 3.2: Eje discriminante
3.3. Manova de una vía

En esta sección desarrollaremos el ejemplo 7 de la página 33 y ampliaremos el concepto de
eje discriminante. El problema supone una generalización del estudiado en la sección anterior,
puesto que trata la relación entre un vector respuesta p-dimensional Y y un factor cualitativo
que, en este caso, distingue entre r ≥ 2 categorías. Por lo tanto, partimos de un diseño,
denominado completamente aleatorizado, similar a (3.7) pero con r muestras independientes
de ni datos cada una. Mantendremos Pr aquí la notación habitual del diseño de experimentos.
En particular, n denotará la suma i=1 ni . La matriz de datos Y sigue entonces un modelo de
distribución Nn (µ, Id, Σ) con µ ∈ V = hv1 , . . . , vr i y Σ > 0. La hipótesis inicial a contrastar
en este caso es H0 : µ1 = . . . = µr , que se corresponde con µ ∈ W = h1n i. Si r > 2 y p > 1 se
verifica, a diferencia de los dos estudios anteriores, que b, según se define en (2.48), es mayor
que 1.
A pesar de que, desde un punto de vista práctico, la comparación de 2 medias es un pro-
blema semejante a la comparaciones de r ≥ 3 medias, el último estudio comporta una mayor
complicación formal dado que no puede resolverse en términos de una distancia T 2 entre un
único par de elementos. Por eso nos limitamos a aplicar la solución general del contraste ex-
puesta en el capítulo anterior a este caso concreto: se obtienen t1 > . . . > tb > 0, los autovalores
42
Capítulo 3 Manova de una vía
positivos2 de S−13 S2 , donde S2 y S3 se calculan según (2.46) y, a partir de los mismos, obtenemos
el valor del estadístico λ de Wilks definido según (2.50), cuya distribución exacta depende de
θ1 ≥ . . . ≥ θb ≥ 0, definidos en el ejercicio 44; por último, se confronta con la distribución
χ2p(r−1) el valor −(n − r) log λ(Y).
En el caso p = 1 el test obtenido es el anova de una vía; en el caso r = 2 es el test (3.9); en
general se denomina manova de una vía, que será asintóticamente válido aunque no se verifique
el supuesto de normalidad si n1 , . . . , nr tienden a infinito.
Desde este punto de vista, el problema de contrastar una hipótesis tipo H0 : µ ∈ W se
reduce a obtener las matrices S2 y S3 adecuadas. En este caso particular, pueden obtenerse
trivialmente de manera similar a SCE y SCH en (2.31).
   
SCH11 . . . SCH1p SCE11 . . . SCE1p
S2 =  .. .. S3 =  .. ..
, (3.13)
   
. . . . 
SCH1p . . . SCHpp SCE1p . . . SCEpp
donde, para h, k = 1, . . . , p,
r
X
SCHhk = ni yi· [h] − y·· [h] · yi· [k] − y·· [k] (3.14)
i=1
ni
r X
X
SCEhk = Yij [h] − yi· [h] · Yij [k] − yi· [k] (3.15)
i=1 j=1
Aunque no vamos a estudiar aquí diseños de experimentos multivariantes con dos o más
factores, el lector debe percatarse de que, si es capaz de resolver el problema en el caso univa-
riante, basta con proceder de manera análoga a (3.14) y (3.15) para obtener la solución general
para el caso multivariante.
El interés de estas dos últimas secciones radica en la vinculación existente entre el manova
de una vía (y también el test (3.9), entendido como caso particular) con el LDA (análisis
discriminate lineal) de Fisher. Por otra parte, el problema de comparación de medias en un
diseño completamente aleatorizado puede entenderse como un problema de regresión lineal,
multivariante en este caso, respecto a r − 1 variables dummys de asignación a categorías, lo
cual justifica a su vez el estudio del problema de regresión lineal multivariante que desarrollamos
en la siguiente sección.
3.3.1. Ejes discriminantes

El concepto de eje discriminante introducido en la sección anterior puede ampliarse cuando
el número de muestras es mayor que 2. Dado un eje hai podemos considerar el estadístico de
contraste Fhai (Y) para la hipótesis inicial de igualdad de medias a partir de los datos proyectados
sobre dicho eje. Nuestro primer objetivo es encontrar el eje ha1 i que lo maximiza. En el caso
r = 2 la solución es (3.11).
Ejercicio 53. Probar que la solución general es el eje ha1 i con
a1 = arg máx{a0 S2 a : a0 S3 a = 1} (3.16)
2
Nótese que descartamos la posibilidad de que dos autovalores sean iguales y de que alguno sea nulo. Se
debe a que puede probarse que, bajo los supuestos del modelo, ello puede ocurrir con probabilidad nula.
43
Capítulo 3 Regresión multivariante
Ejercicio 54. Utilizando el lema 3.0.1 y teniendo en cuenta (1.22), probar que Fha1 i (Y) =
n−r
r−1
· t1 , siendo t1 el primer autovalor de S−1 0
3 S2 y a1 un autovector asociado tal que a1 S3 a1 = 1.
De esta forma construimos el primer vector de puntuaciones discriminantes W[1] = Ya1 . El

proceso puede continuar en principio hasta completar p ejes discriminantes con sus respectivas
puntuaciones: el segundo eje discriminante ha2 i se define como aquél sobre el que debemos
proyectar Y para obtener un vector de puntuaciones W[2] = Ya2 incorrelado con W[1] y con
Fha2 i (Y) máximo, y así sucesivamente hasta obtener ap y el vector de puntuaciones W[p] = Yap .
Los ejes discriminantes son los p autovectores de S−1 3 S2 y los valores máximos del estadístico
F son, salvo el escalar (n − r)/(r − 1), sus respectivos autovalores t1 , . . . , tp . Dado que los
p − b últimos son necesariamente nulos, sólo se contemplan en la práctica los b primeros, de
ahí que en el caso r = 2 consideremos un único eje discriminante. En definitiva, si A denota la
matriz p×p cuyas columnas son los vectores a1 , . . . , ap , podemos transformar la matriz de datos
originales Y en una matiz de idénticas dimensiones con todas las puntuaciones discriminantes
W=YA (3.17)
donde A verifica
t1 0 0 0
 
.. ..

 . 
 .
 0 tb 0 0 
 
A0 S3 A = Id, A 0 S2 A =  (3.18)
 0 0 0 0 

 .. .. 
 . . 
0 0 0 0
El siguiente resultado puede demostrarse a partir de (3.18) y (3.14) y es la clave definitiva

para entender los ejes discriminantes y el significado de los autovalores t1 , . . . , tb :
Ejercicio 55. Para todo k = 1, . . . , p, se verifica:
r
X 2
ni Wi· [k] − W·· [k] = tk (3.19)
i=1
Por otra parte, tal y como se adelantaba en el ejericio 44, los autovalores ti pueden entenderse
respectivamente como estimadores de los autovalores probabilísticos θ1 , . . . , θp de la matriz
Σ−1 µ0 PV |W µ. La hipótesis inicial H0 (1) : θ1 = 0 equivale a H0 : µ1 = . . . = µr = 0, y se
contrasta mediante el manova de una vía a partir de t1 , . . . , tb , tomando como referencia la
distribución χ2p(r−1) . Sin embargo, la veracidad de la hipótesis inicial H0 (2) : θ2 = 0 equivale
en términos intuitivos a que toda la discriminación entre las medias recaiga exclusivamente
en el primer eje discriminante. La hipótesis H0 (2) puede contrastarse a partir de t2 , . . . , tp y
tomando como referencia la distribución χ2(p−1)(r−2) . De esta forma puede evaluarse la capacidad
de discriminación de sucesivos ejes, aunque en la práctica la valoraremos directamente en
términos muestrales ponderando los autovalores t1 , . . . , tb .
Ejercicio 56. Interpretar en los términos de la teoría los cuadros 3.5 y 3.6, correspondientes
a la comparación multivariante de medias entre las tres especies de flores de irisdata.
44
Autovalores
Correlación
Función Autovalor % de varianza % acumulado canónica
a
1 32,192 99,1 99,1 ,985
a
2 ,285 ,9 100,0 ,471
a. Se han empleado las 2 primeras funciones discriminantes
canónicas en el análisis.
Cuadro 3.5: Autovalores y correlaciones canónicas
Lambda de Wilks
Lambda de
Contraste de las funciones Wilks Chi-cuadrado gl Sig.
1 a la 2 ,023 546,115 8 ,000
2 ,778 36,530 3 ,000
Cuadro 3.6: Test de Wilks
3.4. Regresión multivariante

Desarrollamos aquí el ejemplo 8, lo cual da pie al análisis de correlación canónica. El pro-
blema se expresa formalmente así: Y = Xβ + E, donde E ∼ Nn,p (0, Id, Σ) con Σ > 0 y siendo β
una matriz de dimensiones (q + 1) × p del tipo (2.42). El problema de estimación de β queda
resuelto en (2.45). En lo referente al problema de contraste de hipótesis, consideraremos dos
casos de especial interés.
3.4.1. Contraste total: análisis de correlación canónica

Estudiamos primeramente el contraste de la hipótesis inicial H0 : β = 0 que, en términos
de la media µ = Xβ, se expresa mediante H0 : µ ∈ W = h1n i. Por lo tanto, dim V |W = q y
b = mı́n{p, q}. Se denotará por λ(Y)(Z) el estadístico λ de Wilks para el contraste total.
Ejercicio 57. Probar que, en este caso, se verifica
−1
S2 = nSyz Szz Szy (3.20)
−1
S3 = n[Syy − Syz Szz Szy ] (3.21)
El test de Wilks consiste en confrontar −[n − (q + 1)] log λ(Y)(Z) con con la distribución χ2pq ,
donde
b
Y
λ(Y)(Z) = (1 + ti )−1 , t1 > . . . > tb > 0 autovalores positivos de S−1
3 S2 (3.22)
Página 1
i=1
Ejercicio 58. En el caso p = 1, que se corresponde con el problema de regresión múltiple,

tenemos un único número
R2
t1 = (3.23)
1 − R2
45
Es decir que, si p = 1, el test total puede expresarse en función del coeficiente de correlación
múltiple (al cuadrado) definido en (1.21), según (3.23). En el caso multivariante p ≥ 1 podemos
generalizar la relación anterior si definimos r12 > . . . > rb2 > 0 como los autovalores positivos
−1 −1
de Syy Syz Szz Szy .
ti ri2
ri2 = , es decir, ti = , i = 1, . . . , b (3.24)
1 + ti 1 − ri2
Los autovalores r12 > . . . > rb2 > 0 se denominan coeficientes de correlación canónica
muestrales (al cuadrado) y, según hemos visto, contienen información relevante en el contraste
de la hipótesis H0 : β = 0. De hecho, podemos considerarlos como estimadores de los coeficientes
de correlación canónica probabilísticos ρ21 ≥ . . . ≥ ρ2b , que se definen según el teorema 3.4.1, de
manera, que para todo i, θi = ρ2i /(1 − ρ2i ). Es decir, que la hipótesis inicial puede expresarse en
términos de correlaciones canónicas mediante H0 : ρ21 = . . . = ρ2b = 0. No obstante, podemos
ofrecer una interpretación más clara de los coeficientes de correlación canónica.
En lenguaje probabilístico, si Y y Z son vectores aleatorios de dimensiones p y q, respec-
tivamente, buscamos α1 ∈ Rp y β1 ∈ Rq tales que las variables U1 = α10 Y y V1 = β10 Z tengan
varianza 1 y su correlación sea máxima entre todas las proyecciones de Y y Z sobre sendos
ejes de Rp y Rq . En ese caso, los ejes obtenidos, hα1 i y hβ1 i, se denominan primer par de ejes
canónicos, y (U1 , V1 ), el primer par de variables canónicas. La correlación(al cuadrado) entre
ambas se denota por ρ21 y se denomina primer coeficiente de correlación canónica. El siguiente
paso es determinar otro par de ejes y, por lo tanto, otro par de proyecciones (U2 , V2 ), incorrela-
das con (U1 , V1 ) y con una correlación entre sí, ρ22 , máxima, y así sucesivamente hasta llegar a
b = mı́n{p, q}. Consideremos las siguientes matrices cuadradas de orden p y q, respectivamente,
y rango b = mı́n{p, q}:
Σ−1 −1
yy Σyz Σzz Σzy (3.25)
Σ−1 −1
zz Σzy Σyy Σyz (3.26)
Ejercicio 60. Probar que los b primeros autovalores de las matrices (3.25) y (3.26) coinciden
(no así sus respectivos autovectores).
La demostración del siguiente resultado, que se recoge en el manual 59 de la UEx, se basa
fundamentalmente en el lema 3.0.1:
Teorema 3.4.1. Con las notaciones precedentes se verifica:
(i) Los coeficientes de correlación canónicas ρ21 . . . , ρ2b son los b primeros autovalores de la
matriz (3.25).
(ii) Los vectores α1 , . . . , αb que determinan los ejes canónicos asociados a Y pueden obtenerse
como autovectores de la matriz (3.25) asociados a ρ21 . . . , ρ2b , respectivamente. Análoga-
mente, los vectores β1 , . . . , βb que determinan los ejes canónicos para Z pueden obtenerse
como autovectores de la matriz (3.26) asociados a ρ21 . . . , ρ2b , respectivamente.
En definitiva, los ejes canónicos permiten entender de manera más natural la correlación
lineal entre las variables respuestas y las explicativas.
    ρ    
Z1 V1 ←→
1
U1 Y 1
 ..   .   ..   . 
 .  −→  ..   .  ←−  .. 
ρb
Zq Vb ←→ Ub Yp
46
Ejercicio 61. Expresar la definición y el teorema anteriores en términos muestrales.

Ejercicio 62. Probar que, dada una variable aleatoria real Y y un vector aleatorio Z de
dimensión q, la máxima correlación lineal simple entre Y y una combinación lineal de las
componentes de Z, β 0 Z, es el coeficiente (1.18), y se obtiene con β según (1.26).
Sabemos que la hipótesis inicial H0 : µ1 = . . . = µr en un diseño completamente aleatoriza-
do equivale a H0 : β = 0 si parametrizamos el modelo como una regresión lineal multivariante
respecto a r − 1 variables dummys. En ese caso, los autovalores t1 , . . . , tb correspondientes al
manova de una vía, que expresan la capacidad de discriminación de los ejes discriminantes,
pueden calcularse a partir de S2 y S3 definidas según (3.20) y (3.21), siendo Z el vector de
variables dummys. No obstante, dichos autovalores se relacionan con los coeficientes de co-
rrelación canónica según (3.24). Por lo tanto, el propio manova de una vía puede expresarse
en términos de los coeficientes de correlación canónica, calculados a partir de las variables
dummys, de la misma forma que el anova de una vía se expresa en términos del coeficiente de
correlación múltiple R2 . Además, ri expresa, al igual que ti , el poder de discriminación del eje
hai i, con la ventaja a la hora de interpretarlo de que está acotado entre 0 y 1.
Ejercicio 63. Probar que los ejes discriminantes son los propios ejes canónicos para Y que se
obtienen considerando como Z el vector de variables dummys (ver figura 4.2). ¿Cómo podrían
obtenerse mediante SPSS los ejes canónicos corerspondientes a las variables dummys?
Ejercicio 64. Interpretar en los términos de la teoría los coeficientes de correlación canónica
que aparecen en el cuadro 3.5.
Estamos ya en condiciones de interpretar los cuatro tests del manova propuestos en la
página 35: dado que los autovalores t1 , . . . , tb o, mejor, los coeficientes de correlación canónica
asociados, r12 , . . . , rb2 , resumen el grado de correlación entre el vector numérico respuesta Y
y el vector explicativo Z (que en la práctica suele corresponderse con las variables dummys
asociadas a un factor cualitativo) falta por determinar cómo se evalúan exactamente. Pues
bien, el test de Roy considera Pbúnicamente el primer autovalor t1 ; el test de Lawley-Hotelling
considera la suma de todos, i=1 ti ; el de Pillay hace lo mismo pero tomando las correlaciones
canónicas, es decir, bi=1 ri2 ; el test de Wilks, que es el que utilizamos aquí, los multiplica tras
P
sumarles el elemento neutro del producto, es decir, bi=1 (1 + ti )−1 (se invierte porque el TRV
Q
se ha expresado así tradicionalmente).
3.4.2. Contrastes parciales: método Lambda de Wilks

El otro tipo de contraste de interés está relacionado con la depuración del modelo mediante
los algoritmos de selección de variables. Se trata en esta ocasión de contrastar hipótesis iniciales
del tipo H0 : βj1 = . . . = βjk = 0 para k < q y j1 , . . . , jk ∈ {1, . . . , q}. La veracidad de esa
hipótesis conllevaría suprimir del modelo un parte de la matriz Z que se denota por ZD y está
compuesta por las columnas j1 , . . . , jk , dando lugar a un modelo reducido con un nueva matriz
ZR ∈ Mn×(q−k) .
Ejercicio 65. Probar que, si k = 1, el problema puede resolverse haciendo uso de la distribución
2
Tp,n−(q+1) que, salvo una constante, coincide con Fp,n−p−q .
En todo caso, se denota por λ(Y )(ZR |ZD ) el estadístico de Wilks que resuelve este contraste.
El método de Wilks oferece una ventaja a la hora de elaborar un algoritmo de selección de
variables, pues los estadísticos de contraste para los tests parciales pueden obtenerse a partir
de los correspondientes a los test totales para los diferentes modelos reducidos.
47
Capítulo 3 Análisis de perfiles

λ(Y)(Z)
λ(Y)(ZR |ZD ) = (3.27)
λ(Y)(ZR )
h1n ZR i
λ(Y )(ZR |ZD )

*

h1n ZR ZD i H λ(Y )(ZR )
HH
H
λ(Y )(Z) HH
j?
H
h1n i
Figura 3.3: Test de Wilks parcial
En el caso de la regresión lineal multivariante, pomos considerar, además de los conocidos

algoritmos de selección de variables explicativas (hacia adelante, hacia atrás, pasos sucesivos),
otros para la selección de variables respuesta: dado cualquier j = 1, . . . , p, entendemos que el
vector Y[j] es prescindible en el modelo cuando, si consideramos un modelo de regresión lineal
múltiple con Y[j] como variable respuesta y Z, YR como explicativas, Z debería ser eliminada
según el test parcial. Este criterio a la hora de seleccionar variables se relaciona con el concepto
probabilístico de independencia condicional.
Ejercicio 67. Probar que el contraste para Y[j] puede resolverse haciendo uso de la distribución
Fq,n−(p+q) .
Se denota no obstante mediante λ(YR |Y[j])(Z) el estadístico de Wilks que resuelve este
contraste .
Ejercicio 68. Teniendo en cuenta (3.27), probar que
λ(Y)(Z)
λ(YR |Y[j])(Z) = (3.28)
λ(YR )(Z)
Si estamos relacionando un vector numérico Y con un factor cualitativo que distingue r
categorías y parametrizamos el modelo mediante r−1 variables dummys recogidas en una matriz
Z, podemos aplicar una selección de variables respuesta para determinar qué componentes de
Y guardan una relación esencial con el factor. El método Lambda de Wilks se define como
el algoritmo de selección hacia adelante de variables respuestas según el test (3.28), y será de
utilidad en el capítulo 4.
Ejercicio 69. Probar que, en la fase j-ésima del algoritmo Lambda de Wilks, se introduce
la variable que, añadida a las j − 1 ya incluidas anteriormente, aporta una resultado más
significativo en el manova de una vía, es decir, un valor mínimo en el estadístico lambda
de Wilks, siempre y cuando resulte significativo el test parcial que se resuelve mediante la
distribución del ejercicio 67.
El cuadro 3.7 refleja en qué orden se van introduciendo las variables numéricas (todas)
de irisdata según el algoritmo Lambda de Wilks para la discriminación entre especies. En el
cuadro de diálogo de SPSS 3.8 se indica cómo ejecutar el análisis discriminate.
48
Variables no incluidas en el análisis

Lambda de
Paso F para entrar Wilks
0 sepleng 119,265 ,381
sepwidt 49,160 ,599
petleng 1180,161 ,059
petwidt 960,007 ,071
1 sepleng 34,323 ,040
sepwidt 43,035 ,037
petwidt 24,766 ,044
2 sepleng 12,268 ,032
petwidt 34,569 ,025
3 sepleng 4,721 ,023
Cuadro 3.7: Método Lambda de Wilks para irisdata
Cuadro 3.8: Cuadro de diálogos Lambda de Wilks
3.5. Análisis de perfiles

Se trata de una técnica que generaliza el test de Student para muestras relacionadas y da
sentido al contraste H0 : ν = 0 estudiado en la primera sección del capítulo. Este método puede
considerarse una alternativa más robusta al análisis de medidas repetidas (ver Arnold (1981)
y Hair et al. (1999)).
En ocasiones resulta interesante estudiar la evolución de una carácter numérico a lo largo
de una secuencia temporal con p mediciones. En ese caso, contaremos con un vector Y p-
dimensional, de manera que la hipótesis inicial H0 : µ[1] = . . . = µ[p] se interpreta como una
ausencia de evolución, al menos por término medio. También puede ser interesante comparar
las evoluciones en distintas categorías de un factor cualitativo, como en el ejemplo que se recoge
en la figura 3.4, que corresponde del dolor durante seis meses distinguiendo tres tratamientos3
En este caso, que los tres tratamientos tengan efectos idénticos por término medio equivale Pág
a la hipótesis inicial H0 : µ1 = µ2 = µ3 del diseño completamente aleatorizado, que se contrasta
3
J. Rodríguez Mansilla et al. Clinical Rehabilitation (2014).
49
Figura 3.4: Perfiles dolor por tratamientos
mediante el manova de una vía. No obstante, también puede resultar de interés contrastar, por
ejemplo, el paralelismo de los perfiles, que se interpreta como una evolución similar desde la
fase inicial. Si contamos con sólo p = 2 mediciones, una inicial y otra final, estaremos ante un
diseño conocido como de muestras relacionadas. Se resuelve calculando la diferencia D, con
media ν, entre las dos fases. De esta forma, la hipótesis inicial H0 : µ[1] = µ[2] equivale a ν = 0
y se contrasta mediante el test de Student para una muestra aplicado a D. La hipótesis inicial
de paralelismo entre los r perfiles equivale a ν1 = . . . = νr y se contrasta mediante el anova de
una vía.
Sin embargo, cuando consideramos más de 2 fases debemos calcular la diferencia entre cada
variable y la anterior, dando lugar a un vector D en dimensión p − 1. La hipótesis inicial
H0 : µ[1] = . . . = µ[p] se contrasta mediante el test (3.3) aplicado a D, y la de paralelismo
entre los r perfiles, mediante el manova de una vía.
Abordar un análisis de perfiles mediante un manova es sólo una de las posibles opciones
y, seguramente, no la más popular. Los supuestos es los que basa son la normalidad multiva-
riante y la igualdad de matrices de covarianzas (en el caso de incluir un factor intersujeto en
el modelo, como es el tratamiento en el estudio del dolor). Del primero sabemos que puede
obviarse asintóticamente, lo cual justifica la robustez del modelo. Como principal alternativa
podemos destacar4 el modelo de medidas repetidas que, en principio, supone además dos con-
diciones adicionales sobre la matriz o matrices de covarianzas: la igualdad de las varianzas de
las componentes, por un lado, y la igualdad de las covarianzas por otro. Un caso particular de
esta hipótesis es el supuesto de esfericidad (homocedasticidad y covarianzas nulas), que con-
duciría a aplicar un test F , pudiendo aplicarse correcciones en los grados de libertad tanto del
4
Rencher (1996), sección 6.9.
50
numerador como del denominador en dichos test en función del grado de desviación respecto
al modelo esférico. En eso consiste en la práctica el análisis de medidas repetidas. Si no existe
un factor intergrupo y no estamos dispuestos a asumir hipótesis relativas a la distribución del
vector (salvo la continuidad del mismo) contamos con la alternativa de Friedman basada en
rangos.
51
52
Capítulo 4
Problema de clasificación
En este capítulo vamos a abordar el problema de clasificación de una unidad experimental

respecto a r categorías posibles a partir de la medición de p variables numéricas. Por ejemplo,
mediante los datos recogidos en el archivo irisdata podemos elaborar una estrategia para deter-
minar a qué especie (setosa, virgínica o vesicolor) pertenece un lirio a partir de la observación
de sus cuatro medidas morfológicas (longitud y anchura de pétalo y sépalo). Desde el punto
de vista formal y al igual que sucede en el diseño completamente aleatorizado, contamos con
un vector numérico Y p-dimensional y un factor cualitativo I con r categorías que se identi-
fican con sendas distrubuciones en Rp . Debería existir una fuerte relación entre Y y el factor
cualitativo para que pudiéramos discriminar entre las categorías a partir de Y . En ese sentido,
el problema de clasificación puede entenderse como el problema inverso a la comparación de
medias que vimos en el capítulo anterior. La diferencia entre ambos problemas estriba en los
roles que desempeñan el vector numérico y el factor cualitativo en cada caso, como se ilustra
en la figura 4.1.
El problema de clasificación se enmarca en el contexto teórico de la Teoría de la Decisión,
pues buscamos una estrategia adecuada para decidir a qué distribución, pertenece una obser-
vación en Rp . Puede tratarse tanto a partir de un modelo frecuentista como Bayesiano. El
punto de vista Bayesiano se caracteriza por la posibilidad de contemplar preferencias a priori
hacia ciertas categorías. El modelo clásico o frecuentista puede entenderse como un modelo
condicional dado I del Bayesiano, dadas las distribuciones concretas a las que pertenecen las
observaciones numéricas. En un contexto Bayesiano la solución óptima al problema de clasifica-
ción es la denominada estrategia de Bayes. A pesar de las discrepancias a nivel filosófico entre
los marcos teóricos Bayesiano y frecuentista, la estrategia de Bayes podría ser también asumida
según el segundo punto de vista si consideramos una distribución a priori uniforme. Precisa-
mente, dicha estrategia coincide con la estrategia óptima según el criterio minimax, propia del
marco frecuentista. Bajo ciertas condiciones, la estrategia minimax consiste, como veremos, en
el Análisis Discriminate Lineal (LDA) de Fisher, que constituye el núcleo del capítulo. Veremos
también que la estrategia LDA puede entenderse como el reverso de una moneda cuyo anverso
es el manova de una vía, lo cual da pleno sentido al estudio del modelo lineal multivariante.
En general, los elementos básicos de la Teoría de la Decisión e Inferencia Bayesiana pueden
encontrarse en Nogales (1998). Para un desarrollo más detallado de este problema concreto
remitimos al lector a Anderson (1958).
53
Capítulo 4 Planteamiento general
Manova
Factor cualitativo Variables numéricas
Clasificación
Figura 4.1: Manova y clasificación
4.1. Planteamiento general

Nuestro propósito inicial es entender cómo se afronta este problema en el caso sencillo de
que existan sólo dos categorías para extenderlo después de manera natural al caso general de r
Problema Estrategia Regresión Correlación
categorías. No consideraremos aquí ninguna Manova
función de pérdida o costes para nuestro problema
Clasificación Fisher Multivariante Canónica
de decisión y supondremos inicialmente que las distribuciones de probabilidad del modelo son
continuas.
Si tenemos que decidir si una observación y ∈ Rp proviene de un modelo de distribución
P1 , con densidad p1 o, por el contrario, se explica por un modelo P2 , con densidad p2 , podemos
Reducción=maximización
distinguir entre dos tipos de estrategias: una estrategia aleatoria S consisten en calcular la
probabilidad de que la observación provenga de distribución P1 (I = 1) en función de la
observación Y ; una estrategia no Ejes aleatoria
PrincipalesS sólo probabilidad 1 ó 0, es decir, se decanta
directamente por P1 o P2 en función de Y . Los riesgos o probabilidades de error RS (1) y RS (2)
asociados a una estrategia no aleatoria se pueden calcular entonces a partir de p1 y p2 de forma
trivial. Obviamente, cualquier estrategia aleatoria S lleva asociada otra no aleatoria S basada
en un punto de corte en la probabilidad condicional, que por defecto es 0.5.
Podemos considerar un preorden en Varianza intragrupo
el conjunto de las estrategias no aleatorias, de ma-
nera que S1 S2 cuando RS1 (i) ≤ RS2 (i), para i = 1, 2. Se dice S1 S2 cuando alguna de
las desigualdades es estricta. Encontrar un elemento maximal es posible sólo en ciertos casos
triviales. Dada una distribución distribuciónReducción de la
a priori Q concreta, que se identifica con el nú-
mero q = Q({P1 }) ∈ [0, 1], definimos su riesgo de Bayes para q como el error cuadrático medio
dimensión
respecto a la función indicador para {I = 1}
Z Z
q 2
Varianza q · (S(y) − 1) · p1 (y) dy + (1 −Correlación
RS =intergrupo q) · S 2 (y) · p2 (y) dy (4.1)
En el caso particular de que la estrategia sea no aleatoria, se trata de la combinación convexa

de sus riesgos qRS (1) + (1 − q)RS (2). En todo caso, el riesgo se minimiza mediante la esperanza
condicional, es decir, la distribución a posteriori determinada por la Regla de Bayes , que se
denominaEjesestrategia
Discriminantes
de Bayes para q denotándose por Sq : Ejes Canónicos
Variables dummy
q · p1 (y)
Sq (y) = P (I = 1|Y = y) = (4.2)
q · p1 (y) + (1 − q) · p2 (y)
La estrategia no aleatoria asociada Sq consistirá en asignar y ∈ Rp a P1 cuando
p1 (y) 1−q
≥ (4.3)
p2 (y) q
Ejercicio 70. Probar que Sq minimiza el riesgo de Bayes para q en el subconjunto de estrategias
no aleatorias.
54
Capítulo 4 Planteamiento general
Puede probarse también que la clase estrategias no aleatorias {Sq : q ∈ [0, 1]} constituyen
una familia completa maximal, es decir: cualquier estrategia no aleatoria que no pertenece a
la famila es mejorada estrictamente por alguna de la familia, y no existe ninguna dentro de la
familia mejorada estrictamente por alguna otra. Así pues, deberíamos seleccionar una estrategia
de este tipo, pero la elección depende del valor de q que queramos considerar y ello nos situaría
en un marco Bayesiano. Si no estamos en condiciones de proponer una distribución a priori,
podemos optar por escoger la estrategia no aleatoria minimax, que es el elemento maximal
para el orden definido a partir del máximo de los riesgos. Puede probarse que es la estrategia
no aleatoria asociada a la estrategia Bayes para una distribución a priori uniforme, es decir,
S0.5 , y que RS0.5 (1) = RS0.5 (2). Ésta es la que adoptaremos por defecto si no estamos dispuestos
a asumir el marco teórico Bayesiano, teniendo en cuenta que, en todo caso, la introducción de
una distribución a priori en el modelo se traduciría simplemente en una corrección trivial de
la estrategia minimax según (4.3). En definitiva, la estrategia minimax consiste en asignar y a
P1 cuando se verifica
p1 (y) ≥ p2 (y) (4.4)
Es decir, se asigna la observación a la distribución que la hace más verosímil. Se trata pues
de una aplicación directa del Principio de Máxima Verosimilitud y ésta es la idea fundamental
que debe prevalecer. En el caso general de r categorías se procede de forma idéntica, asignando
y a Pi cuando
pi (y) ≥ pj (y), ∀j 6= i (4.5)
Método núcleo de estimación de densidades: Teniendo en cuenta (4.5), es obvio cómo

deberíamos resolver el problema si dispusiéramos de adecuadas estimaciones de las funciones
de densidad: asignando y a Pi cuando
p̂i (y) ≥ p̂j (y) ∀j 6= i (4.6)
Describiremos aquí heurísticamente el denominado método del núcleo de estimación de densi-

dades, empezando por el caso univariante. Para un estudio más detallado remitimos al lector
a Silverman (1986).
Supongamos que contamos con una muestra aleatoria y1 , . . . , yn , correspondiente a una
determinada distribución continua con función de densidad p y queremos estimar el valor de p
en y, que se denota p̂(y). Para ello escogemos un número δ > 0, que denominaremos ancho de
banda, y consideramos el intervalo [y − δ, y + δ], de amplitud 2δ. Si N (y) denota la cantidad de
datos de la muestra en el anterior intervalo y n es suficientemente grande, se sigue de la Ley
Débil de los Grandes Números que P [y − δ, y + δ] ' N (y)/n. Por otra parte, si δ es pequeño
se verifica por el Teorema Fundamental del Cálculo que P [y − δ, y + δ] ' p(y) · 2δ, lo cual
nos induce a definir para cada y ∈ Rp el estimador p̂(y) = N (y)/(2nδ). Si queremos expresar
p̂ en función de los datos de la muestra, hemos de tener en cuenta que un dato yi pertenece
al intervalo anterior si, y sólo si, δ −1 |yi − y| ≤ 1. Definimos entonces la función (denominada
núcleo) 1
si |u| ≤ 1
K(u) = 2 , u ∈ R. (4.7)
0 si |u| > 1
De esta forma,
n
1 X y − yi
p̂(y) = K , x∈R (4.8)
nδ i=1 δ
55
Capítulo 4 Análisis Discriminante Lineal o LDA
En el caso multivariante (dimensión p) no consideramos intervalos de amplitud 2δ centrados

en y sino cubos de volumen 2p δ p , y el núcleo K p asigna el valor 2−p a un punto u cuando
kuk∞ ≤ 1. De esta forma, la función de densidad se estima reemplazando en (4.8) K por
K p y δ por δ p . No obstante, la función de densidad estimada será de tipo escalonado. Un
procedimiento comúnmente utilizado para suavizarla es considerar, en vez del núcleo anterior,
el siguiente:
1 1 0 −1
K̃(u) = p/2
exp − u S u , u ∈ Rp , (4.9)
(2πS) 2
donde S es la matriz de covarianzas muestral. Así, la función de densidad se estima mediante
n
1 X 1 0 −1
p̂(y) = exp − 2 (y − yi ) S (y − yi ) (4.10)
nδ p (2πS)p/2 i=1
2δ
Podemos comprobar que la función anterior se trata, efectivamente, de una densidad. Una vez
estimadas las densidades de las distintas categorías procederemos a establecer las regiones de
clasificación según (4.6). En la literatura estadística encontramos núcleos diferentes a (4.9),
denominado gaussiano, como el triangular, el del coseno o de Epanechnikov, entre otros. Hay
que tener en cuenta que la estimación de las densidades, y por tanto la estrategia de clasifi-
cación, depende de la elección del núcleo K y del ancho de banda δ. Diversos trabajos vienen
a convencernos de que la elección del núcleo es poco determinante. Sin embargo, la elección
del ancho de banda sí lo es. No podemos hablar, desde luego, de un ancho de banda universal,
sino que debe depender del problema considerado. La selección de un ancho de banda excesiva-
mente grande tenderá a estimar la densidad demasiado plana, mientras que uno excesivamente
pequeño la estimará de manera excisivamente abrupta.
Otro inconveniente a tener en cuenta es la denominada “maldición de la dimensión”, que
consiste en que el número de datos requerido para lograr una estimación satisfactoria de la
densidad crece exponencialmente en relación con la dimensión considerada. Por lo tanto, cuan-
do tengamos un amplio número de variables precisaremos de una cantidad ingente de datos
para obtener una estimación fiable de la densidad. Eso explica el hecho de que sigamos hacien-
do hincapié aquí en el método tradicional para clasificar observaciones, denominado Análisis
Discriminante Lineal (LDA), debido a Fisher.
4.2. Análisis Discriminate Lineal

En la Estadística Paramétrica es muy frecuente partir del supuesto de normalidad a la
hora de formular un modelo cuyos estimadores y tests de hipótesis podrán tener distintos
comportamiento ante casos reales en los que este supuesto no se cumpla. En el contexto del
modelo lineal normal, al supuesto de normalidad se le añade el de igualdad de varianzas o de
matrices de covarianzas. Eso es precisamente, lo que haremos a continuación.
Sabemos que para determinar la estrategia aleatoria óptima en un marco Bayesiano se pre-
cisa conocer la distribución a priori y las diferentes funciones de densidad. Si no somos capaces
de aportar una estimación con garantías de las densidades, podemos optar por suponer que
las r distribuciones en juego siguen modelos p-normales con idéntica matriz de covarianzas, es
decir, Pi = Np (µi , Σ), para i = 1, . . . , r. En ese caso, la estrategia dependerá únicamente de
la distribución a priori y de la distancia DΣ2 definida en (1.38). En un problema de clasifica-
ción binaria con una distribución a priori dada por q ∈ [0, 1], la estrategia aleatoria óptima
56
consistiría, según (4.2), en la distribución a posteriori dada por
q · exp{− 12 · DΣ2 (y, µ1 )}

P (I = 1|Y = y) = , (4.11)
q · exp{− 21 · DΣ2 (y, µ1 )} + (1 − q) · exp{− 12 · DΣ2 (y, µ2 )}
de lo que se deduce que

P (I = 1|Y = y) = L − (β0 + y 0 β)

(4.12)
donde
1 − q 1 0 −1
+ µ1 Σ µ1 − µ00 Σ−1 µ0 ,

β0 = log (4.13)
q 2
−1
β = Σ (µ0 − µ1 ) (4.14)
y L denota la denominada función logística, representada en la figura 4.6, que se define mediante
1
L(x) = , x∈R (4.15)
1 + ex
De esta forma, al condicionar sobre los valores numéricos de las observaciones, Yi1 , . . . , Yini , i =
1, 2, en el modelo original Bayesiano, obtenemos el conocido como modelo de regresión logística,
que es del tipo lineal generalizado. Por contra, si condicionamos respecto a las categorías
de las observaciones tendremos dos muestras aleatorias simples e independientes de sendas
distribuciones Np (µi , Σ), lo cual supone un modelo lineal normal multivariante.
Tanto (4.11) como (4.12) expresan la estrategia en función de parámetros poblacionales
que en la práctica serán desconocidos. En la regresión logística, los parámetros se estimarán
siguiendo la técnica propia de los modelo lineales generalizados. Sin embargo, el método LDA se
deriva de la estrategia no aleatoria minimax, es decir, la asociada a la estrategia aleatoria óptima
(4.11) para q = 0.51 , o equivalentente (4.4), pero sustituyendo los parámetros probabilísticos µ1 ,
µ2 y Σ por estimadores óptimos de los mismos, según el modelo lineal normal multivariante. Es
decir, se considera como estimador de µi a la media muestral Yi· de la categoría i-ésima; como
estimador de Σ tomaremos (n1 + n2 )−1 (S1 + S2 ), siendo Si la matriz de covarianzas muestral de
la i-ésima categoría. En ese caso, el método asigna Y a la media más cercana según la distancia
DΣ̂2 . En general, si se trata de un problema de clasificación respecto a r categorías, la estrategia
minimax consiste en asignar Y a la categoría i-ésima cuando se verifica
(Y − yi· )0 Σ̂−1 (Y − yi· ) ≤ (Y − yj· )0 Σ̂−1 (Y − yj· ), ∀j 6= i (4.16)
Si contemplamos una distribución a priori diferente de la uniforme podemos estimarla mediante

las proporciones de la muestra y corregir las desigualdades de (4.16) mediante sumandos que
dependen de las proporciones por categorías. Es lo que podría denominarse estrategia LDA
Bayesiana.
Cada desigualdad en (4.16) da lugar a la división del Rp en dos semiespacios cuya frontera
es una subvariedad afín (p − 1)-dimensional, de ahí que esta estrategia se denomine lineal (de
Fisher) para diferenciarse de la cuadrática (de Fisher también), que veremos más adelante, en
la cual Rp estará fragmentado por cuádricas.
Como ejemplo, utilizaremos el archivo irisdata de Fisher para intentar clasificar una flor
entre las tres especies consideradas en función de sus cuatro medidas morfológicas. El programa
1
Que es la que mejor se asimila a un marco Frecuentista.
57
SPSS diseña la estrategia LDA a partir de los datos ya asignados a categorías y es capaz de
clasificar en función de la misma cualquier otro dato que aparezca desagrupado. También recla-
sifica según la estrategia los propios datos agrupados, como el caso que vemos a continuación.
La reclasificación aporta una estimación de los riesgos de la estrategia, que son, según el cuadro
4.2, del 0 % para setosa, del 2 % para virginica y 4 % para vesicolor.
Discriminante
Grupo mayor
Número Grupo P(G=g | Distancia de

de caso Grupo real pronosticado p D=d) Mahalanobis
1 3 3 ,503 1,000 1,376
Cuadro 4.1: Reclasificación según LDA
Resultados de la clasificacióna
Grupo de pertenencia pronosticado

species setosa vesicolor virginica Total
Original Recuento setosa 50 0 0 50
vesicolor 0 48 2 50
virginica 0 1 49 50
% setosa 100,0 ,0 ,0 100,0
vesicolor ,0 96,0 4,0 100,0
virginica ,0 2,0 98,0 100,0
a. Clasificados correctamente el 98,0% de los casos agrupados originales.
Cuadro 4.2: Estimaciones de los riesgos LDA
Según se indica en el cuadro 4.1, el dato uno se ha clasificado en el grupo 3 porque

P la media de
éste minimiza la distancia de Mahalanobis. Tanto es así que el cociente p̂3 (y)/ i p̂i (y) ' 1, con
las densidades estimadas sustituyendo las medias y matrices de covarianzas por sus estimadores.
Por lo tanto, podemos considerarla una clasificación clara. De hecho, sabemos que es correcta.
La estrategia LDA de Fisher, definida en (4.16), posee buenas propiedades asintóticas.
Concretamente, puede probarse que, en el caso binario (r = 2) los riesgos de la misma convergen
asintóticamente al valor Z ∞
f (x) dx (4.17)
1
2
θ
siendo f la densidad de la distribución N (0, 1), y θ = DΣ2 (µ1 , µ2 ). Se trata del parámetro θ que
aparece en (3.8) en relación con el contraste de la hipótesis inicial H0 : µ1 = µ2 , que se identifica
con θ = 0. Por lo tanto, si µ1 = µ2 , la estrategia de Fisher se comportaría asintóticamente como
un sorteo a cara o cruz. Sin embargo, a medida que las medias se alejan según la métrica de
Mahalanobis, los riesgos asintóticos tienden a 0. En la práctica, que las distribuciones estén bien
diferenciadas suele ser mucho más importante que el cumplimiento de los supuestos del modelo
de cara a lograr una estrategia con riesgos bajos, que es lo que a la postre nos interesa. Eso es
lo que ocurre con irisdata: no estamos en condiciones de asumir la normalidad ni la igualdad
de matrices de covarianzas, pero las tres especies consideradas se diferencian claramente según
Página 1
58
sus medidas morfológicas, de ahí el éxito de la estrategia de Fisher, que queda patente en el
cuadro 4.2.
En definitiva, el manova de una vía Manova
y la estrategia de clasificación lineal de Fisher son
métodos que parten del mismo modelo, aunque en el primer caso es el factor el que desempeña
Factor cualitativo Variables
la función explicativa, mientras que en el segundo es el vector numérico. Un resultado poco
significativo a la hora de comparar las medias no ofrece expectativas de éxito en la clasificación,
Clasificación
justo al contrario que un resultado significativo. Por eso decimos que el manova y la clasificación
son el anverso y el reverso de una misma moneda. De hecho, es el problema de clasificación el
que da pleno sentido al estudio del manova y, dado que este último puede entenderse como una
regresión multivariante respecto a las variables dummys, da sentido al estudio de los coeficientes
de correlación canónica, pues el contraste de igualdad de medias puede expresarse en términos
de los mismos según (3.24).
Problema Estrategia Regresión Correlación

Manova
Clasificación Fisher Multivariante Canónica
Figura 4.2: Esquema general
Una vez hemos entendido el problema de clasificación como un problema de relación entre
un vector aleatorio p-dimensional y un factor con r categorías, cobra especial interés el método
de selección de variables Lambda de Wilks, estudiado en el capítulo 3, pues permite desechar
aquellas componentes del vector que no aportan información particular en el problema de
clasificación.
4.2.1. LDA y ejes discriminantes

Ejercicio 71. Como caso particular probar que, si r = 2, la estrategia (4.16) consiste en
asignar Y a y1 cuando 0
1
Y − (y1 + y2 ) Sc−1 (y1 − y2 ) > 0 (4.18)
2
lo cual sugiere una partición del espacio en función del eje discriminante definido en (3.11) y
que se representa en la figura 3.2.
Este resultado establece una primera conexión entre los ejes discriminantes y el problema
de clasificación, que desarrollaremos en el caso general r ≥ 2. Las observaciones originales com-
ponen en principio un matriz Y ∈ Mn×p , pero pueden expresarse también a través de la matriz
W, definida en (3.17) mediante W = YA, cuyas columnas recogen las diferentes puntuaciones
obtenidas al proyectar sobre los ejes discriminantes, determinados por las correspondientes co-
lumnas de A. Proyectar sobre los ejes discriminantes puede entenderse como un cambio de
coordenadas. Podemos obrar de igual forma con cualquier vector aleatorio Y proyectándolo
sobre los ejes discriminantes para obtener W = A0 Y . Dado que la proyección es una aplicación
lineal, la media aritmética de los datos proyectados coincide con la proyección de la media
aritmética de los datos originales. De esta forma, la estrategia (4.16), expresada en términos
de las puntuaciones discriminantes, asigna W a Wi· cuando, para j = i se alcanza el siguiente
mínimo: 0
mı́n (A0 )−1 (W − Wj· ) Σ̂−1 (A0 )−1( W − Wj· )

(4.19)
1≤j≤r
59
Dado que Σ̂ = (n − r)−1 S3 y en virtud (3.18), la estrategia de Fisher asigna W a la categoría

de Wi· cuando en j = i se alcanza mı́nj kW − Wj· k2 . Es decir, cuando expresamos los datos en
términos de las puntuaciones discriminantes se trata de minimizar la distancia Euclídea. Por
lo tanto, se trata de buscar el siguiente mínimo
p
X 2
mı́n W [k] − Wj· [k] (4.20)
1≤j≤r
k=1
Se sigue entonces de (3.19) que, si k > b, dado que tk = 0, Wj· [k] = W·· [k] para todo j, es
decir, que el sumando k-ésimo en (4.20) será constante en j y, por lo tanto, no interviene en la
búsqueda del mínimo en j. Por eso podemos ignorar las puntuaciones discriminantes asociadas a
autovalores nulos, de manera que el problema queda reducido a minimizar distancias Euclídeas
en dimensión b. Si el valor de b es bajo podemos pues visualizar el problema de clasificación
mediante un gráfico b-dimensional. Por ejemplo, en el caso de irisdata, tenemos el diagrama de
dispersión de la figura 4.3.
Para valores altos de b podemos visualizar igualmente el problema desechando las puntua-
ciones discriminantes asociadas a autovalores pequeños pues, razonando de manera análoga,
los correspondientes sumandos en (3.19) serán casi constantes y, por lo tanto, tendrán escasa
influencia en el problema de minimización. Por ejemplo, en la figura 4.4 representamos las tres
primera puntuaciones discriminantes en un problema de clasificación respecto a 7 categorías,
desechando pues la información aportada por los tres últimos ejes discriminantes. Para deter-
minar si el eje discriminante k-ésimo puede ser despreciado podríamos en principio resolver un
contraste de hipótesis del tipo H0 (k) : θk = 0, según se ve en el apartado 3.3.1. No obstante,
este método requiere de supuestos teóricos relativos a la distribución de los datos y, además,
es muy conservador. Lo habitual es ponderar desde un punto de vista puramente muestral los
autovalores t1 , . . . , tb que, a su vez, se asocian según (3.24) a los coeficientes de correlación
canónica r1 , . . . , rb .
En el caso de irisdata (figura 4.3), podemos apreciar que el peso de la discriminación recae
casi exclusivamente en la primera puntuación discriminante, según sabíamos ya por el cuadro
3.5. En la figura 4.4 (izquierda) se aprecia cierta confusión entre algunas de las variedades de
aceituna a partir de 17 variables numéricas medidas2 al representar las dos primeras puntua-
ciones discriminantes. Sin embargo, la confusión se resuelve en parte al introducir la tercera
puntuación, como se aprecia en la figura de la derecha.
4.2.2. Estrategia cuadrática de Fisher

Se trata de una generalización inmediata de la estrategia lineal. Se asume igualmente la
hipótesis de p-normalidad pero no se asume la igualdad de las r matrices de covarianzas. En
consecuencia, la estrategia consiste en modificar (4.16) reemplazando la estimación conjunta
de la matriz Σ por las diferentes estimaciones Si de cada una de las matrices Σi . Así pues, la
estrategia consiste en asignar Y a la media yi· cuando, para j 6= i, se verifica
(Y − yi· )0 Si−1 (Y − yi· ) ≤ (Y − yj· )0 Sj−1 (Y − yj· ) (4.21)

2
Datos recogidos en el INTAEX de Badajoz.
60
Capítulo 4 Métodos alternativos
Figura 4.3: LDA en irisdata
Al contrario de lo que sucedía en la estrategia lineal, los términos cuadráticos no se anulan en la

inecuación, de ahí el nombre. Para una mejor comprensión de la estrategia, podemos proyectar
los datos sobre los primeros ejes discriminantes, aunque el gráfico no podrá interpretarse en los
términos anteriores.
En el caso de la clasificación de aceitunas según sus variedades, la estrategia cuadrática
de Fisher disminuye ligeramente los riesgos de la lineal, al menos según es estima mediante la
reclasificación. No obstante, precisamos de métodos alternativos de clasificación que presenten
diferencias más radicales respecto al LDA. Estudiaremos muy brevemente algunos de ellos en
la siguiente sección. Las diferentes variantes del método de Fisher se ejecutan en SPSS a través
del cuadro de diálogos 3.8.
4.3. Métodos alternativos

En esta sección expondremos esquemáticamente los tres métodos de clasificación alternati-
vos al LDA más populares, posiblemente, al margen de aquéllos relacionados con la estimación
de densidades.
4.3.1. Regresión logística

Hemos visto que el método LDA puede considerarse el reverso del manova, el test principal
del modelo lineal. Sin embargo, la regresión logística binaria corresponde a un modelo lineal
61
nvar nvar
6,00000
CARRASQUEÑA CARRASQUEÑA
CACEREÑA CACEREÑA
5,00000
CORNICHE CORNICHE
CORNEZUELO CORNEZUELO
MORISCA
Segunda puntuación discriminante
MORISCA
Tercera puntuación discriminante

PICUAL 4,00000
PICUAL
VERDIAL BADAJOZ VERDIAL BADAJOZ
2,50000
2,00000
0,00000
0,00000
-2,50000
-2,00000
-5,00000
-4,00000
-5,00000 -2,50000 0,00000 2,50000 5,00000 -5,00000 -2,50000 0,00000 2,50000 5,00000
Primera puntuación discriminante Primera puntuación discriminante
Figura 4.4: Variedades de aceituna: puntuaciones discriminantes 1, 2 y 3
generalizado. En el primero el parámetro principal se estima mediante (2.45)como solución

a un sistema de ecuaciones lineales y, en el segundo, como a aproximación a un sistema de
ecuaciones no lineales. Vistos con perspectivas, no deberíamos hablar de un modelo LDA y un
modelo de reresión logística binaria, sino entenderlos como métodos diferentes para estimar la
probabilidad a posteriori (4.11) partiendo de un mismo modelo origen. En todo caso, no es de
extrañar pues que la regresión logística binaria ofrezca resultados muy similares a la estrategia
LDA Bayesiana, tal y como se aprecia en la figura 4.5. De ahí que la regresión logística no sea,
en propiedad, una alternativa al método de Fisher.
R2 Lineal R2 Lineal
= 0,986 = 0,993
Página 1
1,0 1,0
Página 1
Probabilidad diabetes Reg Log
,8 ,8
Probabilidad CHD Reg Log
,6 ,6
y=7,12E-3+0,95*x y=1,25E-3+0,99*x
,4 ,4
,2 ,2
,0 ,0
,0 ,2 ,4 ,6 ,8 1,0 ,0 ,2 ,4 ,6 ,8 1,0
Probabilidad diabetes LDA Probabilidad CHD LDA
Figura 4.5: Comparativa entre las probabilides de dabetes tipo II (izquierda) y enfermedad coronaria
(derecha) predichas mediante LDA Bayesiano y regresión logística, a partir de los datos de Diabetes
Data Study (Schoderling) y South African Heart Disease Study (Rousseauw), respectivamen-
te.
62
Página 1 Página 1
Si el factor cualitativo distingue r > 2 categorías podemos aplicar el método de regresión logís-
tica multinomial. A grandes rasgos, consiste en una composición de r − 1 regresiones logísticas
tomando una categoría como referencia. Cada una de estas regresiones permite estimar la pro-
babilidad de que un dato concreto pertenezca a una categoría dada, dividida por la probabilidad
de que pertenezca a la categoría de referencia. Si los r − 1 cocientes resultan ser inferiores a 1,
el dato se asigna a la categoría de referencia; en caso contrario, se asigna a la que aporte un
cociente máximo.
Figura 4.6: Función logística
A pesar de lo dicho anteriomente, en los programas estadísticos habituales se observan

claras diferencias en las salidas ofrecidas según el método empleado (LDA o regresión logística
binaria). La regresión logística aporta directamente estimaciones de Odds Ratios, un concepto
de gran interés en Epidemiología. Además, el método viene tradicionalemente acompañado
del test de Hosmer-Lemeshov, que contrasta la veracidad de la estimación de la distribución
condicional de la respuesta Y dado el vector explicativo. Dicho test consiste en dividir la muestra
en 10 subconjuntos en función de los deciles de la probabilidad estimada de Y = 1, y contrastar
mediante un estadístico χ2 con 9 grados de libertad si las proporciones reales de datos con Y = 1
en dichos subconjuntos se ajustan a lo que marcan las probabilidades medias de los mismos. Un
resultado significativo en este test se traduce en una ineptitud del modelo de regresión logística
para estimar la distribución condicional de la variable categórica y, en particular, para llevar a
cabo la clasificación, lo cual nos conduciría a ensayar con estrategias alternativas de diferente
naturaleza, como las que veremos a continuación. Realmente, nada impediría implementar
este test en el método de Fisher, entre otros. Por otra parte, la abundancia de resultados no
significativos en el test de Hosmer-Lemeshow, incluso con grandes muestras, revela la relativa
intrascendencia en la elección del modelo estadístico si la compramos con el grado de correlación
entre las variables observadas y la respuesta Y .
4.3.2. Vecino más próximo

Dado un valor k = 1, 2 . . ., el método del vecino más próximo para k (K-NN) es un proce-
dimiento de clasificación no paramétrico pues no impone supuesto alguno sobre la distribución
de los datos, salvo que las variables medidas deben ser numéricas. El método se estudia con
detalle en Hastie et al. (2008). Resumidamente, se trata de asignar una observación Y ∈ Rp a
la categoría que tenga mayor presencia en el entorno de Y constituido por las k observaciones
63
de la muestra más próximas. La cercanía se evalúa en principio en función de la métrica Eu-

clídea en Rp , aunque pueden considerarse alternativas. La estrategia depende en gran medida
del valor de k seleccionado, de ahí que, como mínimo, se precise tantear con diferentes valores
y seleccionar aquél cuya estrategia ofrezca menores riesgos estimados.
Grupo
A
B
14,00
12,00
x1
10,00
8,00
5,00 6,00 7,00 8,00 9,00 10,00
x2
Figura 4.7: K-NN 94 %
En la figura 4.7 persentamos el diagrama de dispersión relativo a un vector con p = 2

variables en el cual distinguimos 2 categorías. Si pretendemos determinar una estrategia de
asignación a categorías, los riesgos estimados por reclasificación de la muestra son del 43 %
para LDA, 24 % para la alternativa cuadrática, 45 % para la regresión logística binaria y 6 %
para K-NN con k = 3. En casos como éste se precisa pues una alternativa radicalmente diferente
a LDA.
Ejercicio 72. ¿Cómo se explica que la alternativa cuadrática de Fisher mejore sustancialmente
la estrategia lineal en el ejemplo de la figura 4.7?
Ejercicio 73. ¿Influye la escala en que se miden las componentes de Y en la clasificación
según el método LDA? ¿Influye en la regresión logística? ¿Influye en el método K-NN? Caso
de influir, ¿cómo podríamos compensar ese hecho?
Página 1
4.3.3. Árbol de decisión

Se trata de otro método no paramétrico (aunque lo más apropiado sería calificarlo, al igual
que K-NN, como un procedimiento típico de Minería de Datos) que se puede resumir como
una categorización óptima en la toma de decisión. Se fundamenta, como veremos a continua-
ción, en una idea bastante sencilla. Sin embargo, ante la ausencia de supuestos estadísticos, la
64
Figura 4.8: Vécino más próximio y árbol de decisión para irisdata
optimización se realiza mediante una gran cantidad de ensayos, por lo que no se concibe sin el
uso de un ordenador. De ahí que métodos como éste o como el del vecino más próximo sean
relativamente recientes. Se estudia con detalle en Hastie et al. (2008).
Predicción categórica de una variable numérica: El algoritmo que explicaremos a con-

tinuación fue descrito en Breiman et al. (1984), y se denomina CART o CRT3 .. Se diseña
originalmente para una variable respuesta Y de tipo numérico y un vector p-dimensional X
explicativo y se fundamenta intuitivamente en la propia definición de integral de Riemann
como límite de sumas finitas de áreas de rectángulos. Más generalmente, tenemos en cuenta
que cualquier función en L2 puede aproximarse tanto como se quiera, según la métrica (1.4),
mediante una suma de funciones constantes sobre intervalos, si el espacio origen es R o, en
general, sobre rectángulos medibles, si es Rp .
Dividir el espacio origen Rp en rectángulos medibles se traduce en la práctica en una ca-
tegorización conjunta de las p variables explicativas numéricas, ideada para ajustarse lo mejor
posible (en sentido mínimo cuadrático) a la variable numérica respuesta Y . El método CRT es
de tipo binario pues el vector X será categorizado mediante dicotomizaciones o biparticiones
sucesivas4 . En la figura 4.9 podemos apreciar cómo funciona una dicotomización de este tipo
para el caso p = 1. Ésta se entiende a su vez como una ramificación de un nodo único (nodo
parental) en dos ramas o categorías (nodos filiales), a las cuales se les asigna un valor constante
como aproximación a la variable respuesta Y .
Sabemos (ejercicio 3) que, si queremos aproximar en L2 la variable Y mediante una constante,

ésta debe ser su media respecto a la probabilidad P (en términos muestrales, su media aritmé-
tica), siendo su varianza la medida del error cuadrático. Podemos conseguir un mejor ajuste si
consideramos dos rectángulos, es decir, un punto de corte X1 . En tal caso consideraríamos como
constantes a cada lado del mismo las respectivas medias. La mejoría en términos de L2 sería
la diferencia entre el error cuadrático, es decir, la varianza total por la medida del segmento
completo, y la suma de las nuevas varianzas multiplicadas por las medidas de los respectivos
3
Classification and Regressión Trees
4
Cosa que no ocurre con otros métodos como CHAID.
65
X≤X1 X1 X>X1
Figura 4.9: Dicotomización de una variable explicativa
subsegmentos delimitados por X1 . El punto de corte óptimo X1 será aquél que maximice esa
ganancia. Desde un punto de vista computacional, el problema se reduce a ensayar diferentes
puntos de cortes a lo largo del eje X en busca de la mínima media ponderara entre las va-
rianzas resultantes. Si tenemos p variables debemos ensayar a lo largo de cada uno de los ejes
explicativos.
Una vez encontrado el primer punto de corte óptimo, debemos proceder de igual forma para
obtener un segundo punto de corte X2 en alguna de las variables explicativas, lo cual supondrá
un nuevo nivel de ramificación en algún nodo ya existente, y así sucesivamente. Las medias de
los nodos terminales resultantes nos darán la aproximación categórica a Y . En la figura 4.10
se representa esquemáticamente cuatro iteraciones de ramificación en un problema con p = 2
variables (ejes) explicativas. Ello da lugar a cinco nodos terminales o categorías diferentes. La
mayor proximidad entre colores simboliza un separación más tardía entre los respectivos nodos.
Dejando a un lado la idea que inspira el método, que es integral de Riemann, entenderemos
en lo sucesivo el procedimiento como un algoritmo iterativo de dicotomización, que busca en
cada paso la mínima media ponderada de las varianzas ponderada resultantes al reemplazar
el valor medio del nodo parental por los valores medios de los nodos filiales. Este algoritmo
puede aplicarse pues sin problema en presencia de variables explicativas de tipo cualitativo.
Por lo tanto, no implica la asunción de supuesto alguno respecto al conjunto de p variables
explicativas.
Sobreajuste: Dado que, en la práctica, no contamos con una variable Y ∈ L2 sino con una
muestra Y de tamaño n de la misma, bastaría considerar un árbol en el que cada elemento de la
muestra constituyera un nodo terminal para conseguir un ajuste perfecto con error cuadrático
nulo. Las situaciones de ese tipo se denominan sobreajustes y debe evitarse a toda costa.
Efectivamente, en el contexto de la Minería de Datos, no valoramos una estrategia S concreta
sino el algoritmo que la genera. Es decir, nuestra intención es programar lo mejor posible un
algoritmo que, dada una muestra ξ de tamaño n de (X, Y ) y una observación independiente
X0 de X, determine a partir de ξ una estrategia Sξ para predecir el valor de Y correspondiente
a X = X0 . Así pues, una vez establecido el algoritmo, la predicción Ŷ = Sξ (X0 ) se entiende
66
X4 X1
X2
X3
Figura 4.10: Esquema de árbol con cuatro iteraciones en la ramificación
como una función real de ξ y X0 . Podemos entender intuitivamente que la varianza de Ŷ |X=x0
crece según aumenta el grado de sobreajuste en el algoritmo programado, porque aquello que
está diseñado para explicar demasiado bien una muestra concreta ξ puede diferir mucho de lo
que está diseñado para explicar otra distinta ξ 0 . En ese sentido, el sobreajuste es el enemigo
oculto en la Minería de Datos al igual que la violación de los supuestos del modelo lo era en la
Inferencia Paramétrica.
Para evitar esta patología podemos optar por imponer una cantidad mínima de datos por
nodo. De esta manera renunciamos a árboles complejos a menos que trabajemos con grandes
muestras. También pueden limitarse los niveles de ramificación aplicando mecanismos automá-
ticos de “poda” a partir de un árbol con error 0. La intensidad de dicha poda puede modularse
a partir de cierto parámetro α cuyo valor mínimo es 0 (consultar los detalles en Hastie et al.
(2008)).
Predicción categórica de una variable cualitativa: Pero en este capítulo estamos inten-
tando aproximarnos al valor de una variables respuesta Y que no es numérica sino cualitativa,
con r posibles categorías. El algoritmo del árbol de decisión puede aplicarse también en ese
caso con una serie de modificaciones bastante naturales. Si Y es una variable cualitativa el
algoritmo consiste en asignar a cada nodo filial la categoría de Y predominante en el mismo
(en lugar de la media). La dicotomización se realiza en cada paso de manera que se logre la
máxima reducción en una medida de error de clasificación expresada en función de cierto ín-
dice. Existen diferentes formas de medir dicho error 5 . En este manual y en cada nodo nos
decantamos por el denominado índice de Gini que se define como sigue: si en cierto nodo N se
denota por p̂N (i) la proporción de datos pertenecientes de la categoría i de Y en N , se define
la medida de Gini de N mediante
r
X
G(N ) = p̂N (i) 1 − p̂N (i) (4.22)
i=1
Nótese que, si todos los datos del nodo N pertenecieran a una cierta categoría, los r sumandos
de G(N ) serían nulos. Lo mismo ocurriría si ninguno perteneciera a i. Es más, para cada i,
5
Hastie et al. (2008)
67
p̂N (i) es la media de la función indicador que asigna 1 a los datos en i y 0, al resto, siendo
p̂N (i) 1 − p̂N (i) su varianza. Así pues, el índice de Gini puede interpretarse como una suma
de varianzas o, más concretamente, como una varianza en el sentido más amplio definido en
(1.9), lo cual refuerza la analogía con el caso numérico. Así pues, la ganancia en términos de
Gini al pasar de un nodo parental a dos filiales se define como la diferencia ente la medida de
Gini (varianza) inicial y la media ponderada (en función de las proporciones de datos) de las
varianzas filiales. A lo largo de las diferentes iteraciones podemos ir contabilizando, para cada
variable, la ganancia que supondría en la misma una dicotomización óptima según el índice de
Gini (aunque en la práctica sólo se generará una ramificación en la variable que logre la máxima
reducción). Este contador acumulado da lugar a la denominada función de importancia.
Importancia normalizada
0 25 50 75 100
petleng
Variable independiente
petwidt
sepleng
sepwidt
0,0% 0,2% 0,4% 0,6%
Importancia
Método de crecimiento:CRT
Variable dependiente:species
Figura 4.11: Función de importancia para iris data
4.3.4. Validación de estrategias:

Como ya hemos comentado, los métodos propios de la Minería de Datos, como K-NN y el
árbol de decisión, no se basan en suposiciones sobre las distribuciones de los datos, sino que
elaboran una explicación ad hoc de la clasificación observada en la propia muestra. Luego, si
estimáramos los riesgos inherentes a la estrategia reclasificando directamente las observaciones
de la muestra que hemos utilizado para diseñarla, estaríamos incentivando claramente el sobre-
ajuste, lo cual conduciría a un ingenuo optimismo sobre la validez del método. En este contexto
resulta más apropiado estimar los riesgos a partir de datos cuya categoría de procedencia se co-
nozca, pero que no se hayan utilizado para diseñar la estrategia. Por ejemplo, podemos dividir Página 1
la muestra en dos partes, no necesariamente iguales: una para diseñar la estrategia (muestra
de entrenamiento) y otra para validarla (muestra de validación). Es decir, los datos que se
utilizan para estimar los riesgos de la estrategia no se han usado en el diseño de la misma. Si
finalmente estos riesgos se consideran apropiados, podemos reconstruir la estrategia contando
con la muestra completa.
Por ejemplo, en la figura 4.12 podemos apreciar el árbol de decisión que ofrece el programa
SPSS para determinar la especie a la que pertenece un lirio a partir de sus medidas de pétalo
y sépalo. A partir de (aproximadamente) la mitad la muestra y mediante dos iteraciones (ya
que en el algoritmo hemos impuesto un tamaño mínimo de los nodos de 10) se ha diseñado un
árbol de decisión con tres nodos terminales. Para este árbol se estiman, mediante el resto de la
68
muestra, los riesgos de clasificación errónea, que son concretamente del del 0 % para setosa, 9 %
para vesicolor y del 15 % para virgínica. En la figura 4.11 se expresa la función de importancia
asociada al árbol.
species
Nodo 0
Categoría % n
setosa 33,3 50
setosa vesicolor 33,3 50
vesicolor virginica 33,3 50
virginica Total 100,0 150
petleng
Mejora=0,333
<= 2,450 > 2,450
Nodo 1 Nodo 2
Categoría % n Categoría % n
setosa 100,0 50 setosa 0,0 0
vesicolor 0,0 0 vesicolor 50,0 50
virginica 0,0 0 virginica 50,0 50
Total 33,3 50 Total 66,7 100
petwidt
Mejora=0,260
<= 1,750 > 1,750
Nodo 3 Nodo 4
Categoría % n Categoría % n
setosa 0,0 0 setosa 0,0 0
vesicolor 90,7 49 vesicolor 2,2 1
virginica 9,3 5 virginica 97,8 45
Total 36,0 54 Total 30,7 46
Figura 4.12: Árbol de decisión para irisdata
Técnicas derivadas del árbol de decisión: El árbol de decisión 4.12 se ejecuta en SPSS
a través del cuadro de diálogos 4.8. Existen otros métodos más sofisticados basados en los
mismos principios que son los preferidos por numerosos autores (ver Hastie et al. (2008)) para
afrontar un problema de clasificación: Random Forest y Adaboost, ejecutables en R mediante las
librerías randomForest y Adabag, repectivamente. El primero consiste básicamente en diseñar
diferentes árboles simples a partir de submuestras escogidas al azar (bootstrap) y con un
determinado número de variables escogidas también al azar. Cada árbol simple consiste en una
única partición binaria escogiendo la variable y el punto óptimo; en el árbol complejo final se
clasifica la observación en la categoría mayoritaria entre todos los árboles simples diseñados.
Por su parte, el método adaboost itera un árbol binario incrementando el peso de los datos mal
clasificados; cada árbol binario se pondera a su vez en función de su capacidad predictiva,
Página 1 de
manera que se establece como estrategia de clasificación la suma ponderada de los diferentes
árboles binarios.
4.3.5. Indicaciones sobre Redes Neuronales

Por último, mencionamos otro conjunto de técnicas de regresión y clasificación alternativa
a las anteriores denominada genéricamente por Redes Neuronales, que consisten en explicar
69
la variable respuesta en función de una o varias capas sucesivas de variables ocultas que se
calculan a partir de las observadas, combinando funciones lineales con otras no lineales muy
específicas.
El modelo más simple, con una única capa, consiste en explicar un vector respuesta P =
(P1 , . . . , Pk )0 a partir de las variables observadas Y = (Y1 , . . . , Yp0 ) como P = g ◦ T ◦ Z(Y ),
donde g = (g1 , . . . , gk )0 , Z = (Z1 , . . . , ZM )0 y T = (T1 , . . . , Tk )0 . Las opciones más habituales
para un problema de clasificación respecto a k + 1 categorías son las siguientes:
0
Zm (Y ) = L(α0m + α1m Y ), m = 1, . . . , M (4.23)
0
Tk (Z) = β0k + β1k Z, k = 1, . . . , K (4.24)
eTk
gk (Tk ) = , k = 1, . . . , K (4.25)
1 + eTk
En ese caso, Pk (Y ) se entiende como la probabilidad estimada de que la observación Y perte-
nezca a la categoría k-ésima del factor respuesta. Nótese que, si K = M = 1 y, además, g = Id
(que es la opción que se suele considerar por defecto), el método consistiría e una regresión
logística binaria, luego podríamos considerar las redes neuronales como una generalización en
sentido muy amplio de dicha técnica.
Peso sináptico > 0

Peso sináptico < 0
Sesgo
famhist=0
famhist=1
Sesgo
sbp
H(1:1)
tobacco
H(1:2) chd=0
ldl
H(1:3) chd=1
adiposity
H(1:4)
typea
H(1:5)
IMC
alcohol
age
Función de activación de capa oculta: Sigmoide
Función de activación de capa de salida: Softmax
Página 1
Figura 4.13: Probabilidad de enfermedad coronaria según red neuronal
Los detalles sobre la estimación de estos parámetros podemos encontrarlos en la sección

11.4 de Hastie et al. (2008). Se trata posiblemente del método de clasificación más flexible y se
70
aplica con frecuencia al reconocimiento de imágenes. Por contra, genera un efecto denominado
de “caja negra” que impide entender la influencia de las variables observadas en la respuesta.
En la figura 4.13 podemos apreciar la red neuronal utilizada para pronosticar la probabilidad
de enfermedad coronaria a partir de los datos de South African Heart Disease Study.
Si efectuamos una comparativa entre los métodos estudiados (LDA, regresión logística, ve-
cino más próximo, árbol de decisión y redes neuronales) a través de los datos de Southafrica
Heart Disease, observamos en el cuadro 4.3 muy altas correlaciones entre las probabilidades
pronosticadas de LDA Bayesiano, regresión logística binaria y red neuronal (muy especialmente
entre los dos primeros), como cabría esperar por lo razonado anteriormente. Además, los por-
centajes de clasificación correcta para estos métodos son muy similares (en torno al 85 % para
los sanos y al 55 % para los enfermos coronarios); para el método LDA minimax los porcentajes
de clasificación correcta son del 70 % para los sanos y 73 % para los enfermos; los otros dos
métodos (KNN y árbol) correlacionan menos entre sí y con los demás. Presentan porcentajes
de clasificación correcta en torno al 95 % para los sanos y al 25 % para los enfermos, calculadas
a partir de las muestras de validación.
Reg Log KNN Árbol Red Neuro.

LDA 0.997 0.631 0.735 0.981
Log 0.639 0.738 0.980
KNN 0.524 0.622
Árbol 0.725
Cuadro 4.3: Correlaciones entre probabilidades pronosticadas de enfermedad coronaria
71
72
Capítulo 5
Reducción dimensional
En este capítulo, el único exclusivamente multivariante, se recogen una serie de técnicas

que tienen por objetivo simplificar un conjunto de datos multidimensional, aunque centraremos
mayormente nuestro interés en el análisis de componentes principales. El denominador común
de estas técnicas es que, en última instancia, se fundamentan en el teorema 1.5.1. Hay que
destacar que un estudio de este tipo se enmarca en una fase inductiva, pues el producto final
no es la aplicación de un test de hipótesis sino un gráfico en dimensiones reducidas que permita
una primera visión global de nuestra muestra, a partir de la cual podamos formular distintas
hipótesis. Por ello se ha optado por una redacción de la mayor parte del capítulo en lenguaje
muestral.
Ya hemos comentado que el núcleo del capítulo lo constituye el análisis de componentes
principales, aunque esta técnica puede considerarse un opción particular del denominado aná-
lisis factorial. No parece existir unanimidad de criterios a la hora de catalogar y distinguir los
diferentes métodos. En este volumen se ha optado por presentarlos como un única técnica de
reducción dimensional de carácter dual, que se denomina análisis de componentes principales
desde el punto de vista de la observaciones (filas) y análisis factorial desde el punto de vista de
la variables (columnas). Aunque existen variedades del análisis factorial al margen del análisis
de componentes principales, nos hemos limitado a introducirlas brevemente. Así mismo, el aná-
lisis de correspondencias, que en cierta forma puede entenderse como una extensión del análisis
de componentes principales, se presenta de forma muy escueta, y más aún su generalización,
que es el análisis de correspondencias múltiple. Por último, se ha incluido una sección en la
que se ilustra brevemente la utilidad del uso de componentes principales en un problema de
regresión lineal múltiple.
Desde el punto de vista teórico, debemos entender en esencia que el teorema 1.5.1 viene
a proponer un cambio de la base vectorial inicial, dada por las variables medidas, a una base
ortonormal constituida por los autovectores de la matriz de covarianzas. De esta transformación
obtenemos nuevas variables incorreladas entre sí que se denominarán componentes principales,
lo cual puede permitir una comprensión más fácil de nuestros datos.
5.1. Componentes principales

Antes de abordar un estudio de carácter muestral y con un propósito meramente didáctico,
empezaremos con una definición probabilística de las componentes principales.
Consideremos un vector aleatorio p-dimensional X de media nula y matriz de covarianzas Σ,
73
Capítulo 5 Justificación geométrica
la cual descompone según el teorema 1.5.1 en Γ∆Γ0 , con Γ = (γ1 . . . γp ) y ∆ = diag(δ1 , . . . , δp ).

Para cada j = 1, . . . , p se define la j-ésima componente principal mediante
Uj = γj0 X (5.1)
De esta forma, las p componentes principales ordenadas componen un vector aleatorio p-
dimensional U que se obtiene como transformación lineal (rotación en sentido amplio) del
vector original X, concretamente mediante U = Γ0 X, y que verifica, por lo tanto, Cov(U ) = ∆.
De ello se deduce que las componentes principales son incorreladas entre sí. También se deduce
del teorema 1.5.1 que, de todas las posibles proyecciones de X sobre algún eje de Rp , U1 es la
que presenta una máxima varianza. U1 es concretamente la proyección sobre el eje hγ1 i y su
varianza es δ1 .
Ejercicio 74. Probar que, dada una proyección de U sobre un determinado eje hγi, se verifica
que es incorrelada con U1 si, y sólo si, γ ⊥ γ1 .
Ejercicio 75. Probar que, de todas las posibles proyecciones de X sobre algún eje de Rp
incorreladas con U1 , la que presenta una varianza máxima es U2 .
Así se van obteniendo sucesivamente las componentes principales hasta Up . Así pues, “pasar
a componentes principales” se entiende como considerar, en lugar del vector original X, otro
vector transformado o “rotado” U cuyas componentes son incorreladas y están ordenadas en
función de sus varianzas. ¿Qué ventajas comportan ambas cualidades? Por el momento nos li-
mitaremos a aportar una explicación intuitiva, si bien el el resto del capítulo presentaremos un
razonamiento más formal: primeramente, las ventajas que se derivan de la propiedad probabi-
lística o estadística de incorrelación entre las componentes principales (lo cual bajo el supuesto
de p-normalidad equivaldría a independencia) son las mismas que se derivan de la propiedad
vectorial de la ortogonalidad. Más específicamente, dado que nos proporciona una clara des-
composición de la varianza total, permite conocer con precisión que porcentaje de la misma
perdemos si eliminamos una componente principal concreta. Segundo, y en lo que respecta
precisamente a la eliminación de componentes principales, serán las últimas, es decir, las de
menor varianza, las candidatas a ser suprimidas. Estamos afirmando pues que la importancia
de las diferentes componentes principales viene dada por la magnitud de su varianza y que,
por lo tanto, el orden en el que se obtienen establece una jerarquía real desde el punto de vista
de la reducción dimensional. Intentaremos justificar tal afirmación desde un punto de vista
geométrico en la siguiente sección.
5.2. Justificación geométrica

En esta sección se trata el problema de simplificar la dimensión de un conjunto de datos,
tanto desde el punto de vista de las observaciones, es decir, de las filas de la matriz (1.1), como
de las variables, es decir, de las columnas. Debemos tener presente la notación introducida en el
capítulo preliminar y, en especial, la distancia (1.36) definida en el conjunto de las matrices de
dimensión n × p. Por otra parte, dado k ≤ p, se denota por Hk el conjunto de las subvariedades
afines de dimensión k en Rp . Dada una matriz X ∈ Mn×p , con matriz de covarianzas muestral
S, consideraremos su descomposición tipo (1.40) en función de sus autovalores d1 , . . . , dp y sus
respectivos autovectores, g1 , . . . , gp , pero distinguiendo entre los k primeros y los p−k restantes
mediante 0
D1 0 G1
S = (G1 |G2 ) (5.2)
0 D2 G02
74
El análisis de componentes principales (PCA) se basa en el siguiente resultado, que es conse-

cuencia del lema 1.5.2:
Teorema 5.2.1. Dados 0 ≤ k ≤ p y X ∈ Mn×p con matriz de covarianzas S que descompone
según (5.2), se verifica
mı́n d2n,p X, Xk : Xki ∈ H ∀i ≤ n para algún H ∈ Hk = tr(D2 ),

(5.3)
y se alcanza con Xki = x + PhG1 i (Xi − x) ∀i ≤ n.

La figura 5.1 ilustra el teorema anterior para un problema bidimensional.
6 X[2]

s s
@ @
@
s @
@ @ @ @
s @ @ s @ @s
@ @
@ @s
@ @
@ @s
@s X[1]
-
Figura 5.1: Proyección de observaciones
En el caso trivial k = 0, el teorema afirma que el vector de Rp constante por el que debemos
reemplazar las observaciones Xi con el menor error cuadrático posible es la media aritmética x,
siendo la varianza total muestral, definida en (1.37), la medida de dicho error.
Ejercicio 76. Probar que s2T = pj=1 dj
P
A medida que aumenta el valor de k, la distancia respecto a la simplificación Xk disminuye.

Expresándolo de manera inversa diríamos que al reemplazar Pk las observaciones originales se
explica sólo una parte de la varianza total, concretamente j=1 di . Esta cantidad es máxima
respecto al conjuno de subvariedades afines k-dimensionales y se alcanza proyectando sobre el
subespacio generado por los k primeros autovectores de S, salvo traslaciones. La proporción de
varianza total explicada por los mismos es pues
Pk
j=1 di
(5.4)
tr(S)
Que esta proporción sea próxima a 1 para un valor de k pequeño se traduce en que las n
observaciones reales se encuentran muy próximas a una subvariedad afín de baja dimensión, lo
75
cual equivale a un fuerte grado de correlación lineal (afín) entre las variables medidas para los
datos observados. Ésa es la cualidad que permite obtener una reducción dimensional profunda
mediante la técnica PCA, cosa que parece más factible bajo el supuesto de p-normalidad.
En todo caso, el método puede generalizarse mediante el uso de las denominadas curvas y
superficies principales, que pueden aproximarse mediante iteraciones descritas en Hastie et al.
(2008), si bien es en Gifi (1990) donde se le dedica más espacio a este tipo de problema.
En lo que resta supondremos sin pérdida de generalidad (ya veremos por qué) que x = 0.
En ese caso y para cada j = 1, . . . , p, se denota U[j] = X · gj , es decir, el vector de Rn que
recoge las proyecciones de cada observación Xi , i = 1, . . . , n, sobre el eje hgj i determinado por
el j-ésimo autovector de S. Dicho eje se denomina j-ésimo eje principal y U[j] se denomina
j-ésima componente principal. Las p componentes principales ordenadas constituyen una nueva
matriz U ∈ Mn×p definida mediante
U = XG (5.5)
que expresa las coordenadas de X respecto de la base ortonormal canónica de autovectores G.
Dado k ≤ p, se denota por U la matriz n × k compuesta por las k componentes principales,
cuyas filas y columnas se denotarán con el mismo criterio que en (1.1). Por otra parte, se denota
E = (U[k + 1], . . . , U[p])G02 ∈ Mn×p . En ese caso, se sigue de (5.5) que
X = UG01 + E (5.6)
siendo UG01 la matriz en Mn×p que permite alcanzar las distancia mínima a X en el teorema
5.2.1.
Ejercicio 77. Probar que las componentes principales son incorreladas entre sí. Probar que
el primer eje principal es aquél sobre el que hay que proyectar las observaciones para obtener
una máxima varianza, que vale d1 .
Ejercicio 78. Probar que el segundo eje principal es aquél sobre el que hay que proyectar
las observaciones para obtener la máxima varianza de entre todas la variables incorreladas con
la primera componente principal, que vale d2 , y así sucesivamente. Probar que el último eje
principal es aquél sobre el que hay que proyectar para obtener una mínima varianza.
Ejercicio 79. ¿Cómo se interpreta dp = 0? ¿Cómo se interpreta |Σ| = 0 para un vector
aleatorio con distribución Np (µ, Σ)?
Ejercicio 80. Probar que la matriz de covarianzas de U es SU = D1 .
Así pues, los ejes principales resuelven el problema de maximización de la varianza, mientras
que los ejes discriminantes, estudiados en los capítulos 3 y 4, solucionan el problema de maxi-
mización relativo a la discriminación entre categorías, que a su vez puede entenderse como una
maximización de correlaciones lineales. Esto puede llevar a una cierta confusión en problemas
multivariantes con presencia de un factor cualitativo (como, por ejemplo, el caso del irisdata de
Fisher). La figura 5.2 trata de clarificar la relación que se da en tales casos entre los distintos
tipos de ejes, teniendo en cuenta que, según la notación utilizada en 3.13, la varianza total s2
de una proyección dada, que también podemos denotar mediante SCT , descompone según la
figura 2.5 mediante SCT = SCH + SCE.
Al contrario que en el caso de los ejes discriminantes, el problema de maximización de la

varianza es sensible ante cambios de escala en las variables medidas. En la figura 5.3 se ilustra
el efecto sobre el cálculo de los ejes principales de un cambio de escala en la variable del eje OY.
76
Capítulo 5 Representación de observaciones y variables
Ejes principales
SCT=SCH+SCE
Maximizar
F α SCH/SCE R2=SCH/SCT
Ejes discriminantes Dummys Ejes canónicos
Figura 5.2: Ejes principales, canónicos y discriminantes
Sin embargo, podemos lograr artificialmente un método de reducción dimensional invariante

ante cambios de escala si trabajamos con las variables tipificadas, que es lo que consideraremos
por defecto. Ello equivale a trabajar en todo momento con la matriz de correlaciones muestral
R en lugar de la matriz de covarianzas original S. Dado que tr(R) = p, (5.4) es igual a
Pk
j=1 di
d(k) := (5.7)
p
Al trabajar con variables tipificadas podemos suponer que la media es nula, como indicábamos
antes.
5.3. Representación de observaciones y variables

−1/2
Definimos la matriz F = UD1 de dimensiones n × k que se compone de los siguientes
elementos1  
F1 [1] . . . F1 [k]
F =  ... ..  (5.8)

. 
Fn [1] . . . Fn [k]
Para cada i, la transposición de la fila i-ésima de F es un vector de Rk que se denota2 por
→
−
F i y decimos que está compuesto por las k puntuaciones factoriales de la observación Xi .
→
−
Los diferentes vectores columnas de F, que se denotan por F [j] ∈ Rn , para j = 1, . . . , p, se
1
El paso de la matriz U a F responde, como veremos más adelante, a la intención de facilitar la representación
de las variables, aunque ello vaya en detrimento de la representación de observaciones.
2
Con carácter excepcional distinguiremos en esta sección los vectores de los escalares mediante el uso de
flechas en la notación.
77
Figura 5.3: Efecto de un cambio de escala en los ejes principales
√
denominan ejes factoriales y, salvo una constante (dado que todos tienen norma Euclídea n),
constituyen un sistema ortonormal3 .
Por otra parte, para que se siga verificando la ecuación (5.6), transformamos de manera
1/2
inversa G1 definiendo Λ = G1 D1 ∈ Mp×k . La matriz Λ se expresará así
 
λ1 [1] . . . λk [1]
Λ =  ... ..  (5.9)

. 
λ1 [p] . . . λk [p]
Para cada l = 1, . . . , p, la transposición de la fila l-ésima de Λ es un vector de Rk que se denota

→
− →
−
por λ [l] y se denomina vector de las k componentes de la variable X [l]. Los diferentes vectores
→
−
columnas de Λ se denotan por λ j ∈ Rp , para j = 1, . . . , k, y determinan los k primeros ejes
principales. Por lo tanto, constituyen un sistema ortogonal en Rp . La ecuación (5.6) se expresa
en estos nuevos términos mediante
X = FΛ0 + E (5.10)
O, lo que es lo mismo,
     
X1 [1] . . . X1 [p] F1 [1] . . . F1 [k] λ1 [1] . . . λ1 [p]
 .. ..  =  .. ..  ·  .. ..  + E
 . .   . .   . .  (5.11)
Xn [1] . . . Xn [p] Fn [1] . . . Fn [k] λk [1] . . . λk [p]

n−1 F0 F = Id, F0 E = 0 (5.12)
3
Lo cual facilita la representación de las variables a partir de sus coordenadas respecto a este sistema.
78
5.3.1. Representación de observaciones

Si reproducimos la igualdad (5.10) desde el punto de vista de las filas (observaciones)
tenemos que, para cada i = 1, . . . , n,
→
− Pk →
− →
−
Xi = j=1 Fi [j] · λ j + E i (5.13)
Además, se sigue del teorema 5.2.1 que

n
−1
X →
−
n k E i k2Rp = d2n,p (X, FΛ0 ) (5.14)
i=1
= d2n,p (X, UG01 ) (5.15)
Pp
= j=k+1 dj (5.16)
En consecuencia, al reemplazar las observaciones originales por los puntos expresados me-
diante
Pp las puntuaciones factoriales cometemos unos error global cuya cuya media cuadrática
es j=k+1 dj . Si este término es pequeño podemos considerar que una representación de las
observaciones en dimensión k, como lo que tenemos en la figura 5.5, se aproxima a la realidad.
5.3.2. Representación de variables

Si reproducimos la igualdad (5.10) desde el punto de vista de las columnas (variables)
tenemos que, para cada l = 1, . . . , p,
→
− Pk →
− →
−
X [l] = j=1 λj [l] · F [j] + E [l] (5.17)
Por otra parte, dado que estamos trabajando con los datos tipificados, se verifica que R =
n−1 X0 X. Se deduce nuevamente de (5.10) junto con (5.12) que
R = ΛΛ0 + n−1 E0 E (5.18)
En particular, el coeficiente de correlación entre las variables l-ésima y s-ésima es

→
− →
− →
− →
−
rls = h λ [l], λ [s]i + n−1 h E [l], E [s]i (5.19)
Luego, podemos reproducir la correlación entre las variables originales multiplicando sus com-
ponentes en dimensión k, salvo un error determinado por el segundo sumando que acotaremos
a continuación: en primer lugar, se deduce de la desigualdad de Cauchy-Schwarz que
−1 →− →
− →
− →
−
q
n h E [l], E [s]i ≤ n−1 k E [l]k2Rn · n−1 k E [s]k2Rn (5.20)
Esta desigualdad puede expresarse también en función de las componentes de la variable te-
niendo en cuenta la igualdad (5.18) aplicada a los elementos de la diagonal. Efectivamente,
para todo l = 1, . . . , p,
→
− →
−
n−1 k E [l]k2Rn = 1 − k λ [l]k2Rk (5.21)
79
→
−
Si se denota el término k λ [l]k2Rk por h2l (que estará comprendido en todo caso entre 0 y 1), se
deduce entonces de (5.19), (5.20) y (5.21) que
→
− →
− q
rls − h λ [l], λ [s]i ≤ (1 − h2l )(1 − h2s ) (5.22)

En lo sucesivo, los términos h21 , . . . , h2p se denominarán comunalidades de las respectivas va-
riables. La desigualdad (5.22) se interpreta pues como sigue: si tanto la comunalidad de la
variables l-ésima como de la de la s-ésima son próximas a 1, podremos reproducir la correla-
ción entre ambas variables, salvo un pequeño error, mediante el producto en Rk de sus vectores
de componentes, como los que se representan, por ejemplo, en la figura 5.5. En tal caso, si dichos
vectores están aproximadamente en la misma dirección cabrá pensar en una fuerte correlación
lineal, que será positiva si están en el mismo sentido y negativa en caso contrario; por contra,
si están direcciones aproximadamente perpendiculares, cabrá pensar en una escasa correlación
lineal.
Ejercicio 82. Probar que, si se denota h2 = p1 pj=1 h2j (es decir, la media aritmética de las
P
comunalidades), se verifica
d(k) = h2 (5.23)
Por lo tanto, la media aritmética de las comunalidades equivale a laP proporción de varianza
total explicada por las k primeras componentes principales. Luego, si pj=k+1 dj es próximo a
0 tendremos un buen comportamiento global de las comunalidades y, por lo tanto, el análisis
de los vectores componentes en Rk nos dará una idea bastante aproximada de la correlación
real entre las variables estudiadas.
→
− →
−
Ejercicio 83. Probar que λj [l] es el coeficiente de correlación lineal entre X [l] y U [j]. Por
lo tanto, la matriz Λ recoge las correlaciones entre las variables originales y las componentes
principales.
Podemos apreciar en las ecuaciones (5.13) y (5.17) que los papeles que desempeñan las
matrices F y Λ se permutan tal y como se indica en el cuadro 5.1 según representemos las
observaciones o las variables.
k ejes k coordenadas
n observaciones Λ F
p variables F Λ
Cuadro 5.1: Dualidad observaciones-variables
También hemos de advertir que la representación de las variables mediante sus componentes
constituye tan sólo una simplificación de la información que nos aporta ya de por sí la matriz
de correlaciones R. Concretamente, consideramos en (5.18) la aproximación R ' ΛΛ0 , donde la
diferencia entre ambas se denota por Ψ, denominándose varianzas especificas los elementos de
su diagonal, ψjj , j = 1, . . . , p. De esta forma, la igualdad (5.21) puede expresarse así
1 = h2j + ψjj , 1≤j≤p (5.24)
5.3.3. Representación conjunta de observaciones y variables

Dado que tanto las observaciones como las variables pueden identificarse de manera apro-
ximada con vectores de Rk , según (5.13) y (5.17), respectivamente, podemos representarlas
80
conjuntamente mediante un gráfico k-dimensional que debemos interpretar según la siguien-

te igualdad, que se deduce indistintamente de cualquiera de las dos igualdades mencionadas.
Concretamente, se verifica para cada i = 1, . . . , n y para cada l = 1, . . . , p:
− →
→ −
Xi [l] = h F i , λ [l]i + Ei [l] (5.25)
𝜆𝜆⃗[j]
B C
Figura 5.4: Relación observaciones-variables
En la figura 5.4 se representan conjuntamente en dimensión k = 2 una variable aleatoria

→
− →
−
X [j] medida en n individuos, identificada con el vector λ [j], y tres observaciones multidimen-
→
− → − → − →
−
sionales A , B y C . En este caso, la observación A se caracteriza por un valor de la variable
→
− →
−
j-ésima por encima de la media, la observación B , por un valor por debajo de la media y C ,
por un valor en torno a la media.
Por ejemplo, en un estudio (Bravo et al. (2011)) realizado en el CENSYRA de Badajoz
sobre p = 8 variables que caracterizan la motilidad de los espermatozoides en carneros, a partir
de una muestra de n = 383 observaciones (ver tabla 1.1), se obtuvieron los resultados que se
muestran en los cuadros 5.2 y 5.3. Del cuadro 5.2 se deduce que podemos explicar un h2 = 82 %
de la varianza total proyectando sobre los dos primeros ejes principales, es decir, calculando
→
− →
−
F [1] y F [2]. Según el cuadro 5.3, hay variables como vcl o vap que quedan explicadas casi
perfectamente de esta forma, mientras que bcf y wob quedan deficientemente representadas.
En la parte derecha se recoge la matriz de componentes que se representará en la figura 5.5
junto con las puntuaciones factoriales y que permite simplificar la matriz de correlaciones R.
En la figura 5.5 se aprecia claramente cómo las variables vcl, vsl y vap correlacionan fuerte
y positivamente entre sí, y correlacionan débilmente con el resto; por otra parte, las variables
wob, lin y str correlacionan fuerte y positivamente entre sí y negativamente con alh y bcf, si
81
bien bcf y wob no quedan satisfactoriamente representada por este gráfico, como hemos dicho
anteriormente. Además, podemos apreciar qué espermatozoides presentan valores altos, medios
o bajos para las distintas variables. Se ha superpuesto en el gráfico la circunferencia unidad
para que podamos apreciar qué variables presentan una comunalidad próxima a 1.
Varianza total explicada
Sumas de las saturaciones al cuadrado

Autovalores iniciales de la extracción
% de la % de la
Componente Total varianza % acumulado Total varianza % acumulado
1 3,834 47,929 47,929 3,834 47,929 47,929
2 2,743 34,283 82,213 2,743 34,283 82,213
3 ,860 10,747 92,960
4 ,366 4,578 97,538
5 ,161 2,014 99,552
6 ,033 ,410 99,962
7 ,002 ,030 99,992
8 ,001 ,008 100,000
Método de extracción: Análisis de Componentes principales.
Cuadro 5.2: Autovalores de R
Comunalidades Matriz de componentesa
Inicial Extracción Componente

vcl 1,000 ,991 1 2
vsl vcl -,010 ,995
1,000 ,971
vsl ,632 ,756
vap 1,000 ,993
vap ,455 ,887
LIN% 1,000 ,943 LIN% ,956 -,169
STR% 1,000 ,654 STR% ,678 -,441
WOB% 1,000 ,704 WOB% ,839 ,029
alh 1,000 ,881 alh -,862 ,372
bcf 1,000 ,440 bcf -,639 ,180
Método de extracción:
Método de extracción: Análisis de Componentes principales. Análisis de componentes principales.
Cuadro 5.3: Comunalidades y matriz de componentes
a. 2 componentes extraídos
5.3.4. Rotación de ejes

Dado que nuestra interpretación de un gráfico como el de la figura 5.5 viene dada exclusiva-
mente en términos del producto escalar en Rk , permanecerá invariante ante cualquier rotación
que apliquemos solidariamente a las puntuaciones factoriales y los vectores de componentes
de las variables. En particular, podemos intentar aplicar una rotación que aproxime lo mejor
→
−
posible los diferentes vectores λ [j] a cualquiera de los k ejes de coordenadas. Para conseguir
tal propósito existen diversos métodos iterativos como varimax o equamax, que se describen
con mayor detalle en Rencher (1995). Si finalmente se consigue el objetivo, el análisis de las
variables será más sencillo, pues aquéllas cuyos vectores estén en la dirección de un mismo eje
de coordenadas presentarán una fuerte correlación lineal entre ellas, y guardarán poca corre-
lación con las que estén sobre un eje diferente. En el gráfico 5.5 no se precisaría una rotación
porque las variables ya se agrupan aproximadamente en las direcciones de los ejes OX y OY.
82
Página 1
Capítulo 5 Análisis factorial
Figura 5.5: Puntuaciones factoriales y matriz de componentes
Esto resulta especialmente útil cuando k > 2, porque en ese caso no podemos representar
satisfactoriamente las variables. Además, el análisis numérico de las correlaciones a través de la
matriz de componentes puede resultar en principio complicado. No obstante, tras una rotación
adecuada, el análisis de dicha matriz resulta más sencillo porque deja más claro con qué eje se
identifica cada variable (aquél en el que se obtenga una coordenada próxima a ±1).
En el cuadro de diálogos 5.6 se indica a grandes rasgos cómo ejecutar un analisis de com-
ponentes princiapales (PCA) con SPSS. La ejecución con R se realiza a través del comando
princomp.
Figura 5.6: Cuadro de diálogos PCA
83
Capítulo 5 Análisis factorial
5.4. Análisis factorial

De la figura 5.5 4 se desprende la existencia de dos conglomerados de variables si utilizamos
como criterio de afinidad la correlación lineal: por un lado tenemos vcl, vsp y vap, que se
agrupan en torno al eje OY tras una rotación varimax y, por otro, el resto, que se agrupan en
torno al eje OX (aunque en el segundo grupo distinguimos lin, str y wob que correlacionan
positivamente de alh y bcf, que correlacionan negativamente con las anteriores). Desde un
punto de vista formal podríamos definir factor como una clase de equivalencia en el conjunto de
variables si consideramos una relación basada en la correlación lineal. Desde un punto de vista
práctico, es tarea del investigador experimental dar sentido a dichos factores. En el ejemplo
que nos ocupa el primer factor (eje OY) se identificaría aproximadamente con el concepto
biológico de velocidad espermática, mientras que el segundo (eje OX) se identificaría de
manera aproximada con el de progresividad. Queremos decir a grandes rasgos que esos son
en esencia los dos factores a tener en cuenta en un espermatozoide de este tipo.
Nada nos garantiza en principio una agrupación clara de nuestras variables en torno a un
número reducido de factores, aunque existen parámetros que pueden darnos algunas pistas
antes de proceder con el análisis, como el coeficiente de Kaiser-Meyer-Olkin, definido en (5.26)
2
P
i6=j rij
KM O = P 2
P 2
(5.26)
i6=j rij + i6=j aij
Se basa en la idea de que, de darse una buena agrupación de muchas variables en pocos
factores, deberían observarse fuertes correlaciones simples pero bajas correlaciones parciales.
Concretamente, la suma de coeficientes de correlación simple al cuadrado entre cada posible par
de variables diferentes debe ser grande en relación con la suma de los coeficientes de correlación
2
P
parcial al cuadrado entre cada par, conocidas el resto, que se denota por i6=j aij en (5.26). En
la práctica, valores de KMO inferiores a 0.6 suelen considerase un mal indicio de cara a una
clara simplificación del problema de correlación.
Ejercicio 84. Razonar por qué un valor de KMO próximo a 1 se asocia a una reducción
profunda en el análisis factorial.
En otros contextos se utilizan también métodos como oblimin5 que aplican una transfor-
→
−
mación lineal, no necesariamente ortogonal, alos diferentes vectores λ [j], de manera que las
puntuaciones factoriales transformadas dejan de ser necesariamente incorreladas. Este tipo de
transformación permite una agrupación más clara de las variables en factores aunque éstos no
sean incorrelados y pierdan, por lo tanto, su sentido original.
5.4.1. Modelos con factores latentes

Se entiende en general el análisis factorial (AF) como la búsqueda de los factores subyacentes
en las variables consideradas. Desde ese punto de vista, el PCA podría entenderse como una de
las posibles técnicas de extracción de factores, aunque esta visión es fuente de viejas polémicas.
Muchos especialistas prefieren distinguir claramente entre el PCA y el AF, distinguiendo a su
4
La figura 5.5 no se ha obtenido tal cual de SPSS, sino que es el resultado de una superposición de dos
gráficos distintos. Sin embargo, el programa R puede producir directamente una figura de ese tipo (biplot) a
través de los comandos prcomp y biplot.
5
Rencher (1995).
84
Capítulo 5 Introducción al Análisis de Correspondencias
vez entre un AF exploratorio, orientado a la búsqueda de dichos factores, y otro confirmatorio,

que pretende validar el modelo obtenido.
Los otros métodos más conocidos del análisis factorial (o los métodos del análisis factorial
propiamente dicho, según se vea) nos los estudiaremos aquí, aunque adelantamos que son el de
máxima verosimilitud y el del eje principal. A diferencia de lo explicado hasta el momento en
este capítulo, estamos entendiendo los factores como variables reales pero latentes (no obser-
vadas) a partir de las cuales pueden obtenerse las variables observadas mediante una ecuación
lineal, salvo errores incorrelados. Estos métodos parten de restrictivas suposiciones formales
que, además, no pueden cuestionarse, pues los factores no son observables6 Todos los métodos
tienen como objetivo común una descomposición de la matriz de correlaciones tipo (5.18), es
decir
R = ΛΛ0 + Ψ (5.27)
lo cual puede conseguirse mediante un PCA o bien mediante otras técnicas, algunas de ellas,
como la del eje principal, estrechamente relacionadas al mismo. Las puntuaciones factoriales
pueden obtenerse mediante una regresión lineal donde los coeficientes se estiman como en (2.47)
pero sustituyendo SY F por la estimación de Λ. En este contexto tienen sentido técnicas como la
rotación oblicua introducido anteriormente. Para más detalle consultar Rencher (1995) o Uriel
y Aldás (2005).
Escalamiento multidimensional Se trata de otra forma de abordar el mismo problema de

representación tratado hasta ahora en el capítulo, pero intentando explicar en una baja dimen-
sión, en lugar de la matriz de observaciones X, una matriz de distancias o, más generalmente,
disimilaridades entre datos o variables; también podemos estar interesados en explicar sus si-
milaridades. En ambos casos, el método parte de una matriz de distancias o bien similaridades
que podría ser calculada a partir de las propias observaciones, lo cual puede acabar siendo
equivalente a un análisis factorial. Su interés estriba posiblemente en que ofrece un método
para representar variables de muy distinta naturaleza (de escala, ordinal, cualitativa), puesto
que lo que importa en el análisis son las distancias o similaridades que pueden definirse con
una importante componente subjetiva. Para más detalles consultar Hair et al. (1999) y Hastie
et al. (2008).
5.5. Indicaciones sobre Análisis de Correspondencias

Mientras el objetivo del PCA es representar en un espacio de baja dimensión una matriz X
correspondiente a la medición de varias variables numéricas, el del análisis de correspondencias
es, así mismo, representar en un espacio sencillo una matriz O asociadas a la medición de varias
variables categóricas.
Cuando contamos únicamente con dos variables cualitativas con r y s categorías, respecti-
vamente, la denominada tabla de contingencia o tabla de frecuencias cruzadas O puede iden-
tificarse con r − 1 vectores filas en Rs−1 o, equivalentemente, con s − 1 vectores columnas en
Rr−1 . Si r y s son grandes, puede ser interesante representar las filas o las columnas (o ambas)
en un espacio Rk de dimensión reducida, de manera que se conserve una buena parte de, en
lugar de la varianza total, de la distancia χ2 . Nótese que la varianza total expresa la distancia
6
Lo cual induce, como el lector intuirá, a numerosas polémicas al respecto y a diferentes ópticas a la hora
de abordar este tipo de problemas.
85
Capítulo 5 Multicolinealidad y PCA
en el sentido (1.5) entre la matriz de datos aleatorios X y su proyección sobre el subespacio

h1n i, que identificamos con el determinismo. Es decir, X es la matriz de columnas constantes
más próxima a X según las distancia (1.5). De similar forma el valor χ2 puede interpretarse en
términos de distancia entre la matriz valores observados O y la de valores esperados E.
r X s
X (Oij − Eij )2
χ2 = (5.28)
i=1 j=1
Eij
Nótese que se entiende (5.28) como una distancia entre la tabla de contingencia observada O y
su proyección E sobre el subespacio que identificamos con la independencia entre las variables
cualitativas. Para ser precisos, la matriz de valores esperados E es, de entre aquellas cuyas pro-
porciones condicionadas equivalen a las proporciones marginales, la más próxima a la matriz
de valores observados según cierta distancia, que en tal caso vale χ2 . Los detalles pueden estu-
diarse en Greenacre (1984) y también se desarrollan en el Manual de Análisis Multivariante7 ,
se basa pues en una generalización del teorema 5.2.1, y tiene como producto final un gráfico
denominado biplot, como el de la figura 5.7, que se interpreta de manera parecida al caso
numérico.
Cuando tenemos más de dos variables categóricas podemos optar por varios métodos: el
primero consiste en agrupar las diferentes variables categóricas en dos (con gran cantidad de
categorías) con la intención de aplicar un análisis de correspondencias simple. Por ejemplo, en
un estudio8 de relación entre la especie germinada (se distinguen 11 categorías de leguminosas)
y tres variables que caracterizan las condiciones del terreno, con 4, 3 y 4 categorías anidadas,
respectivamente, se optó por agrupar las tres variables del terreno en una única variable cuali-
tativa denominada grupo que distingue entre 28 categorías diferentes. El biplot de la izquierda
en la figura 5.7, que recoge el 63 % de la distancia χ2 9 , ilustra las asociaciones entre las especies
y las diferentes condiciones del terreno. El cuadro de diálogos 5.8 ilustra cómo se ejecuta esta
técnica con SPSS.
El segundo método consiste en una generalización de la técnica anterior denominada análisis
de correspondencias múltiples. Esta técnica, que también se estudia con detalle en Greenacre
(1984), se basa en la aplicación del análisis simple a una supermatriz denominada de índices.
Su mayor inconveniente radica en la gran pérdida en la explicación de la distancia χ2 (o de la
inercia) que suele conllevar en estos casos la representación gráfica en baja dimensión. Téngase
en cuenta que en estas circunstancias podemos llegar a barajar dimensiones por encima de 1000,
por lo que resultaría poco verosímil explicar más de un 25 % de la inercia en un simple plano.
No obstante, a pesar de su tendencia al simplismo, ofrece una visión global del problema de
correlación que puede resultar muy útil como referencia inicial. Por ejemplo, en la parte derecha
de la figura 5.7 se esquematiza en un biplot la relación conjunta entre 11 variables cualitativas,
algunas de ellas con muchas categorías, relacionadas con un estudio sobre duración da bajas
laborales10
7
Montanero (2008)
8
Pérez-Fernández et al. (2006).
9
Se suele trabajar realmente con la inercia, que es χ2 /n.
10
González-Ramírez et al. (2017).
86
Age
Children
1,0 Civil status
Bronchitis Widow/er Co-diagnostic
T. villosa
O-4 Bronchitis-AC >52 Diagnostic
Renal-colic
Divorced Gender
O-3 Operators 43 Inguinal hernia
D. glomerata Cold Occupation
O-1 Giddiness Children-yes
O-2 Flu-others Payment
0,5 Flu Married
Ax4 Gastr.-NI UNDER 15 Regime
especie Male
Gastroent. Diarrhea Headache Craft Agric-forest-fish Duration
Ax1 Axx2 Zone3
grupo
Ac-pharingitisMigraine Zone4 Meniscus-injury
Bxx3 Elementary
Zone
Axx1 Axx4 38 Sciatica Carpian tunnel
C. cristatus Indirect-payment Zone2
Dimensión 2
Bx3 Bx4 Ax3 Cxx4 Employed

0,0 D. carota Cx4 Administrative Codiag-no Low- back-pain
Professional Dorsalgia
Ax2 Bxx4 F. vulgare Cx3 Tonsilitis-strpt.
E. hirsutum M. arabica Cxx1 Adjustment-disorders
Cxx2 Tonsilitis Technician Depression
Female Zone1
Bx2 Manager
Bxx2 Bxx1 Cxx3 Couple-not-married 32 Service&sales OVER 15
Bx1 Cervical pain Self-employed
Axx3
-0,5 Children-no Anxiety Codiag-yes
Cx1
Direct-payment
Ankle-strain
Single
-1,0 Cx2 <32

Threatenet-abortion
Cervical-strain
R. crispus
Others-pregnancy
-1,5
-1,5 -1,0 -0,5 0,0 0,5 1,0
Dimensión 1
Figura 5.7: Biplots de Análisis de Correspondencias Simple (izda.) y Múltiple (dcha.)
Página 1
Figura 5.8: Cuadro de diálogos Análisis de Correspondencias Simple
5.6. Multicolinealidad y PCA

Página 1
La técnica PCA se utiliza en ocasiones para resolver el problema de multicolinealidad en
una regresión lineal múltiple. Esta afirmación hay que matizarla porque la multicolinealidad
no tiene por qué constituir un problema ni el PCA llega jamás a resolverlo. En todo caso, nos
situamos en el contexto del ejemplo 4 en la página 33, suponiendo que las variables explicativas
Z[1] . . . , Z[q] estén tipificadas. En tal caso, se dice que hay multicolinealidad cuando existe un
alto grado de correlación lineal entre las variables explicativas, lo cual no supondrá perjuicio
alguno a la hora de efectuar predicciones. La multicolinealidad da lugar a un incremento en la
la varianza de los estimadores de la ecuación. Concretamente, se verifica para todo j = 1, . . . , q:
1 1 1
var β̂ j = σ 2 · · 2 ·

2
(5.29)
n sz[j] 1 − Rj∗
donde Rj∗ denota el coeficiente de correlación múltiple de Z[j] respecto al resto de variables
explicativas. Como podemos observar, dado que las variables explicativas están tipificadas, la
87
varianza del estimador depende de la varianza del modelo σ 2 , del tamaño de la muestra n
y del grado de correlación respecto al resto de variables explicativas según el FIV, es decir,
2 −1
(1 − Rj∗ ) . En virtud de (2.37), queda claro que una fuerte dependencia lineal de Z[j] respecto
al resto de variables explicativas dificulta la significación en el contrate parcial H0 : βj = 0.
Hablando en términos más intuitivos, las redundancias lineales entre las variables explicativas
se traducen en una menor confianza en los signos de sus coeficientes estimados, como ya se
comentó en el capítulo 2.
En la figura 5.9 se ilustran los efectos de la multicolinealidad (derecha) sobre la varianza
de los estimadores de regresión a partir de variables tipificadas, en contraposición con una
situación en las que las variables explicativas son incorreladas (izquierda). En ambos casos se
calculan los coeficientes de regresión para dos muestras diferentes, Y1 e Y2 .
Z[2] 6
2 Z[2]
2 r
P Y
u hzi
β̂ 2 Phzi Y 2
1 u
1 r
P Y
u hzi 2
r

β̂ 2 β̂ 2 Phzi Y 1

u

*
Z[1]

r

1 r

1
β̂ 2 β̂ 1
r r - r
2 1 Z[1]
2
β̂ 1 β̂ 1 β̂ 1
Figura 5.9: Interpretación geométrica en Rn de la multicolinealidad
Ejercicio 85. Probar (5.29).

Ejercicio 86. Simular una muestra de tamaño n = 50 de una par de variables Z[1] y Z[2],
ambas con media 0 y desviación típica 1, y con coeficiente de correlación ρ = 0.9. A continua-
ción, simular una variable Y que se obtenga mediante la ecuación Y = 2Z[1] − Z[2] + E, con
E ∼ N (0, σ 2 ) y σ 2 = 1. Proceder a estimar los coeficientes de regresión β1 y β2 mediante una
regresión lineal. Simular de nuevo Y en las mismas condiciones y estimar de nuevo los coeficien-
tes, comparando los resultados. Repetir el procedimiento para n = 100; repetirlo igualmente
para σ 2 = 0.5.
Que los estimadores de los coeficientes de regresión estén sometidos a una fuerte variabi-
lidad sólo tiene trascendencia a la hora de optimizar el modelo, ya que da lugar a resultados
no significativos en los tests parciales. El hecho de que una variable permanezca o salga del
modelo en función del resultado del test parcial es pues muy cuestionable por el efecto de la
multicolinealidad. Para compensarlo se han ideado los diferentes algoritmos de selección de
variables que permiten retener algunas variables en detrimento de otras fuertemente correla-
cionadas con las primeras. No obstante, la selección puede depender de pequeños detalles de la
muestra que podrían llegar a ser arbitrarios desde un punto de vista estadístico. Esta posible
88
arbitrariedad no supone un problema si nuestra intención es predecir con la mayor precisión

posible la variable respuesta y mediante la menor cantidad posible de variables explicativas.
Tan sólo puede tener trascendencia si nuestro objetivo es determinar la influencia real de cada
una de las variables explicativas en la respuesta.
En tal caso, podemos optar por ejecutar la regresión respecto a las componentes principales
de las variables explicativas porque, al ser éstas incorreladas, los tests parciales no pueden verse
contaminados por la multicolinealidad.
Figura 5.10: Interpretación geométrica en Rq+1 de la multicolinealidad
Ejercicio 87. Probar que el coeficiente de correlación múltiple R2 de Y respecto a las varia-
bles explicativas Z permanece invariante si reemplazamos estas últimas por sus componentes
principales U. Probar que R2 puede obtenerse como la suma de los coeficientes de correlación
simple al cuadrado entre Y y cada una de las componentes principales (cosa que no sucede en
general con las variables explicativas originales).
Ejercicio 88. ¿Por qué una fuerte variabilidad de los estimadores se asocia a resultados no
significativos en los tests parciales?
89
Una vez estimado el vector η̂ con los coeficientes de regresión respecto de U, debemos
deshacer el cambio teniendo en cuenta (5.5), obteniendo así la estimación
ˆ
β̂ = Gη̂
Si hemos eliminado las últimas componentes principales en los tests parciales, esta nueva es-
timación de β estará sometida a tantas restricciones lineales como componentes eliminadas,
y será sesgada pero con menor varianza que el EIMV β̂. En las condiciones de la simulación
propuesta en el ejercicio 86, el primer eje principal es h(1, 1)0 i. Luego, si se desecha la segunda
componente principal, la ecuación estimada consistirá en multiplicar Z[1] y Z[2] por un mismo
coeficiente.
Desde un punto de vista práctico, distinguimos pues dos posibles circunstancias: que se
eliminen componentes principales en la regresión lineal, lo cual conduce a considerar una ecua-
ción más estable que puede entenderse como una especie compromiso entre las distintas varia-
bles correlacionadas, como en el ejemplo comentado anteriormente; o bien que no se eliminen
componentes principales, lo cual debe entenderse como que la muestra consta de información
suficiente para determinar qué variables poseen influencia real en la respuesta, en cuyo caso
debemos acatar el resultado que aporten los tests parciales.
En la figura 5.10 representa una ecuación (plano) de regresión construida a partir de dos
variables explicativas fuertemente correlacionadas. En ese caso, el plano estimado se apoya
fundamentalmente, por así decirlo, en el primer eje principal, siendo muy sensible a pequeños
cambios en el segundo (linea verde), ante los cuales bascularía fuertemente. Precisamente, el
contraste parcial de la segunda componente principal determina si la posición respecto al segun-
do eje es estable (significativa), es decir, si la segunda componente principal recoge información
suficiente respecto a Y como para otorgar fiabilidad a la ecuación obtenida. En caso contrario,
se ignora y la ecuación se reduce a la recta que determina la primera componente (roja).
Una opción intermedia puede ser la denominada regresión Ridge que, en última instancia,
se traduce (ver Hastie et al. (2008) sec. 3.4.1) en una ponderación del peso de las diferentes
componentes principales en función de la varianza de las mismas en una estimación sesgada de
la media. Concretamente, se define el estimador Ridge de β para λ ≥ 0 mediante
ridge
β̂ λ = argmin{kY − Y − Zbk2 + λ · kbk2 : b ∈ Rq } (5.30)
El caso λ = 0 proporciona el EIMV β̂. Puede probarse a través del cálculo de derivadas parciales
que el estimador Ridge de β verifica, en términos de las puntuaciones factoriales F [1], . . . , F [q]
de Z (tipificada),
q
ridge X dj
Zβ̂ λ = F[j] F[j]0 Y (5.31)
j=1
d j + λ
Es decir, dado que uno de los efectos de la multicolinealidad puede ser la inflación en el tamaño
del estimador β̂, puede compensarse imponiendo una cota al mismo mediante una penalización
(no eliminación) de las componentes principales con baja varianza.
90
Capítulo 6
Análisis de conglomerados
Recibe el nombre de análisis de conglomerados o análisis clúster un conjunto de técnicas

destinadas a agrupar observaciones por similitud. En otras palabras, son métodos de identi-
ficación de variables categóricas latentes. Cada observación consistirá en p valores numéricos
correspondientes a la medición de sendas variables y, por lo tanto, constituirán puntos de Rp .
Ésa es la razón por la que esta técnica haya sido considerada tradicionalmente como parte de la
Estadística Multivariante, aunque actualmente tiende a catalogarse como Minería de Datos, de
ahí que le dediquemos poco espacio. Para obtener una información más detallada se consultar
Hastie et al. (2008), Mardia et al. (1979) y Hair et al. (1999).
90
80
waiting
70
60
50
1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
eruption
Figura 6.1: Datos geyser Old Faithful
En las dos primeras secciones abordaremos un breve estudio de los dos métodos tradicionales
del análisis de conglomerados: el de k-medias y el jerárquico. En la tercera sección introducire-
mos escuetamente el algoritmo de agrupación EM, basado en un modelo de mezclas. Este tipo
de técnica va más allá de la mera agrupación de observaciones pues tiene el ambicioso objeto
de determinar de manera sencilla y precisa la distribución probabilística que las explica.
91
Capítulo 6 Método jerárquico
En todo caso, para hablar de similitud entre observaciones es preciso definirla previamente
en el espacio Rp . La opción más utilizada es considerar la distancia Euclídea, aunque puede
optarse por medidas de similitud alternativas, algunas bastantes parecidas y otras, como la
correlación entre datos que definiremos más tarde, no tanto. Si optamos por cualquiera de las
dos medidas antes mencionadas debemos tener presente que no son invariantes ante cambios
de escala en cualquiera de las variables consideradas, lo cual afecta de manera decisiva a la
agrupación. De ahí que, en la mayoría de las situaciones, el análisis de conglomerados deba
ir precedido de la tipificación de los datos. No sería el caso si optásemos, por ejemplo, por la
distancia de Mahalanobis, dada la matriz de covarianzas muestral, que es sí invariante.
6.1. Método de k-medias

También conocido como quick-cluster, se utiliza para agrupar los datos en un número k de
conglomerados determinado a priori. La elección de k puede basarse en argumentos formales,
como los que se mencionan en la tercera sección, o bien en argumentos gráficos y, por lo
tanto, intuitivos, como los que se desprenden de la figura 6.1, correspondientes a datos del
geyser Olf Fatithful, de Yellowstone, donde parecen apreciarse con cierta claridad dos grandes
conglomerados.
La técnica consiste en aglomerar todos los datos en torno a k puntos (que se denominan
semillas) en función de la proximidad a éstos, según la distancia considerada (SPSS sólo con-
templa la distancia Euclídea). En ocasiones, estas semillas son establecidas de antemano en
función de conocimientos previos, en cuyo caso el método es trivial. Si queremos formar k
conglomerados pero no contamos con semillas, puede procederse de la siguiente forma: se se-
leccionan k datos, bien aleatoriamente o bien los k primeros, que serán las semillas iniciales.
Los datos restantes se irán aglomerando en torno a ellos. No obstante, si la semilla más cercana
a un dato dista del mismo más que que la semilla más cercana a ésta, dicho dato reemplaza
como semilla a la más cercana y usurpa en lo sucesivo, por así decirlo, su conglomerado. Al
final del proceso, se reconstruyen las semillas como centroides de los conglomerados finales y el
procedimiento se repite sucesivamente hasta conseguir suficiente estabilidad en los centroides
finales.
6.2. Método jerárquico

El método consiste en evolucionar en tantos pasos como datos tengamos desde un estado
inicial, en la que cada individuo constituye un conglomerado unitario, hasta un estado final, en
el que todos los individuos forman parte del mismo conglomerado global. En cada paso se unirán
entre sí los dos conglomerados más similares para formar en lo sucesivo un conglomerado común
indisoluble. El proceso de formación de los conglomerados queda registrado, de manera que se
puede analizar el estado más interesante, que será aquél en el que queden patentes grandes
diferencias entre los conglomerados y pequeñas diferencias dentro de los conglomerados. Eso
querrá decir que en todos los pasos anteriores se unieron conglomerados próximos, pero en
el inmediatamente posterior se unen dos conglomerados distantes, lo cual puede detectarse
gráficamente mediante el dendrograma, que deberá interpretarse subjetivamente.
Así pues, el procedimiento jerárquico proporciona una método con una fuerte componente
subjetiva que permite decidir el número k de conglomerados que subyacen en nuestros da-
92
Capítulo 6 Método jerárquico
Figura 6.2: clúster k-medias para Faithful y jerarquizado para irisdata
tos. No obstante, la interpretación de dendogramas y, por lo tanto, la elección de k, sólo es

aconsejable cuando se manejan muestras pequeñas o moderadas. Si la muestra es grande, el
análisis jerárquico sigue siendo una herramienta igualmente válida para agrupar en un número
k preestablecido de conglomerados.
Hemos dicho anteriormente que cada paso consistirá en la fusión de los dos conglomerados
más similares entre sí. Obviamente, la similitud se determinará en virtud de la medida que
hayamos escogido. En ese sentido y tal y como adelantábamos en la introducción, suele optarse
por una de las siguientes opciones, radicalmente distintas, que son la distancia Euclídea y la
correlación entre datos. Esta última se define como sigue: si X0i , X0k ∈ M1×p son dos datos (filas)
de la muestra con p componentes, se define la correlación entre ambos como la correlación
entre sus transposiciones Xi , Xk ∈ Mp×1 . De esta forma, si los puntos obedecen a una misma
tendencia entre las variables tendrán una similitud próxima a 1. Esta medida de similitud,
que no es una distancia, se utiliza a veces para agrupar en conglomerados cuando se aprecian
diferentes tendencias entre los datos. No obstante, este tipo de agrupación puede llevarse a
cabo de una forma más sofisticada mediante el algoritmo EM que veremos a continuación.
En todo caso, debemos tener presente que cualquier medida de similitud entre datos m se
aplica a cada par de puntos en Rp , mientras que los conglomerados son conjuntos (unitarios
o no). Por ello, queda aún pendiente determinar una medida m̃ de similitud entre conjuntos
partiendo de la medida m de similitud entre puntos seleccionada. En ese sentido contamos con
varias opciones. El SPSS utiliza por defecto la vinculación inter-grupos, que consiste en definir
la medida entre dos conglomerados A y B mediante
X
m̃(A, B) = [card(A × B)]−1 m(a, b) (6.1)
a∈A, b∈B
En la figura 6.3 presentamos el dendrogama correspondiente a 25 flores de irisdata aglome-

radas en función de sus cuatro medidas morfológicas. Se han utilizados las opciones que SPSS
ofrece por defecto: distancia Euclídea y vinculación intergrupos.
En un análisis de este tipo hay que tener muy presente que los datos extremos constituyen
conglomerados unitarios hasta fases muy avanzadas del análisis, como es el caso de la flor 66, lo
cual supone uno de los mayores inconvenientes del método jerárquico. Haciendo caso omiso de
la misma, se perfilan, de manera subjetiva, entre dos y tres conglomerados de datos. Si optamos
93
Capítulo 6 Algoritmo EM
por dos, podremos comprobar que el más pequeño está compuesto exclusivamente por flores
tipo setosa, mientras que el más grande está compuesto por flores tipo vesicolor y virgínica.
De lo dicho hasta ahora se deduce que el mayor problema a la hora de formar conglomerados
es determinar de una manera mínimamente objetiva el número k de clústers. No obstante,
podemos optar por diferentes combinaciones entre ambas técnicas para lograr una solución
aproximada: por ejemplo, podemos seleccionar a partir de la muestra original una pequeña
muestra piloto y determinar k a partir del dendrograma de la segunda. También puede invertirse
el orden agrupando primeramente respecto a un número elevado m de semillas, que da lugar
a m centroides finales. Éstos se someten entonces a un análisis jerárquico, de manera que los
grupos correspondientes a centroides próximos se unirán dando lugar a un número menor de
conglomerados homogéneos.
Cluster Dendrogram for Solution HClust.2
40
30
Height
20
10
0
17
9
33
3
22
16
24
29
32
12
13
18
19
25
23
27
20
30
15
26
4
31
10
1
34
2
35
8
14
21
28
5
11
Observation Number in Data Set Z

Method=ward; Distance=euclidian
Figura 6.3: Dendrograma para irisdata
6.3. Algoritmo EM
En la sección anterior destacamos lo conflictivo que resulta determinar el número k de con-
glomerados a configurar a partir de la observación de la muestra. Aunque ya propusimos dos
métodos (bastante subjetivos) para tal fin, existen diversos procedimientos semiautomáticos
basados en principios bastante intuitivos, como el método gráfico del codo y el de Calinsky-
Harabasz. El método Bayesiano que describimos a continuación está basado en un modelo de
mezclas, es decir, en la aproximación de distribucionesde probabilidad p-dimensional mediante
94
combinaciones convexas de k distribuciones p-normales. Las condiciones de partidas son muy

similares a las del modelo de regresión logística, con la salvedad de que la variable cualitativa
I no tiene que ser necesariamente binaria y, además, es latente, es decir, no es observable. En
definitiva, sobre un cierto espacio de probabilidad Bayesiano contamos con una variable alea-
toria no observada I con valores en {1, . . . , k}, siendo k en principio desconocido, y un vector
aleatorio observado Y con valores en Rp . De la primera suponemos que sigue una distribución
multinomial con probabilidades a priori {q1 , . . . , qk } y, de la segunda, que la distribución con-
dicional de Y supuesto que I = j es p-normal con media µj y matriz de covarianzas Σj , para
j = 1, . . . , k. Se denotará por pj la correspondiente función de densidad.
Obviamente, cuanto mayor sea el número de componentes k que integren la mezcla y menos
restricciones lineales impongamos a las respectivas matrices de covarianzas, es decir, cuanto
mayor sea el número de componentes d del parámetro θ del modelo, mayor será el logaritmo de
la función de verosimilitud, L(θ; Y). De ahí que, para evitar sobreajustes, se utilice el criterio
de información bayesiano (BIC, ver Schwarz (1978)) a la hora de valorar la aptitud del modelo
a seleccionar.
d
BIC = L(θ; Y) − log n (6.2)
2
El algoritmo Bayesiano denominado EM (esperanza-maximización), que se describe con mayor
detalle en Hastie et al. (2008), tiene un carácter general, pues proponer una alternativa a la
maximización del logaritmo de la versosimilitud en circunstancias como la mezcla de normales,
en las que la distribución (marginal) de las variables observadas conducen a un logaritmo
de sumas, mientras que la distribución (producto) de las variables observadas y las latentes
condicionada a las primeras conduce a una suma de logaritmos, mucho más manejable desde
el punto de vista del Cálculo Diferencial.
Así pues, dado un valor de k determinado, unas restricciones lineales concretas sobre las
matrices de covarianzas y una estimación inicial θ̂0 del parámetro, que procuraremos que sea
lo más razonable posible, el algoritmo EM permite mejorarla en términos de la verosimilitud
de las observaciones Y = (Y1 , . . . , Yn )0 . Se basa en lo siguiente: dado Y ∈ Mn×p , si L, L̃ y LY
denotan respectivamente los logaritmos de las funciones de verosimilitud de las distribuciones
marginales, producto y condicional dado Y (a posteriori), se verifica para todo θ que

L(θ; Y) = Eθ̂0 L(θ; Y)|Y (6.3)

= Eθ̂0 L̃(θ; Y, I1 , . . . , In )|Y − Eθ̂0 LY (θ; I1 , . . . , In ) (6.4)
El método esperanza-maximización consiste en reemplazar θ̂0 por el parámetro θ̂1 que maximice
la esperanza condicional dado Y expresada en primer término en (6.4). Es decir, se trata de
encontrar mediante las técnicas de Calculo Diferencial habituales los valores de (qj , µj , Σj ), con
j = 1, . . . , k, que maximicen la siguiente expresión
n X
X k
Pθ̂0 [Ii = j|Yi ] · [log qj + log pj (µj , Σj ; Yi )] (6.5)
i=1 j=1
Las probabilidades a posteriori Pθ̂0 [Ii = j|Yi ] se calculan mediante la regla de Bayes según se
indica en (6.7). Dado que, en virtud de la desigualdad de Jensen, el término que figura restando
en (6.4) alcanza su máximo en θ = θ̂0 , se deduce que
L(θ̂1 ; Y) ≥ L(θ̂0 ; Y) (6.6)
95
Si el método se aplica de manera iterativa se puede ir mejorando la estimación hasta alcanzar

un grado de estabilidad prefijado, llegando pues a una solución particular para el valor de k y
las restricciones determinadas. La solución general será aquella que maximize el BIC de entre
todas las soluciones particulares para diferentes valores de k y niveles de restricción respecto a
la matriz de covarianzas. En definitiva, los pasos del algoritmo EM-clúster son los siguientes:
1. Se fijan k y las restricciones particulares sobre las matrices de covarianzas.
2. En la fase s = 0 se efectúa una estimación inicial θ̂0 mediante los parámetros muestrales
q̂j0 , µ̂j0 y Σ̂j0 , j = 1, . . . , k, obtenidos con un algoritmo de k-medias. La forma de estimar
las matrices de covarianzas, tanto en este paso como los sucesivos, dependerá de las
restricciones lineales que les impongamos.
3. En la fase s-ésima, con s ≥ 0, se calcula mediante la regla de Bayes las probabilidades a

posteriori para θ̂s (de esta forma se obtiene la esperanza condicional (6.5)):
q̂js · pj (µ̂js , Σ̂js ; Yi )

Pθ̂s (I = j|Yi ) = Pk (6.7)
m=1 q̂ms · pm (µ̂ms , Σ̂ms ; Yi )
4. Se calcula el parámetro θ̂s+1 mejorado cuyas componentes maximizan (6.5) mediante1

n
1X
q̂j,s+1 = P (I = j|Yi ) (6.8)
n i=1 θ̂s
n
1 X
µ̂j,s+1 = Pθ̂s (I = j|Yi ) · Yi (6.9)
nq̂j,s+1 i=1
n
1 X
Σ̂j,s+1 = Pθ̂s (I = j|Yi ) · (Yi − µ̂j,s+1 )(Yi − µ̂j,s+1 )0 (6.10)
nq̂j,s+1 i=1
5. θ̂s+1 reemplaza a θ̂s y vuelven a aplicarse los pasos 3 y 4 sucesivamente hasta alcanzar
suficiente estabilidad en L(θ̂, Y).
6. Este procedimiento se lleva a cabo con diferentes valores de k y diversos grados de res-
tricción para la matriz de covarianzas, seleccionando el modelo que maximize el BIC. En
dicho modelo, cada observación Yi se asigna al clúster (índice) que maximiza la probabi-
lidad a posteriori (6.7) en la última iteración.
Este método puede ejecutase haciendo uso del paquete mclust del programa R, que consi-
dera 14 tipo de restricciones sobre las matrices de covarianzas. Si, por ejemplo, lo aplicamos a
los datos de Old Faithful, el método proporciona un valor máximo del BIC para k = 3 com-
ponentes, con matrices de covarianzas asociadas a elipses con el mismo volumen, excentricidad
y orientación (EEE). En la parte izquierda de la figura 6.4 se muestran la comparativa entre
los 126 diferentes modelos considerados; en la parte derecha, se muestran los diferentes clústers
con las tres distribuciones 2-normales superpuestas; por último, en la figura 6.5 se muestra la
densidad estimada como mezcla de dichas distribuciones. En ocasiones como ésta y a la vista
1
En este caso mostramos las estimaciones de las matrices de covarianzas que corresponderían al modelo sin
restricciones.
96
del gráfico, puede procederse a agrupar clusters (verde y azul) cuya separación no resulte na-
tural. Es decir, que este método aparentemente automático también puede precisar en última
instancia de decisiones subjetivas, por lo que, en definitiva, no difiere tanto de los comentados
en la sección anterior.
Classification
-2500
90
80
-3000
waiting
BIC
70
-3500
60
EII EVE
VII VEE
EEI VVE
VEI EEV
50
EVI VEV
VVI EVV
-4000
EEE VVV
1 2 3 4 5 6 7 8 9 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
Number of components eruption
Figura 6.4: Comparativa de BIC y EM-clúster par Old Faithful

90
0.015
0.025
0.0
35 4
0.0
80
0.03
waiting
0.02
70
0.01
0.005
0.005
0.015
60
0.025
3
0.0
50
0.02
0.01
1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
eruption
Figura 6.5: Densidad estimada para Old Faithful
6.4. Análisis de conglomerados bietápico

En este última sección apuntaremos algunas indicaciones acerca de otra técnica automática
de selección de conglomerados implementada en el programa SPSS (de hecho, para entender
o, mejor dicho, intuir su funcionamiento, hay que remitirse a día de hoy al propio manual del
programa). Se trata de un método similar aunque más sofisticado que el propuesto al final de la
sección6.2. Supone que todas las variables son independientes, siendo normales tipificadas las
97
numéricas y multinomiales las cualitativas (no obstante, el procedimiento resulta ser bastante
robusto ante la violación de este supuesto).
En un primer paso se generan por similaridad un cierto número de conglomerados que no
puede superar un máximo establecido. Los conglomerados se van formando conforme se leen
los datos en el orden en el que aparecen enel archivo (lo cual puede resultar trascendental). La
medida de similaridad por defecto es el logaritmo de la verosimilitud, calculada a partir de las
p variables numéricas. De los conglomerados formados sólo interesará, en la siguiente fase, sus
centroides (medias en dimensión p). El siguiente paso es un procedimiento de aglomeración de
dichos centroides hasta llegar a un conglomerado (centroide) único. En esta fase se consideran
las proporciones de las categorías correspondientes a las variables cualitativas a la hora de
calcular el logaritmo de la verosimilitud y los centroides. El árbol resultante debe podarse
para obtener una solución óptima en función del BIC y teniendo además en cuenta unas cotas
máximas de complejidad a niveles horizontal y vertical predeterminadas.
Este método posee dos desventajas, al menos aparentes: su complejidad y los fuertes su-
puestos sobre los que descansa teóricamente. Por contra, tiene como ventaja principal, a parte
de su automatismo, que contempla el uso de variables cualitativas.
98
Bibliografía
Anderson, T.W. (1958), “An Introduction to Multivariate Statistical Analysis”, Wiley.

Arnold, S.F. (1981), “The Theory of Linear Models and Multivariate Analysis”, Wiley.
Bilodeau, M. y Brenner, D. (1999), “Theory of Multivariate Statistics”, Springer.
Bravo, J.A., Montanero, J., Calero, R., Roy, T.J., “Identification of sperm moyility characte-
ristics in ejaculates from Ile de France rams”, Animal Reproduction Science, vol 129, 22-29
(2011).
Breiman, L., Friedman, J., Stone, C.J. y Olsen, R.A. (1984) “Classification and Regression
Trees”, Taylor & Francis.
Dillon, W.R. y Goldstein, M. (1984), “Multivariate Analysis. Methods and Aplications”, Wiley.
Dobson, A.J. (1990), “An Introduction to Generalized Linear Models”, Chapman & Hall.
Flury, B. (1997), “A First Course in Multivariate Statistics”, Springer.
Gifi, A. (1990), “Nonlinear Multivariante Analysis”, Wiley.
González-Ramirez, C., Montanero-Fernández, J., Peral-Pacheco, D., “A multifactorial study on
duration of temporary disabilities in Spain”, Archives of Environmental & Occupational Health,
vol 72(6), 328-335 (2017).
Greenacre, M.J. (1984), “Theory and Applications of Correspondence Analysis”, Academic
Press.
Hair, J.F., Anderson, R.E., Tatham, R.L., y Black, C.B. (1999), “Análisis Multivariante”, Pren-
tice Hall.
Hastie, T., Tibshirani, R. y Friedman, J. (2008), “ The Elements of Statistical Learning”,
Springer.
Mardia, K.V., Kent, J.T. y Bibby, J.M. (1979), “ Multivariate Analysis”, Academic Press.
Montanero, J. (2008), “Manual 56: Modelos Lineales”, Servicio de Publicaciones UEx.
http://hdl.handle.net/10662/2443
Montanero, J. (2008), “Manual 59: Análisis Multivariante”, Servicio de Publicaciones UEx.
http://hdl.handle.net/10662/2444
Nogales, A.G. (1998), “Estadística Matemática”, Servicio de publicaciones UEx.
Peña, D. (2010), “Regresión y Diseño de Experimentos”, Alianza editorial.
Pérez-Fernández, M.A., Calvo Magro, E., Montanero-Fernández, J., Oyola-Velasco, J.A. “Seed
Germination in response to chemicals: Effect of nitrogen and pH in media”, Journal of Envi-
99
ronmental Biology, vol 27(1), 13-20 (2006).

Rencher, A.C. (1995), “Methods of Multivariate Analysis”, Wiley.
Rodríguez-Mansilla, J., González López-Arza, M.V,. Varela-Donoso, E., Montanero-Fernández,
J., González-Sánchez, B., Garrido-Ardila, E.M., “The effects of ear acupressure massage therapy
and no therapy on symptons of dementia: a randomized controlled trial”, Clinical Rehabilitation,
vol 29(7) 683-693 (2015).
Silverman, B. W. (1986), “Density Estimation for Statistics and Data Analysis”, Chapman &
Hall.
Schwarz, G. (1978), “Estimating de dimension of a model”, The Annals os Statistics, 6(2),
461-464.
Uriel, E. y Aldás, J. (2005), “Análisis Multivariante Aplicado”, Thomson.
100
Índice alfabético
árbol de decisión, 64 correlación entre datos, 93

índice de Gini, 67 correlaciones parciales, 84
covarianza, 10
algoritmo CRT, 65 covarianzas parciales, 14
algoritmo EM, 94 criterio de información Bayesiano BIC, 95
análisis clúster k medias, 92 curvas principales, 76
análisis clúster jerárquico, 92
análisis de componentes principales PCA, 73 dendrograma, 92
análisis de conglomerados bietápico, 97 desigualdad de Bonferroni, 39
análisis de conglomerados o de clusters, 91 desigualdad de Cauchy-Schwarz, 79
análisis de correspondencias, 85 diagonalización de una matriz simétrica, 17,
análisis de la varianza o anova, 27 19, 74
análisis de medidas repetidas, 51 diseño completamente aleatorizado, 25
análisis de perfiles, 49 diseño completamente aleatorizado multivarian-
análisis discriminante cuadrático de Fisher, 57, te, 33, 42, 56
60 distancia χ2 , 86
análisis discriminante lineal LDA de Fisher, distancia dn,p , 15
56, 60 distancia de Mahalanobis, 16, 20, 22, 56, 58,
análisis factorial, 84 92
autovalor, 16, 17, 35, 42, 44, 46, 60 distancia en L2 , 8
autovector, 16, 17, 44, 46, 73 distancia Euclídea, 8, 22, 60, 92, 93
distribución χ2 , 23, 32, 33, 36, 39, 43, 45
base ortonormal, 17, 73 distribución F -Snedecor, 23, 28, 33
distribución t-Student, 24
centroide, 92 distribución T 2 -Hotelling, 32, 35, 38, 40
coeficiente de correlación lineal , 10 distribución condicional, 21, 95
coeficiente de correlación lineal múltiple, 11, distribución de Wishart, 32
12, 28, 30, 46 distribución multinomial, 95
coeficiente de correlación parcial, 30 distribución normal matricial, 31
coeficiente de correlación simple, 12 distribución normal multivariante, 15, 19
coeficiente de Kaiser-Meyer-Olkin KMO, 84 distribución normal multivariante esférica, 22,
coeficientes de correlación canónica, 14, 46, 47, 24
59 distribución a posteriori, 54, 57, 62, 96
coeficientes de correlación parcial, 14 distribución a priori, 54–57, 95
componentes principales, 73, 74, 76, 80, 90 distribuciones de Wilks, Lawley-Hotelling, Roy
comunalidades, 80, 82 y Pillay, 33
contraste de hipótesis, 27, 34
contraste total en regresión, 45 ecuación de regresión lineal, 12
contrastes para la media, 27 ejes canónicos, 46
contrastes parciales en regresión, 47, 89 ejes discriminantes, 42, 43, 47, 59
101
ejes factoriales, 77 matriz de componentes, 81
ejes principales, 78, 81 matriz de correlaciones, 12, 77, 81
elipsoide de confianza, 37 matriz de covarianzas, 12, 20, 22
elipsoides, 22 matriz idempotente, 16
escalamiento multidimensional, 85 matriz ortogonal, 16, 23
espacio L2 , 7, 65 matriz positiva, 16
espacio Euclídeo, 7 media, 9
esperanza, 9 medidas de similitud, 93
esperanza condicional, 15, 21 minería de datos, 64
esperanza-maximización EM, 95 modelo lineal multivariante, 37
estimación, 26, 34 modelo lineal normal, 22, 24
estimación sesgada, 90 modelo lineal normal multivariante, 31
estimador de máxima verosimilitud EMV, 26, muestra aleatoria, 7
34 muestra aleatoria simple de una distribución
estimador insesgado de mínima varianza EIMV, normal, 25, 33, 37
26, 34, 37, 40 muestras independientes de distribuciones nor-
estimador Ridge, 90 males con idéntica varianza, 25
estrategia de Bayes, 53, 55 muestras independientes de distribuciones nor-
estrategia LDA, 57 males con idénticas matrices de cova-
estrategia LDA Bayesiana, 57, 62 rianzas, 33, 40
estrategia minimax, 55 multicolinealidad, 87
estrategia no aleatoria, 53
nodo filial, 65
factor, 84 nodo paranetal, 65
factor de inflación de la varianza FIV, 30, 88 nodos terminales, 66
factores latentes, 84 norma en L2 , 8
familia completa maximal, 55 norma Euclídea, 8, 23
función característica, 19
función de densidad, 20, 21, 55 odds ratio, 63
función logística, 57, 63 ortogonalidad en Rn , 8
ortogonalidad en L2 , 8
grados de libertad, 23, 28
parametrización de un modelo, 26
heterocedasticidad, 27 preorden en estrategias, 54
principio de Invarianza, 23, 24, 27, 35
importancia, 68
principio de Máxima Verosimilitud, 55
incorrelación, 8, 10, 14, 20
independencia, 8, 14, 20 problema de clasificación, 53
inercia, 86 producto escalar en Rn , 8
producto interior en L2 , 7
lema fundamental de Neyman-Pearson, 28 proporción de varianza total explicada, 80
leyes de los Grandes Números, 9, 38, 55 proyección ortogonal, 8, 16
puntuaciones discriminantes, 42, 44, 59
método de Mínimos Cuadrados, 8 puntuaciones factoriales, 79, 81
método del eje principal, 85
método del vecino más próximo KNN, 63 redes neuronales, 69
método Lambda de Wilks, 48, 59 regla de Bayes, 54, 96
método núcleo, 55 regresión lineal múltiple, 25
manova, 34, 42, 47, 50, 53 regresión lineal multivariante, 34, 45
regresión logística, 61
riesgo, 54, 58
rotación de ejes factoriales, 82
rotación oblicua, 84, 85
rotación varimax, 82
semilla, 92
sobreajuste, 66, 95
subespacios h1i y h1n i, 9
suficiencia, 27
tabla de contingencia, 85
teoría de la Decisión, 53
teorema Central del Límite, 38
teorema de Lehmann-Scheffé, 26
teorema de los multiplicadores finitos de La-
grange TMFL, 37, 41
teorema Fundamental del Cálculo, 55
teorema Límite Central, 9
test F , 28, 35, 43
test de Friedman, 51
test de Hosmer-Lemeshov, 63
test de Lawley-Hotelling, 35, 47
test de Pillay, 35, 47
test de razón de verosimilitudes TRV, 28, 35,
38, 40
test de Roy, 35, 47
test de Student, 28, 35, 40
test de Wilks, 35, 43, 45, 47
test UMP-invariante, 28, 38, 40
tipificación, 77, 92
tolerancia, 30
variabilidad, 10
variabilidad parcial, 13
variables dummys o ficticias, 26, 30, 47, 48, 59
varianza, 9
varianza parcial, 13
varianza total, 10, 16, 74, 75
varianzas específicas, 80
vector aleatorio, 7, 19
vector de componentes de la variable, 78, 80
vinculación inter-grupos, 93

Manual de Estadística Multivariante

Cargado por

Copyright:

Formatos disponibles

Manual de Estadística Multivariante

Cargado por

Información del documento

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Manual de Estadística Multivariante

Cargado por

Copyright:

Formatos disponibles

MANUAL ABREVIADO DE

Jesús Montanero Fernández

Badajoz, diciembre de 2018

2. Modelo lineal multivariante 19

3. Aplicaciones del modelo 37

4.3.4. Validación de estrategias: . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Índice alfabético 101

Proyección ortogonal: La noción de perpendicularidad se relacionará bajo ciertas condi-

1.2. Principales parámetros probabilísticos

Media: Primeramente definimos la media de una variable aleatoria X como su esperanza, es

Varianza: Dada una variable aleatoria X ∈ L2 de media µ, se define su varianza mediante

var[X] = E[(X − µ)2 ] (1.6)

denotándose también por la letra σ 2 . La raíz cuadrada positiva de la varianza se denomina

var[X] = E[X 2 ] − E[X]2 (1.7)

Ejercicio 2. Probar (1.7).

Ejercicio 3. Probar que µ = Ph1i X.

var[X] = d2 (X, E[X]) (1.8)

Varianza total: Si X es un vector aleatorio p-dimensional de componentes X[1], . . . , X[p], se

Este parámetro puede interpretarse en términos de la distancia d2[p] definida en el espacio de

d2[p] (X, Y ) = EP kX − Y k2Rp

varT [X] = d2[p] (X, E[X]) (1.11)

Covarianza: Dado un vector aleatorio p-dimensional X, se define la covarianza entre dos

Coeficiente de correlación: La desigualdad anterior invita a definir el denominado coefi-

Cov[X] = E[(X − µ)(X − µ)0 ] (1.14)

E[X̃] = A0 E[X] + b, Cov[X̃] = A0 Cov[X]A (1.15)

1.3. Principales parámetros muestrales

Covarianza muestral: Dados X, Y ∈ Rn , se denota

S = n−1 (X − X)0 (Y − Y) ∈ Mp×p (1.20)

Coeficientes de correlación simple y múltiple: Dados X1 ∈ Rn y X2 ∈ Mn×q se define

1.4. Regresión lineal

mı́n kX1 − E[X1 ] − b0 (X2 − E[X2 ])k2 : b ∈ Rp2

Ecuación de regresión lineal: La solución se obtiene proyectando ortogonalmente el vector

X1 − E[X1 ] − β 0 (X2 − E[X2 ]) ⊥ hX2 − E[X2 ]i (1.24)

Figura 1.2: Ecuación de regresión lineal

En el caso p2 = 1 (regresión simple), el problema se reduce a encontrar el escalar β ∈ R tal

Ejercicio 9. Probar que, en general, la ortogonalidad (1.24) se alcanza con el vector

Por otra parte, si se define

se verifica que X1 = α + β 0 X2 + E, donde E, denominado variabilidad parcial de X1 dado X2 ,

var[E] = σ12 − Σ12 Σ−1

Figura 1.3: Descomposición de la varianza

Coeficientes de correlación parcial Razonando de manera análoga en el caso general

Incorrelación e independencia: Se dice que X1 y X2 son incorreladas cuando Σ12 = 0,

explicar parte alguna de la variabilidad de X1 como función lineal de la variabilidad de X2 .

X1 y X2 incorreladas ⇔ hX1 − E[X1 ]i ⊥ hX2 − E[X2 ]i (1.31)

X1 y X2 independientes ⇔ M(X1 ) ⊥ M(X2 ) (1.32)

1.5. Nociones básicas de Álgebra Lineal

hA, Bi = tr(A0 B) (1.33)

d2n,p (A, B) = n−1 tr[(A − B)0 (A − B)] (1.35)

Matriz positiva: En el conjunto de matrices cuadradas m × m, podemos definir el siguiente

Distancia de Mahalanobis: Dada una matriz A ∈ Mm×m simétrica y positiva podemos

Se trata de una generalización de la métrica Euclídea, que se obtendría en el caso A = Id.

Matriz idempotente: Se dice que una matriz A ∈ Mm×m es idempotente cuando A2 = A.

Autovalores y autovectores: Dada una matriz A ∈ Mm×m , se dice que δ ∈ R es un

(i) Si A es simétrica, su rango coincide con el número de autovalores no nulos.

(v) La traza de una matriz simétrica es la suma de sus autovalores y el determinante, el

(vi) La inversa de una matriz simétrica positiva también es positiva.

(vii) La matriz PV de la proyección ortogonal sobre el subespacio lineal V ⊂ Rn consta de

Cuadro 1.1: Matriz de datos muestra tamaño n = 38 y dimensión p = 8

Modelo lineal multivariante

2.1. Distribución normal multivariante

PE2 Y 6 PE1 ⊕E2 Y