0% encontró este documento útil (0 votos)
138 vistas123 páginas

Muñoz Armayones, Sandra TFG PDF

Descargar como pdf o txt
Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1/ 123

FACULTAD DE MATEMÁTICAS

DEPARTAMENTO DE ESTADÍSTICA E INVESTIGACIÓN OPERATIVA

Trabajo Fin de Grado

Técnicas Multivariantes para el Análisis de Datos Ómicos

Sandra Muñoz Armayones

Dirigido por:
Dña. Inmaculada Barranco Chamorro
2016
Resumen

El constante aumento en la generación de datos ómicos y el desarrollo de tecnologías que


permiten su análisis han hecho que el interés por estudiar simultáneamente datos procedentes
de distintas técnicas ómicas sea cada vez mayor, con el propósito de conocer las relaciones
subyacentes entre ellos.

Para ello se requieren herramientas matemáticas que puedan adaptarse al análisis de una gran
cantidad de datos, reduciendo a su vez la complejidad de los mismos y facilitando así su inter-
pretación. Esto nos lleva a considerar métodos de proyección, que serán descritos en el trabajo
además de aplicarse en diferentes conjuntos de datos reales.

En la introducción al trabajo, se han presentado las principales técnicas ómicas así como la tec-
nología de microarrays, además de explicar qué se entiende por integración de datos. Planteamos
la necesidad de utilizar técnicas de Análisis Multivariante para analizar el tipo de datos que nos
ocupa.

En el primer capítulo del trabajo se explica con detalle el Análisis de Componentes Principales
(ACP), una técnica capaz de crear un pequeño conjunto de variables que resuman la información
de las originales y que permitan posteriores análisis más profundos de los datos.
Comenzamos introduciendo la notación, además de definir componentes principales y compo-
nentes principales muestrales. Se presentan ejemplos para explicar el círculo de correlación y la
utilidad del ACP para solucionar problemas de colinealidad en la regresión lineal múltiple. Se
explican tres criterios para escoger las CP significativas.
Por último, se ha realizado una aplicación a un conjunto de datos de expresión génica, en la que
utilizamos ACP como una técnica exploratoria de datos.

En el segundo capítulo vemos la Descomposición en Valores Singulares (DVS). Este será un


paso intermedio en varias de las técnicas estadísticas expuestas a lo largo del trabajo y nos
permite descomponer matrices rectangulares como producto de otras. Se detallan propiedades de
la descomposición, la utilidad de la técnica para aproximar matrices y su representación gráfica
como un biplot.

En el tercer capítulo nos centramos en el Análisis de Correspondencias (AC). Esta técnica es


aplicable tan sólo a variables categóricas y define unos índices. A estos índices se les denomina
coordenadas principales y estándar y se obtienen a partir de la descomposición del estadístico

3
4

chi-cuadrado, χ2 . Se realiza una aplicación a datos de expresión génica utilizando el paquete


made4 de R.

El cuarto capítulo trata sobre Análisis de Coinercia (ACoi). Se presenta un coeficiente que
permite medir la correlación entre conjuntos de datos donde se trabaja con las mismas muestras.
Con el propósito de realizar un análisis integrado, consideramos dos conjuntos de datos de
expresión génica. Se realiza un Análisis de Correspondencias como paso previo al Análisis de
Coinercia, el cuál nos permite cuantificar y visualizar la relación existente entre ambos conjuntos.

La última técnica que trataremos se describe en el quinto capítulo y se conoce como Análisis de
Correlación Canónica, que explora las relaciones de dependencia entre conjuntos de variables.Se
ha resuelto detalladamente el problema que supone encontrar los dos primeros vectores de corre-
lación canónica haciendo uso de la Descomposición en Valores Singulares. A partir de esto se
han definido las variables de correlación canónica y se han planteado contrastes de significación
para elegir las más relevantes.
Se ilustra con un estudio de datos de expresión génica en grupos de ratones sometidos a distintas
dietas.

Para finalizar, se han detallado todos los paquetes de R empleados, describiendo cada paquete
así como todas las funciones y argumentos que hemos usado.
Abstract

Constant growth in omics data generation and development of technology that allows this data
analysis has resulted in an increasing interest in studying different kinds of omics techniques, in
order to know the mutual interactions between these data sets.

Mathematical tools to analyse large data sets are required. They must be able to reduce comple-
xity and make the interpretation of these data easier. So we consider projection methods which
will be described along this work. As illustrations, applications to different real data sets are
included.

In the Introduction, the main omics techniques are presented. Microarray technology and data in-
tegration are explained. The need for using multivariate analysis techniques is also contemplated.

In Chapter 1, we focus on Principal Components Analysis (PCA). This technique is able to create
a little set of variables which summarize information and permit deeper analysis of data. We
introduce the appropriate notation, and define principal and sample principal components. Exam-
ples to explain the correlation circle are given. We show how useful this method can be to deal
with highly correlated variables in linear regression. Three different options to choose important
components are described. Finally, we apply PCA to explore microarray gene expression data.

In Chapter 2, we study Singular Value Decomposition (SVD). This is a common tool in multiva-
riate analysis used to decompose a rectangular matrix as product of other matrices. We highlight
properties and biplot representation of this technique.

In Chapter 3, Correspondence Analysis (CA) is presented as a technique applicable to categorical


variables. Indexes, called principal and standard coordinates, are obtained from the decomposi-
tion of a χ2 statistic. An application is carried out by using made4 package of R.

Chapter 4 is devoted to Coinertia Analysis (CIA). This technique allows us to obtain a coefficient
which explains the existing correlation between two data sets containing the same samples. An
application in which we perform an integrated analysis is given. Quantification and visualization
of the relationships between the two data sets, under consideration, is possible thanks to Corres-
pondence Analysis, which is a previous step in order to apply CIA.

In Chapter 5, Canonical Correlation Analysis (CCA) is proposed as a technique to explore

5
6

dependence between variable sets. A method to find the two first canonical correlation vectors
is studied in detail. Canonical correlation variables are presented, and significance tests are
proposed to choose the most relevant ones. CCA is applied to a nutritional study in mice.

Finally, an Appendix is given with the R packages and functions used in this work.
Índice general

1. Análisis de Componentes Principales 17


1.1. Notación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2. Componentes Principales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.2.1. Transformaciones lineales de las variables originales . . . . . . . . . . 19
1.3. Cálculo y propiedades de las Componentes Principales Muestrales . . . . . . . 24
1.4. Componentes principales para datos estandarizados . . . . . . . . . . . . . . . 26
1.4.1. Círculo de correlación . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.5. Componentes Principales de una matriz de datos bivariante . . . . . . . . . . . 33
1.6. Análisis de Componentes Principales y Regresión Lineal . . . . . . . . . . . . 34
1.6.1. Aplicación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1.7. Escogiendo el número de componentes . . . . . . . . . . . . . . . . . . . . . . 42
1.8. Aplicación práctica del ACP . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
1.8.1. Exploración de datos . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
1.8.2. Detección de efecto lote en experimentos ómicos . . . . . . . . . . . . 50

2. Descomposición en Valores Singulares 56


2.1. Descomposición . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.1.1. Propiedades de la descomposición . . . . . . . . . . . . . . . . . . . . 57
2.1.2. Aproximación matricial . . . . . . . . . . . . . . . . . . . . . . . . . 59
2.2. Biplots y matriz de aproximación . . . . . . . . . . . . . . . . . . . . . . . . . 60

3. Análisis de Correspondencias 61
3.1. El método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.2. Descomposición Chi-Cuadrado . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2.1. Descomposición y propiedades . . . . . . . . . . . . . . . . . . . . . . 65
3.3. Coordenadas Principales y Coordenadas Standard . . . . . . . . . . . . . . . . 66
3.4. Análisis de Correspondencias en la Práctica . . . . . . . . . . . . . . . . . . . 68
3.5. Aplicación práctica del AC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.5.1. Exploración de los datos . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5.2. Análisis de Correspondencias . . . . . . . . . . . . . . . . . . . . . . 74

7
8 ÍNDICE GENERAL

4. Análisis de Coinercia 78
4.1. El Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.2. Aplicación práctica del ACoi . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.2.1. Análisis de Correspondencias . . . . . . . . . . . . . . . . . . . . . . 82
4.2.2. Análisis de Coinercia . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

5. Análisis de Correlación Canónica 95


5.1. El método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.1.1. Primer par de vectores de correlación canónica . . . . . . . . . . . . . 97
5.1.2. Correlación Canónica y Descomposición Singular . . . . . . . . . . . 99
5.2. Variables de correlación canónica . . . . . . . . . . . . . . . . . . . . . . . . 100
5.3. Contrastes de significación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.4. Aplicación práctica del ACC . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.4.1. ACC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

6. Paquetes de R 115
Introducción y Contexto

De la biología a la ómica
Se conoce como Ómica al estudio de las diferentes componentes que participan y/o regulan
procesos biológicos complejos.

Las técnicas ómicas se basan en el análisis de una gran cantidad de datos. La principales técnicas
ómicas son:
Genómica: Se encarga de secuenciar, montar y analizar la estructura y función de los
genomas; siendo un genoma el conjunto completo de ADN dentro de una sola célula en
un organismo.

Transcriptómica: Mide los niveles de ARN de una población celular.

Proteómica: Estudia la estructura y función de las proteínas.

Inómica: Estudia el perfilado de elementos, las interacciones entre ellos y su regulación


bioquímica.

Metabolómica: Mide niveles de metabolitos, siendo éstas las moléculas que intervienen en
las reacciones de las células de un organismo.

Fenómica. Evalúa rasgos. Relaciona genética, epigenética (cambios en el fenotipo o


expresión génica) y factores ambientales.
Nos centraremos en la tecnología de microarrays (o chips de ADN) para el estudio de la trans-
criptómica.

Esta tecnología se basa en el Dogma Central de la Biología Molecular, una hipótesis formulada
en 1953 para referirse a los procesos llevados a cabo en la transmisión y expresión de la herencia
genética. Estos procesos se conocen como replicación (ADN se duplica), transcripción (ADN se
sintetiza) y por último; síntesis de proteínas, que se lleva a cabo en los ribosomas a partir del
ADN transcrito.

Describiremos el procedimiento para obtener datos mediante este tipo de tecnología con detalle.

9
10 ÍNDICE GENERAL

Podemos describir el proceso analítico de los microarrays ayudándonos del siguiente esquema:

Como se observa, tenemos que partir de una cuestión biológica y orientar el experimento para
encontrar respuesta. Una vez que se ha diseñado el experimento, podemos proceder a medir los
niveles de ARN.

El ácido ribonucleico que contiene la información del ADN original y la traspasa al ribosoma,
donde tienen lugar la síntesis de las proteínas, recibe el nombre de ARNm o ARN mensajero.

Se representan en un chip las secuencias biológicas, de manera que se pueda cuantificar el nivel
de transcripción en una matriz numérica. En cada una de las celdas del chip se almacenan copias
de un segmento de ARNm, de manera que todas las celdas tienen el mismo número de copias
ÍNDICE GENERAL 11

pero en cada una de ellas aparece una secuencia distinta. El número de ARNm de cada segmento
será variable.

Este ARNm extraído de las muestras se identifica con unos marcadores fluorescentes. Luego se
hibrida, de manera que cadenas de ARNm complementarias se combinan.
Tras esto se lava el microarray para eliminar aquellas muestras que no se hayan hibridado y se
escanea. En el escaneo se ilumina con un láser que revela el color del marcador fluorescente en
función de la cantidad de ARNm que haya en cada celda del chip. El resultado final se denomina
nivel de expresión.

Posteriormente, se lleva a cabo un análisis de imagen; donde se convierte la cantidad de se-


cuencias hibridadas (medidas por fluorescencia) en una intensidad de luz (número), para luego
normalizar los datos. Esto forma parte de lo que se conoce como preprocesamiento.

En el proceso, existen varios elementos que pueden interferir en la medida de la expresión génica.
Estos son:

La fluorescencia, debida al rendimiento del escáner o a distinta eficiencia de las etiquetas.

La densidad de impresión.

El experimento biológico, por la pureza de las muestras o la manipulación de las mismas.

Por ello se hace el preprocesamiento, donde se busca disminuir tendencias erróneas, así como la
varianza de los datos. Dentro del preprocesamiento se llevan a cabo varios pasos: cuantificación
de la imagen, exploración de los datos, corrección de fondo; normalización; sumarización y, por
último, determinación de la calidad.

En esta etapa del análisis se pasa de microarrays a una matriz numérica analizable, conocida
como matriz de expresión; donde las columnas son las condiciones y las filas los genes. Esta
matriz cuantifica la abundancia de la variable i en la muestra j. Los valores de expresión para
cada muestra son independientes, mientras que para cada variable no.

De cada muestra se obtiene información que la describe. Este tipo de información se conoce
como metadatos o variables fenotípicas.

Una vez que hemos llegado a esta etapa, sólo queda hacer el análisis estadístico correspondiente
en función del objetivo que persigamos y luego interpretarlo.
12 ÍNDICE GENERAL

El análisis de alto nivel que se realiza antes de la verificación biológica depende del estudio
que se desee llevar a cabo. Existen tres tipos: comparación de clases, descubrimiento de clase o
predicción de clase.
Es por ello por lo que se usan diferentes técnicas o bien de estadística descriptiva, o bien de
estadística inferencial.

Los test de hipótesis son útiles para los estudios de comparación de clases, mientras que la
parte de la estadística descriptiva que define distancias entre genes y muestras se usa para el
descubrimiento de clases. El modelo lineal discriminante y las máquinas de vector soporte son
efectivas para la predicción de clase.

En los problemas de “Análisis de grupos de genes” se mide la asociación entre grupos de varia-
bles y una variable fenotípica de interés.

Los resultados de la mayoría de los análisis ómicos consisten en una o más listas de identificado-
res. Estos identificadores son nombres o números únicos que necesitan estar:

asociados con una o más bases de datos (anotación)

interpretados en un contexto biológico (anotación funcional)

Se conoce como “GO” (Gene Ontology) a la base de datos de anotaciones creada para proveer
un vocabulario apropiado para describir genes y atributos del producto génico, en otras palabras,
ÍNDICE GENERAL 13

del material bioquímico que resulta de la expresión de un gen.

Esta base de datos está organizada por función molecular, proceso biológico y componente
celular.

Análisis interactivo e integración de datos


Existe una verdadera necesidad de analizar e integrar datos para cualquier investigación biológica.

La primera década del siglo XXI fue significativa, pues comenzaron a generarse una gran canti-
dad de datos biológicos debido a la creciente disponibilidad de secuencias de genomas además
del desarrollo de tecnologías de alto rendimiento.

Todas estas tecnologías se caracterizan porque miden un sólo tipo de información en muchas
variables, simultáneamente. En muchos casos, éstas pertenecen a un tipo específico de tecnología
ómica; como por ejemplo la transcriptómica o la proteómica.

En un principio, cada una de las aproximaciones anteriores se utilizaba por separado ya que
las tecnologías eran muy costosas al estar poco avanzadas. No obstante, a medida que fueron
mejorando y siendo más asequibles iba creciendo el interés por considerar simultáneamente
datos de distinto tipo, a fin de entender mejor los procesos biológicos que intervienen en un
mismo problema.

Mientras más datos para trabajar y más posibilidades de obtenerlos había, mayor era el interés
de orientar esta parte de la biología hacia una modelización y análisis de organismos como un
conjunto.

Para el mismo estudio se emplearían datos de distinta clase (expresión, proteínas, metabolitos,...),
y por tanto, se necesitaban métodos y herramientas para analizarlos de manera conjunta.
Es por ello por lo que se han desarrollado muchos métodos en los últimos años que permiten
alcanzar nuestro objetivo. Existen, o bien métodos basados en machine-learning (aprendizaje
automático), redes bayesianas, máquinas de vector soporte y métodos basados en gráficos; o
bien métodos de estadística multivariante, que lleva usándose mucho tiempo para combinar y
visualizar datos multivariantes.

El objetivo de la integración de datos como la entendemos es combinar diferentes recursos


de datos que contribuyan a una mejor comprensión del fenómeno general de estudio. Pero el
concepto “integración de datos” no tiene siempre el mismo significado.

Por una parte, el término puede usarse para describir herramientas y métodos que combinan y
analizan múltiples recursos de datos; un significado puramente informático. Esto conduce a que
14 ÍNDICE GENERAL

el usuario en ocasiones no sepa en concreto de qué base de datos procede la información con que
trabaja, pues esta está dispersa entre distintas bases.

Por otra parte, se pueden combinar estudios relacionados a fin de obtener una conclusión de
más peso. Esto se describiría como una aproximación mixta que se centra en combinar estudios
uniendo los datos (meta-análisis) a la vez que los reprocesa y reanaliza (integración de los datos).

Sin embargo, pueden considerarse distintos tipos de datos (medidos o no en los mismos indivi-
duos) y tratar de combinarlos de manera que ayuden a interpretar los procesos biológicos que
interfieren. Éste es el tipo de estudio de datos en el que nos centraremos.

Hay muchos tipos de análisis de datos ómicos en función del problema biológico (mapa genético,
clasificación, extracción de características,...) o estadístico, el tipo de datos (similar o heterogé-
neo) y la etapa de integración (comparación de datos, normalización, filtrado de calidad,...).

El método que usemos para nuestro análisis tiene que:

Reducir dimensión eficientemente

Representar simultáneamente muestras y variables de cada conjunto de datos

Reducir el problema que surge al haber muchas variables y pocos individuos

Integrar datos suplementarios en un espacio común con los datos originales.

Técnicas Clásicas para el Análisis de Datos Ómicos


Los métodos clásicos para el análisis de datos multivariantes son Regresión Múltiple, que deter-
mina si existe relación de dependencia entre dos ó más variables; Análisis Discriminante, cuyo
objetivo es describir diferencias significativas entre grupos sobre los que se observan p variables;
y ANOVA o Análisis de la Varianza, que compara las medias de dos o más variables.

Estos métodos clásicos se aplican cuando el número de individuos es mucho mayor que el de
variables.
Se asume:

Variables independientes

Más observaciones que variables

Normalidad multivariante

Existencia de una variable dependiente

Pocos datos perdidos


ÍNDICE GENERAL 15

En los problemas que nos plantearemos no se verifican la mayoría de asunciones necesarias


para aplicar técnicas clásicas. En conjuntos de datos ómicos, por lo general; se miden muchas
variables simultáneamente y hay muy pocas muestras analizadas, por lo que se necesitan otro
tipo de técnicas para su análisis.

La mejor forma de solucionar el problema es reduciendo la dimensión de los datos mediante


métodos de proyección, donde se evalúan todas las variables de forma conjunta proporcionando
así modelos más estables, mostrando relaciones subyacentes (denominadas variables latentes) y
además asegurando poca pérdida de información.

Estos métodos son:

Análisis de Componentes Principales (ACP)

Descomposición en Valores Singulares (DVS)

Análisis de Correspondencias (AC)

Análisis de Coinercia (ACoi)

Análisis de Correlación Canónica (ACC)


Capítulo 1

Análisis de Componentes Principales

El objetivo del Análisis de Componentes Principales (ACP) es reducir la dimensionalidad de un


conjunto de datos multivariante mediante la sustitución de las variables originales por transfor-
maciones lineales de las mismas, obteniendo componentes incorreladas y con poca pérdida de
información.

De este modo, el método permite hacer una representación gráfica que podamos interpretar y
también adecuar los datos para realizar cualquier otro tipo de análisis.

Esta técnica, iniciada por Pearson en 1901 y desarrollada por Hotelling en 1933, no requiere
supuesto de normalidad; sólo que exista el vector de medias y la matriz de varianzas y covarianzas
del vector aleatorio objeto de estudio.

Mostraremos que el ACP es, además, útil para visualizar clusters en los datos e identificar outliers.

Expresaremos la variación en el conjunto de variables originales X = (X1 , . . . , Xp )0 en términos


de un conjunto de variables incorreladas Y = (Y1 , . . . , Yp )0 , en orden decreciente de “importan-
cia”.

Las primeras componentes explicarán una gran proporción de la variabilidad total, lo que nos
será útil para resumir los datos en una dimensión menor a la inicial.

1.1. Notación
Sea un vector aleatorio p-dimensional X = (X1 , . . . , Xp )0

Se supone que existe µ, vector de esperanzas de X:

E(X) = (µ1 , . . . , µp )0 = µ

17
18 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Sea Σ la matriz de varianzas y covarianzas de las variables originales:


 
σ11 . . . σ1p
Σ = Cov(X) =  ... . . . ... 
 
σp1 . . . σpp
con
σij = Cov(Xi , Xj ) = E((Xi − µi )(Xj − µj )0 )
σii = Cov(Xi , Xi ) = V ar(Xi ) = E(Xi2 ) − µ2i
Entonces, X ∼ (µ, Σ)

Consideramos una m.a.s. de X de tamaño n. Construimos la matriz de datos muestrales siguiente:


 
x11 . . . x1p
X =  ... . . . ...  = (x1 , . . . , xn )0 = (x(1) , . . . , x(p) )
 
xn1 . . . xnp
donde
xi = (xi1 , . . . , xip )0 es un vector p × 1 con el valor de cada variable para el individuo
i-ésimo.
x(j) = (x1j , . . . , xnj )0 es un vector n×1 que contiene la muestra de tamaño n de la variable
Xj .
Sea Σ
b la matriz de varianzas y covarianzas muestrales de X
n
1 X
Σ
b = (xi − x)(xi − x)0
n − 1 i=1

donde x es el vector de medias muestrales.

Para j, s ∈ {1, . . . , p} se denota la varianza muestral de x(j) como σ̂x2(j)


n
1 X
σ̂x2(j) = σ̂jj = (xij − xj )2
n − 1 i=1
La covarianza muestral entre x(j) y x(s) se denota por σ̂x(j) ,x(s)
n
1 X
σ̂x(j) ,x(s) = σ̂js = (xij − xj )(xis − xs )
n − 1 i=1

Por último, denotamos por λ1 , λ2 , ..., λp los autovalores de la matriz Σ y por e1 , e2 , . . . , ep sus
autovectores.
Del mismo modo, denotamos sus análogos muestrales como λ̂1 , λ̂2 , ..., λ̂p y ê1 , ê2 , . . . , êp ; res-
pectivamente.
1.2. COMPONENTES PRINCIPALES 19

1.2. Componentes Principales


La primera componente principal se define como la combinación lineal de las variables originales
con varianza máxima de entre todas las combinaciones posibles.

La segunda componente principal se define como la combinación que maximiza la varianza


acumulada y que; además, es incorrelada con la primera componente principal. Se sigue este
mismo razonamiento para definir el resto de componentes.

A continuación introducimos notación adecuada para tratar con transformaciones lineales de


las variables. Así mismo, se recogen propiedades de dichas transformaciones lineales que se
utilizarán en las secciones siguientes.

1.2.1. Transformaciones lineales de las variables originales


Sea tj = (tj1 , . . . , tjp )0 ∈ Rp , con j ∈ {1, . . . , p}.

Consideremos la siguiente transformación:


      
x11 . . . x1p tj1 tj1 x11 + . . . + tjp x1p y1j
Xtj =  ... . . . ...   ...  =  ..
=
.. 
     
. . 
xn1 . . . xnp tjp tj1 xn1 + . . . + tjp xnp ynj
El resultado obtenido por esta transformación se denota como y (j) y es la muestra de tamaño n
de la nueva variable Yj = t0j X.

Estas transformaciones de las variables originales presentan las propiedades que se detallan en la
siguiente proposición.

Proposición 1.2.1 Para cada una de las transformaciones Yj (j = 1, . . . , p), definimos las
siguientes características poblacionales:

Media poblacional:

E(Yj ) = E(t0j X) = t0j E(X) = t0j µ = tj1 µ1 + . . . + tjp µp

Varianza poblacional:

V ar(Yj ) = V ar(t0j X) = t0j Σtj

Covarianza entre variables:

Cov(Yi , Yj ) = Cov(t0i X, t0j X) = t0i Cov(X, X)tj = t0i V ar(X)tj = t0i Σtj
20 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Además, definimos los siguientes estadísticos muestrales análogos:


Media muestral: y (j) = t0j x

Varianza muestral: σ̂y2 = t0j Σt


b j
(j)

Dados tj , dj ∈ Rp , la covarianza muestral de las transformaciones lineales t0j X y d0j X es:


σ̂t0j X,d0j X = t0j Σd
b j

Proposición 1.2.2 La matriz de varianzas y covarianzas Σ presenta las siguientes propiedades:


Σ es simétrica y semidefinida positiva. Por tanto, los autovalores son reales y positivos
cumpliendo:
λ1 ≥ λ2 . . . ≥ λp ≥ 0

Sea ej el autovector unitario asociado al autovalor λj , con k ej k2 = 1. Los autovectores


de esta matriz son ortogonales, es decir

e0i ej = 0, si i 6= j

Disponemos estos autovectores por columnas en una matriz a la que denotamos por E.

