Muñoz Armayones, Sandra TFG PDF
Muñoz Armayones, Sandra TFG PDF
Muñoz Armayones, Sandra TFG PDF
Dirigido por:
Dña. Inmaculada Barranco Chamorro
2016
Resumen
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.
3
4
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.
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.
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
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
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.
Metabolómica: Mide niveles de metabolitos, siendo éstas las moléculas que intervienen en
las reacciones de las células de un organismo.
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.
En el proceso, existen varios elementos que pueden interferir en la medida de la expresión génica.
Estos son:
La densidad de impresión.
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:
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
Esta base de datos está organizada por función molecular, proceso biológico y componente
celular.
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.
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,...).
Estos métodos clásicos se aplican cuando el número de individuos es mucho mayor que el de
variables.
Se asume:
Variables independientes
Normalidad multivariante
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.
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
E(X) = (µ1 , . . . , µp )0 = µ
17
18 CAPÍTULO 1. ANÁLISIS DE COMPONENTES PRINCIPALES
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
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:
Varianza poblacional:
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
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
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.
p p
X X
3. T raza(Σ) = σjj = T raza(Λ) = λj
j=1 j=1
1.2. COMPONENTES PRINCIPALES 21
Demostración:
ΣE = EΛ
E0 ΣE = Λ
λ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:
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}.
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.
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
Σt1 = λIt1
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:
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:
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.
Yj = e0j X, j = 1, . . . , p
λ
Pp j
s=1 λs
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
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
b 0 xi ,
yi = E i = 1, . . . , n
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
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
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.
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.
Gracias a este gráfico podremos visualizar las variables más correladas con las componentes
principales que recojan más información.
= Σ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.
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.
library(mclust)
library(FactoMineR)
data(banknote)
str(banknote)
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.
summary(banknote)
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.
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)
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.
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
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.
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
Con la columna “Dist” se indica la distancia de cada categoría al centro de los ejes del gráfico.
Una vez vistos los resultados obtenidos con el comando summary, podríamos; por lo tanto,
definir las dos primeras CP como:
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).
plot(res.pca,choix = "var")
Length
0.5
Left
Diagonal Right
Dim 2 (21.30%)
0.0
Top
Bottom
−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.
|R − λI| = 0 ⇔ (1 − λ)2 − ρ2 = 0
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
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.
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.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.
data(USairpollution)
names(USairpollution)
str(USairpollution)
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)
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.
USairpollution$negtemp<-USairpollution$temp*(-1)
USairpollution$temp<-NULL
cor(USairpollution[,-1])
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)
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
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.
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
−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
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
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
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.
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
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.
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.
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
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.
summary(res.pca)
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
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
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
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 )).
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.
superior a 0,8.
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
Dim 1 (19.63%)
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)
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
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.
## [1] "data.frame"
dim(bcData)
## [1] 18 12630
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.
bcData<- bcData[-(1:2),]
head(names(bcData))
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)
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.
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
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 |
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
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
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
aa=cbind.data.frame(bcData[,3],res.pca$ind$coord)
bb=coord.ellipse(aa, bary=TRUE)
plot.PCA(res.pca, habillage=3, ellipse=bb)
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
2.1. Descomposición
Comenzaremos definiendo la descomposición en valores singulares de una matriz X.
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
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.
α1 0 ··· 0
0 α2 ··· 0
D= , siendo α1 ≥ α2 · · · ≥ αp los valores singulares de X.
.. .. .. ..
. . . .
0 0 · · · αp
2
c) Los cuadrados de los valores singulares de X, α12 , α22 , · · · , αp , son los autovalores
de X0 X y XX0 .
Demostración:
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.
d) X es simétrica si y sólo si X0 = X
Queda probado por tanto que en este caso los vectores singulares izquierdos y
derechos para X son los mismos.
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.
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.
min(tr[(X − X)(X
b b 0 ])
− X)
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.
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).
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
H0 : P1 (Bj ) = · · · = Pn (Bj ), j = 1, · · · , p
B1 ··· Bj ··· Bp
•Independencia de caracteres
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:
gráfico.
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:
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
(Nij −Eij )
Sea C una matriz de elementos cij = 1/2 , calcularemos su descomposición en valores
Eij
singulares.
C = ΓD∆0 ,
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)}.
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 .
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
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.
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.
A partir de las ecuaciones que definen las coordenadas principales rk y sk , sabemos que:
r0k Ark = λk , s0k Bsk = λk
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
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.
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).
N.j s2kj
Ca (j, sk ) =
λk
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
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 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 .
Datos
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
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
summary(khan)
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.
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
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.
k.coa<-ord(k.data, type="coa")
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).
k.coa$ord$li[1:2,1:3]
k.coa$ord$co[1:2,1:3]
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)
head(cumsum(k.coa$ord$eig*100/sum(k.coa$ord$eig)))
Por tanto, el 30.63076 % de la varianza está explicada por los dos primeros ejes.
• 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.
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
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
Análisis de Coinercia
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 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
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.
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.
(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.
inX = traza(XQX X0 D)
inY = traza(YQY Y0 D)
80 CAPÍTULO 4. ANÁLISIS DE COINERCIA
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.
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.
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
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.
library(made4)
data(NCI60)
class(NCI60)
## [1] "list"
names(NCI60)
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)
Veamos los datos de ambos grupos para las 4 primeras variables y las 4 primeras muestras.
NCI60$Ross[1:4,1:4]
NCI60$Affy[1:4,1:4]
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,]
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.
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]
data.coa1$ord$co[1:4,1:3]
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)))
Por tanto, el 29.47025 % de la varianza está explicada por los dos primeros ejes.
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)))
Por tanto, el 27.38023 % de la varianza está explicada por los dos primeros ejes
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
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
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
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
Aunque ambos conjuntos tengan genes diferentes, se espera que muestren patrones y tendencias
similares.
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
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
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
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
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
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:
Si suponemos que :
X µ ΣXX ΣXY
∼ ,
Y ν ΣYX ΣYY
V ar(X) = ΣXX(p×p)
V ar(Y ) = ΣYY(q×q)
máx{a01 ΣXY b1 }
a1 ,b1
Puesto que ΣXY 0 es equivalente a ΣYX , podemos afirmar que λ = µ, y por tanto el sistema de
ecuaciones resultante es:
Teorema 5.1.1 Los primeros vectores canónicos verifican las siguientes relaciones
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
K = ΓD∆0
−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:
−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.
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.
−1/2
bi = ΣYY δ i ,
con i = 1, . . . , m (m = min{p, q})
Una vez que conocemos el valor exacto de nuestros índices, podemos definir las variables de
correlación canónica
ηi = a0i X
ϕi = b0i Y
1/2
siendo ρi = λi el correspondiente coeficiente de correlación canónica, con i = 1, . . . , m.
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:
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
•Test de Bartlett-Lawley:
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
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.
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
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
Puede establecerse esta igualdad porque el estadístico es invariante por transformaciones lineales.
|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 ϕ.
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
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.
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.
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"
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.
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
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)
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
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.
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)
A continuación, se muestran algunos de los valores que toman los individuos cuando introducimos
los coeficientes obtenidos.
head(res.cc$scores$xscores)
head(res.cc$scores$yscores)
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)
## [,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
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.
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
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
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.
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.
[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.
[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
[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.
[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/.
[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