Ep×p = (e1 , e2 , . . . , ep )

La matriz E es ortogonal:
E0 E = EE0 = Ip

Todas estas propiedades también se cumplen para su equivalente muestral Σ


b

El siguiente teorema nos será muy útil para definir las combinaciones lineales de las variables
originales que permiten obtener las componentes principales. Dichas combinaciones lineales
se obtendrán a partir de la descomposición espectral de la matriz de varianzas y covarianzas Σ,
como veremos más adelante.

Teorema 1.2.1 Teorema de la descomposición espectral

En las condiciones anteriores, se tiene que:


1. E0 ΣE = Λ = diag{λ1 , λ2 , . . . , λp }
p
X
0
2. Σ = EΛE = λj ej e0j
j=1

p p
X X
3. T raza(Σ) = σjj = T raza(Λ) = λj
j=1 j=1
1.2. COMPONENTES PRINCIPALES 21

Demostración:

1. Sea Σ la matriz de varianzas y covarianzas de X, {λ1 , λ2 , . . . , λp } sus autovalores y


E = (e1 , e2 , . . . , ep ) la matriz que contiene sus autovectores.

Sea ej el j-ésimo autovector de Σ. Por la definición de autovector, sabemos que se satisface


la siguiente igualdad:
Σej = λj ej

Sea Λ = diag{λ1 , λ2 , . . . , λp }, podemos expresar la ecuación anterior en forma matricial:

ΣE = EΛ

Por la Proposición 1.2.2 sabemos que E es ortogonal, luego multiplicando por E0 a


izquierda en ambos lados de la igualdad, se obtiene:

E0 ΣE = Λ

2. Partiendo de ΣE = EΛ y multiplicando a derecha por E0 en ambos lados de la igualdad,


se obtiene Σ = EΛE0 .

 
λ1 0 ··· 0  
e01 p
0
 0 λ2 ··· 0 
..  = X λ e e0
Σ = EΛE = (e1 , e2 , . . . , ep ) 
 
.. .. .. ..  .  j j j
 . . . .  0 j=1
ep
0 0 · · · λp

3. Por el primer apartado sabemos que Λ = E0 ΣE y por lo tanto T r(Λ) = T r(E0 ΣE)

Para probar esta parte del teorema hay que tener en cuenta que la matriz E es ortogonal y
que T r(AB) = T r(BA) para dos matrices A y B para las que se puedan definir ambos
productos.
Por consiguiente, se satisface:

T r(E0 ΣE) = T r(E0 EΣ) = T r(Σ)

Queda así demostrado que T r(Λ) = T r(Σ)


22 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Para definir las transformaciones de las variables que se llevan a cabo en la determinación de cada
componente principal es necesario calcular el valor de los coeficientes tj = (tj1 , . . . , tjp )0 ∈ Rp ,
con j ∈ {1, . . . , p}.

Comenzaremos calculando los coeficientes para la primera componente principal.

Si la expresión de esta componente es Y1 = t11 X1 + t12 X2 + ... + t1p Xp , y su varianza pobla-


cional es t01 Σt1 ; esta varianza podría aumentar sin límite elevando el valor de los coeficientes
t01 = (t11 , t12 , ..., t1p ).

Por consiguiente, para encontrar la varianza máxima es necesario introducir la siguiente restric-
ción en nuestro problema: t01 t1 = 1

Si queremos maximizar una función dependiente de un gran número de variables y sujeta a una,
o más de una restricción, usamos el método de Los multiplicadores de Lagrange.

La técnica de los multiplicadores de Lagrange prueba que, si λ1 es el mayor autovalor de Σ;


entonces requiriendo que t01 t1 = 1, puede verse que t1 es el autovector de la matriz Σ correspon-
diente a λ1 y que este autovalor es la varianza de la primera CP.

Del mismo modo, con los estadísticos muestrales correspondientes se obtiene la primera compo-
nente principal muestral, que definiremos en la próxima sección.

Veamos más detenidamente el razonamiento que se sigue para obtener las transformaciones de
las variables.

Como hemos planteado anteriormente, para la primera CP tenemos que resolver el problema:
max(t01 Σt1 )
s.a. : t01 t1 = 1
Consideramos la siguiente función y buscamos el máximo, derivando su expresión e igualándola
a cero:
L(t1 ) = t01 Σt1 − λ(t01 t1 − 1)
(1.2)
∂(L(t1 ))
= 2Σt1 − 2λIt1 = 0
∂(t1 )
(1.3)
(Σ − λI)t1 = 0 (1.4)
Por el Teorema de Rouché-Frobenius, la matriz (Σ − λI) ha de ser singular para que el sistema
tenga solución distinta de cero.
1.2. COMPONENTES PRINCIPALES 23

Luego: |Σ − λI| = 0 ⇔ λ es autovalor de Σ.

A partir de la igualdad dada en 1.4 se obtiene:

Σt1 = λIt1

Sustituimos esta igualdad en la expresión de la varianza:

V ar(Y1 ) = t01 Σt1 = t01 λIt1 = λt01 t1 = λ

Luego, para maximizar la varianza de la combinación lineal que define la primera componente
principal hay que tomar el mayor autovalor (λ1 ) y el correspondiente autovector (e1 ) de Σ.

Para la segunda componente; además, se requiere que sea incorrelada con la primera, es decir:

Cov(Y2 , Y1 ) = Cov(t02 X, t01 X) = t02 Σt1 = 0

Por la condición anterior se verifica Σt1 = λIt1 , de forma que sustituyendo esta igualdad en la
expresión de la covarianza obtenemos lo siguiente:

t02 Σt1 = t02 λIt1 = 0 ⇔ t02 t1 = 0

Como sabemos, los autovectores de Σ son ortonormales, así que cumplen ambas condiciones.

Se sigue el mismo razonamiento con todas las componentes, resolviéndose así el problema de
maximizar la varianza de la combinación lineal correspondiente que además es incorrelada con
las combinaciones anteriores.

Obtenemos, por tanto, la transformación de las variables para cada componente.

Definición 1.2.1 Se define la j-ésima componente principal como:

Yj = e0j X, j = 1, . . . , p

con ej autovectores unitarios de Σ.

Teorema 1.2.2 La variabilidad explicada por la j-ésima CP es

λ
Pp j
s=1 λs

A su vez, la proporción de variabilidad total no explicada por las k primeras CP es


Pp
λj
Pj=k+1
p
s=1 λs
24 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

1.3. Cálculo y propiedades de las Componentes Principales


Muestrales
Una vez resuelto el problema que supone calcular los coeficientes para las transformaciones de
las variables originales X, podemos definir las CP muestrales.

Definición 1.3.1 Se define la j-ésima componente principal muestral como:


Yj = ê0j X, j = 1, . . . , p.

con êj autovectores unitarios de Σ.


b

A continuación introducimos notación adecuada para trabajar con las CP muestrales.


Puntuaciones correspondientes a la j-ésima componente principal

Estos son los valores muestrales de la j-ésima componente principal. Se obtienen multipli-
cando el autovector ê0j por la muestra del individuo correspondiente.
   0 
y1j êj x1
y (j) =  ...  =  ...  = Xêj , j = 1, . . . , p
   
0
ynj êj xn

Puntuaciones de las p componentes principales se recogen en la siguiente matriz


 
x01 ê1 . . . x01 êp
Y =  ... .. ..  = XE
b = (y , . . . , y )

. .  (1) (p)
x0n ê1 . . . x0n êp

Estos son los valores muestrales de las p componentes principales, recogidos en una matriz
de dimensión n × p.

Cada columna es una muestra de tamaño n de la correspondiente CP. Cada fila muestra los
valores de todas las CP para cada uno de los individuos de la muestra.
Transformación de los datos muestrales de las variables originales en los datos muestrales
de las CP

Podemos escribir la matriz anterior por filas, siendo cada fila y 0i los valores muestrales de
las p componentes principales para el individuo i-ésimo.
   0 
x01 ê1 . . . x01 êp y1
 .. . .. .
..  = XE b =  ... 
Y= .
 

x0n ê1 . . . x0n êp y 0n
1.3. CÁLCULO Y PROPIEDADES DE LAS COMPONENTES PRINCIPALES MUESTRALES25

Los valores muestrales de las CP para el individuo i-ésimo son

b 0 xi ,
yi = E i = 1, . . . , n

A continuación introducimos varios resultados que caracterizan a las componentes principales


muestrales.

Proposición 1.3.1 Se mantienen las distancias entre los datos originales y los datos transfor-
mados:
d2 (y h , y i ) = d2 (xh , xi )

Demostración:
La prueba se basa en la ortogonalidad de la matriz E.
b

d2 (y h , y i ) = (y i − y h )0 (y i − y h ) = (xi − xh )0 E
bEb 0 (xi − xh ) = (xi − xh )0 (xi − xh ) = d2 (xh , xi )

Proposición 1.3.2 Las componentes principales muestrales tienen las siguientes propiedades:

Matriz de varianzas-covarianzas:

Σ b 0Σ
by = E bE b = diag{λ̂1 , λ̂2 , . . . , λ̂p }
b =Λ

Demostración:

Σ b 0 X) = E
b y = Cov(E b 0 Cov(X)E b 0Σ
b =E bEb =Λ
b

Esta última igualdad se obtiene por el Teorema de la Descomposición Espectral (Teorema


1.2.1) .

La varianza muestral de la j-ésima CP es igual al autovalor correspondiente, λ̂j

σ̂y2 = λ̂j , ∀j = 1, ..., p


(j)

Covarianza muestral entre dos CPs: σ̂y(j) y(s) = 0, ∀j 6= s

Pp Pp Pp
j=1 σ̂x2(j) = tr(Σ)
b =
j=1 λ̂j = j=1 σ̂y2
(j)
26 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

1.4. Componentes principales para datos estandarizados


En ocasiones, las variables que caracterizan los datos con los que trabajamos están medidas en
escalas muy distintas.

La técnica de Análisis de Componentes Principales es sensible a cambios escalares, es decir, si


intentamos igualar escalas multiplicando una variable por un escalar obtendremos diferentes
autovalores y autovectores. Esto se debe a que la descomposición de autovalores se hace en
la matriz de covarianzas y no en la de correlación; por tanto, los resultados que obtendremos
usando ambas matrices serán distintos.

Es por ello por lo que a veces es preferible trabajar con los datos estandarizados, porque de esta
forma las magnitudes de las variables son similares. Se resuelve así el problema que supone que
las variables con mayor varianza influyan más en la determinación de la primera componente
principal.

Definición 1.4.1 La matriz de datos estandarizados es la siguiente:


. . . x1pσ̂−x
 x11 −x1 p

σ̂1 p

X(s) = 
 .. .. 
. . 
xn1 −x1 xnp −xp
σ̂1
... σ̂p

donde xj es la media muestral y σ̂j la desviación típica de x(j) , siendo x(j) un vector de dimensión
n × 1 conteniendo la muestra de tamaño n para la variable Xj (j = 1, . . . , p).
Si calculamos las componentes principales a partir de la matriz de datos estandarizados, esto
equivale a calcular las componentes a partir de la matriz de correlaciones en lugar de la matriz
de covarianzas.

En la práctica, sólo debemos utilizar la matriz Σb para obtener las componentes principales
cuando las variables originales tengan la misma escala de medida.

Si no se cumple esta condición emplearemos la matriz de datos estandarizados, de forma que la


suma total de los autovalores será p (pues se asumen varianzas iguales a 1 para todas las varia-
λ̂
bles) y la proporción de variabilidad recogida por la j-ésima componente principal muestral es pj .

1.4.1. Círculo de correlación


El círculo de correlación es una herramienta que permite conocer la contribución de cada variable
en la determinación de una CP.
Esta contribución, si es que existe, puede ser tanto positiva como negativa. Depende de la proxi-
midad de la variable en su representación a la periferia del círculo y al eje definido por la CP.
1.4. COMPONENTES PRINCIPALES PARA DATOS ESTANDARIZADOS 27

A su vez, permite conocer relaciones entre las variables.

Gracias a este gráfico podremos visualizar las variables más correladas con las componentes
principales que recojan más información.

Para comenzar, definimos la covarianza entre dos conjuntos de variables como:

Cov(X, Y ) = E(XY 0 ) − E(X)E(Y 0 ) = E(XX 0 E) − µµ0 E = V ar(X)E =

= ΣE = EΛE0 E = EΛ
Por tanto, la correlación lineal entre la variable original Xj y la CP Yk se expresa como:

Cov(Xj , Yk ) λk ekj λk ekj λk ekj
ρXj ,Yk = p =p = √ =
V ar(Xj )V ar(Yk ) V ar(Xj )V ar(Yk ) σ j λk σj

cumpliéndose:
p Pp
X λk ekj 2 σj2
ρ2Xj ,Yk = k=1
= =1
k=1
σj2 σj2
Es obvio que si estudiamos la correlación de la variable Xj con las dos primeras componentes
principales, ρ2Xj ,Y1 + ρ2Xj ,Y2 <= 1, por tanto los puntos estarán siempre dentro de un círculo de
radio 1.

A continuación, mostraremos con un ejemplo la interpretación del círculo de correlación y su


utilidad.

Aplicación
El fichero de datos contiene 6 medidas distintas de billetes de 1000 francos suizos, de manera
que 100 de ellos son verdaderos (los primeros) y los restantes 100 falsos.
Los datos se encuentran en el paquete mclust.

Llevaremos a cabo un Análisis de Componentes Principales ayudándonos del paquete


FactoMineR.

library(mclust)
library(FactoMineR)
data(banknote)
str(banknote)

'data.frame': 200 obs. of 7 variables:


$ Status : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 ...
28 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

$ Length : num 215 215 215 215 215 ...


$ Left : num 131 130 130 130 130 ...
$ Right : num 131 130 130 130 130 ...
$ Bottom : num 9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
$ Top : num 9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 ...
$ Diagonal: num 141 142 142 142 142 ...

Con la función str vemos la estructura de los datos con los que trabajamos. Tenemos un
data.frame con 200 observaciones para las que se miden las siguientes variables:

1. Status: factor que representa el estado del billete. “genuine” indica que es verdadero,
mientras que “counterfeit” indica que no.

2. Length: longitud del billete en mm. (X1 )

3. Lef t: anchura del borde izquierdo en mm. (X2 )

4. Right: anchura del borde derecho en mm. (X3 )

5. Bottom: anchura del margen inferior en mm. (X4 )

6. T op: anchura del margen superior en mm. (X5 )

7. Diagonal: longitud de la diagonal en mm. (X6 )

Por tanto, sólo existe una variable categórica (Status).

summary(banknote)

Status Length Left Right


counterfeit:100 Min. :213.8 Min. :129.0 Min. :129.0
genuine :100 1st Qu.:214.6 1st Qu.:129.9 1st Qu.:129.7
Median :214.9 Median :130.2 Median :130.0
Mean :214.9 Mean :130.1 Mean :130.0
3rd Qu.:215.1 3rd Qu.:130.4 3rd Qu.:130.2
Max. :216.3 Max. :131.0 Max. :131.1
Bottom Top Diagonal
Min. : 7.200 Min. : 7.70 Min. :137.8
1st Qu.: 8.200 1st Qu.:10.10 1st Qu.:139.5
Median : 9.100 Median :10.60 Median :140.4
Mean : 9.418 Mean :10.65 Mean :140.5
3rd Qu.:10.600 3rd Qu.:11.20 3rd Qu.:141.5
Max. :12.700 Max. :12.30 Max. :142.4
1.4. COMPONENTES PRINCIPALES PARA DATOS ESTANDARIZADOS 29

Una vez visto un resumen de cada variable con summary, llevamos a cabo un Análisis de
Componentes Principales con la función PCA de FactoMineR.

Además de especificar el data.frame de los datos con los que trabajamos, con el argumento
“graph = FALSE” indicamos que no queremos ninguna representación gráfica, y con “quali.sup=1”
indicamos que la primera variable es categórica y debe tratarse de forma diferente al resto.

res.pca=PCA(banknote, graph = FALSE,quali.sup=1 )

Con el comando summary podemos ver un resumen del ACP llevado a cabo. Muestra los
autovalores y autovectores de la matriz de varianzas y covarianzas muestrales Σ,
b correlaciones
con componentes principales, así como otra información de utilidad; diferenciando individuos,
variables y variables categóricas.

summary(res.pca)

Veamos este resumen con detalle.


En primer lugar aparecen los autovalores, porcentaje de varianza explicada y porcentaje de
varianza acumulado para las 6 CP.

Eigenvalues
Dim.1 Dim.2 Dim.3 Dim.4
Variance 2.946 1.278 0.869 0.450
% of var. 49.093 21.301 14.484 7.496
Cumulative % of var. 49.093 70.394 84.878 92.374
Dim.5 Dim.6
Variance 0.269 0.189
% of var. 4.478 3.148
Cumulative % of var. 96.852 100.000

Podemos observar que las 2 primeras componentes explican un 70.394 % de la variabilidad total.

En la primera columna del apartado “Individuals” se indica la distancia en la representación de


cada individuo a la categoría a la que pertenece. Tan sólo se indica información de los 5 primeros
individuos.

Individuals
Dist Dim.1 ctr cos2 Dim.2 ctr
1 | 3.970 | 1.747 0.518 0.194 | 1.651 1.066
2 | 2.533 | -2.274 0.878 0.806 | -0.539 0.114
3 | 2.457 | -2.277 0.880 0.859 | -0.108 0.005
4 | 2.414 | -2.284 0.885 0.895 | -0.088 0.003
5 | 4.234 | -2.632 1.176 0.386 | 0.039 0.001
30 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

cos2 Dim.3 ctr cos2


1 0.173 | 1.424 1.166 0.129 |
2 0.045 | 0.533 0.163 0.044 |
3 0.002 | 0.717 0.296 0.085 |
4 0.001 | -0.606 0.211 0.063 |
5 0.000 | 3.196 5.878 0.570 |

En las columnas “Dim.1”, “Dim.2” y “Dim.3” se muestran los valores y (j) para j = 1, 2, 3; en
otras palabras, los valores de las 3 primeras CP para los primeros individuos de la muestra.

En la columna “ctr” se representa la contribución de cada individuo en la determinación de cada


CP.
La calidad de representación de un elemento (ya sea individuo o variable) en los ejes de un
gráfico se mide por el coseno al cuadrado del ángulo existente entre el vector del elemento y
su proyección en el eje correspondiente. Estos valores, que a su vez se corresponden con la
correlación al cuadrado entre estos elementos y cada CP, aparecen en la columna “cos2”.

Además de la contribución y la calidad de representación, en el apartado “Variables” aparecen los


valores de (ê1 , ê2 , ê3 ) necesarios para la determinación de cada CP. Se muestran en las columnas
“Dim.1”, “Dim.2” y “Dim.3”.

Variables
Dim.1 ctr cos2 Dim.2 ctr cos2
Length | -0.012 0.005 0.000 | 0.922 66.503 0.850 |
Left | 0.803 21.880 0.644 | 0.387 11.694 0.149 |
Right | 0.835 23.686 0.698 | 0.285 6.374 0.081 |
Bottom | 0.698 16.545 0.487 | -0.301 7.088 0.091 |
Top | 0.631 13.534 0.399 | -0.103 0.837 0.011 |
Diagonal | -0.847 24.350 0.717 | 0.310 7.504 0.096 |
Dim.3 ctr cos2
Length -0.016 0.031 0.000 |
Left 0.096 1.069 0.009 |
Right 0.115 1.525 0.013 |
Bottom 0.544 34.052 0.296 |
Top -0.734 62.027 0.539 |
Diagonal 0.106 1.297 0.011 |

Por último, aparece información sobre las categorías del factor Status.

Supplementary categories
Dist Dim.1 cos2 v.test Dim.2
counterfeit | 1.548 | 1.498 0.937 12.314 | -0.348
genuine | 1.548 | -1.498 0.937 -12.314 | 0.348
1.4. COMPONENTES PRINCIPALES PARA DATOS ESTANDARIZADOS 31

cos2 v.test Dim.3 cos2 v.test


counterfeit 0.051 -4.340 | 0.001 0.000 0.012 |
genuine 0.051 4.340 | -0.001 0.000 -0.012 |

Con la columna “Dist” se indica la distancia de cada categoría al centro de los ejes del gráfico.

Se indican las coordenadas y calidad de la representación (“Dim.” y “cos2”, respectivamente);


además de un test “v.test” para los individuos de cada categoría. Veremos este test con más
detenimiento en otras aplicaciones prácticas del ACP.

Una vez vistos los resultados obtenidos con el comando summary, podríamos; por lo tanto,
definir las dos primeras CP como:

Y1 = −0,012X1 + 0,803X2 + 0,835X3 + 0,698X4 + 0,631X5 − 0,847X6

Y2 = 0,922X1 + 0,387X2 + 0,285X3 − 0,301X4 − 0,103X5 + 0,310X6


Se observa que la primera CP es esencialmente la suma de las variables Lef t (X2 ) y Right (X3 )
y la resta de la variable Diagonal (X6 ).

La segunda CP, sin embargo, depende en su mayor parte de la variable Length(X1 ).

Representaremos las correlaciones muestrales entre las variables y estas dos CP, que serán los
ejes de abcisas y ordenadas; respectivamente. La suma de estas correlaciones al cuadrado para
cada variable dará como resultado la longitud del vector en la representación.

Como hemos especificado con anterioridad, las correlaciones al cuadrado aparecen en la columna
“cos2” del resumen.

Hacemos uso de la función plot y con el comando “choix” especificamos que sólo queremos
una representación de las variables.

Puede verse en el gráfico que, efectivamente, la primera CP está muy bien explicada por las
variables X2 (Lef t), X3 (Right) y X6 (Diagonal).

Consideremos la variable Lef t (X2 ), para la que tenemos que:

rX2 ,Y1 = (cos(X2 , Y1 )) = 0,80249611

rX2 ,Y2 = (cos(X2 , Y2 )) = 0,386


Es sencillo comprobar que se corresponde con la respresentación gráfica. La suma de estas
correlaciones muestrales al cuadrado es igual al módulo del vector representado, que en este
caso es 0,793.
32 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

plot(res.pca,choix = "var")

Variables factor map (PCA)


1.0

Length
0.5

Left
Diagonal Right
Dim 2 (21.30%)

0.0

Top

Bottom
−0.5
−1.0

−1.0 −0.5 0.0 0.5 1.0

Dim 1 (49.09%)
1.5. COMPONENTES PRINCIPALES DE UNA MATRIZ DE DATOS BIVARIANTE 33

Algo muy parecido sucede con la variable Right (X3 ), que está un poco más correlacionada con
la primera CP y menos con la segunda CP; aunque a efectos de la representación la diferencia es
muy pequeña. El módulo del vector es 0,779 en este caso.

Por último, la variable Diagonal (X6 ) tiene correlaciones rX6 ,Y1 = (cos(X6 , Y1 )) = −0,84675
y rX6 ,Y2 = (cos(X6 , Y2 )) = 0,309838. Los signos de estas correlaciones los hemos obtenido
a partir de la representación gráfica, ya que con la información que obtenemos de la función
“res.pca” aparecen las correlaciones al cuadrado.

El valor del módulo del vector es 0,813; lo que significa que es bastante cercano a la periferia
del círculo, como se puede observar en el gráfico.

Por otra parte, la segunda CP, Y2 , está explicada en su mayor parte por la variable Length(X1 ).
p
X1 tiene muy poca correlación con la primera CP (pues (cos2 (X1 , Y1 )) apenas es cero),
mientras que rX1 ,Y2 = (cos(X1 , Y2 )) = 0,92195. Además, si sumamos ambas correlaciones a
cuadrado, obtenemos que el módulo del vector en el círculo es 0,85.

Variables cercanas a la periferia indican que la suma de sus correlaciones con las CP son muy
cercanas a 1, lo que indica que las CP explican una alta proporción de las variables.

También podemos deducir del gráfico que variables como Lef t y Right o T op y Bottom están
muy correladas entre sí, debido a que están muy próximas entre sí en la representación. Por otra
parte, las variables Diagonal y Bottom están negativamente correladas de forma significativa;
ya que están representadas en lados opuestos. Si existen variables ortogonales, eso indica que no
están correladas.

1.5. Componentes Principales de una matriz de datos biva-


riante
Consideremos una matriz de datos bivariante con dos variables X1 y X2 con coeficiente de
correlación lineal ρ. La matriz de correlación sería:
 
1 ρ
R=
ρ 1

Necesitamos encontrar los autovalores y autovectores de R.


En primer lugar, calculamos los autovalores:

|R − λI| = 0 ⇔ (1 − λ)2 − ρ2 = 0

El resultado de esta última igualdad es λ1 = 1 + ρ, λ2 = 1 − ρ.


34 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Para el cálculo de los autovectores, se resuelve en primer lugar la ecuación Re1 = λ1 e1 .


Obtenemos el siguiente sistema:

e11 + ρe12 = (1 + ρ)e11 ,
ρe11 + e12 = (1 + ρ)e12

Como las dos ecuaciones son idénticas, tomamos e11 = e12 e introducimos la restricción de
normalización (e01 e1 = 1).

Obtenemos que las coordenadas para el primer autovector son e11 = e12 = √1
2

Siguiendo el mismo razonamiento para resolver la ecuación Re2 = λ2 e2 , obtenemos: e21 =


√1 = −e22 .
2

Por lo tanto, las dos componentes principales se expresan como sigue:


1
Y1 = √ (X1 + X2 )
2
1
Y2 = √ (X1 − X2 )
2
El cálculo de las varianzas de las componentes es sencillo, pues como hemos definido con
anterioridad, se corresponde con el autovalor λj y tiene como resultado:

V ar(Y1 ) = 1 + ρ, V ar(Y2 ) = 1 − ρ

Si ρ fuera negativo, el orden de los autovalores cambiaría y por tanto las componentes principales
se expresarían al contrario. Si ρ fuera nulo, los autovalores serían igual a uno.

Un par de aclaraciones necesarias:

Hay una forma arbitraria de escoger el signo de los elementos de ej . Se suele tomar ej1
positivo.

Los coeficientes que definen las dos componentes no dependen de ρ a pesar de que la
proporción de varianza explicada por cada uno si lo haga.

1.6. Análisis de Componentes Principales y Regresión Lineal


Cuando queremos llevar a cabo una regresión lineal, es común que se nos planteen dificultades
tales como un número demasiado elevado de predictores para la variable objetivo o bien colinea-
lidad entre estos predictores; en otras palabras, variables con alta correlación entre sí.
1.6. ANÁLISIS DE COMPONENTES PRINCIPALES Y REGRESIÓN LINEAL 35

Gracias a la técnica de Análisis de Componentes Principales pueden abordarse estos problemas,


permitiendo simultáneamente reducir dimensión y explicar un gran porcentaje de variabilidad de
la variable respuesta.

Veamos un ejemplo práctico que lo ilustre.

1.6.1. Aplicación
En este ejemplo, trabajamos con un conjunto de datos que recoge la contaminación del aire en
varias ciudades de EEUU.

Los datos que manejamos se encuentran en el paquete HSAUR3, que se detalla en el último
capítulo de esta memoria junto con el resto de paquetes utilizados en las aplicaciones prácticas
de las técnicas.

Con el paquete MVA podemos representar boxplot bivariantes, por lo que nos será de gran utilidad.

oldopt <- options(digits=3)


options(width=60)
on.exit( {options(oldopt)} )
require(HSAUR3)
require(MVA)

data(USairpollution)
names(USairpollution)

## [1] "SO2" "temp" "manu" "popul" "wind"


## [6] "precip" "predays"

En nuestros datos se miden 7 variables. Estas variables son:

SO2: contenido de SO2 en el aire (microgramos por metro cúbico)


temp: temperatura media anual en grados Fahrenheit
manu: número de empresas de producción con 20 trabajadores o más
popul: población según el censo de 1970 (miles de personas)
wind: velocidad media del viento por año (millas por hora)
precip: precipitación media anual (pulgadas)
predays: número medio de días con precipitaciones por año
36 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

str(USairpollution)

## 'data.frame': 41 obs. of 7 variables:


## $ SO2 : int 46 11 24 47 11 31 110 23 65 26 ...
## $ temp : num 47.6 56.8 61.5 55 47.1 55.2 50.6 54 49.7 51.5 ...
## $ manu : int 44 46 368 625 391 35 3344 462 1007 266 ...
## $ popul : int 116 244 497 905 463 71 3369 453 751 540 ...
## $ wind : num 8.8 8.9 9.1 9.6 12.4 6.5 10.4 7.1 10.9 8.6 ...
## $ precip : num 33.36 7.77 48.34 41.31 36.11 ...
## $ predays: int 135 58 115 111 166 148 122 132 155 134 ...

Observamos que trabajamos con un data.frame de dimensión 41 × 7, y por lo tanto, medimos las
7 variables anteriores en 41 individuos.
Veamos un resumen de los datos.

summary(USairpollution)

## SO2 temp manu


## Min. : 8.00 Min. :43.50 Min. : 35.0
## 1st Qu.: 13.00 1st Qu.:50.60 1st Qu.: 181.0
## Median : 26.00 Median :54.60 Median : 347.0
## Mean : 30.05 Mean :55.76 Mean : 463.1
## 3rd Qu.: 35.00 3rd Qu.:59.30 3rd Qu.: 462.0
## Max. :110.00 Max. :75.50 Max. :3344.0
## popul wind precip
## Min. : 71.0 Min. : 6.000 Min. : 7.05
## 1st Qu.: 299.0 1st Qu.: 8.700 1st Qu.:30.96
## Median : 515.0 Median : 9.300 Median :38.74
## Mean : 608.6 Mean : 9.444 Mean :36.77
## 3rd Qu.: 717.0 3rd Qu.:10.600 3rd Qu.:43.11
## Max. :3369.0 Max. :12.700 Max. :59.80
## predays
## Min. : 36.0
## 1st Qu.:103.0
## Median :115.0
## Mean :113.9
## 3rd Qu.:128.0
## Max. :166.0

No tendremos en cuenta la variable SO2, y nos centraremos en dos variables basadas en la


ecología humana (popul, manu) y cuatro basadas en el clima (temp, wind, precip, predays).

Dado que las variables están medidas en escalas muy distintas, obtenemos las componentes
a partir de la matriz de correlaciones en lugar de usar la matriz de varianzas y covarianzas.
1.6. ANÁLISIS DE COMPONENTES PRINCIPALES Y REGRESIÓN LINEAL 37

Como especificamos en la sección anterior, esto es equivalente a trabajar con la matriz de datos
estandarizados.

Con la siguiente instrucción, transformamos la variable temp por su equivalente negativa,


ganando así en interpretación. La llamaremos negtemp.

USairpollution$negtemp<-USairpollution$temp*(-1)
USairpollution$temp<-NULL

Calculamos la matriz de correlaciones.

cor(USairpollution[,-1])

## manu popul wind precip


## manu 1.00000000 0.95526935 0.23794683 -0.03241688
## popul 0.95526935 1.00000000 0.21264375 -0.02611873
## wind 0.23794683 0.21264375 1.00000000 -0.01299438
## precip -0.03241688 -0.02611873 -0.01299438 1.00000000
## predays 0.13182930 0.04208319 0.16410559 0.49609671
## negtemp 0.19004216 0.06267813 0.34973963 -0.38625342
## predays negtemp
## manu 0.13182930 0.19004216
## popul 0.04208319 0.06267813
## wind 0.16410559 0.34973963
## precip 0.49609671 -0.38625342
## predays 1.00000000 0.43024212
## negtemp 0.43024212 1.00000000

Cabe destacar que las variables popul y manu están muy relacionadas, pues su coeficiente de
correlación es 0,95527.

Construimos una matriz scatterplot de las 6 variables, incluyendo además el histograma para cada
variable en la diagonal. Esta matriz esta formada por una serie de gráficos de puntos dispersos,
donde cada eje representa una variable y cada punto un individuo de la muestra. Puede sernos
útil para visualizar outliers.
38 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

data("USairpollution",package="HSAUR3")
panel.hist<-function(x,...){
usr<-par("usr");on.exit(par(usr))
par(usr=c(usr[1:2],0,1.5))
h<-hist(x,plot=FALSE)
breaks<-h$breaks;nB<-length(breaks)
y<-h$counts; y<-y/max(y)
rect(breaks[-nB],0,breaks[-1],y,col="grey",...)
}
USairpollution$negtemp<-USairpollution$temp*(-1)
USairpollution$temp<-NULL
pairs(USairpollution[,-1],diag.panel = panel.hist,pch=".",cex=1.5)

0 1500 3000 10 30 50 −75 −60 −45

3000
manu

1500
0
3000

popul
1500
0

12
wind

10
8
6

precip
50
30
10

predays
140
80
40
−45

negtemp
−60
−75

0 1500 3000 6 8 10 12 40 80 140

Figura 1.1: Matriz scatterplot de las 6 variables consideradas en los datos.


1.6. ANÁLISIS DE COMPONENTES PRINCIPALES Y REGRESIÓN LINEAL 39

Ahora, llevamos a cabo el Análisis de Componentes Principales.

usair_pca<-princomp(USairpollution[,-1],cor=TRUE)
summary(usair_pca,loadings=TRUE)

## Importance of components:
## Comp.1 Comp.2 Comp.3
## Standard deviation 1.4819456 1.2247218 1.1809526
## Proportion of Variance 0.3660271 0.2499906 0.2324415
## Cumulative Proportion 0.3660271 0.6160177 0.8484592
## Comp.4 Comp.5 Comp.6
## Standard deviation 0.8719099 0.33848287 0.185599752
## Proportion of Variance 0.1267045 0.01909511 0.005741211
## Cumulative Proportion 0.9751637 0.99425879 1.000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## manu -0.612 0.168 -0.273 -0.137 0.102 0.703
## popul -0.578 0.222 -0.350 -0.695
## wind -0.354 -0.131 0.297 0.869 -0.113
## precip -0.623 -0.505 0.171 0.568
## predays -0.238 -0.708 -0.311 -0.580
## negtemp -0.330 -0.128 0.672 -0.306 0.558 -0.136

Gracias al comando summary se observa que las tres primeras componentes explican el
0.8484592 % de la varianza.

Si queremos interpretar los datos usando etiquetas para las componentes, tendremos que exami-
nar los coeficientes que definen cada una de ellas y que están recogidos en las columnas “Comp.”
del resumen.
Por ejemplo, la primera componente podría etiquetarse como “calidad de vida”, la segunda
como “tiempo húmedo” (dado que los coeficientes más altos acompañan a las variables precip
y predays), y la tercera como “tipo de clima” (las variables con más peso son precip y negtemp).

Podemos representar los individuos mediante un boxplot bivariante con ejes las 3 primeras
componentes principales, ya que es un método más objetivo que las matrices scatterplot para
etiquetar outliers.
Está basado en la construcción de un par de elipses concéntricas, de forma que en la de menor
radio se incluye el 50 % de los datos y en la otra se incluyen los datos que puedan ser outliers.

Podemos deducir que “Chicago” es un outlier; así como “Phoenix” y “Philadelphia”.


“Phoenix” ofrece la mejor calidad de vida y “Buffalo” sería la ciudad con el ambiente más
húmedo.
40 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

pairs(usair_pca$scores[,1:3],ylim=c(-6,4),xlim=c(-6,4),
panel=function(x,y,...){
text(x,y,abbreviate(row.names(USairpollution)),
cex=0.6)
bvbox(cbind(x,y),add=TRUE)
})

−6 −4 −2 0 2 4

4
Phnx Phnx

2
Miam
NwOr
Chrl LttR Albq NwOr LttR Chrl Albq
Miam
Jcks
Rchm Jcks Rchm
Wlmn
Nshv SlLC Nshv Wlmn SlLC
Mmph
Atln
Nrfl
Albn
Cncn
Lsvl SnFr Mmph
AtlnCncn
Nrfl
SnFr Albn
Lsvl
Hrtf
Clmb DsMnWcht
Omah
KnsC DllsDnvr KnsCHrtf
Clmb Dnvr
Wcht
DsMn
Omah

0
Prvd Wshn Wshn
Dlls Prvd
SttlPttsIndn St.L St.L Sttl
Hstn
Bltm Hstn BltmIndn
Ptts

Comp.1 Bffl Mlwk


Mnnp
Clvl Clvl
Mlwk
Bffl
Mnnp

−2
Dtrt Dtrt
Phld Phld

−4
−6
Chcg Chcg

Phnx Phnx
4

Albq Albq
SnFr SnFr
2

Dnvr Dnvr
Chcg SlLC Chcg SlLC
Dlls Dlls
Phld Wcht Phld Wcht
Dtrt Omah
St.L
KnsC Dtrt Omah
St.L
KnsC
Hstn
DsMn
Wshn Hstn Wshn DsMn
0

Mlwk Bltm
Mnnp IndnMmph Rchm
LttR Mmph
Bltm
Rchm
LttR Indn Mlwk
Mnnp
Cncn
Wlmn
Lsvl
Nshv NshvCncn
Wlmn
Lsvl
Clvl
Clmb
Ptts
Prvd
Sttl
Atln
NrflJcks
Albn
Hrtf NwOr
Chrl
Miam
Comp.2 Miam
Jcks
NwOr
Atln Nrfl
Clmb Albn
Chrl
Hrtf
Ptts
Prvd
Clvl
Sttl
−2

Bffl Bffl
−4
−6

4
2

Mnnp
Bffl
Mlwk DsMnSlLC Bffl Mnnp
Mlwk
DsMn
Albn Albq Albn Omah SlLC
Dnvr
Omah
Wcht Wcht Dnvr Albq
Clvl Prvd
Ptts Prvd
Ptts
Sttl HrtfWlmn
Clmb SttlClvl
Hrtf
Wlmn
Clmb
Dtrt Indn
KnsCSnFr Indn
KnsC
Dtrt SnFr
0

St.L Nrfl Chrl


Cncn ChrlNrfl
Cncn St.L
Wshn
Bltm LsvlRchm Wshn
Bltm
Lsvl
Rchm

Phld
Dlls NshvLttRPhnx
Atln
Mmph
LttR
Nshv
Atln
Mmph
Phld
Dlls Phnx
Comp.3
−2

Hstn Jcks
NwOr Jcks Hstn
NwOr
Chcg Chcg
Miam Miam
−4
−6

−6 −4 −2 0 2 4 −6 −4 −2 0 2 4

Figura 1.2: Boxplots bivariantes de las tres primeras componentes principales.


1.6. ANÁLISIS DE COMPONENTES PRINCIPALES Y REGRESIÓN LINEAL 41

A continuación, determinaremos qué variables son las mejores predictoras del grado de contami-
nación en el aire de una ciudad (variable SO2).

Esta pregunta suele responderse con regresión lineal múltiple, pero no puede aplicarse en este
caso debido a un problema de colinealidad en los datos (como apuntamos anteriormente, las
variables popul y manu están altamente correlacionadas).

Podríamos prescindir de estas dos variables, aunque es mejor aplicar regresión lineal múltiple a
las componentes principales de las variables originales. Se sigue este paso porque las componen-
tes principales son incorreladas entre sí.

Pero, para resolver el problema de esta forma necesitamos preguntarnos: ¿Cuántas componentes
principales necesitaría usar como variables explicativas en la regresión?

Llevamos a cabo la regresión lineal en R con las 3 primeras CP que, como indicamos antes,
explican casi el 85 % de la variabilidad.

usair_reg<-lm(SO2~usair_pca$scores[,1:3],
data=USairpollution)
summary(usair_reg)

Call:
lm(formula = SO2 ~ usair_pca$scores[, 1:3], data = USairpollution)

Residuals:
Min 1Q Median 3Q Max
-36.42 -10.98 -3.18 12.09 61.27

Coefficients:
Estimate Std. Error t value
(Intercept) 30.049 2.907 10.34
usair_pca$scores[, 1:3]Comp.1 -9.942 1.962 -5.07
usair_pca$scores[, 1:3]Comp.2 -2.240 2.374 -0.94
usair_pca$scores[, 1:3]Comp.3 -0.375 2.462 -0.15
Pr(>|t|)
(Intercept) 1.8e-12 ***
usair_pca$scores[, 1:3]Comp.1 1.1e-05 ***
usair_pca$scores[, 1:3]Comp.2 0.35
usair_pca$scores[, 1:3]Comp.3 0.88
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
42 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Residual standard error: 18.6 on 37 degrees of freedom


Multiple R-squared: 0.418,Adjusted R-squared: 0.371
F-statistic: 8.87 on 3 and 37 DF, p-value: 0.000147

Claramente, la puntuación de la primera componente principal es la más predictiva, pues es la


que tiene menor error.
Etiquetada como “calidad de vida”, esta combinación de las variables originales es la que mejor
predice el grado de contaminación en el aire (SO2 ).

1.7. Escogiendo el número de componentes


Hay distintos criterios para decidir cuantas componentes principales necesitamos. Estos son:
Criterio del porcentaje.
Tomamos las componentes suficientes para explicar un porcentaje alto de la variabilidad
de las variables originales. La variabilidad explicada por cada componente (autovalores)
suele estabilizarse en su representación a partir de un cierto número. Éste índica la última
componente que tendremos en cuenta.

Normalmente, tomamos componentes hasta que la varianza explicada esté en torno a un


70 % o 90 %; aunque el punto óptimo de corte suele ser λ∗ = 00 7 ∗ v, siendo v = tr(Σ)
p

Criterio de Kaiser.
Cuando trabajamos con la matriz R, dado que se asumen varianzas iguales a uno no
tendremos en cuenta los autovalores menores que esa cantidad. Estudios de Montecarlo
afirman que el valor óptimo de λ∗ para el punto de corte en la gráfica es 0’7.
Excluir las componentes cuyos autovalores estén por debajo de la media de los autovalores.
Todos estos métodos dan lugar a diferentes conclusiones.

1.8. Aplicación práctica del ACP


Cuando el tipo de estudio que llevamos a cabo es de descubrimiento de clases, necesitamos
definir distancias entre genes o muestras para reconocer patrones o identificar efecto lote.

Es por ello por lo que surge la necesidad de reducir la dimensión de la matriz de datos, para
facilitar la interpretación mediante la observación de tendencias en la visualización gráfica.
La forma de hacerlo es mediante técnicas de estadística descriptiva como ACP o clustering.

En estos ejemplos veremos el Análisis de Componentes Principales. Para hacer una reducción
efectiva nos basamos en el porcentaje de variabilidad explicada, ya que puede servirnos para
1.8. APLICACIÓN PRÁCTICA DEL ACP 43

alcanzar los objetivos anteriores.

Este método suele usarse como una reducción previa; es decir, como una técnica exploratoria
necesaria para interpretar luego los datos ayudándonos de otras técnicas.

Datos
Esta sección introduce los datos que usaremos en los siguientes ejercicios prácticos.

1. El conjunto de datos chicken contiene datos de expresión génica de 7406 genes (varia-
bles) medidos en 43 pollos (muestras). También hay una variable categórica que contiene
6 dietas en función de distintos niveles.

2. breastCancer es un conjunto de datos de expresión génica (microarrays, tipo ’hgu95-


A’) obtenidos en un estudio donde se ha analizado el efecto de diferentes tratamientos y
tiempos de exposición en niveles de expresión génica de un grupo de mujeres afectadas
por cáncer de pecho. El estudio está disponible en la base de datos GEO (Gene Expression
Omnibus). Para este estudio nos hemos basado tan sólo en 18 muestras, guardadas en el
fichero Breast-Cancer.txt.

Paquetes
1. FactoMineR
Es un paquete de R que permite realizar Análisis Exploratorio de Datos Multivariante
y Minería de Datos. En concreto, tiene implementados métodos para realizar Análi-
sis de Componentes Principales, Análisis de Correspondencias, y Análisis de Corres-
pondencias Múltiples, así como otros métodos estadísticos multivariantes avanzados.
Destacamos que con él se pueden obtener gráficos de gran calidad. Este paquete ha
sido desarrollado por Husson; Josse, Lê, y Mazet (2007). La última versión es Facto-
mineR 1.31.5 publicada el 07/01/2016, y que se encuentra disponible en https://cran
.rproject.org/web/packages/FactoMineR/index.html.

1.8.1. Exploración de datos


Cargamos el paquete FactoMineR y construimos el data.frame “poulet” a partir del fichero
“chicken.txt”.

También calculamos la dimensión de la matriz y obtenemos como resultado 7406 filas y


43 columnas.

Dado que en la estructura de un data.frame las columnas se corresponden con las


variables y las filas con los individuos, trasponemos la matriz de partida y la renombramos;
obteniendo 43 observaciones de 7406 variables.
44 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

require(FactoMineR)
poulet <- read.table(file.path("data", "chicken.txt"),
header=T,sep="\t", dec=".", row.names=1)
dim(poulet)

## [1] 7406 43

poulet=as.data.frame(t(poulet))
dim(poulet)

## [1] 43 7406

Creamos una nueva variable categórica definiendo 6 dietas y, convirtiéndola previamente


en un factor, la añadimos como columna con la instrucción cbind.data.frame a
nuestro data frame “poulet”.
Con la función colnames especificamos el nombre de la variable.

poulet =cbind.data.frame(as.factor(c(rep("N",6),
rep("J16",5), rep("J16R5",8), rep("J16R16",9),
rep("J48",6), rep("J48R24",9))), poulet)
colnames(poulet)[1] = "Diet"

Ahora, procedemos a hacer el ACP, indicando que representamos los individuos con
diferente color basándonos en nuestra variable categórica “Dieta”; que es la primera, y que
no introducimos gráficos.
Hemos convertido la matriz de datos al formato data.frame para poder ahora aplicar la
función PCA.

res.pca = PCA(poulet, quali.sup=1, graph=F)

El comando summary nos dará información sobre el ACP.

summary(res.pca)

Mostraremos los resultados más significativos.


En primer lugar aparecen los autovalores (varianza, porcentaje de varianza explicada y
porcentaje de varianza acumulado) para las 42 CP que considera.
1.8. APLICACIÓN PRÁCTICA DEL ACP 45

Eigenvalues
Dim.1 Dim.2 Dim.3 Dim.4
Variance 1453.572 692.788 536.208 434.453
% of var. 19.627 9.354 7.240 5.866
Cumulative % of var. 19.627 28.981 36.222 42.088

Dim.41 Dim.42
Variance 38.259 33.786
% of var. 0.517 0.456
Cumulative % of var. 99.544 100.000

Para los primeros individuos se calcula la distancia a la categoría a la que pertenecen,


además de los valores de y (j) , j = 1, 2, 3, que aparecen en las columnas “Dim.1”, “Dim.2”
y “Dim.3”; respectivamente.

Individuals
Dist Dim.1 ctr cos2 Dim.2
N_1 | 72.888 | 6.168 0.061 0.007 | -17.356
N_2 | 69.634 | 15.472 0.383 0.049 | -14.354
N_3 | 70.814 | 16.951 0.460 0.057 | -14.533
j16_3 | 82.817 | 8.092 0.105 0.010 | 9.506
j16_4 | 74.501 | 28.090 1.262 0.142 | 4.704

ctr cos2 Dim.3 ctr cos2


N_1 1.011 0.057 | -47.024 9.590 0.416 |
N_2 0.692 0.042 | -37.041 5.951 0.283 |
N_3 0.709 0.042 | -41.914 7.619 0.350 |
j16_3 0.303 0.013 | 43.142 8.072 0.271 |
j16_4 0.074 0.004 | 41.583 7.499 0.312 |

También se indica la contribución de cada individuo en la determinación de cada CP


(“ctr”), así como la calidad de la representación de estos (“cos2”).

Para las variables se obtienen los primeros autovectores (las columnas “Dim.1” y “Dim.2”
se corresponden con ê1 y ê2 ). También aparecen la contribución de cada variable en una
CP determinada y la calidad de representación o correlación al cuadrado entre ellas; como
hemos comentado en ejemplos anteriores.
46 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

Variables
Dim.1 ctr cos2 Dim.2 ctr cos2
A4GALT | -0.526 0.019 0.276 | -0.060 0.001 0.004
A4GNT | -0.124 0.001 0.015 | 0.217 0.007 0.047
AACS | 0.332 0.008 0.110 | 0.036 0.000 0.001
AADACL1 | 0.253 0.004 0.064 | 0.284 0.012 0.081
AADACL2 | 0.510 0.018 0.260 | -0.486 0.034 0.236

Por último, para cada categoría de la variable “Dieta” aparece la distancia al centro de
los ejes de la representación gráfica. Cabe señalar que cada categoría se proyecta en el
baricentro de los individuos que pertenecen a la misma.

Aparecen también las coordenadas de las categorías para las primeras CP, la correlación
entre ellas y un test de normalidad “v.test”; considerando como variable de agrupación
la variable categórica. En este test se basan las elipses de confianza que se mostrarán
posteriormente.

Supplementary categories
Dist Dim.1 cos2 v.test
J16 | 55.795 | 26.069 0.218 1.607
J16R16 | 42.365 | 34.124 0.649 2.984
J16R5 | 41.933 | 18.534 0.195 1.506
J48 | 79.310 | -73.356 0.855 -5.021
J48R24 | 48.622 | -25.591 0.277 -2.238
N | 46.001 | 14.122 0.094 0.967

Dim.2 cos2 v.test


J16 | 7.307 0.017 0.653
J16R16 | -7.524 0.032 -0.953
J16R5 | -13.613 0.105 -1.602
J48 | -22.706 0.082 -2.251
J48R24 | 36.262 0.556 4.594
N | -8.340 0.033 -0.827

En el siguiente gráfico aparecen las proyecciones de los individuos que resultan del ACP
así como las distintas categorías del factor “Dieta”. Cada individuo se representa con un
color distinto en función de la dieta que siga.
1.8. APLICACIÓN PRÁCTICA DEL ACP 47

Se representan las proyecciones porque, como sabemos por teoría, las distancias al cuadra-
do entre elementos originales y transformados se conservan (d2 (y h , y i ) = d2 (xh , xi )).

plot(res.pca, choix="ind", habillage=1)

Individuals factor map (PCA)

J16
J16R16
100

J16R5
J48 j48r24_9
J48R24
N
j48r24_5

j48r24_8
j48r24_7
50
Dim 2 (9.35%)

J48R24

j48r24_3
j16_5
j48r24_2 N_7
j48r24_6 j16_3 j16_6
j48r24_1 j16_4
N_6 J16
j16r16_1
j16r16_3 j16_7
j16r5_7
0

j16r16_8
j16r16_6
j48_7 j16r5_1 j16r16_2
N j16r5_8 J16R16
j48r24_4 J16R5 j16r16_5 j16r16_7
j48_1 j16r5_2 j16r5_6 j16r16_9
N_2 N_3 j16r5_5
J48 N_1
j48_6 j16r5_3 N_4 j16r16_4
j48_3 j16r5_4
j48_4

j48_2
−50

−100 −50 0 50

Dim 1 (19.63%)

Podemos ver que los ejes “Dim 1” y “Dim 2” son las CP1 y CP2 respectivamente, y
aparece indicado el porcentaje de variabilidad explicado por cada una de ellas.
La primera CP distingue los individuos con las dietas “j48” y “j48R24” del resto; mientras
que la segunda CP distingue los individuos que siguen las dietas “j16” y “j48R24” de los
restantes individuos de la muestra.

A continuación, representamos en un círculo de correlación la relación entre las variables


iniciales y las dos primeras CP, teniendo en cuenta sólo aquellas con correlación igual o
48 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

superior a 0,8.

plot(res.pca, choix="var", lim.cos2.var=0.8)

Variables factor map (PCA)


1.0

GTF2H1
TTLL5
MAS1
0.5

ELL2 GSPT1
F7 MAP1LC3C
MLF1
Dim 2 (9.35%)

0.0

C16orf48 GRIA3
TTC151
PRC1 TTC8
KLHL101
C6orf66 C7orf30
HS2ST1HSPA8
UBXD7 PKD1
CHRNA4 KCNJ15
CD72 C6orf601
ITPA NDUFA2 TCTE3
KIAA1524 MRPS17
TNNI1 STK38
PCNXL21
CCNH CNNM11
SOCS7 CHMP71 MAGI11
KCNAB1
TRIM24
−0.5

BAI3 BTBD121
−1.0

−1.0 −0.5 0.0 0.5 1.0

Dim 1 (19.63%)

Como ya vimos en este capítulo, variables próximas en el círculo de correlación están


muy correladas entre sí, si son ortogonales no están correladas y si están en lados opuestos
están negativamente correladas de forma significativa. Además, cuando las variables están
cercanas al centro, la información está en otros ejes (componentes principales).

Por último, visualizamos las elipses de confianza de las clases de la variale categórica
tras el ACP con el comando plotellipses(). Se calculan por defecto a un nivel de
confianza 0,95 suponiendo la normalidad de los datos originales.
1.8. APLICACIÓN PRÁCTICA DEL ACP 49

plotellipses(res.pca)

Confidence ellipses around the categories of Diet

J16
J16R16
100

J16R5
J48 j48r24_9
J48R24
N
j48r24_5

j48r24_8
j48r24_7
50
Dim 2 (9.35%)

J48R24

j48r24_3
j16_5
j48r24_2 N_7
j48r24_6 j16_3 j16_6
j48r24_1 J16 j16_4
N_6 j16r16_1
j16r16_3 j16_7 j16r5_7
0

j16r16_8j16r16_7
j48_7 j16r5_1 j16r16_2
N j16r5_6 J16R16j16r16_6
j48r24_4 j16r16_5 j16r5_8
N_2 j16r5_2 j16r16_9
J16R5
j48_1 N_1 N_3 j16r5_5
j48_6
J48
j48_3 j16r5_3 N_4 j16r16_4
j48_4 j16r5_4

j48_2
−50

−100 −50 0 50

Dim 1 (19.63%)

Estas elipses son útiles para comparar categorías, decidir si una observación procede o
no de una determinada población con distribución normal bivariante y además, detectar
puntos anómalos.
En este caso, las elipses de las categorías “J48R24” y “J48” no se solapan, lo que significa
que estas son significativamente diferentes. Sin embargo, las elipses correspondientes a las
4 categorías restantes muestran que no existen diferencias significativas entre ellas.
50 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

1.8.2. Detección de efecto lote en experimentos ómicos


Los datos producidos en estudios de tipo microarrays pueden confundirse por el conocido
como efecto lote. Este efecto se debe al procesamiento de los datos en distintos lotes (días,
platos, operarios,...) y puede crear confusión si oculta el efecto de los tratamientos objeto
de estudio.

Cuando hay diferentes condiciones experimentales, si individuos de cada grupo están bien
equilibrados entre los lotes el problema no es grave (o bien se compensa o bien se puede
eliminar con técnicas ANOVA).
Sin embargo, si lotes y grupos están mezclados no se sabe si las diferencias que detectemos
se deben a condiciones experimentales diferentes o a haber sido procesados en diferentes lo-
tes. Esto conlleva que los datos puedan perder su valor y el estudio estadístico no sea válido.

A continuación, se recoge un ejemplo en el que utilizamos ACP para detectar el efecto


lote. Trabajaremos con el conjunto “Breast_Cancer.txt” que, tal y como especificamos en
el apartado datos de esta sección 1.8, contiene datos de expresión génica medidos en 18
mujeres con cáncer de pecho.

bcData <- read.table(file.path("data",


"Breast_Cancer.txt"), head=T, sep="\t", row.names=1,
as.is=TRUE)
class(bcData)

## [1] "data.frame"

dim(bcData)

## [1] 18 12630

Podemos observar que “bcData” es un data.frame con 18 filas y 12630 columnas.

Se ha utilizado la función substr, que permite seleccionar determinadas posiciones de


una cadena de caracteres; y hemos recortado los nombres de las filas de nuestra matriz
para facilitar la lectura de las salidas.

rownames(bcData)<- substr(rownames(bcData),4,8)

No vamos a considerar dos de las primeras filas de la matriz, pues son datos de control
que carecen de interés en nuestro estudio. Por tanto, nuestro data.frame será de dimensión
1.8. APLICACIÓN PRÁCTICA DEL ACP 51

16 × 12630.

Además, vamos a considerar las 3 primeras variables como categóricas (Tratamiento,


Tiempo y Lote).

bcData<- bcData[-(1:2),]
head(names(bcData))

## [1] "Treatment" "Time"


## [3] "Batch" "Treatment.Combination"
## [5] "X100_g_at" "X101_at"

for(i in 1:3) bcData[,i]<- as.factor(bcData[,i])

Como hemos convertido las variables Tratamiento, Tiempo y Lote en factores, la función
summary nos muestra una tabla de frecuencia con cada una de sus categorías.
La variable “Batch” que indica el lote en que se procesaron los datos, es la que nos
interesará para detectar si los datos son válidos o no.

summary(bcData$Treatment)

## E2 E2+ICI E2+Ral E2+TOT


## 4 4 4 4

summary(bcData$Time)

## 8 48
## 8 8

summary(bcData$Batch)

## A B
## 8 8

Una vez visto esto, llevamos a cabo el ACP sin incluir la variable “Treatment.Combination”,
pues aporta información redundante.

res.pca = PCA(bcData[,-4], quali.sup=1:3, graph=F)


summary(res.pca)

Con el comando summary, hemos obtenido un resumen del ACP.


52 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

En particular, información sobre los autovalores, los individuos y las variables. Además;
información sobre las 3 variables categóricas y cada clase (E2, E2+ICI, E2+Ral, E2+TOT,
Time8, Time48, A y B).

Mostraremos tan sólo los resultados más relevantes sobre algunos de los individuos y
variables. La información obtenida es la misma que en la aplicación anterior de ACP (datos
“poulet”).

Eigenvalues
Dim.1 Dim.2 Dim.3 | Dim.15
Variance 2145.25 1426.41 1249.83 | 452.18
% of var. 17.03 11.32 9.92 | 3.59
Cumulative % of var. 17.03 28.35 38.27 | 100.00

Individuals
Dist Dim.1 ctr cos2 Dim.2
99 | 124.565 | 60.994 10.839 0.240 | 10.537
38 | 116.159 | -41.266 4.961 0.126 | -59.426
39 | 133.028 | 89.899 23.545 0.457 | -57.436

ctr cos2 Dim.3 ctr cos2


99 0.486 0.007 | -20.336 2.068 0.027 |
38 15.473 0.262 | -32.974 5.437 0.081 |
39 14.454 0.186 | -21.334 2.276 0.026 |

Variables
Dim.1 ctr cos2 Dim.2 ctr cos2
X100_g_at | 0.364 0.006 0.132 | -0.018 0.000 0.000 |
X101_at | -0.199 0.002 0.040 | -0.459 0.015 0.210 |
X102_at | 0.768 0.027 0.589 | 0.004 0.000 0.000 |

Dim.3 ctr cos2


X100_g_at -0.343 0.009 0.118 |
X101_at 0.286 0.007 0.082 |
X102_at -0.182 0.003 0.033 |

Supplementary categories
Dist Dim.1 cos2 v.test Dim.2
E2 | 63.713 | 24.946 0.153 1.204 | -47.038
Time 8 | 33.218 | -17.948 0.292 -1.501 | 10.387
A | 41.098 | 31.590 0.591 2.642 | 19.788
B | 41.098 | -31.590 0.591 -2.642 | -19.788
1.8. APLICACIÓN PRÁCTICA DEL ACP 53

cos2 v.test Dim.3 cos2 v.test


E2 0.545 -2.785 | -25.005 0.154 -1.582 |
Time 8 0.098 1.065 | -12.985 0.153 -1.423 |
A 0.232 2.029 | -16.493 0.161 -1.807 |
B 0.232 -2.029 | 16.493 0.161 1.807 |

Es sencillo observar que las 2 primeras CP explican un 28.347 % de la variabilidad total.

Representamos los resultados del ACP obtenidos para los individuos, dando diferentes
colores para aquellos que fueron procesados en el lote A de los que fueron procesados en
el B.
54 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES

plot(res.pca, choix="ind", habillage=3)

Individuals factor map (PCA)

A
B
15904
50

15902
15900 15908

E2+ICI A
E2+Ral
Time 8 15906
15909
E2+TOT 13099 15910
15903
Dim 2 (11.32%)

15905
15901 15907 Time 48
15911
B
−50

E2

13138 13139

13140
−100

−50 0 50 100

Dim 1 (17.03%)

Los individuos que se corresponden con el lote A están representados en color negro, y los
que se corresponden con B en color rojo.
Además, las clases de cada una de las 3 variables categóricas que hemos considerado están
representadas en color rosa.
1.8. APLICACIÓN PRÁCTICA DEL ACP 55

Añadimos ahora las elipses de confianza:

aa=cbind.data.frame(bcData[,3],res.pca$ind$coord)
bb=coord.ellipse(aa, bary=TRUE)
plot.PCA(res.pca, habillage=3, ellipse=bb)

Individuals factor map (PCA)

A
B
15904 15902
50

15908

15900
E2+ICI
E2+Ral 15906 A
Time 8 13099
15903 15909 E2+TOT 15910
Dim 2 (11.32%)

15901 15905
15907 Time 48
15911
B
−50

E2

13138 13139

13140
−100

−50 0 50 100

Dim 1 (17.03%)

Podemos ver que las elipses no se solapan, lo que indica que existen diferencias signi-
ficativas entre los datos del lote A y el lote B. Esto significa, por tanto; que hay efecto
lote.
Capítulo 2

Descomposición en Valores Singulares

La Descomposición en Valores Singulares (DVS) es un método algebraico que permite


descomponer una matriz rectangular como producto de otras tres matrices. Es una de las
aproximaciones que se utiliza en Análisis Multivariante cuando la dimensión del conjunto
de variables con el que se trabaja supone un problema.
Por consiguiente, puede decirse que es un paso intermedio que se utiliza en multitud de
técnicas estadísticas.

Así lo hemos utilizado en Análisis de Correspondencias, para descomponer el estadístico


χ2 y poder obtener las coordenadas principales y standard de los elementos de la tabla de
contingencia.
Del mismo modo, este método se ha empleado en Análisis de Correlación Canónica para
obtener las variables de correlación canónica a partir de la matriz ΣXX −1/2 ΣXY ΣYY −1/2 .

2.1. Descomposición
Comenzaremos definiendo la descomposición en valores singulares de una matriz X.

Definición 2.1.1 Sea X una matriz de orden n × p y rango p ≤ n, entonces la descompo-


sición en valores singulares de X es:

X = UDV0
donde U y V son matrices de dimensión n × p y p × p; respectivamente, y D es una matriz
de dimensión p × p conteniendo los valores singulares de X.

56
2.1. DESCOMPOSICIÓN 57

En U se recogen los denominados vectores singulares izquierdos de X y en V los vectores


singulares derechos. Ambas matrices son ortogonales, es decir:

U0 U = UU0 = I

V0 V = VV0 = I

Esta ortogonalidad que caracteriza a las matrices U y V nos será de gran utilidad para
explicar otras propiedades de la descomposición en valores singulares para una matriz
rectangular.

La matriz D se define como:

 
α1 0 ··· 0
 0 α2 ··· 0 
D= , siendo α1 ≥ α2 · · · ≥ αp los valores singulares de X.
 
.. .. .. ..
 . . . . 
0 0 · · · αp

2.1.1. Propiedades de la descomposición


En el siguiente Lema se recogen las principales propiedades de las matrices que resultan
de la descomposición. Nos serán de gran utilidad para las secciones siguientes.

Lema 2.1.1 En las condiciones anteriores, se tiene que:

a) Vp×p es la matriz de autovectores de X0 X

b) Un×p es la matriz de autovectores de XX0

2
c) Los cuadrados de los valores singulares de X, α12 , α22 , · · · , αp , son los autovalores
de X0 X y XX0 .

d) Si X es simétrica y n=p, las matrices U y V son iguales.

Demostración:

a) Sea X = UDV0 la descomposición en valores singulares de X. Por la ortogonalidad


de la matriz U, podemos expresar el producto X0 X como:

X0 X = VDU0 UDV0 = VD2 V0


58 CAPÍTULO 2. DESCOMPOSICIÓN EN VALORES SINGULARES

Multiplicando a derecha en ambos lados de la igualdad por V y por la ortogonalidad


de esta matriz, obtenemos:

X0 XV = VD2
Esta expresión satisface la definición de autovector para X0 X, con V el conjunto de
autovectores de esta matriz y D2 sus autovalores.

Luego queda probado que V es el conjunto de autovectores de X0 X.

b) Análoga al caso anterior.

c) Siguiendo el razonamiento anterior, las matrices X0 X y XX0 pueden expresarse


como:

X0 X = VDU0 UDV0 = VD2 V0

XX0 = UDV0 VDU0 = UD2 U0


Teniendo de nuevo en cuenta la definición de autovector, es sencillo ver que D2 es la
matriz de autovalores para ambas matrices.

Por tanto, si denotamos por λ1 , . . . , λp los autovalores de XX0 y X0 X, puede verse


que los valores singulares de X recogidos
p en la matriz D son equivalentes la raíz de
estos autovalores. Es decir, αj = + λj , con j = 1, . . . , p.

d) X es simétrica si y sólo si X0 = X

Si es así, también se cumplirá lo siguiente:

X0 X = XX0 ⇔ VD2 V0 = UD2 U0

Queda probado por tanto que en este caso los vectores singulares izquierdos y
derechos para X son los mismos.

En el siguiente Corolario se especifica la relación entre la Descomposición en Valores


Singulares y el Análisis de Componentes Principales.
2.1. DESCOMPOSICIÓN 59

Corolario 2.1.1 a) Las Componentes Principales de X0 X están relacionadas con los


autovectores de XX0

b) Las Componentes Principales de XX0 están relacionadas con los autovectores de


X0 X.

Demostración:

a) Las componentes principales de X0 X están dadas por Z = XV, pues como hemos
visto en el primer capítulo de la memoria, las combinaciones lineales de máxima
varianza de las variables originales se obtienen con los autovectores y V es el
conjunto de autovectores de X0 X.

Volviendo a hacer uso de la ortogonalidad de V, podemos expresar Z como:

Z = UDV0 V = UD

Estas CP son, por tanto, una versión reescalada de U; vectores singulares izquierdos
de X y autovectores de XX’.
b) Análoga al caso anterior.

2.1.2. Aproximación matricial


En este apartado explicaremos cómo aproximar la matriz inicial por una de menor dimen-
sión sin perder información relevante.

Si queremos trabajar con una matriz Xb aproximación de menor rango de la matriz X,


podemos calcularla mediante la aproximación de mínimos cuadrados.

El problema a resolver sería:

min(tr[(X − X)(X
b b 0 ])
− X)

La matriz que buscamos es X b = UDr V0 = Pr αj uj v 0 , siendo Dr la matriz diagonal


j=1 j
con los r ≤ p primeros valores singulares de X.

Por las propiedades de la descomposición especificadas con anterioridad, la solución al


problema puede expresarse también como X b = Zr V0 , con Zr las r primeras componentes
r
0 0
principales de X X y Vr sus r primeros autovectores.
60 CAPÍTULO 2. DESCOMPOSICIÓN EN VALORES SINGULARES

2.2. Biplots y matriz de aproximación


En este apartado, mostramos cómo la aproximación de la matriz X puede representarse en
un biplot.
Recordemos que un biplot es una representación gráfica que combina filas y columnas
de una matriz de datos, pudiendo así interpretar la distancia entre los mismos y obtener
conclusiones.

Sea X̂ = Ur Dr Vr0 , las columnas de Vr dan información sobre las columnas o variables,
mientras que las columnas de Ur dan información sobre las filas o muestras de X.

Esto se debe a la relación entre componentes principales y autovectores de una matriz que
vimos en la sección anterior.
Si Z = UD son las componentes principales de la matriz X0 X, entonces la matriz X se
expresará como:
X = UDV0 = ZV0
y por tanto podremos obtener información sobre las filas de la matriz X.
Del mismo modo, si W = VD son las componentes principales de la matriz XX0 ,
entonces a partir de X0 = WU0 podremos obtener información sobre las columnas de X.
Cada componente principal definirá un eje y los inviduos serán puntos dispersos.

Para poder visualizar el problema, normalmente sólo se consideran dos CP.


Capítulo 3

Análisis de Correspondencias

El Análisis de Correspondencias (AC) es una técnica iniciada por J.P. Benzécri en 1930 que
analiza las asociaciones entre filas y columnas de una tabla de contingencia. Ésta se define
como una tabla con dos entradas donde se muestran las frecuencias de cada una de las
clases de dos variables cualitativas. En general, consideraremos tablas de dimensión (n×p).

El Análisis de Correspondencias tiene como objetivo principal desarrollar índices simples


que muestren las relaciones entre las categorías de las filas y las columnas, indicando qué
categorías columnas tienen más peso en una categoría fila y viceversa.
A su vez, es útil para reducir la dimensión de la tabla de contingencia. Para ello, se calcu-
larán los índices en orden decreciente de importancia, de manera que podamos resumir
esta tabla sin perder información relevante. Si sólo usamos dos índices, podremos mostrar
los resultados en gráficos de dos dimensiones.

3.1. El método
Este método fue desarrollado para variables cualitativas exclusivamente.
Cada entrada de la tabla es el número de observaciones de la muestra que se correspon-
de simultáneamente con la i-ésima categoría fila y la j-ésima categoría columna (para
i = 1, ..., n y j = 1, ..., p).

Una vez que tenemos establecido el tipo de variables que medimos, podemos plantearnos
dos tipos de análisis para los datos: homogeneidad de poblaciones o independencia de
caracteres.

61
62 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

•Homogeneidad de poblaciones

Se considera una variable categórica B con p modalidades.


El problema se estudia para n poblaciones, siendo la hipótesis de homogeneidad:

H0 : P1 (Bj ) = · · · = Pn (Bj ), j = 1, · · · , p

Seleccionamos una muestra Ni de cada población y construimos la tabla de contingencia.

Definición 3.1.1 Sea Nij el número de elementos de la muestra tomada de la población


Pi que presentan la categoría Bj . Se tiene que Nij ∼ Bi(Ni , Pi (Bj )).

Se construye la tabla de contingencia como:

B1 ··· Bj ··· Bp

P1 N11 · · · N1j · · · N1p N1 = N1.


.. .. .. .. ..
. . . . .
Pi Ni1 · · · Nij · · · Nip Ni = Ni.
.. .. .. ... .. ..
. . . . .
Pn Nn1 · · · Nnj · · · Nnp Nn = Nn.

N·1 · · · N.j · · · N.p N


Pp Pn
donde Ni = Ni· = j=1 Nij y N·j = i=1 Nij

•Independencia de caracteres

En este problema se consideran dos variables categóricas A y B con n modalidades y p


modalidades respectivamente, siendo la hipótesis de independencia:

H0 : P (Ai ∩ Bj ) = P (Ai )P (Bj ), ∀i, j

Para definir la tabla de contingencia, seleccionamos una muestra de tamaño N y observa-


mos (A, B) en cada elemento.
3.1. EL MÉTODO 63

Definición 3.1.2 Sea Nij el número de elementos de la muestra que presentan Ai ∩ Bj .


Se tiene que Nij ∼ Bi(N, P (Ai ∩ Bj )) y se define la tabla de contingencia como:

B1 ··· Bj ··· Bp Marginal de A

A1 N11 · · · N1j · · · N1p N1.


.. .. .. .. ..
. . . . .
Ai Ni1 · · · Nij · · · Nip Ni.
.. .. .. ... .. ..
. . . . .
An Nn1 · · · Nnj · · · Nnp Nn.

Marginal de B N·1 · · · N.j · · · N.p N


Pp Pn
donde Ni· = j=1 Nij y N·j = i=1 Nij

A continuación, haremos uso del concepto de tabla de contingencia para obtener la que
se conoce como tabla de frecuencias relativas, matriz con la que trabajamos y que es
equivalente para los dos tipos de estudio.
Se mostrarán las frecuencias relativas conjuntas y marginales.
Nij
Definición 3.1.3 Sea fij = N
, se puede definir la siguiente tabla de frecuencias relati-
vas:

f11 · · · f1j · · · f1p f1.


 
 .. .. .. ..
 . . . .

 P p
F =  fi1 · · · fij · · · fip  fi. = j=1 fij
 
 . .. . . .  ..
 .. . . ..  .
fn1 · · · fnj · · · fnp fn.

f·1 · · · f.j · · · f.p


 
Estas frecuencias marginales filas f1. . . . fi. . . . fn. ’ y columnas f·1 · · · f.j · · · f.p
también se conocen como pesos fila y columna, respectivamente.

Como comentamos en la introducción al capítulo, el objetivo de esta técnica es definir la


relación entre filas y columnas de la tabla mediante unos índices. ¿Cómo se miden estas
relaciones entre categorías filas y columnas?.

Podemos estudiarlas gráficamente, representando las categorías e interpretando las posicio-


nes relativas de los puntos en función de los pesos correspondientes a cada fila o columna.
Más adelante en este capítulo veremos cómo interpretar las distancias entre elementos del
64 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

gráfico.

Obtendremos las coordenadas para la representación de cada clase gracias a un sistema de


índices similares al ACP, pero descomponiendo el valor del estadístico “chi-cuadrado” en
lugar de la varianza total.

En la siguiente sección se mostrará cómo obtener este estadístico y su descomposición,


además de algunas propiedades.

3.2. Descomposición Chi-Cuadrado

Para calcular el valor del estadístico χ2 , en primer lugar se estima el valor esperado de
cada elemento de la tabla de contingencia Eij , para luego comparar estos valores con los
valores observados correspondientes Nij .

El valor esperado de cada elemento Nij dependerá del tipo de análisis de los datos.

• Homogeneidad de poblaciones:

Ni· N·j
ÊH0 (Nij ) = Ni· P (Bj ) = N

• Independencia de caracteres:

N·j Ni· N·j


ÊH0 (Nij ) = N P (Ai ∩ Bj ) = N ( NNi· )( N
) = N

Hemos comprobado que el valor de Eij es el mismo para ambos casos. Así, el estadístico
propuesto es:

n X p
X
2 (Nij − Eij )2
χ = (3.1)
i=1 j=1
Eij

Ni. N.j
siendo Eij = ÊH0 (Nij ) = N
.

La razón por la que se divide por Eij en la expresión del estadístico es porque ha de
considerarse la distinta precisión de cada coordenada fila o columna. Debemos tener en
cuenta que cada fila o columna tiene un peso distinto.
3.2. DESCOMPOSICIÓN CHI-CUADRADO 65

3.2.1. Descomposición y propiedades


Para calcular los índices que permiten interpretar las relaciones entre clases, comenzaremos
definiendo una matriz a partir del estadístico χ2 .

(Nij −Eij )
Sea C una matriz de elementos cij = 1/2 , calcularemos su descomposición en valores
Eij
singulares.

Esta descomposición se define como

C = ΓD∆0 ,

siendo Γ(n×p) la matriz de autovectores de CC0 , ∆(p×p) los de C0 C y D(p×p) la matriz


diagonal de valores singulares αj , j ∈ 1, . . . , p.

Estos p valores singulares cumplen αj = + λj , con λj los autovalores de CC0 y C0 C,


p

como vimos en el Lema 2.1.1 del segundo capítulo de esta memoria.

Tras la descomposición de la matriz C, podemos deducir:

R
X 1/2
cij = λk γik δjk (3.2)
k=1

con γik elemento de la matriz Γ , δjk elemento de ∆ y R = rango(C) ≤ mı́n{(n − 1), (p − 1)}.

Por la normalidad de los autovectores contenidos en Γ y ∆ y la ecuación (3.2), se satisface


la siguiente igualdad:
R p
n X
X X
0
tr(CC ) = λk = c2ij = χ2
k=1 i=1 j=1

Se demuestra así que el estadístico χ2 definido en (3.1) puede expresarse como la suma de
los R autovalores no nulos de la matriz CC0 .

Denotemos por A = diag(N1· , . . . , Nn· ) y B = diag(N·1 , . . . , N·p ). Es sencillo compro-


bar que se cumple
( p p
C (N·1 , ..., N.p )0 = C B1p = 0
p √
C0 (N1. , ..., Nn. )0 = C A1n = 0
66 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

Una vez descompuesto el estadístico χ2 , los autovectores δ k y γ k son objeto de interés


para calcular los índices que permiten analizar la correspondencia entre filas y columnas.
Por las propiedades de la descomposición en valores singulares, sabemos que los autovec-
tores γ k proporcionan información sobre las filas de la tabla de contingencia y δ k sobre las
columnas.
A partir de esto, obtenemos el siguiente resultado.

Proposición 3.2.1 Las relaciones de dualidad entre el espacio fila y columna (para
k = 1, ..., R) están dadas por:
(
δ k = √1λk C0 γ k ,
γk = √1 Cδ
λk k

mientras que las proyecciones de las filas y las columnas de C están dadas por:
( √
Cδ k = λk γ k ,

C0 γ k = λk δ k

Las aproximaciones de χ2 que se propongan dependerán de la información que podamos


recoger con los mayores autovalores.

Si suponemos que el primer autovalor es dominante, eso significa que podemos hacer la
1/2
aproximación cij ≈ λ1 γi1 δj1 ; de manera que si δj1 y γi1 toman valores grandes y del
mismo signo; esto indica una asociación positiva entre las categorías fila i y columna j.
Si por el contrario toman valores de signo contrario entonces la asociación será negativa.

3.3. Coordenadas Principales y Coordenadas Standard


Hay dos posibles representaciones de las entradas de la tabla de contingencia, en función
de si asumimos o no que existe relación entre los datos.

Cuando partimos de que filas y columnas en la tabla están relacionadas, representaremos


estas filas y columnas utilizando las llamadas coordenadas principales. La solución que
obtenemos se conoce como una solución simétrica del problema.

En el caso contrario, usaremos las denominadas coordenadas standard para las columnas,
obteniéndose así lo que se conoce como una solución asimétrica.
Comenzaremos definiendo las Coordenadas Principales.
3.3. COORDENADAS PRINCIPALES Y COORDENADAS STANDARD 67

Definición 3.3.1 Sean rk(n×1) las proyecciones de A−1/2 C en δ k y sk(p×1) las proyecciones
de B−1/2 C0 en γ k .
Los vectores rk y sk se conocen como coordenadas principales de filas y columnas,
respectivamente. ( √
rk = A−1/2 Cδ k = λk A−1/2 γ k

sk = B−1/2 C0 γ k = λk B−1/2 δ k
con k = 1, . . . , R.

Proposición 3.3.1 Estos vectores satisfacen las siguientes igualdades:


 0
rk (N1. , ..., Nn. )0 = 0
s0k (N·1 , ..., N.p )0 = 0

A partir de las ecuaciones que definen las coordenadas principales rk y sk , sabemos que:
r0k Ark = λk , s0k Bsk = λk

De la relación de dualidad entre los autovectores de la Descomposición singular obtenemos:

r
1 N −1
rk = √ A−1/2 CB1/2 sk = A Xsk
λk λk
r
1 −1/2 0 1/2 N −1 0
sk = √ B C A rk = B X rk
λk λk

Proposición 3.3.2 Se satisfacen las relaciones:


p
1 X Nij
ri = s
λi j=1 j Ni.
n
1 X Nij
sj = r
λj i=1 i N.j

En la siguiente proposición se detallan la media y varianza de las Coordenadas Principales.

Proposición 3.3.3 La media de los vectores rk y sk es:


1 0
rk = r (N1. , ..., Nn. )0 = 0
N k
1
sk = s0k (N·1 , ..., N.p )0 = 0
N
La varianza de estos es:
68 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

n
1 X r0 Ar λk
V ar(rk ) = 2
Ni. rki = k k = , (3.5)
N i=1 N N
p
1 X 2 s0k Bsk λk
V ar(sk ) = N.j skj = = (3.6)
N j=1 N N

Para la solución asimétrica del problema necesitamos definir las Coordenadas Standard.

Definición 3.3.2 Se definen las Coordenadas Standard como:


(
rks = A−1/2 γ k
sks = B−1/2 δ k

Por último, haremos uso de las contribuciones para evaluar el peso de cada fila (o columna)
en las varianzas de los factores dadas en (3.5) y (3.6).

Definición 3.3.3 La contribución de la fila i a la varianza del factor rk se define como:


2
Ni. rki
Ca (i, rk ) =
λk
mientras que la contribución de la columna j a la varianza del factor sk es:

N.j s2kj
Ca (j, sk ) =
λk

3.4. Análisis de Correspondencias en la Práctica


Como hemos mencionado anteriormente, la representación gráfica en los k ejes de las n
filas y las p columnas de X viene dada por los elementos rk y sk ; respectivamente.

Las coordenadas filas y columnas pueden ser representados en un gráfico bidimensional. Si


están representados próximos el uno al otro (lejanos al origen), esto indicaría que la i-ésima
categoría fila tiene una alta frecuencia condicional Nij /N.j en la expresión de la j-ésima
categoría columna, y la j-ésima categoría columna tiene una alta frecuencia condicional
Nij /Ni. en la expresión de la i-ésima categoría fila; lo que indicaría una asociación positiva
entre ambas.

Razonando del mismo modo, si ambas categorías estuvieran representadas muy lejos la una
de la otra, la asociación entre ellas sería negativa y por tanto la contribución condicional
3.5. APLICACIÓN PRÁCTICA DEL AC 69

de frecuencia sería muy pequeña.

Normalmente, un gráfico bidimensional explica satisfactoriamente los datos.

Dependiendo del gráfico, la interpretación de los datos se puede resumir como sigue:

• Gráfico para filas de la tabla de contingencia: Se representan los elementos en los ejes
r1 y r2 . La proximidad entre clases de una misma variable categórica indica frecuencias
similares.

• Gráfico para columnas de la tabla de contingencia: Se representan los elementos en los


ejes s1 y s2 . La interpretación de este gráfico es similar al anterior. Proximidad entre clases
de una misma variable categórica indica frecuencias similares.

• Gráfico conjunto: La proximidad entre filas y columnas indica el peso que tiene un
elemento de la matriz sobre el otro. Una fila que está bastante distante de una columna
particular indica que no hay casi observaciones en la columna para esta fila (y viceversa).
Filas y columnas con fuerte asociación se proyectan en la misma dirección desde el origen.
El origen es la media de rk y sk .

Todas las interpretaciones anteriores dependerán de la calidad de la representación gráfica


que se evalua, como en ACP, usando el porcentaje de varianza acumulada.

3.5. Aplicación práctica del AC


Esta es una aplicación del Análisis de Correspondencias simétrico, por lo que asumiremos
que existe relación entre filas y columnas y calcularemos las coordenadas principales para
ambas.

Nos basamos en el porcentaje de variabilidad explicada para elegir el número k de ejes


que recogen la mayor cantidad de información posible.

Datos

Se introducen los datos que usaremos como ilustración de la técnica expuesta.


La base de datos khan está disponible en el paquete made4. El conjunto de datos Khan
contiene 2308 perfiles de expresión génica (individuos) y valores de expresión para 6567
clones, de los cuales 3789 eran genes conocidos y 2778 eran EST (“expressed sequence
70 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

tag”); que se introducen para estudiar la expresión de genes en los 4 tipos de tumores de
células pequeñas azules de la infancia:

a) Neuroblastoma (NB)
b) Rabdomiosarcoma (RMS)
c) Un conjunto de linfomas No Hodgkin, llamado Linfoma Burkitt (BL)
d) Familia de tumores Ewing (EWS)

En este caso, el subconjunto de datos que analizamos está formando por 306 genes para
muestras de 64 pacientes.

Paquetes

El paquete made4 facilita el análisis multivariante de datos de microarrays de expresión


génica. Proporciona un conjunto de funciones que utilizan y extienden funciones gráficas
y estadísticas multivariantes disponibles en ade4. Además, acepta datos de expresión
génica en una amplia variedad de formatos de entrada, incluyendo los de Bioconductor,
AffyBatch, ExpressionSet, marrayRaw y data.frame o matrix.
made4 requiere instalar ade4 y además scatterplot3d. El paquete se encuentra en biocon-
ductor: “http://bioconductor.org/biocLite.R”

En primer lugar realizamos un análisis exploratorio de los datos y posteriormente un


Análisis de Correspondencias Simétrico.

3.5.1. Exploración de los datos


Comenzamos cargando los paquetes necesarios para ejecutar el Análisis de Corresponden-
cias.

library(ade4)
library(made4)

Este subconjunto del conjunto de datos original es una lista con 7 elementos. Con la
función names, vemos los nombres de cada uno de ellos.

data(khan)
class(khan)

## [1] "list"

names(khan)
3.5. APLICACIÓN PRÁCTICA DEL AC 71

## [1] "train" "test"


## [3] "train.classes" "test.classes"
## [5] "annotation" "gene.labels.imagesID"
## [7] "cellType"

Veamos con detenimiendo cada uno de estos elementos.

a) “khan$train” es un data.frame de dimensión 306 × 64, lo que significa que hay 64


muestras para el conjunto de entrenamiento.
b) “khan$test” es un data.frame de dimensión 306 × 25, lo que indica que hay 25
muestras para el conjunto test.
c) “khan$train.classes” es un factor de longitud 64 con diferentes niveles, mostrando los
distintos tipos de cáncer que presentan los individuos del conjunto de entrenamiento.
d) “khan$test.classes” es un factor de longitud 25 con diferentes niveles que indican los
tipos de cáncer de cada uno de los individuos del conjunto test.
e) “khan$annotation” es un data.frame (306 × 8) con 8 anotaciones distintas para los
306 genes que estudiaremos.
f ) “khan$gene.labels.imagesID” es un objeto de tipo carácter de longitud 306 con
información sobre imágenes etiquetadas de los genes que estudiaremos.
g) “khan$cellType” es un objeto de tipo carácter de longitud 64 con información sobre
el tipo de célula de cada elemento del conjunto de entrenamiento. Veremos de qué se
trata esto más adelante.

summary(khan)

## Length Class Mode


## train 64 data.frame list
## test 25 data.frame list
## train.classes 64 factor numeric
## test.classes 25 factor numeric
## annotation 8 data.frame list
## gene.labels.imagesID 306 -none- character
## cellType 64 -none- character

Como hemos comentado anteriormente, trabajamos con conjunto entrenamiento y test en


los que medimos 306 genes, de manera que los individuos de cada conjunto son distintos y
en ningún momento se solapan.

Ahora, veamos “khan$cellType” con más profundidad.


72 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

khan$cellType

[1] "T" "T" "T" "T" "T" "T" "T" "T" "T" "T" "T" "T" "T" "C"
[15]"C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"
[29]"C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C"
[43]"C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "C" "T" "T" "T"
[57]"T" "T" "T" "T" "T" "T" "T" "T"

Al ejecutar esta orden, aparecen los tipos de células de las 64 muestras del conjunto de
entrenamiento. “T” indica que la muestra ha sido obtenida por una biopsia del tumor,
mientras que “C” indica que se ha obtenido a partir de una línea celular.
Es fácil comprobar que en el estudio de los tumores de tipo “RMS” y “EWS” se han usado
ambos métodos, mientas que en los otros dos la muestra se ha obtenido a partir de la línea
celular.

En las siguientes instrucciones hemos especificado que sólo trabajaremos con el conjunto
de entrenamiento y que “k.class” son los tipos de cáncer que presentan las 64 muestras.

k.data<-khan$train
dim(k.data)

## [1] 306 64

k.class<-khan$train.classes

La función overview() da una visión general de los datos. Dibuja un diagrama de cajas,
un histograma y un dendograma de análisis jerárquico.

El dendograma de análisis jerárquico es una herramienta de análisis cluster, método que


permite construir grupos a partir de objetos multivariantes. Para ello, se construye una
matriz que contenga las medidas de similitud; aunque cabe destacar que normalmente las
variables categóricas conducen a valores próximos mientras que las continuas dan lugar a
valores distantes.

Este análisis se puede plantear o bien para clasificar filas o bien para clasificar columnas
de la tabla de contingencia. En este caso se clasifican columnas, que contienen las 64
muestras. Muestras que a priori son indistinguibles pueden tener alguna subclase que
gracias a este análisis podemos identificar.
3.5. APLICACIÓN PRÁCTICA DEL AC 73

Se lleva a cabo un clustering jerárquico aglomerativo, donde se empieza por la partición


más fina posible y se va agrupando. Se usa la correlación de Pearson como medida de
similitud.

overview(k.data,classvec=k.class, labels=k.class)

Histogram
0.0 0.2 0.4 0.6 0.8
Height

RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
NB
RMS
RMS
RMS
RMS
RMS
RMS
RMS
EWS
EWS
RMS
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
NB
NB
NB
NB
NB
NB
NB
NB
NB
NB
NB
boxplot Histogram

30
15000
25

20
Frequency

10000

15

10
5000

0 0
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
EWS
BL−NHL
EWS
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
BL−NHL
NB
NB
NB
NB
NB
NB
NB
NB
NB
NB
NB
RMS
NB
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS
RMS

10

15

20

25

30

35

El dendograma representa las clases de cáncer de cada individuo, agrupadas por similitud.
Se establecen dos grupos, uno de ellos conteniendo en su mayor parte individuos con el
tipo de tumor “RMS” y el otro conteniendo los 3 tipos restantes.

El diagrama de cajas y bigotes representa los valores que toman las variables para cada
individuo en función del primer y tercer cuantil.
74 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

El histograma representa en el eje de abcisas el valor que pueden tomar las variables y en
el eje de ordenadas la frecuencia con que se toma este valor para los 64 individuos.

3.5.2. Análisis de Correspondencias


Procedemos al Análisis de Correspondencias de los datos. Emplearemos la función ord
del paquete made4 para especificar que ha de llevarse a cabo este análisis y el argumento
“type=coa” para indicar que es un análisis de tipo simétrico.

k.coa<-ord(k.data, type="coa")

Analicemos los resultados de esta instrucción más detenidamente.

k.coa$ord

## Duality diagramm
## class: coa dudi
## $call: dudi.coa(df = data.tr, scannf = FALSE, nf = ord.nf)
##
## $nf: 63 axis-components saved
## $rank: 63
## eigen values: 0.1713 0.1383 0.1032 0.05995 0.04965 ...
## vector length mode content
## 1 $cw 64 numeric column weights
## 2 $lw 306 numeric row weights
## 3 $eig 63 numeric eigen values
##
## data.frame nrow ncol content
## 1 $tab 306 64 modified array
## 2 $li 306 63 row coordinates
## 3 $l1 306 63 row normed scores
## 4 $co 64 63 column coordinates
## 5 $c1 64 63 column normed scores
## other elements: N

Los resultados de ordenación contienen una lista de longitud 14. Esta lista incluye 63
autovalores (eig), el vector de pesos columnas de longitud 64 (cw) y el vector de pesos
filas de longitud 306 (lw). Además, las Coordenadas Principales filas rj (li) y columnas
sj (co), así como sus equivalentes normalizadas (l1) y (c1).

En la siguiente instrucción, hemos ejecutado los valores de las 3 primeras Coordenadas


Principales filas y columnas para los dos primeros genes e individuos; respectivamente.
3.5. APLICACIÓN PRÁCTICA DEL AC 75

k.coa$ord$li[1:2,1:3]

## Axis1 Axis2 Axis3


## 25725 0.7001400 -0.08693461 -0.1278178
## 193913 -0.2486831 0.64857220 -0.2381108

k.coa$ord$co[1:2,1:3]

## Comp1 Comp2 Comp3


## EWS.T1 -0.3923287 -0.3255404 -0.02589332
## EWS.T2 -0.3461686 -0.2956390 0.05687684

La suma de los autovalores será equivalente al estadístico chi-cuadrado de la tabla de


contingencia, tal y como vimos en el desarrollo teórico de la técnica.

Veamos el porcentaje de varianza explicada por cada eje j de los 63 que hay en total.

k.coa$ord$eig*100/sum(k.coa$ord$eig)

## [1] 16.947461298 13.683297888 10.207338717 5.931017086


## [5] 4.911896717 3.810337299 3.027912076 3.004793647
## [9] 2.419837082 2.319143755 2.042227963 1.886343953
## [13] 1.841744206 1.695622185 1.585943486 1.437215463
## [17] 1.342418983 1.251040188 1.167146174 1.114077872
## [21] 1.076647536 0.989724516 0.919956542 0.865911339
## [25] 0.812460108 0.789851377 0.755875491 0.736710674
## [29] 0.694142462 0.652011530 0.614248981 0.582822698
## [33] 0.548383426 0.513845412 0.493681883 0.446536643
## [37] 0.440410135 0.425398824 0.416579264 0.398563871
## [41] 0.376038015 0.367570772 0.356010253 0.332082395
## [45] 0.327998743 0.318501534 0.295054657 0.283405148
## [49] 0.257741008 0.241029558 0.227885174 0.225137815
## [53] 0.214519150 0.203081716 0.179324545 0.172844643
## [57] 0.168069405 0.151730905 0.142298857 0.129021644
## [61] 0.115421971 0.106023909 0.008629433

Con cumsum calculamos la varianza acumulada por los primeros ejes:

head(cumsum(k.coa$ord$eig*100/sum(k.coa$ord$eig)))

## [1] 16.94746 30.63076 40.83810 46.76911 51.68101 55.49135


76 CAPÍTULO 3. ANÁLISIS DE CORRESPONDENCIAS

Por tanto, el 30.63076 % de la varianza está explicada por los dos primeros ejes.

Análisis de Correspondencias- Visualización de resultados

Hay muchas funciones en el paquete made4 que permiten visualizar resultados de un


análisis de ordenación. La forma más simple de hacerlo es usando la función plot. Re-
presentaremos los autovalores, variables y casos (microarrays).

• Gráfico de autovalores: en el gráfico superior izquierdo se muestran los autovalores. Sólo


tendremos en cuenta aquellos hasta que comienza a estabilizarse la gráfica, en este caso
los dos primeros.

• Gráfico para muestras: Se representan los individuos en los ejes s1 y s2 , pues son las
columnas de la matriz inicial. Perfiles próximos en la representación indican que existe
similitud entre ellos.

El eje s1 distingue la familia de tumores de tipo RMS de los demás perfiles de expresión
génica, mientras que el eje s2 separa las muestras pertenecientes al tipo Neuroblastoma
(NB) y Linfoma Burkitt (BL) frente al resto.

• Gráfico para variables: Se representan los genes (filas de la matriz inicial) en los ejes r1
y r2 . Genes próximos en la representación indican que existe similitud entre ellos.
Podemos observar los grupos de genes que separa cada una de los dos ejes.

• Gráfico conjunto: Se proyectan filas y columnas de la tabla de contingencia en el mismo


gráfico. Genes y muestras con fuerte asociación se proyectan en la misma dirección desde
el origen. Se revela así la relación existente entre los genes y los tipos de cáncer que se
asocian.
3.5. APLICACIÓN PRÁCTICA DEL AC 77

plot(k.coa, classvec = k.class)


d = 0.5

BL.C1
BL.C7
BL.C5
BL.C8
BL.C6
BL.C3
BL.C2
BL.C4
Eigenvalues
NB.C12
NB.C9
NB.C6
NB.C8
NB.C5
NB.C4
NB.C10
NB.C11
NB.C2
NB.C1
NB.C7 RMS.C10
RMS.C6
NB.C3
RMS.C5
RMS.C7
RMS.C4
RMS.C8
RMS.C3
RMS.C9
RMS.C11
EWS.C2 RMS.T5
RMS.T1
RMS.C2
EWS.C4
EWS.C10
EWS.C11
EWS.C1
EWS.C6 RMS.T2
EWS.C9
EWS.C3 RMS.T9 RMS.T4
RMS.T6 RMS.T11
EWS.C7 RMS.T3
RMS.T8
RMS.T10
EWS.C8
EWS.T2 RMS.T7
EWS.T1
EWS.T4
EWS.T6 EWS.T13
EWS.T19
EWS.T9

EWS.T14 EWS.T3
EWS.T7
EWS.T11
EWS.T15
EWS.T12

d = 0.5 d = 0.5

609663
241412 609663
241412
80109
183337
740604 80109
183337
740604
767183781047
193913 BL−NHL
767183781047
193913
700792
200814 700792
200814

417226 417226 NB

812965 296448 812965 RMS 296448


298062
207274 298062
207274
244618
1409509
882506461425 244618
1409509
882506461425
245330 245330
291756 343867 291756 343867
EWS

122159 122159
814260 814260
770394
377461 754600 39093
345553 770394
377461 754600 39093
345553
52076 1473131
866702 52076 1473131
866702
377731 377731
357031 357031
43733 43733

Figura 3.1: Plot del Análisis de Correspondencias.


A. autovalores
B. proyección de muestras microarrays de pacientes con varios tipos de tumores
C. proyección de genes
D.biplot mostrando genes y muestras
Capítulo 4

Análisis de Coinercia

El Análisis de Coinercia (ACoi) es una técnica de análisis multivariante que permite


identificar tendencias en conjuntos de datos que contienen las mismas muestras, es decir,
que se han recogido sobre los mismos individuos. A su vez, permite reducir la dimensión
de los conjuntos más similares.

Esta herramienta no se limita al análisis de conjuntos que contienen el mismo número de


variables y puede aplicarse en caso de que haya más variables que muestras.

El Análisis de Correspondencias es un paso previo y necesario para su aplicación.

4.1. El Método
Sean Xn×p e Yn×q dos matrices de datos con los mismos n individuos a los que se miden
p y q variables, respectivamente. Nos interesa estudiar la relación entre estos dos conjuntos
de variables.

Como ya sabemos, la relación entre dos variables cuantitativas se mide con el coeficiente
de correlación lineal ρ, que toma valores entre -1 y 1 dependiendo de la robustez y el
sentido de la relación entre ellas.

Para estudiar la relación entre conjuntos de variables mediante el Análisis de Coinercia,


necesitamos definir un índice conocido como coinercia.

Para definir este índice, necesitamos definir previamente una coestructura entre las tablas
de datos iniciales con las que trabajamos. Pasaremos de las matrices de datos originales X

78
4.1. EL MÉTODO 79

e Y a considerar la matriz producto X0 Y de dimensiones p × q.


Esquemáticamente, tiene la siguiente forma:

X Y → X’Y

de manera que la matriz con la que trabajaremos tendrá dimensión p × q y por tanto sus
filas recogerán las p variables que caracterizan a X y sus columnas las q variables que
caracterizan a Y.
También es posible trabajar con las q variables de Y como filas y las p variables de X
como columnas de esta nueva matriz, que en ese caso tendrá dimensión q × p.

El cálculo de la coinercia se basará en los tripletes estadísticos (X, QX , D) e (Y, QY , D),


siendo QX p×p una matriz diagonal representando los pesos de las p variables de X y
QY q×q los pesos de las q variables de Y. Estas matrices diagonales representan por tanto
los pesos de las columnas de las matrices iniciales. Dn×n es una matriz diagonal que
recoge los pesos de los n individuos (filas).

Nótese que hemos definido con la misma matriz D los pesos de las filas para ambas ta-
blas de datos pues aunque se midan distintas variables se trabaja con los mismos individuos.

Cada triplete es resultado de un Análisis de Correspondencias, de forma que las matrices


X e Y serán las tablas χ2 derivadas de los conjuntos de datos originales.

Observación 4.1.1 En las matrices X e Y se recogen los denominados residuos de


Pearson obtenidos para cada casilla como:

(Nij − Eij )
p
Eij

siendo Nij el valor observado en la tabla de contingencia y Eij su valor esperado corres-
pondiente, según se vio en la sección 3.2 del tercer capítulo.

A partir de estos elementos, podemos definir los conceptos de inercia y coinercia.

La inercia es una medida de la variabilidad en los datos y se define como la distancia χ2


entre un elemento y su perfil medio, teniendo en cuenta el peso de cada elemento.

Definición 4.1.1 La expresión de la inercia en el contexto que nos ocupa es:

inX = traza(XQX X0 D)


inY = traza(YQY Y0 D)
80 CAPÍTULO 4. ANÁLISIS DE COINERCIA

Definición 4.1.2 Se define la coinercia entre dos espacios como:

Coin(XY) = traza(XQY X0 DYQY Y0 D)

Si los datos están centrados, la inercia será una suma de varianzas mientras que la coinercia
será una suma de covarianzas al cuadrado.

A continuación, definiremos los pasos que se llevan a cabo para la identificación de ten-
dencias.

Estas tablas de datos X e Y se representan en dos hiperespacios de dimensiones p y q,


respectivamente.
Con un análisis previo, se maximiza la inercia (o variabilidad explicada) de ambos hiperes-
pacios. Esto se lleva a cabo mediante un análisis de corespondencias de cada conjunto de
datos. Como resultado se obtienen dos conjuntos de ejes, uno para cada conjunto de datos.

El primer par de ejes en el que se representarán las relaciones entre ambos conjuntos de
variables se elige maximizando la coinercia, de forma que represente la tendencia más
importante entre los conjuntos de datos.
El segundo par se obtiene con el mismo procedimiento, pero requiriendo también que sea
ortogonal al primer par. Así se va repitiendo el proceso con todos los ejes.

¿Cómo podemos medir la similitud entre ordenaciones?

Para ello utilizaremos el conocido como coeficiente de correlación vector, denominado


coeficiente R-V. Este coeficiente es una extensión multivariante del coeficiente de correla-
ción de Pearson, con la diferencia fundamental de que mide la correlación existente entre
tablas de datos en lugar de entre variables.

Definición 4.1.3 El valor del coeficiente de correlación vector R-V es equivalente al


cuadrado del coeficiente de correlación lineal ρ y tiene la siguiente expresión:

Coin(XY)
p
Coin(XX)Coin(YY)
Toma valores entre 0 y 1, pues se define para matrices semidefinidas positivas. Mientras
más cercano a 1 sea el valor de este coeficiente, mayor será la correlación entre los dos
conjuntos de variables.
4.2. APLICACIÓN PRÁCTICA DEL ACOI 81

4.2. Aplicación práctica del ACoi


En esta sección vemos como el Análisis de Coinercia puede revelar patrones o tendencias
que siguen conjuntos de datos donde el número de variables excede en gran medida del
número de muestras, como es el caso de los análisis de microarrays.

Comenzaremos aplicando Análisis de Correspondencias a los conjuntos iniciales para


identificar gráficamente la relación entre variables de cada matriz de datos.

Una vez calculados los ejes con coinercia máxima, podremos representar gráficamente
la divergencia entre perfiles de expresión génica obtenidos desde diferentes plataformas
microarrays. Los genes que definan las tendencias principales en el análisis podrán identi-
ficarse fácilmente.

Datos

Examinamos dos conjuntos de datos de expresión génica con las mismas 60 líneas celu-
lares, que serán las muestras en nuestro estudio. Estas líneas derivan de pacientes con 9
fenotipos tumorales distintos (leucemia (LE), melanomas (ME), cáncer de ovarios (OV),
de pecho (BR), próstata (PR), pulmón (LU), renal (RE), colon (CO) y del sistema central
nervioso (CNS)).

La expresión génica de estas células fue analizada por los grupos Ross2000 y Staunton2001
en diferentes plataformas microarrays.

Compararemos los resultados obtenidos por la compañía Affymetrix (Staunton et al., 2001)
con los resultados obtenidos usando spotted arrays (o microarrays de dos colores) por la
compañía Stanford (Ross et al., 2000).
Los datos completos de Stanford y Affymetrix pueden encontrarse en http://discover.
nci.nih.gov/datasetsNature2000.jsp y en http:\www-genome.wi.mit.
edu/mpr/NCI60/; respectivamente.

1. El conjunto de datos del grupo Ross contiene un subconjunto de los genes originales de
tamaño 1375. Aquellos con más del 15 % de valores perdidos se han quitado y el resto
se han calculado por el método de los k vecinos más cercanos, considerando 16 vecinos
y la distancia euclídea. Además, el conjunto de datos pre-procesado contiene valores
logarítmicos.

2. El conjunto de datos Staunton contiene valores para un subconjunto de 1517 del total
de genes estudiados. Algunos pasos llevados a cabo en el preprocesamiento de los datos
82 CAPÍTULO 4. ANÁLISIS DE COINERCIA

han sido reemplazar los valores que se diferencien de la media menos de 100 unidades por
100, eliminar genes cuya expresión ha sido invariante a lo largo de las 60 líneas celulares y
seleccionar los subconjuntos de genes que conlleven un cambio mínimo de al menos 500
unidades de diferencia con la media. También se han centrado los datos y se ha aplicado
una transformación logarítmica (base 2).

Ambos conjuntos de datos están disponibles en el fichero NCI60 del paquete made4.

En esta aplicación práctica se han seleccionado 144 variables de las 1375 que contiene la
base de datos original del grupo Ross, y 144 de las 1517 originales del grupo Staunton.

Consideramos por tanto 144 variables distintas para cada laboratorio medidas sobre los
mismos individuos, cada uno de los cuales presenta uno de los 9 tipos de cáncer que
estudiamos.

El Análisis de Coinercia permitirá visualizar los genes con patrones de expresión similares
a través de las dos plataformas consideradas.

Paquetes

Para llevar a cabo este análisis en R, hay que ejecutar la instrucción cia en el paquete
ade4 descrito en ejemplos anteriores.

4.2.1. Análisis de Correspondencias


En este caso, no hay asunción previa de la relación entre dos conjuntos de variables. Nos
interesa integrar estos dos conjuntos de datos y escoger variables simultáneamente.

Por tanto, llevaremos a cabo un análisis asimétrico de los datos.

Cargamos la librería made4, y ejecutamos los datos que necesitamos.

library(made4)
data(NCI60)
class(NCI60)

## [1] "list"

Vemos que los datos están almacenados en una lista.


4.2. APLICACIÓN PRÁCTICA DEL ACOI 83

names(NCI60)

## [1] "Ross" "Affy" "classes" "Annot"

Esta lista contiene los siguientes elementos:

a) “NCI60$Ross” es un data.frame de dimensión 144 × 60. Se corresponde con las


medidas de los 144 genes que interesan para el estudio llevadas a cabo por el
laboratorio Ross.
b) “NCI60$Affy” es un data.frame de dimensión 144 × 60. Se corresponde con medidas
de 144 genes llevadas a cabo por el laboratorio Affymetrix.
c) “NCI60$classes” es una matriz de dimensión 60 × 2 que contiene información sobre
el tipo de cáncer de cada individuo.
d) “NCI60$Annot” es un data.frame de dimensión 144 × 4 conteniendo 4 anotaciones
distintas de los genes medidos en cada una de las muestras.

Con la función summary obetenemos información del análisis llevado a cabo por los dos
laboratorios sobre los mismos 60 individuos.

summary(NCI60)

## Length Class Mode


## Ross 60 data.frame list
## Affy 60 data.frame list
## classes 120 -none- character
## Annot 4 data.frame list

Veamos los datos de ambos grupos para las 4 primeras variables y las 4 primeras muestras.

NCI60$Ross[1:4,1:4]

## BREAST_BT549 BREAST_HS578T BREAST_MCF7


## 484963 0.713 0.544 -1.643
## 510395 -2.293 -2.584 -2.471
## 76539 1.130 1.443 -0.062
## 509732 0.837 0.824 -0.057
## BREAST_MCF7ADRr
## 484963 1.162
## 510395 -2.863
## 76539 0.318
## 509732 0.227
84 CAPÍTULO 4. ANÁLISIS DE COINERCIA

NCI60$Affy[1:4,1:4]

## BREAST_BT549 BREAST_HS578T BREAST_MCF7


## V00594_s_at 0.7131644 0.6079392 -4.40948789
## M60854_at -0.2952515 -0.3567663 0.05465057
## X53331_at 2.1143670 0.3103401 0.29865832
## L19686_rna1_at 0.5050693 0.1812558 0.53823619
## BREAST_MCF7ADRr
## V00594_s_at 0.261889367
## M60854_at 0.001557987
## X53331_at 0.000000000
## L19686_rna1_at -0.027256449

Con la función table se obtiene una tabla de frecuencias para cada clase de cáncer.

table(NCI60$classes[,2])

##
## BREAST CNS COLON LEUK MELAN NSCLC
## 8 6 7 6 8 9
## OVAR PROSTATE RENAL
## 6 2 8

También podemos ver las anotaciones de las primeras variables de cada conjunto:

NCI60$Annot[1:4,]

## Ross.ImageID Symbol.Source Affy.Affy ID


## Hs.418241 484963 MT2A V00594_s_at
## Hs.397609 510395 RPS16 M60854_at
## Hs.365706 76539 MGP X53331_at
## Hs.407995 509732 MIF L19686_rna1_at
## Symbol.AnnAffy
## Hs.418241 MT2A
## Hs.397609 RPS16
## Hs.365706 MGP
## Hs.407995 MIF

Procedamos al Análisis de Correspondencias Asimétrico de los datos obtenidos por el


grupo Ross.
4.2. APLICACIÓN PRÁCTICA DEL ACOI 85

La función ord simplifica el funcionamiento de métodos de ordenación canónicos como Análisis


de Componentes Principales, Correspondencias y Correpondencias Asimétrico. Proporciona una
“envoltura” con la cual se puede llamar a cada uno de esos métodos.
Para llevar a cabo el Análisis de Correspondencias Asimétrico hemos especificado “type="nsc’.

data.coa1 <-ord(NCI60$Ross, type="nsc")

Veamos los resultados de la ordenación:

data.coa1$ord

## Duality diagramm
## class: nsc dudi
## $call: dudi.nsc(df = data.tr, scannf = FALSE, nf = ord.nf)
##
## $nf: 59 axis-components saved
## $rank: 59
## eigen values: 0.007109 0.004584 0.003671 0.002675 0.001829 ...
## vector length mode content
## 1 $cw 60 numeric column weights
## 2 $lw 144 numeric row weights
## 3 $eig 59 numeric eigen values
##
## data.frame nrow ncol content
## 1 $tab 144 60 modified array
## 2 $li 144 59 row coordinates
## 3 $l1 144 59 row normed scores
## 4 $co 60 59 column coordinates
## 5 $c1 60 59 column normed scores
## other elements: N

Esta lista de longitud 14 incluye los 59 autovalores que resultan de la DVS de la matriz inicial,
los pesos de las columnas y de las filas, además de las nuevas coordenadas para filas y columnas.
Ya vimos con profundidad la estructura de los resultados en la aplicación para el Análisis de
Correspondencias del capítulo anterior.

Para el Análisis de Correspondencias asimétrico se trabaja con las Coordenadas Principales y


Coordenadas Standard, como vimos en teoría.

Hemos hallado también los valores de coordenadas filas (rj ) y columnas (sj ) para los primeros
genes e individuos; respectivamente.
86 CAPÍTULO 4. ANÁLISIS DE COINERCIA

data.coa1$ord$li[1:4,1:3]

## Axis1 Axis2 Axis3


## 484963 -0.08471002 -0.02782602 0.013984051
## 510395 0.18763603 -0.00994279 -0.008648991
## 76539 -0.03424693 0.01577408 -0.027110886
## 509732 -0.01198044 0.06460366 0.005213023

data.coa1$ord$co[1:4,1:3]

## Comp1 Comp2 Comp3


## BREAST_BT549 -0.10275027 -0.05317834 -0.03500724
## BREAST_HS578T -0.18105947 -0.09919005 -0.03431942
## BREAST_MCF7 0.11365401 0.02823559 -0.08839012
## BREAST_MCF7ADRr -0.03574366 -0.03563459 -0.04106029

La varianza acumulada por los primeros ejes viene dada por el siguiente resultado:

head(cumsum(data.coa1$ord$eig*100/sum(data.coa1$ord$eig)))

## [1] 17.91786 29.47025 38.72162 45.46232 50.07146 54.35127

Por tanto, el 29.47025 % de la varianza está explicada por los dos primeros ejes.

Ahora, llevamos a cabo un Análisis de Correspondencias Asimétrico de los datos; considerando


sólo los resultados obtenidos por el grupo Affymetrix.

data.coa2 <-ord(NCI60$Affy, type="nsc")

Al igual que con el grupo Ross, en los resultados de ordenación se obtienen autovalores, pesos
filas y columnas, coordenadas filas y columnas, coordenadas normalizadas,etc.
La varianza acumulada viene dada por

head(cumsum(data.coa2$ord$eig*100/sum(data.coa2$ord$eig)))

## [1] 17.54391 27.38023 34.10021 39.07819 43.68462 47.91515

Por tanto, el 27.38023 % de la varianza está explicada por los dos primeros ejes

Análisis de Correspondencias- Visualización de resultados


Hay muchas funciones en made4 para visualizar resultados de un análisis de ordenación. La forma
más simple de ver gráficamente los resultados de ord es usando la función plot. Dibujaremos
los autovalores, variables y casos (microarrays).
Ya vimos la interpretación de estos gráficos con detalle en el capítulo anterior.
4.2. APLICACIÓN PRÁCTICA DEL ACOI 87

plot(data.coa1, classvec = NCI60$classes[,2])


d = 0.1

BREAST_MDAMB435
MELAN_SKMEL5
Eigenvalues BREAST_MDAN

MELAN_UACC62
MELAN_SKMEL28
MELAN_UACC257
LEUK_SR
MELAN_M14
LEUK_MOLT4
MELAN_SKMEL2LEUK_CCRFCEM
MELAN_MALME3M
LEUK_K562
LEUK_RPMI8226
RENAL_SN12C NSCLC_H23LEUK_HL60
BREAST_MCF7
NSCLC_H522
NSCLC_HOP92 NSCLC_H460
MELAN_LOXIMVI COLON_SW620
BREAST_MDAMB231
CNS_SNB19
CNS_SF268 COLON_COLO205
COLON_HT29
COLON_HCT15
CNS_SNB75 PROSTATE_PC3
COLON_HCT116COLON_HCC2998
CNS_U251 BREAST_T47D
OVAR_IGROV1
BREAST_MCF7ADRr
RENAL_CAKI1 COLON_KM12
RENAL_TK10
RENAL_ACHN
CNS_SF539
BREAST_BT549
RENAL_UO31
NSCLC_HOP62
OVAR_OVCAR8 PROSTATE_DU145
RENAL_7860
NSCLC_A549
OVAR_OVCAR5
CNS_SF295 RENAL_A498
NSCLC_H322M
RENAL_RXF393 NSCLC_EKVX
OVAR_SKOV3
NSCLC_H226 OVAR_OVCAR3
BREAST_HS578T OVAR_OVCAR4

d = 0.2

248589 241935
365476
428733
365558
366310
488599
377132
486821 323730
471364
487878
429771
510395
510116
510363
470386 376296
470385
487394487513 162479
486700
488370
108837486215
47359 429288
375834 136590 510534

21822

Figura 4.1: Plot del Análisis de Correspondencias(Ross).


A. autovalores
B. proyección de muestras microarrays (proyecciones columnas)
C. proyección de genes (proyecciones filas).
88 CAPÍTULO 4. ANÁLISIS DE COINERCIA

plot(data.coa2, classvec = NCI60$classes[,2])

d = 0.05

COLON_COLO205
COLON_HCC2998
NSCLC_A549
OVAR_OVCAR4OVAR_OVCAR3
Eigenvalues OVAR_SKOV3
BREAST_MDAMB231
RENAL_RXF393
COLON_KM12
NSCLC_H332M
COLON_SW620
RENAL_A498
OVAR_OVCAR8
RENAL_ACHN LEUK_RPMI8266
COLON_HT29
BREAST_MCF7ADRr
NSCLC_H226 BREAST_T47D
COLON_HCT15
RENAL_TK10 NSCLC_H460
RENAL_UO31 BREAST_MCF7
OVAR_IGROV1
COLON_HCT116
NSCLC_EKVX
OVAR_OVCAR5
LEUK_SR
BREAST_HS578T
RENAL_CAKI1PROSTATE_DU145LEUK_HL60
LEUK_CCRFCEM
LEUK_K562
RENAL_7860
BREAST_BT549 LEUK_MOLT4
NSCLC_H522
CNS_SF539 MELAN_LOXIMVI
CNS_SF268 NSCLC_H23
NSCLC_HOP92
RENAL_SN12C
NSCLC_HOP62 PROSTATE_PC3
CNS_SF295
CNS_SNB75 CNS_U251
CNS_SNB19

BREAST_MDAMB435
MELAN_UACC62
MELAN_MALME3M
MELAN_SKMEL2
BREAST_MDN
MELAN_UACC257
MELAN_M14
MELAN_SKMEL28

MELAN_SKMEL5

d = 0.1
U78095_at

M94250_at
M35878_at L33930_s_at
M59807_at D17793_at
Z24727_at M80563_at
M93036_at
X68314_at
U62015_at
D49824_s_at
J03764_at
M55998_s_at

M92934_at
J04469_at
L19686_rna1_at
Z18951_at
V00594_s_at L07548_at
X92896_at
M83088_at D14811_at
M84349_at
M29277_at
U49785_at
M15841_at
M23294_at HG1612−HT1612_at
L06419_at

J03040_at

Plot del Análisis de Correspondencias(Affy).


A. autovalores
B. proyección de muestras microarrays (proyecciones columnas)
C. proyección de genes (proyecciones filas).
En ambas representaciones de las muestras analizadas podemos observar que las líneas celulares
correspondientes a melanomas (ME) se distinguen del resto por el segundo eje.

A veces es útil tener una visión de separación por grupos de cada una de las Coordenadas.
La función heatplot representa en el eje de abcisas las coordenadas (Comp)) y en el eje de
4.2. APLICACIÓN PRÁCTICA DEL ACOI 89

ordenadas los pesos columnas (60 muestras).

heatplot(data.coa1$ord$co[,1:10],dend="none",
classvec=NCI60$classes, scale="none")

Color Key

−2 0 1 2
Value

BREAST_BT549
BREAST_HS578T
BREAST_MCF7
BREAST_MCF7ADRr
BREAST_MDAMB231
BREAST_MDAMB435
BREAST_MDAN
BREAST_T47D
CNS_SF268
CNS_SF295
CNS_SF539
CNS_SNB19
CNS_SNB75
CNS_U251
COLON_COLO205
COLON_HCC2998
COLON_HCT116
COLON_HCT15
COLON_HT29
COLON_KM12
COLON_SW620
LEUK_CCRFCEM
LEUK_HL60
LEUK_K562
LEUK_MOLT4
LEUK_RPMI8226
LEUK_SR
MELAN_LOXIMVI
MELAN_M14
MELAN_MALME3M
MELAN_SKMEL2
MELAN_SKMEL28
MELAN_SKMEL5
MELAN_UACC257
MELAN_UACC62
NSCLC_A549
NSCLC_EKVX
NSCLC_H226
NSCLC_H23
NSCLC_H322M
NSCLC_H460
NSCLC_H522
NSCLC_HOP62
NSCLC_HOP92
OVAR_IGROV1
OVAR_OVCAR3
OVAR_OVCAR4
OVAR_OVCAR5
OVAR_OVCAR8
OVAR_SKOV3
PROSTATE_DU145
PROSTATE_PC3
RENAL_7860
RENAL_A498
RENAL_ACHN
RENAL_CAKI1
RENAL_RXF393
RENAL_SN12C
RENAL_TK10
RENAL_UO31
Comp1

Comp2

Comp3

Comp4

Comp5

Comp6

Comp7

Comp8

Comp9

Comp10

Con este gráfico vemos qué peso tiene cada uno de los individuos en cada uno de los 10 primeros
ejes.
Aquellos que presentan el fenotipo tumoral leucemia (LE) son los que más peso tienen en el
primer eje, mientras que los que presentan fenotipo tumoral colon (CO) contribuyen más a la
determinación del segundo eje. Esto coincide con los resultados obtenidos tras el Análisis de
Correspondencia de ambos conjuntos de datos.
90 CAPÍTULO 4. ANÁLISIS DE COINERCIA

4.2.2. Análisis de Coinercia


Llevaremos a cabo un Análisis de Coinercia para comparar ambos conjuntos. Con la implemen-
tación del Análisis de Coinercia del paquete ade4 (usando la función cia), se supone que los
pesos de las filas de ambos conjuntos de datos son los mismos.

Aunque ambos conjuntos tengan genes diferentes, se espera que muestren patrones y tendencias
similares.

coin <- cia(NCI60$Ross, NCI60$Affy)


names(coin)

## [1] "call" "coinertia" "coa1" "coa2"

coin$coinertia

Coinertia analysis
call: coinertia(dudiX = coa1, dudiY = coa2, scannf =
cia.scan, nf = cia.nf)
class: coinertia dudi

$rank (rank) : 59
$nf (axis saved) : 2
$RV (RV coeff) : 0.7859656

eigenvalues: 2.266e-05 9.904e-06 4.342e-06 2.335e-06 ...

vector length mode content


1 $eig 59 numeric Eigenvalues
2 $lw 144 numeric Row weigths (for coa2 cols)
3 $cw 144 numeric Col weigths (for coa1 cols)

data.frame nrow ncol content


1 $tab 144 144 Crossed Table (CT): cols(coa2) x cols(coa1)
2 $li 144 2 CT row scores (cols of coa2)
3 $l1 144 2 Principal components (loadings for coa2 cols)
4 $co 144 2 CT col scores (cols of coa1)
5 $c1 144 2 Principal axes (loadings for coa1)
6 $lX 60 2 Row scores (rows of coa1 cols)
7 $mX 60 2 Normed row scores (rows of coa1)
8 $lY 60 2 Row scores (rows of coa2)
9 $mY 60 2 Normed row scores (rows of coa2)
4.2. APLICACIÓN PRÁCTICA DEL ACOI 91

10 $aX 2 2 Corr coa1 axes / coinertia axes


11 $aY 2 2 Corr coa2 axes / coinertia axes

CT rows = cols of coa2 (144) / CT cols = cols of coa1 (144)

Obtenemos un lista con los 59 autovalores, pesos filas (lw) y columnas (cw).
Se obtienen también las coordenadas para las variables de los conjuntos del grupo Affymetrix
(li) y Ross (co); así como las coordenadas para las muestras.
Las coordenadas para los individuos para el grupo Ross se recogen en $lX, mientras que las del
grupo Affymetrix se recogen en $lY .

El coeficiente RV es una medida de similitud global entre ambos conjuntos de datos. Mientras
más cercano a 1, mayor será la correlación entre ambos conjuntos.

coin$coinertia$RV

## [1] 0.7859656

Hemos obtenido un valor igual a 0,7859656, que indica una correlación bastante fuerte.

Para visualizar las líneas celulares con perfiles de expresión génica similares o distintas usamos
plotarrays.

Las muestras están representadas en las coordenadas dadas por lX y lY para cada matriz. Cada
línea celular está coloreada por su fenotipo (por ejemplo, colon de verde, pecho de rojo, melano-
ma de rosa, etc).

En el gráfico podemos ver puntos (Ross) y flechas (Affy) unidas por líneas. La distancia entre
estos puntos y flechas indica la similitud entre perfiles.
92 CAPÍTULO 4. ANÁLISIS DE COINERCIA

plotarrays(coin, classvec=NCI60$classes[,2], lab="", cpoint=3)

d=1

Figura 4.2: Comparación de los perfiles de expresión génica de spotted arrays y Affymetrix de
NCI60.
4.2. APLICACIÓN PRÁCTICA DEL ACOI 93

Con la función plot, se representan simultáneamente el gráfico anterior y las proyecciones


génicas para cada conjunto de datos.
En el primer gráfico puede verse como el primer eje separa las líneas celulares correspondientes
a leucemia (LE, color marrón) y colon (CO, color verde) del resto de fenotipos tumorales y que
las líneas celulares correspondientes a melanomas (ME, color rosa) se distinguen del resto por el
segundo eje.

Respecto a los genes representados en los gráficos B y C, podemos afirmar que aquellos localiza-
dos al final de los ejes son los que más intervienen en la definición de los mismos.

La interpretación conjunta de estos tres gráficos revela que genes y líneas celulares proyectadas
en la misma dirección desde el origen tienen una fuerte asociación.
94 CAPÍTULO 4. ANÁLISIS DE COINERCIA

plot(coin, classvec=NCI60$classes[,2])
d=1

OVAR_OVCAR4
OVAR_SKOV3
OVAR_OVCAR3
COLON_COLO205 NSCLC_A549
COLON_HCC2998
NSCLC_H322M
COLON_KM12 NSCLC_EKVX
COLON_HT29 RENAL_A498
RENAL_TK10
RENAL_UO31
PROSTATE_DU145
OVAR_OVCAR5 RENAL_RXF393
OVAR_IGROV1RENAL_ACHN
BREAST_T47D
RENAL_CAKI1 NSCLC_H226
OVAR_OVCAR8
RENAL_7860
BREAST_HS578T
BREAST_MDAMB231
BREAST_MCF7ADRr
COLON_SW620
COLON_HCT15
COLON_HCT116
BREAST_BT549
NSCLC_HOP62
CNS_SF295
BREAST_MCF7 NSCLC_H460
PROSTATE_PC3 CNS_SF539
NSCLC_HOP92
CNS_SF268
NSCLC_H522 MELAN_LOXIMVI
CNS_U251
NSCLC_H23 CNS_SNB19
LEUK_RPMI8226 RENAL_SN12C
LEUK_HL60 CNS_SNB75
LEUK_K562
LEUK_CCRFCEM
LEUK_MOLT4
LEUK_SR

MELAN_MALME3M
MELAN_SKMEL2
MELAN_M14
MELAN_UACC62
MELAN_UACC257
BREAST_MDAMB435
MELAN_SKMEL28
BREAST_MDAN

MELAN_SKMEL5

CIA of df1 NCI60$Ross and df2 NCI60$Affy


d = 0.01 d = 0.005

U78095_at

X510534 X21822
M94250_at Z24727_at
L33930_s_atHG2788.HT2896_at
M35878_at
X162479 X67951_at
M59807_at U62015_at
X429288 X136590
X47359
X486215 M93036_at
D17793_at
X470385
X211995 X302394 M31994_at
M80563_at
X68314_at J03764_at
M55998_s_at
X486700 X65614_at D49824_s_at
V00594_s_at
X510363
X510116 X487394 X375834
X108837 M92934_at
X510395
X429771 X488370
X487513 Z18951_at
X487974
X470386 M11749_at
J04469_at
L19686_rna1_at
X377620
X428733
X366310
X365558
X365476X471855
X486821 D49400_at
D14811_at
X92896_at
X04325_at L06419_at
M23294_at
X241935
X209758 X487878 HG1612.HT1612_at
M29277_at
M15841_at
X248589
J03040_at

variables df1 NCI60$Ross variables df2 NCI60$Affy

Figura 4.3: Análisis de Coinercia.


A. Representa las 60 muestras microarrays proyectadas en un espacio único. Los 60 círculos
representan el conjunto Ross y las 60 flechas representan el conjunto Affy.
B. proyecciones génicas para Ross.
C. proyecciones génicas para Affy.
Capítulo 5

Análisis de Correlación Canónica

El principal propósito del Análisis de Correlación Canónica (ACC) es la exploración de la


correlación entre conjuntos de variables.

Existen tres tipos de correlación:


Correlación simple, si se relacionan dos variables.
Correlación múltiple, si se relacionan una variable y un vector aleatorio.
Correlación canónica, si se relacionan dos vectores aleatorios.
Esta herramienta del análisis estadístico multivariante desarrollada en un principio por Hotelling
se basa en las proyecciones, ya que estructuras de datos multivariantes muy complejas son más
fáciles de comprender si estudiamos proyecciones de menor dimensión.

La proyección de una variable define un índice que tiene correlación máxima con el índice de
otra variable para cada muestra separadamente. El objetivo es maximizar la correlación entre las
proyecciones de menor dimensión de los dos conjuntos de datos.
Los vectores de correlación canónica se obtendrán a partir del análisis de covarianza conjunto de
las dos variables.

5.1. El método
Supongamos que tenemos dos conjuntos de variables aleatorias X ∈ Rp y Y ∈ Rq ,y por
tanto dos matrices X e Y de orden n × p y n × q, respectivamente. Las columnas de éstas
corresponderían a las variables y las filas a las unidades experimentales.

Asumiremos, sin pérdida de generalidad, que las variables de ambas matrices están estanda-
rizadas y que p < q. Supondremos además que ningún conjunto de variables explica el otro,
planteando por tanto el problema de forma simétrica.

95
96 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

En este análisis, el problema es encontrar un índice que describa la posible relación lineal entre
X e Y . Estos índices se conocen como variables canónicas y exactamente podemos contruir
m = mı́n{p, q} parejas de variables.

Estas variables son las combinaciones de las variables originales a0i X y b0i Y , con i = 1, . . . , m. Se
escogen los coeficientes ai y bi , conocidos como vectores canónicos, de forma que la correlación
lineal entre estos índices sea máxima.

El problema por tanto se plantea cómo la búsqueda de ai = (ai1 , · · · , aip )0 y bi = (bi1 , · · · , biq )0
tal que:

ηi = a0i X = ai1 X1 + · · · + aip Xp , ϕi = b0i Y = bi1 Y1 + · · · + biq Yq

y la correlación lineal entre ηi y ϕi sea máxima.

Veamos cómo expresar el coeficiente de correlación entre las proyecciones de X e Y .

Si suponemos que :
     
X µ ΣXX ΣXY
∼ ,
Y ν ΣYX ΣYY

Entonces sabemos que se cumplen las siguientes propiedades:

V ar(X) = ΣXX(p×p)

V ar(Y ) = ΣYY(q×q)

Cov(X, Y ) = E((X − µ)(Y − ν)0 ) = ΣXY = Σ0YX(p×q)

De esta forma, es fácil obtener la expresión de ρ(ηi , ϕi ) que buscamos:

Cov(ηi , ϕi ) Cov(a0i X, b0i Y ) a0i ΣXY bi


ρ(ηi , ϕi ) = p =p =
V ar(a0i X) V ar(b0i Y ) (a0i ΣXX ai )1/2 (b0i ΣYY bi )1/2
p p
V ar(ηi ) V ar(ϕi )

Para cualquier escalar c ∈ R+ se tiene que:

ρ(cηi , ϕi ) = ρ(ηi , cϕi ) = ρ(ηi , ϕi )

Sabido esto, podemos plantear apropiadamente el problema a resolver.


Comenzaremos calculando el primer par de vectores de correlación canónica a1 y b1 con detalle.
5.1. EL MÉTODO 97

5.1.1. Primer par de vectores de correlación canónica


El problema que tratamos resolver es el siguiente:

máx{a01 ΣXY b1 }
a1 ,b1

bajo las restricciones


 0
a1 ΣXX a1 = 1
b01 ΣYY b1 = 1
Las restricciones impuestas sirven para expresar que las varianzas de η1 y ϕ1 están normalizadas,
pues si aumentaramos el valor de a1 y b1 las varianzas correspondientes lo harían también.

Resolvemos este problema por el método de los multiplicadores de Lagrange.

Consideramos la siguiente función L y buscamos el máximo.


λ µ
L = (a01 ΣXY b1 ) − (a01 ΣXX a1 − 1) − (b01 ΣYY b1 − 1)
2 2
Para calcular el máximo igualamos las derivadas parciales de esta expresión a cero.
∂L
= ΣXY b1 − λΣXX a1 = 0
∂a1
∂L
= ΣYX a1 − µΣYY b1 = 0
∂b1
Despejándo obtenemos las siguientes relaciones:

ΣXY b1 = λΣXX a1
ΣYX a1 = µΣYY b1

Multiplicando por a01 en la primera ecuación y por b01 en la segunda se obtiene:

a01 ΣXY b1 = λa01 ΣXX a1




b01 ΣYX a1 = µb01 ΣYY b1

Puesto que ΣXY 0 es equivalente a ΣYX , podemos afirmar que λ = µ, y por tanto el sistema de
ecuaciones resultante es:

a01 ΣXY b1 = λa01 ΣXX a1 (5.2)


b01 ΣYX a1 = λb01 ΣYY b1 (5.3)
Si despejamos b1 de 5.3, se obtiene la relación b1 = λ−1 ΣYY −1 ΣYX a1 .
98 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

Sustituyendo en la ecuación 5.2 se satisface la siguiente igualdad:

ΣXY (λ−1 ΣYY −1 ΣYX a1 ) = λΣXX a1

que equivale a (ΣXX −1 ΣXY ΣYY −1 ΣYX )a1 = λ2 a1

Es sencillo observar que a1 es el autovector correspondiente al valor propio λ2 de la matriz


(ΣXX −1 ΣXY ΣYY −1 ΣYX ).

Operando de la misma forma con la ecuación restante obtenemos la relación correspondiente a


b1 . A partir de esto podemos afirmar el siguiente resultado.

Teorema 5.1.1 Los primeros vectores canónicos verifican las siguientes relaciones

(ΣXX −1 ΣXY ΣYY −1 ΣYX )a1 = λ2 a1

(ΣYY −1 ΣYX ΣXX −1 ΣXY )b1 = λ2 b1 ,


siendo λ uno de los multiplicadores de Lagrange asociado al problema anteriormente planteado.

Teorema 5.1.2 Los vectores canónicos, cumpliendo a01 ΣXX a1 = 1 y b01 ΣYY b1 = 1, se relacio-
nan por
a1 = λ−1 ΣXX −1 ΣXY b1 ,
b1 = λ−1 ΣYY −1 ΣYX a1 .
Definiéndose la primera correlación canónica como ρ1 = λ1 1/2 , donde λ1 es el mayor autovalor
de ΣXX −1 ΣXY ΣYY −1 ΣYX .
Demostración: La demostración puede verse en [1]

El procedimiento para obtener el resto de vectores que definen las variables de correlación
canónica es similar, suponiendo además que cada variable de correlación canónica es incorrelada
con todas las anteriores.
Esto permite obtener cada par de variables ordenadas por importancia.

Aparte de las relaciones entre los vectores de correlación canónica obtenidas mediante el
método de los multiplicadores de Lagrange, podemos abordar el problema planteado usando
Descomposición en Valores Singulares; tal y como veremos a continuación.
5.1. EL MÉTODO 99

5.1.2. Correlación Canónica y Descomposición Singular


Para resolver el problema planteado al comienzo de esta sección, partiremos de la siguiente
matriz:
−1/2 −1/2
K = ΣXX ΣXY ΣYY
siendo ΣXX la varianza del vector X, ΣYY la del vector Y y ΣXY = Σ0YX la covarianza entre
ambos conjuntos de variables X e Y .

Llevando a cabo la descomposición en valores singulares de K obtenemos:

K = ΓD∆0

siendo Γ = (γ 1 , γ 2 , . . . , γ m ) autovectores de KK0 , ∆ = (δ 1 , δ 2 , . . . , δ m ) autovectores de K0 K


y D matriz diagonal con la raíz de los m autovalores no nulos de KK0 .

Teorema 5.1.3 Los vectores de correlación canónica a1 y b1 son:


−1/2
a1 = ΣXX γ 1

−1/2
b1 = ΣYY δ 1
siendo γ 1 y δ 1 los autovectores recogidos en Γ y ∆ que están asociados al mayor autovalor de
KK0 .

Demostración:

Sabemos por las propiedades de la descomposición en valores singulares que Γ = (γ 1 , γ 2 , . . . , γ m )


es el conjunto de autovectores de KK0 cumpliendo KK0 = ΓD2 Γ0 .

Desarrollamos el producto KK0 como sigue


−1/2 −1/2 −1/2 −1/2 −1/2 −1/2
KK0 = ΣXX ΣXY ΣYY ΣYY ΣYX ΣXX = ΣXX ΣXY Σ−1
YY ΣYX ΣXX

Si γ 1 es el autovector de KK0 asociado al mayor autovalor se cumple:

−1/2 −1/2
ΣXX ΣXY Σ−1
YY ΣYX ΣXX γ 1 = λ1 γ 1

−1/2
Multiplicamos a izquierda por ΣXX en ambos lados de la igualdad:
−1/2 −1/2
Σ−1 −1
XX ΣXY ΣYY ΣYX (ΣXX γ 1 ) = λ1 ΣXX γ 1
100 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

Comparando esto con la relación para a1 descrita en el Teorema 5.1.1 queda probada la primera
igualdad.

Siguiendo el mismo razonamiento con K0 K probaremos la segunda igualdad.




Es posible probar que las correlaciones canónicas son invariantes por transformaciones lineales
y, por tanto, pueden calcularse partiendo de las matrices de correlaciones.

En la siguiente sección veremos con detalle la expresión explícita de los m pares de vectores y
variables de correlación canónica.

5.2. Variables de correlación canónica


Definimos estas variables a partir de los vectores de correlación canónica.

Definición 5.2.1 Los m vectores de correlación canónica ai y bi , se definen como:


−1/2
ai = ΣXX γ i

−1/2
bi = ΣYY δ i ,
con i = 1, . . . , m (m = min{p, q})

γ i es el i-ésimo autovector recogido en Γ y δ i el i-ésimo autovector en ∆ asociado al autovalor


λi correspondiente.

Una vez que conocemos el valor exacto de nuestros índices, podemos definir las variables de
correlación canónica

Definición 5.2.2 Las variables de correlación canónica se definen como

ηi = a0i X

ϕi = b0i Y
1/2
siendo ρi = λi el correspondiente coeficiente de correlación canónica, con i = 1, . . . , m.

El siguiente resultado recoge propiedades fundamentales de este método.


5.2. VARIABLES DE CORRELACIÓN CANÓNICA 101

Teorema 5.2.1 Si suponemos ρ1 > ρ2 > · · · > ρm , se cumple:


1. Las variables canónicas η1 , · · · , ηm y ϕ1 , · · · , ϕm están incorreladas.

0 1, i = j
Cov(ηi , ηj ) = ai ΣXX aj =
0, i 6= j


1, i = j
Cov(ϕi , ϕj ) = b0i ΣYY bj =
0, i 6= j
2. La primera correlación canónica ρ1 es la correlación máxima entre una proyección de X
y una de Y .
3. ρ2 es la correlación máxima entre proyecciones de X incorreladas con η1 y proyecciones
de Y incorreladas con ϕ1 .
4. Cov(ηi , ϕj ) = 0 si i 6= j
Demostración:

1. Inmediata por la ortogonalidad que caracteriza a los autovectores.

4. Sea Cov(ηi , ϕj ) = E(a0i XY 0 bj ) = a0i ΣXY bj .

Por el Teorema 5.1.2 obtuvimos una expresión del vector bj . Sustituyendo en la expresión de la
covarianza se muestra:
Cov(ηi , ϕj ) = a0i (λj ΣXX aj) = 0


Teorema 5.2.2 Sean ηi y ϕi las i-ésimas variables de correlación canónica, con i = 1, · · · , m.


Definimos η = (η1 , · · · , ηm ) y ϕ = (ϕ1 , · · · , ϕm ) los vectores m-dimensionales. Entonces la
matriz de varianzas y covarianzas es:
   
η Ik Λ
V ar =
ϕ Λ Im
1/2 1/2
con Λ = diag(λ1 , · · · , λm )
1/2
Este teorema prueba que los coeficientes de correlación canónica ρi = λi son las covarianzas
entre las√variables canónicas ηi y ϕi y los índices η1 = a01 X y ϕ1 = b01 Y tienen covarianza
máxima λ1 = ρ1 .
102 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

5.3. Contrastes de significación


Una vez obtenidas las m relaciones canónicas, nos puede interesar decidir cuáles son significati-
vas.
Suponiendo normalidad multivariante, pueden realizarse contrastes como el test secuencial de
Barlett-Lawley o contrastes de independencia de los conjuntos de variables originales.

•Test de Bartlett-Lawley:

Suponemos X ∼ Np (0, ΣXX ) e Y ∼ Nq (0, ΣYY )

En este test se plantea que las primeras k correlaciones canónicas obtenidas a partir de la matriz
de varianzas y covarianzas de X e Y son no nulas, a diferencia del resto.

Para k = 0, 1, · · · , m; y ρ0 = 1, se plantea

H0k : ρk > ρk+1 = · · · = ρm = 0

Si se acepta la hipótesis nula, se acepta por tanto que las k primeras variables de correlación
canónica expresan satisfactoriamente la dependencia entre variables; mientras que si se acepta la
hipótesis alternativa quiere decir que no podemos prescindir de ninguna variable de correlación
canónica.

Sea Lk es el estadístico del contraste


" k
# " m #
1 X Y
Lk = − n − 1 − k − (p + q + 1) + ρ−2
i log (1 − ρ2i )
2 i=1 i=k+1

se aceptará la hipótesis nula cuando este no sea significativo.

Contrastes de independencia
Este contraste se lleva a cabo suponiendo normalidad multivariante de X e Y y planteando la
hipótesis nula de que X e Y son independientes.

Puesto que se supone normalidad, la hipótesis nula será equivalente a plantear que las variables
son incorreladas (ΣXY = 0). Por tanto, el contraste es el siguiente:

H0 : ΣXY = 0
H1 : ΣXY 6= 0
5.3. CONTRASTES DE SIGNIFICACIÓN 103

Para resolverlo tenemos dos opciones:

• Test de razón de verosimilitudes

Sea x = (x1 , · · · , xn ) una muestra aleatoria de X con función de densidad definida como f (x, θ)
y θ ∈ Θ un parámetro desconocido. Si existe una subregión paramétrica Θ0 de Θ, se puede
plantear el test:

H0 : θ ∈ Θ0
H1 : θ ∈ Θ − Θ0
La razón de verosimilitud es el estadístico:
L(x, θˆ0 )
λR =
L(x, θ̂)
que toma valores entre 0 y 1, siendo θ̂ el estimador de máxima verosimilitud de θ ∈ Θ , θˆ0 el
correspondiente a θ ∈ Θ0 y L(x, θ) la función de verosimilitud definida como:
n
Y
L(x, θ) = f (xi , θ)
i=1

Aplicando este test a nuestro caso, el estadístico que obtendríamos sería:


|Σ|
b |R|
b
λR = =

b XX ||Σ
b YY | |R
b XX ||R
b YY |
que sigue una distribución Lambda de Wilks Λ(p, n − 1 − q, q) y cuando toma valores pequeños
y significativos supone el rechazo de la hipótesis nula H0 .

Puede establecerse esta igualdad porque el estadístico es invariante por transformaciones lineales.

La distribución Lambda de Wilks se define a partir de dos conjuntos de variables independientes


con distribución Wishart: A ∼ Wp (Σ, m), B ∼ Wp (Σ, n)

|A| 1
λR = = ∼ Λ(p, m, n)
|A + B| |(I + A−1 B)|

• Principio de Unión-Intersección
Partimos de que disponemos de los vectores η y ϕ.

El Principio de Unión-Intersección plantea tanto contrastes multivariantes como test univariantes.


Por ejemplo, si nuestra hipótesis nula es del tipo H0 : µ = µ0 , haciendo uso del vector a podemos
transformar el contraste como sigue:
104 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

Sustituimos X por X a = Xa, de manera que la hipótesis nula sería:

H0 = ∩a H0 (a), con H0 (a) : µ(a) = µ0 (a)

Si aplicamos este Principio a nuestro caso, la hipótesis nula sería H0 : ρ(η, ϕ) = 0

Esto nos lleva a estudiar la significación de la primera correlación canónica a través del test:

H0 : ρ1 = maxρ(η, ϕ) = 0
H1 : ρ1 > 0

Sólo se acepta la hipótesis nula en caso de que la correlación no sea significativa.

Observación 5.3.1 Como hemos mencionado anteriormente, ésta técnica de la que posterior-
mente derivarían otras técnicas de Análisis Multivariante fue introducida por Hotelling. El
objetivo era encontrar el trasfondo de las relaciones entre cuerpo y mente a través de test
mentales y medidas biométricas.

Algunas de las aplicaciones de la técnica se han llevado a cabo en psicología, ecología o en


datos génicos.

Los únicos resultados que se conocen sobre la distribución de las correlaciones canónicas son
asintóticos (Muirhead, 1982), debido a la gran complejidad que supone.

5.4. Aplicación práctica del ACC


El Análisis de Correlación Canónica es un análisis exploratorio que permite ilustrar las correla-
ciones entre dos conjuntos de datos en que se miden las mismas unidades experimentales.

Dado que en este caso, el número de variables supera considerablemente el número de muestras
analizadas; tendremos que tener esto en cuenta a la hora de realizar el análisis.

Comprobaremos que con los dos primeros pares de variables de correlación canónica podemos
representar con gran exactitud la relación entre los dos conjuntos de variables.

Datos
En el siguiente ejercicio práctico trabajaremos con la base de datos nutrimouse, disponible
en el paquete CCA.
En esta base de datos se recoge un estudio de la nutrición de ratones. Se estudiaron un total de
40 ratones, clasificados en función de dos factores:
1. Genotipo: El estudio se ha llevado a cabo en ratones de tipo salvaje y P P ARα deficientes.
5.4. APLICACIÓN PRÁCTICA DEL ACC 105

2. Dieta: Se usaron aceites para la preparación experimental de dietas. Acites de colza y maíz
(50/50) fueron usados para una dieta de referencia (REF), aceite de coco hidrogenado para
una dieta de ácidos grasos saturados (COC), aceite de girasol para una dieta rica en ácidos
grasos w6 (SUN), aceite de semilla de lino para una dieta rica en ácidos grasos w3 (LIN) y
aceites de maíz/colza/pescado enriquecido (42.5/42.5/15) para una dieta de pescado (FISH).
Mediremos 141 variables agrupadas en los siguientes conjuntos:

1. Los niveles de expresión de 120 genes medidos en células del hígado. Los genes elegidos
se han seleccionado de un total de 30000, pues se consideran relevantes en el contexto de
nuestro estudio.
2. Concentraciones de 21 ácidos grasos hepáticos medidos por cromatología de gases.

Paquetes
Para este ejemplo trabajaremos con el paquete CCA (disponible en “http://CRAN.R-project.org/”),
que contiene las herramientas necesarias para llevar a cabo un Análisis de Correlación Canónica
en conjuntos de datos con más variables que observaciones. Además, permite trabajar con valores
perdidos.

5.4.1. ACC
En primer lugar, cargamos el paquete CCA previamente instalado.

require(CCA)

Una vez hecho esto, cargamos los datos contenidos en “nutrimouse”. Con la función View
podemos ver la matriz de datos en la consola de Rstudio.

data("nutrimouse")

Con class podemos ver la clase del objeto con el que trabajamos.

class(nutrimouse)

## [1] "list"

Es una lista de 4 elementos, que contiene:


1. “nutrimouse$gene”, un data.frame de dimensión 40 × 120 que contiene los niveles de ex-
presión de 120 genes en 40 muestras. Esta será una de las matrices con las que trabajaremos
más adelante.
2. “nutrimouse$lipid”, un data.frame de dimensión 40 × 21. Esta será la otra matriz utilizada
en el análisis.
106 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

3. “nutrimouse$diet", un factor de longitud 40 indicando la dieta que sigue cada uno de los
individuos.
4. “nutrimouse$genoType”, un factor de longitud 40 con distintos niveles que indican el tipo
de ratón analizado en función de su genotipo.
Dado que los datos están almacenados en una lista, para evitar problemas que puedan surgir
más adelante los convertiremos en una matriz, diferenciando las dos clases existentes entre las
variables especificadas con anterioridad.

Los datos de los ratones a los que se han medido variables relacionadas con expresión génica se
guardan en la matriz X y el resto de datos en la matriz Y.
Con la función dim comprobamos que las dimensiones de las matrices se corresponden con lo
explicado.

X <- as.matrix(nutrimouse$gene)
Y <- as.matrix(nutrimouse$lipid)
dim(X)

## [1] 40 120

dim(Y)

## [1] 40 21

Un paso preliminar para aplicar luego la técnica es conocer las correlaciones entre variables de
ambas matrices. Haremos esto con la función matcor.

correl <- matcor(X, Y)

Podemos ver gráficamente esta matriz con la función img.matcor, donde “type=2” indica que
tenemos en cuenta las 3 matrices generadas anteriormente; de dimensiones p × p, q × q y p × q
respectivamente.
5.4. APLICACIÓN PRÁCTICA DEL ACC 107

img.matcor(correl, type = 2)

X correlation Y correlation

Cross−correlation

−1.0 −0.5 0.0 0.5 1.0

Figura 5.1: Matrices de correlación para las variables de X (lado superior izquierdo), para las
variables de Y (lado superior derecho), para las variables de ambos conjuntos (parte inferior).
Los valores se representan por colores que van desde azul (correlación negativa) hasta rojo
(correlación positiva).
108 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

Para ilustrar la técnica se han seleccionado 10 variables de la matriz X con sample y el argu-
mento “size=10”. Hemos fijado una semilla para que siempre haga el muestreo seleccionando
los mismos individuos.
En cambio, sí seleccionamos todas las variables de Y, ya que tan sólo son 21.

Aparecerán los valores que toman las variables de X seleccionadas al azar, así como las de Y
para los 40 individuos.

Dado que se supone que las variables están estandarizadas cuando se aplica esta técnica, con la
función scale estandarizamos los datos de ambas tablas de datos.

La función cc() se usa para llevar a cabo el análisis en este caso, donde hay que tener en cuenta
que existen muchas más variables que individuos analizados.

set.seed(24)
Xr<- X[,sample(1:120, size=10)]
Xscale<- scale(Xr)
Yscale<- scale(Y)
res.cc <- cc(Xscale, Yscale)

Con la función names podemos ver todos los nombres de los elementos que se ha generado al
aplicar ACC.

names(res.cc)

## [1] "cor" "names" "xcoef" "ycoef" "scores"

Si vemos los resultados de “res.cc” paso a paso, podemos observar que en primer lugar aparecen
los valores de las correlaciones de los 10 índices resultantes (cor(a0i X, b0i Y )), i = 1, . . . , 10. La
primera combinación de variables originales es la que tiene mayor correlación entre índices.

res.cc$cor

## [1] 0.9888384 0.9795369 0.9409709 0.9218689 0.8894283


## [6] 0.8042567 0.7679251 0.6598241 0.6176711 0.4321375

Cabe destacar que se han calculado 10 combinaciones posibles de las variables de cada matriz
porque, como sabemos por teoría, se lleva a cabo el mínimo entre p y q.

Con barplot, respresentamos las correlaciones canónicas en función de las 10 variables de


correlación canónica obtenidas en el análisis.
5.4. APLICACIÓN PRÁCTICA DEL ACC 109

barplot(res.cc$cor, xlab = "Dimension", ylab =


"Canonical correlations", names.arg = 1:10, ylim = c(0,1))
1.0
0.8
Canonical correlations

0.6
0.4
0.2
0.0

1 2 3 4 5 6 7 8 9 10

Dimension

Figura 5.2: Correlación entre las variables canónicas, teniendo en cuenta los 10 índices calcula-
dos.
110 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

También podemos ver los nombres de las variables que se han seleccionado de ambos conjuntos
para llevar a cabo el análisis.

res.cc$names

## $Xnames
## [1] "CYP8b1" "CYP27a1" "RARa" "Lpin3" "PMDCI"
## [6] "apoA.I" "CYP3A11" "RXRb2" "SHP1" "CYP2b10"
##
## $Ynames
## [1] "C14.0" "C16.0" "C18.0" "C16.1n.9" "C16.1n.7"
## [6] "C18.1n.9" "C18.1n.7" "C20.1n.9" "C20.3n.9" "C18.2n.6"
## [11] "C18.3n.6" "C20.2n.6" "C20.3n.6" "C20.4n.6" "C22.4n.6"
## [16] "C22.5n.6" "C18.3n.3" "C20.3n.3" "C20.5n.3" "C22.5n.3"
## [21] "C22.6n.3"
##
## $ind.names
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11"
## [12] "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22"
## [23] "23" "24" "25" "26" "27" "28" "29" "30" "31" "32" "33"
## [34] "34" "35" "36" "37" "38" "39" "40"

Tras esto, aparecen los coeficientes para cada conjunto de variables. Esto en teoría se correspon-
de con los valores de ai y bi para cada uno de los 10 índices que hemos obtenido, respectivamente.

Con la función head hemos reducido la salida de R para mostrar tan sólo los valores para las
primeras variables de cada conjunto.

head(res.cc$xcoef)
head(res.cc$ycoef)

[,1] [,2] [,3] [,4] [,5]


CYP8b1 -0.04317032 -0.13249783 -0.1225827 0.03851034 0.1431669
CYP27a1 0.13837072 -0.06704582 0.6847454 -0.90146905 -0.5080292
RARa -0.06755834 0.36427522 0.6446277 0.74865435 -0.5921778

[,1] [,2] [,3] [,4] [,5]


C14.0 -3.709529 -12.235375 -19.771392 -20.354826 -13.501884
C16.0 -20.181158 -52.688694 -86.390523 -97.289840 -71.384173
C18.0 -14.229198 -37.673860 -63.710750 -70.225616 -53.371176
5.4. APLICACIÓN PRÁCTICA DEL ACC 111

A continuación, se muestran algunos de los valores que toman los individuos cuando introducimos
los coeficientes obtenidos.

head(res.cc$scores$xscores)

[,1] [,2] [,3] [,4] [,5]


1 0.6203762 -1.10043552 -0.8632726 0.28122447 -1.1996885
2 1.3464731 -0.29545361 0.2785762 0.37691778 0.6708977
3 1.0441307 -0.33255241 -0.4895570 0.06214316 1.8445435

head(res.cc$scores$yscores)

[,1] [,2] [,3] [,4] [,5]


1 0.5162834 -1.3123080 -1.1923051 0.4132465 -0.8649199
2 1.3433274 -0.1746114 0.2325537 0.6365109 1.6488218
3 0.9831786 -0.3782187 -0.4989220 0.1365050 1.0234334

También se muestran las correlaciones entre las variables originales y las de correlación canó-
nica. Se muestran algunos valores de cor(X, a0i X); aunque se pueden obtener también los de
cor(Y , a0i X), cor(X, b0i Y ), cor(Y , b0i Y ).

head(res.cc$scores$corr.X.xscores)

## [,1] [,2] [,3] [,4]


## CYP8b1 0.09787891 -0.311267976 -0.02742548 -0.09783216
## CYP27a1 0.52252379 -0.227431346 0.47454092 -0.06563921
## RARa 0.07793188 0.237850465 0.42938647 0.45099119
## Lpin3 -0.03257523 0.262348515 0.10301516 0.40953864
## PMDCI 0.75313911 -0.601752574 0.03463925 0.24570849
## apoA.I 0.69139446 -0.000527109 0.03644501 -0.06229574
## [,5] [,6] [,7] [,8]
## CYP8b1 0.38206052 0.730094529 0.01819970 -0.42071311
## CYP27a1 -0.07135934 0.376644816 0.25243205 0.04629565
## RARa 0.14691865 0.219137273 0.37205274 -0.24975868
## Lpin3 0.19061760 0.114730180 0.22379031 -0.26325679
## PMDCI -0.06771857 0.001600979 0.01355798 0.04560853
## apoA.I 0.35553712 0.137263174 0.58591644 -0.14603883
112 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

## [,9] [,10]
## CYP8b1 -0.155085083 0.053087487
## CYP27a1 -0.005629969 -0.482592885
## RARa 0.260445559 -0.459725158
## Lpin3 -0.092898047 -0.757801903
## PMDCI 0.046948637 -0.006803365
## apoA.I -0.034246900 -0.075619656

res.cc$scores$corr.Y.xscores
res.cc$scores$corr.X.yscores
res.cc$scores$corr.Y.yscores

Con la función plt.cc se representan variables e individuos.

Variables e individuos pueden representarse por el par (ηs , ηt ) o bien por (ϕs , ϕt ); con 1 ≤ s <
t ≤ m.

Dado que las correlaciones entre las dos primeras variables canónicas son cercanas a 1, las
representaciones gráficas en los ejes definidos por ambas (η1 , η2 ) y (ϕ1 , ϕ2 ) son similares y
hemos representado tanto variables como individuos haciendo uso de (η1 , η2 ).

En el gráfico A se representan las correlaciones entre variables, diferenciando dos círculos (uno
de radio 0.5 y otro de radio unidad) que revelen los patrones más destacados en el anillo entre
las dos circunferencias.

Los puntos representados en rojo indican las correlaciones entre las 10 variables de X y las dos
primeras variables de correlación canónica a01 X y a02 X, por tanto las coordenadas de cada punto
son (cor(a01 X, X), cor(a02 X, X)) = (cor(η 1 , X), cor(η 2 , X)).

Sin embargo, los triángulos indican la correlación de las 21 variables de la matriz Y con las
variables de correlación canónica a01 Y y a02 Y , siendo por tanto sus coordenadas:
(cor(a01 X, Y ), cor(a02 X, Y )) = (cor(η 1 , Y ), cor(η 2 , Y ))

Se asume que las variables Xj e Yk tienen varianza unidad, sus proyecciones en el plano están
dentro de un círculo de radio 1 centrado en el origen, llamado el círculo de correlación.
Variables con una fuerte relación se proyectan en la misma dirección desde el origen. Mientras
más distancia hay al origen más fuerte es la relación.

En el gráfico B se representan los individuos caracterizados por el genotipo “ppar” o “wt” y la


dieta que siguen. Son las puntuaciones de los individuos para la combinación de la variable X,
representadas en las dos primeras dimensiones (a01 X, s02 X)=(η1 ,η2 ).
5.4. APLICACIÓN PRÁCTICA DEL ACC 113

plt.cc(res.cc, var.label = TRUE, ind.names =


paste(nutrimouse$genotype,nutrimouse$diet, sep = "-"))
1.0

ppar−sun
2

ppar−lin
ppar−sunppar−coc
C16.1n.9 ppar−refppar−coc
0.5

C18.1n.9 ppar−coc
C18.1n.7
C14.0
1

C20.2n.6C16.1n.7
C22.4n.6 ppar−coc
C20.1n.9
C18.2n.6 Lpin3 ppar−ref
ppar−fish
RARa
RXRb2 wt−coc
Dimension 2

Dimension 2

C22.5n.6 ppar−ref ppar−sun


CYP2b10
C18.3n.6 C20.4n.6 ppar−ref wt−ref
C20.3n.9 apoA.I
0.0

ppar−lin
SHP1 wt−ref
wt−coc
0

C18.3n.3 wt−sunwt−refwt−ref
C20.3n.3 CYP27a1
C20.3n.6 wt−sun
wt−cocwt−sun
wt−sun
CYP8b1 ppar−fish wt−coc
ppar−lin
ppar−lin
−0.5

C16.0
C18.0
ppar−fish
−1

C20.5n.3
C22.5n.3 PMDCI wt−lin
wt−lin
C22.6n.3 wt−lin
wt−lin
wt−fish
wt−fish
ppar−fish
CYP3A11
−1.0

wt−fish

−1.0 −0.5 0.0 0.5 1.0 −2 −1 0 1

Dimension 1 Dimension 1
114 CAPÍTULO 5. ANÁLISIS DE CORRELACIÓN CANÓNICA

Es sencillo observar que la variable de correlación canónica η1 (dimensión 1) separa los indi-
viduos por genotipo ya que mientras que en el lado izquierdo están representados los ratones
P P ARα deficientes, en el derecho aparecen los ratones de tipo salvaje.
Por otra parte, la variable de correlación canónica η2 (dimensión 2) separa los tipos de dietas en
función de su composición. La separación de dietas en los ratones de tipo P P ARα es menos
precisa.

Las relaciones entre los dos gráficos pueden revelar asociaciones entre variables y unidades.
Capítulo 6

Paquetes de R

En este capítulo se hace un listado de todos los paquetes de R utilizados en las aplicaciones
prácticas de la memoria, explicando sus características principales. También se han detallado las
funciones de las que hemos hecho uso y los argumentos que se han considerado relevantes.

Paquete “base”
El paquete base es el que ejecuta las funciones básicas de R.

Funciones
1. as.data.frame: Comprueba si un objeto es data.frame y si no lo es, lo transforma.
2. as.factor: Comprueba si un objeto es un factor y si no lo es, lo transforma.
3. as.matrix: Comprueba si un objeto es una matriz y si no lo es, lo transforma.
4. cbind.data.frame: Toma la secuencia de los argumentos de un data.frame y combina
sus columnas.
5. class: Muestra la clase a la que pertenece un objeto.
6. colnames: Conjunto de nombres de las columnas de un objeto de tipo matriz.
7. cumsum: Calcula la suma acumulada de los valores de su argumento.
8. data: Especifica un conjunto de datos.
9. dim: Dimensión de un objeto.
10. file.path(): Construye el camino a un fichero desde componentes en otra plataforma.
11. getwd(): Devuelve el directorio de trabajo que se está usando.
12. length: Longitud de un objeto.
13. library & require: Se cargan y adjuntan los paquetes especificados.
14. names: Nombres de un objeto.

115
116 CAPÍTULO 6. PAQUETES DE R

15. options: Permite al usuario fijar funciones globales que afectan a la forma en que se
computan los resultados.
digits: controla la cantidad de dígitos cuando se trabaja con valores numéricos.
width: controla el número de columnas máximo por línea para mostrar vectores y
matrices.
16. paste: Concatena vectores tras convertirlos en caracteres.
17. read.table: Lee un fichero en formato tabla y crea un “data.frame” a partir de él, cuyos
casos se corresponden con las filas y las variables con los campos en el fichero.
dec: El caracter a usar con los puntos decimales.
header: Valor lógico indicando si los nombres de las variables se asignan a la primera
fila.
row.names: Puede ser un vector con los nombres de las filas o bien un número con la
columna de la tabla que contiene los nombres de las filas.
sep: Separador de caracteres.
18. rownames: Conjunto de nombres de las filas de un objeto de tipo matriz.
19. sample: Toma una muestra de un tamaño especificado de los elementos de un objeto, con
o sin reemplazamiento.
20. sapply: Devuelve un vector o matriz de la misma dimensión que la dada, resultado de
aplicar una función especificada.
21. scale: Función genérica que centra y escala las columnas de una matriz numérica.
22. str: Alternativa a summary que muestra la estructura de un objeto.
23. sum: Calcula la suma de los valores de su argumento.
24. summary: Resúmenes de resultados de varios modelos.
nbelements: Número de elementos a mostrar.
25. table: Usa factores de clasificación cruzada para construir una tabla de contingencia de
los elementos de cada combinación de los factores.
26. View: Permite visualizar cualquier objeto de R con forma de matriz.

Paquete “CCA”
CCA proporciona un conjunto de funciones que extienden la función cancor() del paquete
stats con nuevas salidas numéricas y gráficas.
La función cancor() calcula las correlaciones canónicas entre dos matrices de datos. El
paquete CCA incluye una extensión que permite trabajar con más variables que observaciones,
así como con valores perdidos.
117

Funciones
1. barplot: Crea una representación gráfica con barras verticales u horizontales.
names.arg: Vector de nombres para ser representado bajo cada barra o grupo de
barras.
xlab: nombre eje X.
xlim: límites eje X.
ylab: nombre eje Y.
ylim: límites eje Y.
2. cc: Ejecuta un Análisis de Correlación Canónica. Completa la función cancor().
3. img.matcor: Muestra imágenes de las matrices de correlación entre 2 conjuntos de
datos.
type: Caracter determinando la clase de plot. “type=1” representa una matriz (p + q) ×
(p + q) y “type=2” 3 matrices (p × q), (q × q) y (p × p).
4. matcor: Matrices de correlación entre dos conjuntos de datos.
5. plt.cc: Utiliza la función plt.var(), plt.ind() o ambas para representar indivi-
duos, variables o ambos en un gráfico con ejes las variables canónicas.
ind.names: Vector conteniendo los nombres de los individuos.
var.label: Valor lógico indicando si el eje debe incluirse en la representación de
variables.

Paquete “FactoMineR”
FactoMineR permite realizar análisis exploratorio de datos multivariantes y minería de datos.

Funciones
1. PCA: LLeva a cabo un Análisis de Componentes Principales de los datos.
graph: Toma valores lógicos. Indica si incluiremos gráficos o no.
qualisup: Variables cualitativas suplementarias.
2. plot: Representa gráficos.
choix: Elegimos qué representar (variables o individuos).
habillage: Damos color a los elementos representados en función del valor que tome
este parámetro.
lim.cos2.var: Indica la calidad mínima de la proyección de los objetos que dibujare-
mos.
3. plotellipses: Dibuja las elipses de confianza para cada clase de una variable categóri-
ca. El nivel por defecto es 0,95.
4. substr: Extrae o sustituye subcadenas de un vector de caracteres.
118 CAPÍTULO 6. PAQUETES DE R

Paquete “graphics”
Graphics reune un conjunto de funciones para gráficos básicos.

Funciones
1. hist: Genera un histograma con los datos proporcionados.
plot: Toma valores lógicos. Con FALSE se devuelve una lista con los puntos para
ambos ejes (“breaks” y “counts”) en lugar de representar el histograma. TRUE es el
valor por defecto.
2. par: Ajusta parámetros gráficos
mf row: Indica cómo se dispondrán los plots por filas.
3. pairs: Produce una matriz de scatterplots.
cex: Valor numérico con la cantidad por la que se rigen el texto y los símbolos
representados. Comienza con el 1.
diag.panel: Función para aplicar en las diagonales.
panel: Función que se usa para representar el contenido de cada panel.
pch: Símbolo para usar en el plot.
xlim: Límites eje X.
ylim: Límites eje Y.
4. plot: Representa gráficos.
5. text: Añade texto a un plot.

Paquete “HSAUR3”
Este paquete proporciona funciones, conjuntos de datos, análisis y ejemplos de la tercera edición
del libro A Handbook of Statistical Analysis using R (Torsten Hothorn and Brian S. Everiitt,
Chapman & Hall, CRC, 2014).

Paquete“knitr”
knitr se usa como herramienta alternativa a Sweave. Mientras que Sweave es un sistema que
proporciona un marco de trabajo para generar documentos que combinan simultáneamente
texto y código R, con el paquete knitr se consigue lo mismo mediante un diseño más flexible y
características nuevas; entre ellas la posibilidad de tener un control más preciso de los gráficos.
119

Funciones
1. opts_chunk: Opciones para los diferentes chunk de R, siendo un chunk cada una de las
partes de nuestro documento en las que introducimos código de R.
2. opts_chunk$set(): Suele utilizarse en el primer chunk del documento, donde se fijan
las opciones globales para los siguientes chunk.
concordance: Toma valores lógicos. Permite escribir ficheros de concordancia, per-
mitiendo así situar los números de las líneas del código con los de las líneas de los
resultados, siendo esto de gran ayuda cuando tenemos un error al generar el PDF.
También podemos usar la sintaxis de este paquete para manejar cómo incluir el código R en
nuestro documento. Las opciones disponibles son:
cache: FALSE por defecto. Guarda los resultados en cache.
comment: Caracter para comentarios.
echo: TRUE por defecto. Muestra el código y los resultados.
error: FALSE por defecto. Muestra errores.
eval: TRUE por defecto. Evalúa el código e incluye los resultados.
f ig.height: 7 por defecto. Alto en pulgadas para las figuras.
f ig.width: 7 por defecto. Ancho en pulgadas para las figuras.
message: TRUE por defecto. Muestra mensajes.
results: “markup” por defecto. Controla la forma en que se muestran los resultados al
compilar el documento a PDF. Puede ser “markup”, “hide”, “asis”, y “hold”.
tidy: FALSE por defecto. Muestra el código de forma organizada.
warning: TRUE por defecto. Muestra advertencias.

Paquetes “ade4” y “made4”


ade4 es un paquete desarrollado por el laboratorio de Biometría y Biología Evolucionaria
(UMR 5558) de la Universidad de Lyon. Contiene funciones para analizar datos ecológicos y
ambientales en el marco de métodos exploratorios euclídeos.
made4 necesita ade4. Además, facilita el análisis multivariante de datos de expresión génica
obtenidos con microarrays.

Funciones
1. cia: Análisis de Coinercia con dos conjuntos de datos.
2. coord.ellipses: Construye elipses de confianza
bary: Toma un valor lógico. Se calculan las coordenadas de la elipse alrededor del
baricentro de los individuos.
120 CAPÍTULO 6. PAQUETES DE R

3. heatplot: Utiliza la función heatmap2, usando por defecto un rango de colores de rojo
a verde. Además, dibuja dendogramas de los casos y las variables, usando como métrica
similitud entre correlaciones y clustering. Es útil para análisis exploratorio de los datos.
classvec: Factor o vector que describe las clases en filas o columnas de los datos.
dend: Indica si los dendogramas se dibujan para filas, columnas o ambos. También
puede no dibujarse.
scale: Indica la escala a utilizar en función de para qué se dibujen los dendogramas.
Por defecto es fila. Puede ser “none”.
4. ord: Lleva a cabo un ACP (“pca”), Análisis de Correspondencias (“coa”) o Análisis de
Correspondencias Asimétrico (“nsc”) en datos de expresión génica.
type: Caracter “coa”, “pca” o “nsc” indicando la transformación de datos requerida.
Por defecto es “coa”.
5. overview: Representa un boxplot, histograma y dendograma de análisis jerárquico.
classvec: Vector o factor que describe las clases en columnas del conjunto de datos.
labels: Vector, etiquetas para las muestras de los plots.
6. plotarrays: Gráfico de proyecciones de las muestras coloreadas por grupos. Útil para
visualizar coordenadas de arrays resultantes de distintos análisis de datos microarray.

Paquete “mclust”
Basado en clustering, clasificación y estimación de la densidad, incluyendo regularización
Bayesiana.

Paquete “MVA”
Proporciona funciones, conjuntos de datos, análisis y ejemplos del libro An Introduction to
Applied Multivariate Analysis with R (Brian S. Everitt & Torsten Hothorn, Springer, 2011)

Funciones
1. bvbox: Representa boxplot en dos dimensiones.
add: Se añade a un plot existente.

Paquete “stats”
El paquete stats proporciona funciones estadísticas.
121

Funciones
1. cor: Matriz de correlaciones entre columnas de dos conjuntos de datos.
2. lm: Ajusta modelos lineales. Puede usarse para modelos de regresión, análisis de varianza
y de covarianza.
3. princomp: Lleva a cabo un Análisis de Componentes Principales. Devuelve un objeto
del tipo princomp.
cor: Valor lógico que indica si se usa o no la matriz de correlación.
Bibliografía

[1] Cuadras CM. (2014). Nuevos Métodos de Análisis Multivariante. CMC Editions.

[2] Culhane A. (2014). Introduction to Multivariate Analysis of Microarray Gene Expression


Data using MADE4.

[3] Culhane A.(adapted by Sánchez A.) (2015) Cross-platform data integration using Coinertia
Analysis.

[4] Culhane A., Perrière G, G. Higgins D. (2003) Cross-platform comparison and visualisation
of gene expression data using co-inertia analysis. BMC Bioinformatics.

[5] Culhane A., Thioulouse J., Perriere G., Higgins DG. (2005) MADE4:an R package for
multivariate analysis of gene expression data Bioinformatics 21(11): 2789-90.

[6] Dufour A.B. (2008) Enseignements de Statistique en Biologie. URL: http://pbil.univ-


lyon1.fr/R/pdf/course6.pdf

[7] Everitt B., Hothorn T. (2015). An Introduction to Applied Multivariate Analysis with R.
Springer.

[8] González I., Déjean S., G.P. Martin P., Baccini A. (2008) CCA: An R Package to Extend
Canonical Correlation Analysis. Journal of Statistical Software.

[9] González I., Déjean S. (2012). CCA: Canonical correlation analysis. R package version 1.2.
http://CRAN.R-project.org/package=CCA

[10] Härdle W.K., Simar L. (2015). Applied Multivariate Statistical Analysis. Springer-Verlag.

[11] Husson F., Josse J., Le S., Mazet J. (2015). FactoMineR: Multivariate Explora-
tory Data Analysis and Data Mining. R package version 1.31.4. http://CRAN.R-
project.org/package=FactoMineR

[12] Jobson J.D.(1992). Applied Multivariate Data Analysis.


Volume II:Categorical and Multivariate Methods. Springer-Verlag.

[13] Johnson R. A., Wichern D.W. (2002). Applied Multivariate Statistical Analysis. Pearson
Education, Prentice Hall.

122
BIBLIOGRAFÍA 123

[14] Lê Cao K., GP Martin P., Christèle-Granié, Besse P. (2009) Sparse canonical methods for
biological data integration: application to a cross-platform study. BMC Bioinformatics.

[15] Lê S., Josse J., Husson F. (2008). FactoMineR: an R package for multivariate analysis.
Journal of Statistical Software.

[16] Mardia K.V., Kent J.T., Bibby J.M. (1979) Multivariate Analysis. Academic Press.

[17] Peña D.(2002) Análisis multivariante de datos. MC GRAW HILL INTERAMERICANA.

[18] R Core Team (2014). R: A language and environment for statistical computing. R Founda-
tion for Statistical Computing, Vienna, Austria. URL: http://www.R-project.org/.

[19] R Studio URL: http://www.rstudio.org/

[20] Sánchez A. Multivariate Methods for the integrate analysis of omics data. (1): Exploratory
Data Analysis. URL: http://eib.stat.ub.edu/Integrative+Analysis+of+Omics+Data

[21] Sánchez A. Multivariate Methods for the Integration and Visualization of Omics Data.
URL: http://eib.stat.ub.edu/Integrative+Analysis+of+Omics+Data

[22] Santamaria R. Análisis de Datos de Microarray. Universidad de Salamanca.

También podría gustarte