Introducción

En este trabajo se presenta una Introducción al lenguaje de estadística R. Inicialmente se específica qué es R y algunas de las propiedades que tiene este lenguaje de programación. Parte fundamental en cualquier proyecto de investigación donde se consideren conjuntos de datos es hacer una exploración preliminar, esto se logra por medio de diferentes tipos de gráficas; pues estás nos brindan una perspectiva de la distribución de las variables e inclusive nos indican alguna posible relación que pueda brindar una hipótesis. Se exploran las gráficas de barras y sectores para variables cualitativas, posteriormente se presentan los histogramas, diagramas de cajas y diagramas de dispersión para variables cuantitativas.

Con R se pueden llevar a cabo múltiples análisis estadísticos, como estimación de parámetros (por lo regular para estimar la media o la varianza) con sus respectivos intervalos de confianza, pruebas de hipótesis, pruebas de normalidad (para verificar si una variable numérica tiene distribución normal), la prueba \(\chi^2\) de Pearson considerando la distribución de una variable cualitativa o la posible relación entre dos variables. En este trabajo se hace énfasis en tres tipos de análisis estadísticos: Análisis de varianza (ANOVA), análisis de regresión lineal y el análisis de componentes principales (ACP).

En la parte final se presenta Bioconductor, un proyecto de código abierto para el análisis de datos en Genómica. Se expone un caso de abundancias de microorganismos al nivel taxonómico Género para determinar por medio de un análisis de coordenadas principales las diferencias entre individuos afectados con cierto tipo de enfermedades.

¿Qué es R?

R es un lenguaje y entorno para cómputo estadístico y gráficos. Es un proyecto GNU que es similar al lenguaje y entorno S que fue desarrollado en Bell Laboratories (anteriormente AT&T, ahora Lucent Technologies) por John Chambers. R proporciona una amplia variedad de técnicas estadísticas (modelado lineal y no lineal, pruebas estadísticas clásicas, análisis de series temporales, clasificación, agrupamiento,…) y gráficas. Una de las fortalezas de R es la facilidad con la que se pueden producir gráficas bien diseñadas con calidad de publicación, que incluyen símbolos matemáticos y fórmulas donde sea necesario. R está disponible como Software Libre bajo los términos de la Licencia Pública General GNU de la Free Software Foundation en forma de código fuente. Compila y se ejecuta en una amplia variedad de plataformas UNIX y sistemas similares (incluidos FreeBSD y Linux), Windows y MacOS.

R, como S, está diseñado en torno a un verdadero lenguaje informático y permite a los usuarios agregar funcionalidad adicional mediante la definición de nuevas funciones. R se puede extender (fácilmente) a través de paquetes . Hay alrededor de ocho paquetes suministrados con la distribución R y muchos más están disponibles a través de la familia de sitios de Internet CRAN que cubren una amplia gama de estadísticas modernas (R-Project 2022).

Exploración de datos.

Utilizaremos el conjunto de datos iris, desarrollado por Ronald Fisher en el artítulo The use of multiple measurements in taxonomic problems (Fisher, 1936). Este conjunto brinda las medidas en centímetros de las variables longitud y ancho del sépalo, longitud y ancho del pétalo; respectivamente para 50 flores de tres especies distintas de iris. Las especies son iris setosa, versicolor y virginica. Por lo tanto contiene 150 observaciones (filas) y 5 variables (columnas), en la siguiente tabla se muestran las primeras observaciones de este conjunto:

Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa
4.6 3.4 1.4 0.3 setosa
5.0 3.4 1.5 0.2 setosa
4.4 2.9 1.4 0.2 setosa
4.9 3.1 1.5 0.1 setosa
5.4 3.7 1.5 0.2 setosa
4.8 3.4 1.6 0.2 setosa
4.8 3.0 1.4 0.1 setosa
4.3 3.0 1.1 0.1 setosa
5.8 4.0 1.2 0.2 setosa
5.7 4.4 1.5 0.4 setosa
5.4 3.9 1.3 0.4 setosa
5.1 3.5 1.4 0.3 setosa
5.7 3.8 1.7 0.3 setosa
5.1 3.8 1.5 0.3 setosa
5.4 3.4 1.7 0.2 setosa
5.1 3.7 1.5 0.4 setosa
4.6 3.6 1.0 0.2 setosa
5.1 3.3 1.7 0.5 setosa
4.8 3.4 1.9 0.2 setosa
5.0 3.0 1.6 0.2 setosa
5.0 3.4 1.6 0.4 setosa
5.2 3.5 1.5 0.2 setosa
5.2 3.4 1.4 0.2 setosa
4.7 3.2 1.6 0.2 setosa
4.8 3.1 1.6 0.2 setosa
5.4 3.4 1.5 0.4 setosa
5.2 4.1 1.5 0.1 setosa
5.5 4.2 1.4 0.2 setosa
4.9 3.1 1.5 0.2 setosa
5.0 3.2 1.2 0.2 setosa
5.5 3.5 1.3 0.2 setosa
4.9 3.6 1.4 0.1 setosa
4.4 3.0 1.3 0.2 setosa
5.1 3.4 1.5 0.2 setosa
5.0 3.5 1.3 0.3 setosa
4.5 2.3 1.3 0.3 setosa
4.4 3.2 1.3 0.2 setosa
5.0 3.5 1.6 0.6 setosa
5.1 3.8 1.9 0.4 setosa
4.8 3.0 1.4 0.3 setosa
5.1 3.8 1.6 0.2 setosa
4.6 3.2 1.4 0.2 setosa
5.3 3.7 1.5 0.2 setosa
5.0 3.3 1.4 0.2 setosa
7.0 3.2 4.7 1.4 versicolor
6.4 3.2 4.5 1.5 versicolor
6.9 3.1 4.9 1.5 versicolor
5.5 2.3 4.0 1.3 versicolor
6.5 2.8 4.6 1.5 versicolor
5.7 2.8 4.5 1.3 versicolor
6.3 3.3 4.7 1.6 versicolor
4.9 2.4 3.3 1.0 versicolor
6.6 2.9 4.6 1.3 versicolor
5.2 2.7 3.9 1.4 versicolor
5.0 2.0 3.5 1.0 versicolor
5.9 3.0 4.2 1.5 versicolor
6.0 2.2 4.0 1.0 versicolor
6.1 2.9 4.7 1.4 versicolor
5.6 2.9 3.6 1.3 versicolor
6.7 3.1 4.4 1.4 versicolor
5.6 3.0 4.5 1.5 versicolor
5.8 2.7 4.1 1.0 versicolor
6.2 2.2 4.5 1.5 versicolor
5.6 2.5 3.9 1.1 versicolor
5.9 3.2 4.8 1.8 versicolor
6.1 2.8 4.0 1.3 versicolor
6.3 2.5 4.9 1.5 versicolor
6.1 2.8 4.7 1.2 versicolor
6.4 2.9 4.3 1.3 versicolor
6.6 3.0 4.4 1.4 versicolor
6.8 2.8 4.8 1.4 versicolor
6.7 3.0 5.0 1.7 versicolor
6.0 2.9 4.5 1.5 versicolor
5.7 2.6 3.5 1.0 versicolor
5.5 2.4 3.8 1.1 versicolor
5.5 2.4 3.7 1.0 versicolor
5.8 2.7 3.9 1.2 versicolor
6.0 2.7 5.1 1.6 versicolor
5.4 3.0 4.5 1.5 versicolor
6.0 3.4 4.5 1.6 versicolor
6.7 3.1 4.7 1.5 versicolor
6.3 2.3 4.4 1.3 versicolor
5.6 3.0 4.1 1.3 versicolor
5.5 2.5 4.0 1.3 versicolor
5.5 2.6 4.4 1.2 versicolor
6.1 3.0 4.6 1.4 versicolor
5.8 2.6 4.0 1.2 versicolor
5.0 2.3 3.3 1.0 versicolor
5.6 2.7 4.2 1.3 versicolor
5.7 3.0 4.2 1.2 versicolor
5.7 2.9 4.2 1.3 versicolor
6.2 2.9 4.3 1.3 versicolor
5.1 2.5 3.0 1.1 versicolor
5.7 2.8 4.1 1.3 versicolor
6.3 3.3 6.0 2.5 virginica
5.8 2.7 5.1 1.9 virginica
7.1 3.0 5.9 2.1 virginica
6.3 2.9 5.6 1.8 virginica
6.5 3.0 5.8 2.2 virginica
7.6 3.0 6.6 2.1 virginica
4.9 2.5 4.5 1.7 virginica
7.3 2.9 6.3 1.8 virginica
6.7 2.5 5.8 1.8 virginica
7.2 3.6 6.1 2.5 virginica
6.5 3.2 5.1 2.0 virginica
6.4 2.7 5.3 1.9 virginica
6.8 3.0 5.5 2.1 virginica
5.7 2.5 5.0 2.0 virginica
5.8 2.8 5.1 2.4 virginica
6.4 3.2 5.3 2.3 virginica
6.5 3.0 5.5 1.8 virginica
7.7 3.8 6.7 2.2 virginica
7.7 2.6 6.9 2.3 virginica
6.0 2.2 5.0 1.5 virginica
6.9 3.2 5.7 2.3 virginica
5.6 2.8 4.9 2.0 virginica
7.7 2.8 6.7 2.0 virginica
6.3 2.7 4.9 1.8 virginica
6.7 3.3 5.7 2.1 virginica
7.2 3.2 6.0 1.8 virginica
6.2 2.8 4.8 1.8 virginica
6.1 3.0 4.9 1.8 virginica
6.4 2.8 5.6 2.1 virginica
7.2 3.0 5.8 1.6 virginica
7.4 2.8 6.1 1.9 virginica
7.9 3.8 6.4 2.0 virginica
6.4 2.8 5.6 2.2 virginica
6.3 2.8 5.1 1.5 virginica
6.1 2.6 5.6 1.4 virginica
7.7 3.0 6.1 2.3 virginica
6.3 3.4 5.6 2.4 virginica
6.4 3.1 5.5 1.8 virginica
6.0 3.0 4.8 1.8 virginica
6.9 3.1 5.4 2.1 virginica
6.7 3.1 5.6 2.4 virginica
6.9 3.1 5.1 2.3 virginica
5.8 2.7 5.1 1.9 virginica
6.8 3.2 5.9 2.3 virginica
6.7 3.3 5.7 2.5 virginica
6.7 3.0 5.2 2.3 virginica
6.3 2.5 5.0 1.9 virginica
6.5 3.0 5.2 2.0 virginica
6.2 3.4 5.4 2.3 virginica
5.9 3.0 5.1 1.8 virginica

Tipos de variables

Las técnicas de visualización y los resúmenes estadísticos dependen del tipo de la variable a analizar. Tomando en cuenta los valores que una variable puede tomar, se clasifican en dos grupos: variables numéricas y variables categóricas (o cualitativas). El conjunto de datos iris contiene cuatro variables numéricas: Sepal.Length, Sepal.Width, Petal.Length y Petal.Width. La variable Species es cualitativa y consta de tres diferentes categorías (descritas previamente).

En algunos conjuntos de datos las variables numéricas son variables de conteo. También es común utilizar números para codificar a las variables categóricas y pueden ser ordinales o nominales dependiendo de la extensión de información que la interpretación numérica proporciona. Para variables ordinales, tales como el grado de una enfermedad, aunque los números no tienen su significado usual, preservan un rango.

Gráficas

En esta sección se mostrarán las gráficas más comunes para la exploración de datos. Tales gráficas dependen del tipo de variable y proporcionan información para la comprensión de datos, averiguando los posibles valores de cada característica y encontrando cómo tales valores varían entre los elementos de la muestra. Las gráficas fueron realizadas (excepto la de sectores) utilizando el package ggplot2 (Wickham, 2016)

Gráfica de barras.

Para variables categóricas, las gráficas de barras son una de las maneras más simples de visualizar datos. Usando una gráfica de barras, podemos visualizar los valores posibles (categorías) que una variable cualitativa puede tomar, así como el número de veces que cada categoría ha sido observada en nuestra muestra. La altura de cada barra en esta gráfica muestra el número de veces que la categoría correspondiente ha sido observada. Las alturas de las barras suman \(n\), el tamaño de la muestra. Las gráficas de barras nos indican cómo están distribuidos los valores observados de una variable categórica en nuestra muestra. En la siguiente Figura se muestra la gráfica de barras para la variable Species del conjunto de dato iris, como se había indicado previamente hay 50 observaciones de cada especie en el conjunto de datos.

species_color <- c('#1627dc', '#ffb600', '#ff2e16')
ggplot(iris)+
    geom_bar(aes(Species, fill=Species), color="black", show.legend = FALSE)+
    labs(y="Frecuencia")+
    scale_fill_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )
Fig 1. a ) Diagrama de barras para la variable Species

Fig 1. a ) Diagrama de barras para la variable Species

Gráfica de sectores.

Podemos usar una gráfica de sectores para visualizar las frecuencias relativas de diferentes categorías para una variable cualitativa. En una gráfica de sectores, el área de un círculo está dividida en sectores donde cada uno representa las posibles categorías de la variable. El área de cada sector es proporcional a su frecuencia. A continuación se muestra la respectiva gráfica de sectores para la variable Species, donde se observa una distribución uniforme de los datos. Esta gráfica fue realizada con el package plotrix (Lemon, 2006)

tab_Species <- table(iris_tibble$Species) %>% as_vector() %>% unname()

pie3D(tab_Species, labels = c("Setosa", "Versicolor", "Virginica"), col=species_color)

Histogramas.

Los histogramas son usados comúnmente para visualizar variables numéricas, son similares a una gráfica de barras después que los valores de la variable son agrupados en un número finito de intervalos. Para cada intervalo, la altura de la barra corresponde a la frecuencia de la observación en ese intervalo, también es común usar la densidad en lugar de la frecuencia relativa o el porcentaje. Además de la localización (media o mediana) y la dispersión de una distribución, la forma de un histograma también nos muestra cómo los valores observados están localizados. A continuación se muestran los histogramas para cada una de las variables numéricas en el conjunto de datos iris, para obtener más información se identifica la especie. Se puede observar que el ancho del sépalo se distribuye de manera uniforme alrededor de la media, además los valores más altos son para la especie setosa. Por otro lado, la longitud del petalo tiene una distribución bimodal, se pueden identificar los subconjuntos definidos por la especie, en este caso las flores de la especie setosa tienen menor longitud de pétalo y son significativamente diferentes a las otras dos especies.


Sepal. Length

ggplot(iris)+
    geom_histogram(aes(Sepal.Length, fill=Species), color="black", show.legend = TRUE)+
    labs(y="Frecuencia", x="Longitud del sépalo")+
    scale_fill_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Width

ggplot(iris)+
    geom_histogram(aes(Sepal.Width, fill=Species), color="black", show.legend = TRUE)+
    labs(y="Frecuencia", x="Ancho del sépalo")+
    scale_fill_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Petal.Length

ggplot(iris)+
    geom_histogram(aes(Petal.Length, fill=Species), color="black", show.legend = TRUE)+
    labs(y="Frecuencia", x="Longitud del petalo")+
    scale_fill_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Petal.Width

ggplot(iris)+
    geom_histogram(aes(Petal.Width, fill=Species), color="black", show.legend = TRUE)+
    labs(y="Frecuencia", x="Ancho del pétalo")+
    scale_fill_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )



Diagrama de cajas

El diagrama de cajas es una forma gráfica de representar de manera resumida un conjunto de datos numéricos \(x_1,..., x_n\). Esta representación está compuesta por una caja y por dos marcas en dos de sus extremos opuestos que asemejan brazos.Para dibujar estos diagramas se necesita determinar el centro de la caja, su altura y los tamaños de los brazos superior e inferior.

Una manera de construir un diagrama de caja y brazos es a través de los cuantiles. La altura de la caja parte del primer cuartil \(Q_{0.25}\) y se extiende hasta el tercer cuartil \(Q_0.75\). El segundo cuartil \(Q_{0.5}\), es decir, la mediana, se encuentra dentro de la caja pero no necesariamente en el centro.

La altura de la caja es entonces el así llamado rango intercuartil

\[RIC=Q_{0.75}-Q_{0.25}\]

El rango intercuartil mide la longitud del intervalo más pequeño que contiene el 50% de los datos centrales alrededor de la mediana.Por su nombre en inglés, el rango intercuartil también puede denotarse por las letras IQR, interquartile range. Las longitudes de los brazos se puede establecer como 1.5 veces el rango intercuartil RIC, y en este caso, los brazos tienen idéntica longitud. A los valores que se encuentren abajo de la marca del brazo inferior o arriba de la marca del brazo superior se les llama valores atípicos. A los valores que se encuentren arriba de \(Q_{0.75} +3RIC\) o abajo de \(Q_{0.25}-3RIC\) se les llama extremadamente atípicos (outliers).

A continuación se muestran los diagramas de caja para las cuatro variables numéricas en el conjunto iris. Una pregunta que surge es si los valores de las variables varían considerablemente entre el tipo de especie, por ejemplo, para el ancho del pétalo las diferencias parecen significativas, pero no es así para el ancho del sépalo. Posteriormente se llevará a cabo un análisis de varianza (ANOVA) para determinar si tales diferencias son significativas. Para esta sección se utilizó como referencia principal el libro (Rincón 2012a).


Sepal. Length

ggplot(iris)+
    geom_boxplot(aes(x=Species, y=Sepal.Length, color=Species), show.legend = FALSE)+
    labs(y="Longitud del sépalo", x="Especie")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Width

ggplot(iris)+
    geom_boxplot(aes(x=Species, y=Sepal.Width, color=Species), show.legend = FALSE)+
    labs(y="Ancho del sépalo", x="Especie")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Petal.Length

ggplot(iris)+
    geom_boxplot(aes(x=Species, y=Petal.Length, color=Species), show.legend = FALSE)+
    labs(y="Longitud del pétalo", x="Especie")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Petal.Width

ggplot(iris)+
    geom_boxplot(aes(x=Species, y=Petal.Width, color=Species), show.legend = FALSE)+
    labs(y="Ancho del pétalo", x="Especie")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Diagramas de dispersión

El diagrama de dispersión es una descripción gráfica al considerar dos variables numéricas a la vez. En particular, estaremos interesados en detectar alguna posible relación estadística entre dos variables, \(x\) y \(y\). Gráficamente los datos \((x_1, y_1), . . . , (x_n, y_n)\) pueden representarse como puntos en un plano.

Los diagramas de dispersión ayudan a identificar visualmente alguna posible relación entre las dos variables en estudio. En general, pueden existir varios tipos de relaciones entre dos variables pero distinguiremos dos casos: relaciones lineales y relaciones no lineales. Particularmente estaremos interesados en detectar tendencias de relación lineal entre las variables. Diremos entonces que: 1. Existe correlación positiva entre las variables cuando los valores de \(y\) tienden a crecer de manera lineal conforme los valores de \(x\) crecen. 2. Existe correlación negativa entre las variables cuando los valores de \(y\) tienden a decrecer de manera lineal conforme los valores de \(x\) crecen 3. No existe correlación entre las variables cuando ninguna de las dos tendencias anteriores se presenta.

A continuación se muestran los 6 diagramas de dispersión considerando las posibles combinaciones de las 4 variables numéricas del conjunto de datos iris, en cada uno se agrega la línea de regresión que indica alguna posible relación lineal. Por ejemplo al graficar la longitud y ancho del sépalo no se aprecia ningún patrón o relación. Por otro lado, en la gráfica de longitud y ancho del pétalo se identifica una relación lineal positiva entre las variables.

Sepal.Length vs Sepal.Width

ggplot(iris)+
    geom_point(aes(Sepal.Length, Sepal.Width, color=Species), size=2)+
    geom_smooth(aes(Sepal.Length, Sepal.Width), method = "lm", color = "blue", se=FALSE)+
    labs(x="Longitud del sépalo", y="Ancho del sépalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Length vs Petal.Length

ggplot(iris)+
    geom_point(aes(Sepal.Length, Petal.Length, color=Species), size=2)+
    geom_smooth(aes(Sepal.Length, Petal.Length), method = "lm", color = "blue", se=FALSE)+
    labs(x="Longitud del sépalo", y="Longitud del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Length vs Petal.Width

ggplot(iris)+
    geom_point(aes(Sepal.Length, Petal.Width, color=Species), size=2)+
    geom_smooth(aes(Sepal.Length, Petal.Width), method = "lm", color = "blue", se=FALSE)+
    labs(x="Longitud del sépalo", y="Ancho del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Width vs Petal.Length

ggplot(iris)+
    geom_point(aes(Sepal.Width, Petal.Length, color=Species), size=2)+
    geom_smooth(aes(Sepal.Width, Petal.Length), method = "lm", color = "blue", se=FALSE)+
    labs(x="Ancho del sépalo", y="Longitud del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Sepal.Width vs Petal.Width

ggplot(iris)+
    geom_point(aes(Sepal.Width, Petal.Width, color=Species), size=2)+
    geom_smooth(aes(Sepal.Width, Petal.Width), method = "lm", color = "blue", se=FALSE)+
    labs(x="Ancho del sépalo", y="Ancho del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Petal.Length vs Petal.Width

ggplot(iris)+
    geom_point(aes(Petal.Length, Petal.Width, color=Species), size=2)+
    geom_smooth(aes(Petal.Length, Petal.Width), method = "lm", color = "blue", se=FALSE)+
    labs(x="Longitud del pétalo", y="Ancho del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Análisis estadísticos con R

De los múltiples análisis estadísticos que se pueden llevar a cabo utilizando R, en esta sección nos enfocamos en tres: Análisis de varianza (ANOVA), análisis de regresión lineal y el análisis de componentes principales. Este último es un análisis de estadística multivariable.

Análisis de varianza (ANOVA)

Los modelos de análisis de varianza (ANOVA), generalizan el t-test y son usados para comparar las medias de múltiples grupos identificados por una variable categórica con más de dos categorías posibles. La variable categórica es llamada factor y es considerada como la variable explicativa o independiente. La variable numérica, denotada como \(Y\), cuya media se compara entre los diferentes grupos, es considerada la variable de respuesta. Nos enfocaremos en los análisis de un solo factor (conocidos como one-way ANOVA).

Para llevar a cabo el análisis de varianza es necesario calcular diferentes variaciones. La variación entre grupos es denotada por \(SS_B\) y calculada por

\[SS_B=\sum_{i=1}^{k}n_i(\bar{y}_i- \bar{y})^2,\]

donde \(k\) es el número de grupos.

La variación en los grupos es denotada por \(SS_W\) y calculada por

\[SS_W=\sum_{i=1}^{k}\sum_{j=1}^{n_i}(y_{ij}-\bar{y}_i)^2.\]

La variación total en \(Y\) es dada por

\[SS=\sum_{i=1}^{k}\sum_{j=1}^{n_i}(y_{ij}-\bar{y})^2.\]

Entonces, tendremos que la variación total satisface

\[SS=SS_B+SS_W.\]

Se denotará la media poblacional total de \(Y\) por \(\mu\) y las medias de los grupos por \(\mu_i\). Entonces podemos expresar la hipótesis nula de no diferencia en las medias entre los grupos por

\[H_0:\mu_1=\mu_2=\mu_3=\mu_4=\mu.\]

La hipótesis alternativa \(H_A\) es que al menos una de las medias de los grupos \(\mu_i\) es diferente de la media total \(\mu\).

La estadística de prueba para examinar la hipótesis nula es llamada F-estadística y es definida por

\[F=\frac{SS_B/(k-1)}{SS_W/(n-k)},\]

donde \(n\) es el tamaño total de la muestra, y \(k\) es el número de grupos. El numerador es llamado la media cuadrada para grupos, y el denominador es llamado el error medio cuadrático (EMC).

Nótese que la estadística de prueba está basada en la comparación de la variación entre grupos (la cual es explicada por el factor) y la variación en los grupos (la cual es aleatoria). Cuando las medias de los grupos son significativamente diferentes, y su variación es relativamente grande comparada con la variable aleatoria en los grupos, el valor de la estadística es grande.

La distribución \(F\) es una distribución de probabilidad continua y está determinada por dos parámetros \(df_1\) y \(df_2\), y es denotada por \(F(df_1,df_2)\). Nos referimos a \(df_1\) y \(df_2\) como los grados de libertad del numerador y grados de libertad del denominador respectivamente.

Para el ANOVA de un factor, la estadística \(F\) tiene distribución \(F(df_1=k-1,df_2=n-k)\) suponiendo que la hipótesis nula es verdadera. Aquí \(df_1=k-1\) es el número de grupos menos 1, son los grados de libertad del numerador y \(df_2=n-k\), el tamaño de la muestra menos el número de grupos, son los grados de libertad del denominador.

Una de las condiciones para usar este análisis es que las observaciones en cada grupo sean independientes e idénticamente distribuidas (con distribución normal). Una de las suposiciones al realizar el ANOVA es que todos los grupos tienen la misma varianza poblacional \(\sigma^2\), la cual es desconocida.

Cuando las medias de los grupos son muy diferentes una de otra, la variación entre grupos \(SS_B\) es alta. Como resultado, la F-estadística es grande, y dichos valores son considerados como extremos si la hipótesis nula es verdadera y proporcionan evidencia sólida contra tal hipótesis. Para calcular el nivel de significancia observado \(p_{obs}\), se encuentra la probabilidad de los valores tan extremos o más que el valor observado de la estadística de prueba.

Ejemplo. A continuación se llevará a cabo el análisis de varianza para la variable Petal.Length respecto a la variable Species, en este caso el número de observaciones es \(n=150\) y son \(k=3\) grupos determinados por las categorías de la variable Species (setosa, versicolor y virginica). En la Figura se puede apreciar una notable diferencia entre los valores para la longitud del pétalo entre cada especie.

ggplot(iris)+
    geom_boxplot(aes(x=Species, y=Petal.Length, color=Species), show.legend = FALSE)+
    labs(y="Longitud del pétalo", x="Especie")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

Llevando a cabo el análisis de varianza con la función aov de R se obtienen los siguientes resultados:

Petal.Length_aov <- aov(Petal.Length~Species, data=iris)
summary(Petal.Length_aov)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  437.1  218.55    1180 <2e-16 ***
## Residuals   147   27.2    0.19                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
numSummary(iris[,c("Petal.Length"), drop=FALSE], 
           groups=iris$Species, statistics=c("mean", "sd", "IQR", "quantiles"),
           quantiles=c(0,.25,.5,.75,1))
##             mean        sd   IQR  0% 25%  50%   75% 100% Petal.Length:n
## setosa     1.462 0.1736640 0.175 1.0 1.4 1.50 1.575  1.9             50
## versicolor 4.260 0.4699110 0.600 3.0 4.0 4.35 4.600  5.1             50
## virginica  5.552 0.5518947 0.775 4.5 5.1 5.55 5.875  6.9             50

La primera columna muestra los grados de libertad (Df), los cuales son \(k-1=2\) y \(n-k=147\), respectivamente. Los valores de la segunda columna (Sum Sq), son las variaciones entre los grupos y en los grupos: \(SS_B=437.1\) y \(SS_W=27.2\). El valor observado de la estadística de la estadística de prueba, \(F\), es 1180 con un valor \(p<2\times 10^{-16}\). Por lo tanto existe suficiente evidencia estadística para rechazar la hipótesis nula y se puede afirmar que la diferencia en las medias es significativamente distinta entre los grupos. También se agrega una tabla con el resumen estadístico (media, desviación estándar, cuantiles…) para la variable Petal.Length.

El contenido de esta sección está basado en el capítulo 9 de (Shahbaba, 2012).

Análisis de regresión lineal simple

La regresión lineal simple es un modelo paramétrico para predecir una variable de respuesta cuantitativa \(Y\) a partir de una sola variable \(X\), suponiendo que existe una relación lineal entre \(X\) y \(Y\)

\[Y\approx \beta_0+\beta_1 X\]

Donde \(\beta_0\) y \(\beta_1\) con constantes (desconocidas) que representan la intersección y la pendiente del modelo lineal; son denominados coeficientes o parámetros del modelo. Se calculan valores estimados de tales parámetros, los cuales son denotados por \(\hat{\beta_0}\) y \(\hat{\beta}_1\), entonces

\[\hat{y}=\hat{\beta}_0+\hat{\beta}_1 x\]

donde \(\hat{y}\) indica el valor estimado de \(Y\).

Estimación de los coeficientes

Para estimar los coeficientes \(\beta_0\), \(\beta_1\), sean

\[(x_1, y_1), (x_2,y_2), \dots, (x_n, y_n)\]

\(n\) parejas de observaciones, cada una consiste en un valor observado de \(X\) y de \(Y\). El objetivo es obtener valores estimados \(\hat{\beta_0}\), \(\hat{\beta_1}\) tales que el modelo lineal ajuste bien los datos disponibles. El criterio más utilizado es el de mínimos cuadrados.

Sea \(\hat{y}_i=\hat{\beta}_0+\hat{\beta}_1x_i\) el valor estimado para \(Y\) considerando el \(i\)-ésimo valor de \(X\). Entonces \(e_i=y_i-\hat{y}_i\) representa el \(i\)-ésimo residuo. Se define la suma del cuadrado de los residuos (RSS) como

\[RSS=e_1^2+e_2^2+\cdots+e_n^2,\]

o de manera equivalente

\[RSS=(y_1-\hat{\beta}_0-\hat{\beta}_1x_1)^2+(y_2-\hat{\beta}_0-\hat{\beta}_1x_2)^2+\cdots+(y_n-\hat{\beta}_0-\hat{\beta}_1x_n)^2.\]

Los valores estimados de los coeficientes \(\hat{\beta}_0\) y \(\hat{\beta}_1\) calculados por medio de mínimos cuadrados minimizan el RSS (multiplicadores de Lagrange) y están dados por

\[\begin{equation} \begin{aligned} &\hat{\beta}_{1}=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}} \\ &\hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1} \bar{x} \end{aligned} \end{equation}\]

donde

\[\bar{y}=\frac{1}{n}\sum_{i=1}^ny_i \quad \text{ y }\quad \bar{x}=\frac{1}{n}\sum_{i=1}^nx_i\]

son las medias muestrales.

Ejemplo. Se lleva a cabo el análisis de regresión lineal considerando como variable de respuesta a Petal.Width y la variable explicativa a Petal.Length. La gráfica de dispersión muestra una relación lineal positiva entre estas dos variables.

ggplot(iris)+
    geom_point(aes(Petal.Length, Petal.Width, color=Species), size=2)+
    geom_smooth(aes(Petal.Length, Petal.Width), method = "lm", color = "blue", se=FALSE)+
    labs(x="Longitud del pétalo", y="Ancho del pétalo")+
    scale_color_manual(values = species_color)+
    theme_bw()+
 theme(
    axis.text.x=element_text(size=13),
    axis.text.y=element_text(size=13),
    axis.title.x = element_text(face="bold", size=14),
    axis.title.y = element_text(face="bold", size=14)
  )

ggsave("Petalos.jpg", width = 95, height = 60, units="mm", scale=1.5, dpi=520)

Al utilizar la función lm para llevar a cabo el análisis de regresión lineal simple, los resultados son los siguientes:

iris_lm <- lm(Petal.Width~ Petal.Length, data=iris)

summary(iris_lm)
## 
## Call:
## lm(formula = Petal.Width ~ Petal.Length, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56515 -0.12358 -0.01898  0.13288  0.64272 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.363076   0.039762  -9.131  4.7e-16 ***
## Petal.Length  0.415755   0.009582  43.387  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2065 on 148 degrees of freedom
## Multiple R-squared:  0.9271, Adjusted R-squared:  0.9266 
## F-statistic:  1882 on 1 and 148 DF,  p-value: < 2.2e-16

El coeficiente constante \(\beta_0=-0.36\) es significativo (\(p_{value}=4.7\times 10^{-16}\)), de igual manera el coeficiente \(\beta_1=0.415\); este valor indica que por al incrementar un centímetro la longitud del pétalo se estima que haya un incremento de 0.415 centímetros en el ancho del pétalo. La estadística \(R^2=0.9271\) está relacionada con la bondad del ajuste, la estimación es mejor conforme tal valor es cercano a 1.

data_lm <- tibble(
  Species=iris$Species,
  residuos=iris_lm$residuals,
  res_std=scale(iris_lm$residuals),
  Petal.Length=iris$Petal.Length,
  Petal.Width=iris$Petal.Width,
  fitted=iris_lm$fitted.values
)


graph_res_total <- ggplot(data_lm)+
  geom_point(aes(fitted, residuos, color=Species), size=2)+
  geom_smooth(aes(fitted, residuos), size=1.2, color="blue", se=FALSE)+
  scale_color_manual(values = species_color)+
  #scale_x_continuous(breaks = seq(-3, 4, by=1))+
  labs(x="Valor ajustado Petal.Width", y="Residuos", color="Species")+
  #geom_vline(xintercept = c(perPHm, perPHM), linetype = "dashed", color="red", size=1.2)+
  theme_bw()+
  theme(
    #legend.position="top",
    legend.title = element_text(size=11, face="bold"),
    legend.text = element_text(size=10),
    plot.title = element_text(size=14),
    axis.text.x=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.title.x = element_text(face="bold", size=13),
    axis.title.y = element_text(face="bold", size=13),
    #strip.text.x = element_text(size = 11)
    #legend.position = c(.95, .95),
    #legend.justification = c("right", "top"),
    #legend.box.just = "right"
   # legend.margin = margin(6, 6, 6, 6)
  )+
  guides(color=guide_legend(title="Species", override.aes = list(size=3)))
graph_res_total
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

graph_fitted_total <- ggplot(data_lm)+       
  geom_point(aes(Petal.Width, fitted, color=Species), size=2)+
  geom_line(aes(Petal.Width, Petal.Width), size=1.2, color="blue", se=FALSE)+
  scale_color_manual(values = species_color)+
  #scale_x_continuous(breaks = seq(-3, 4, by=1))+
  labs(x="Petal.Width observado", y="Petal.Width ajustado", color="Species")+
  #geom_vline(xintercept = c(perPHm, perPHM), linetype = "dashed", color="red", size=1.2)+
  theme_bw()+
  theme(
    #legend.position="top",
    legend.title = element_text(size=11, face="bold"),
    legend.text = element_text(size=10),
    plot.title = element_text(size=14),
    axis.text.x=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.title.x = element_text(face="bold", size=13),
    axis.title.y = element_text(face="bold", size=13),
    #strip.text.x = element_text(size = 11)
    #legend.position = c(.95, .95),
    #legend.justification = c("right", "top"),
    #legend.box.just = "right"
   # legend.margin = margin(6, 6, 6, 6)
  )+
  guides(color=guide_legend(title="Species", override.aes = list(size=3)))  
## Warning: Ignoring unknown parameters: se
graph_fitted_total

graph_res_hist_total <- ggplot(data_lm)+
  geom_histogram(aes(residuos, fill=Species), color="black")+
  #geom_smooth(aes(fitted, residuos), size=1.2, color="blue", se=FALSE)+
  scale_fill_manual(values = species_color)+
  #scale_x_continuous(breaks = seq(-3, 4, by=1))+
  labs(x="Residuos", y="Frecuencias", fill="Species")+
  #geom_vline(xintercept = c(perPHm, perPHM), linetype = "dashed", color="red", size=1.2)+
  theme_bw()+
  theme(
    #legend.position="top",
    legend.title = element_text(size=11, face="bold"),
    legend.text = element_text(size=10),
    plot.title = element_text(size=14),
    axis.text.x=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.title.x = element_text(face="bold", size=13),
    axis.title.y = element_text(face="bold", size=13),
    #strip.text.x = element_text(size = 11)
    #legend.position = c(.95, .95),
    #legend.justification = c("right", "top"),
    #legend.box.just = "right"
   # legend.margin = margin(6, 6, 6, 6)
  )+
  guides(color=guide_legend(title="Species", override.aes = list(size=3)))  
  

graph_res_hist_total
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

graph_qq_total <- ggplot(data_lm, aes(sample=res_std))+   
  stat_qq_line(aes(color=Species), size=1)+
  stat_qq(aes(color=Species), size=2)+
  #stat_qq(aes(Group=Quality), shape=21)+
  #geom_point(aes(ICA, fitted, color=Quality), size=2)+
  #geom_line(aes(ICA, ICA), size=1.2, color="blue", se=FALSE)+
  scale_color_manual(values = species_color)+
  #scale_x_continuous(breaks = seq(-3, 4, by=1))+
  labs(x="Distribución teórica", y="Distribución muestral", color="Species")+
  #geom_vline(xintercept = c(perPHm, perPHM), linetype = "dashed", color="red", size=1.2)+
  theme_bw()+
  theme(
    #legend.position="top",
    legend.title = element_text(size=11, face="bold"),
    legend.text = element_text(size=10),
    plot.title = element_text(size=14),
    axis.text.x=element_text(size=12),
    axis.text.y=element_text(size=12),
    axis.title.x = element_text(face="bold", size=13),
    axis.title.y = element_text(face="bold", size=13),
    #strip.text.x = element_text(size = 11)
    #legend.position = c(.95, .95),
    #legend.justification = c("right", "top"),
    #legend.box.just = "right"
   # legend.margin = margin(6, 6, 6, 6)
  )+
  guides(color=guide_legend(title="Species", override.aes = list(size=1)))  
  

graph_qq_total

Análisis de componentes principales (ACP).

La idea central del análisis de componentes principales (ACP) es reducir la dimensión de un conjunto de datos que consiste en una considerable cantidad de variables relacionadas entre sí, manteniendo tanto como sea posible la variación presente en el conjunto de datos. Esto se consigue transformando las variables originales en un nuevo conjunto de variables, las componentes principales (PC’s), que no están correlacionadas y están ordenadas de tal forma que las primeras describen la mayor proporción de la varianza en el conjunto de datos original.

Definición y desarrollo de las componentes principales

Sea \(\mathbf{x}\) un vector de \(p\) variables aleatorias. Se considera obtener a partir del conjunto original,algunas (\(\ll p\)) variables que preserven la mayor parte de la información dada por las varianzas, correlaciones o covarianzas.

El primer paso es buscar una función lineal \(\mathbf{\alpha_1\cdot x}\) que tenga máxima varianza, donde \(\mathbf{\alpha_1}\) es un vector de \(p\) constantes \(\alpha_{11}, \alpha_{12}, \dots, \alpha_{1p}\), de tal manera que

\[\mathbf{\alpha_1\cdot x}=\alpha_{11} x_1+\alpha_{12} x_2+\cdots +\alpha_{1p} x_p=\sum_{j=1}^p\alpha_{1j}x_j.\]

Posteriormente, se busca una función \(\mathbf{\alpha_2 \cdot x}\), sin correlación con \(\mathbf{\alpha_1 \cdot x}\) y que tenga máxima varianza, y así sucesivamente, tal que en la \(k\)-ésima etapa se encuentra una función \(\mathbf{\alpha_k \cdot x}\) que tiene máxima varianza y con correlación nula con el complemento \(\mathbf{\alpha_1\cdot x}, \mathbf{\alpha_2\cdot x}, \dots, \mathbf{\alpha_{k-1}\cdot x}\). La \(k\)-ésima variable obtenida \(\mathbf{\alpha_k\cdot x}\) es la \(k\)-ésima componente principal. Se obtienen \(p\) componentes principales esperando que la mayor parte de la variación en \(\mathbf{x}\) sea explicada por \(m\) CP’s, donde \(m\ll p\).

Si un conjunto con \(p(>2)\) variables tienen correlaciones considerables, entonces las primeras componentes principales explican la mayor parte de la variación en las variables originales. Las últimas CP’s identifican las direcciones en las cuales hay poca variación.

Consideremos el caso donde el vector de variables aleatorias \(\mathbf{x}\) tienen una matriz de covarianza \(\mathbf{\Sigma}\) conocida.

Resulta que, para \(k=1,2,\cdots, p\), la \(k\)-ésima componente principal está dada por \(z_k=\mathbf{\alpha_k\cdot x}\) donde \(\mathbf{\alpha_k}\) es un eigenvector de \(\mathbf{\Sigma}\) correspondiente al \(k\)-ésimo mayor eigenvalor \(\lambda_k\). Además, si \(\mathbf{\alpha_k}\) se normaliza (\(\mathbf{\alpha_k\cdot \alpha_k}=1\)), entonces \(var(z_k)=\lambda_k\).

Para desarrollar la forma de las CP’s, consideremos primero \(\mathbf{\alpha_1' x}\); el vector \(\mathbf{\alpha_1}\) maximiza \(var(\mathbf{\alpha_1\cdot x})=\mathbf{\alpha_1' \Sigma \alpha_1}\) restringido a \(\mathbf{\alpha_1\cdot \alpha_1}=1\).

Para maximizar \(\mathbf{\alpha_1' \Sigma \alpha_1}\) sujeta a \(\mathbf{\alpha_1\cdot \alpha_1}=1\), la aproximación usual es utilizar multiplicadores de Lagrange. Maximizar

\[\mathbf{\alpha_1' \Sigma \alpha_1}-\lambda(\mathbf{\alpha_1\cdot \alpha_1}-1),\]

donde \(\lambda\in \mathbb{R}\) es un multiplicador de Lagrange. Al derivar respecto a \(\mathbf{\alpha_1}\) resulta

\[\mathbf{\Sigma\alpha_1}-\lambda \mathbf{\alpha_1}=0,\]

o bien

\[(\mathbf{\Sigma}-\lambda \mathbf{I_p})\mathbf{\alpha_1}=\mathbf{0}.\]

Por lo tanto, \(\lambda\) es un eigenvalor de \(\mathbf{\Sigma}\) y \(\mathbf{\alpha_1}\) es el eigenvector correspondiente. Para decidir cuál de los \(p\) eigenvectores brinda \(\mathbf{\alpha_1'x}\) con varianza máxima, nótese que la cantidad a ser maximizada es

\[\mathbf{\alpha_1' \Sigma \alpha_1}=\mathbf{\alpha_1'}\lambda \mathbf{\alpha_1}=\lambda \mathbf{\alpha_1'}\mathbf{\alpha_1}=\lambda,\]

así, \(\lambda\) debe ser tan grande como sea posible. Luego, \(\mathbf{\alpha_1}\) es el eigenvector correspondiente al mayor eigenvalor de \(\mathbf{\Sigma}\), y \(var(\mathbf{\alpha_1\cdot x})=\mathbf{\alpha_1'\Sigma \alpha_1}=\lambda_1\), el eigenvalor mayor.

En general, la \(k\)-ésima componente principal de \(\mathbf{x}\) es \(\mathbf{\alpha_k\cdot x}\) y \(var(\mathbf{\alpha_k\cdot x})=\lambda_k\), donde \(\lambda_k\) es el \(k\)-ésimo eigenvalor mayor de \(\mathbf{\Sigma}\), y \(\mathbf{\alpha_k}\) es el correspondiente eigenvector.

La segunda componente principal \(\mathbf{\alpha_2\cdot x}\), maximiza \(\mathbf{\alpha_2'\Sigma \alpha_2}\) sujeta a tener covarianza nula con \(\mathbf{\alpha_1}\), es decir, \(cov(\mathbf{\alpha_1\cdot x, \alpha_2\cdot x})=0\), pero

\[cov(\mathbf{\alpha_1\cdot x, \alpha_2\cdot x})=\mathbf{\alpha_1' \Sigma \alpha_2}=\mathbf{\alpha_2' \Sigma \alpha_1}=\mathbf{\alpha_2'}\lambda_1 \mathbf{\alpha_1'}=\lambda_1\mathbf{\alpha_2'\alpha_1}=\lambda_1\mathbf{\alpha_2'\alpha_1}\]

Por lo tanto, cualquiera de las ecuaciones,

\[\begin{equation} \begin{array}{rr} \alpha_{1}^{\prime} \Sigma \alpha_{2}=0, & \alpha_{2}^{\prime} \Sigma \alpha_{1}=0 \\ \alpha_{1}\cdot \alpha_{2}=0, & \alpha_{2}\cdot \alpha_{1}=0 \end{array} \end{equation}\]

podría ser usada para especificar la correlación cero entre \(\mathbf{\alpha_1\cdot x}\) y \(\mathbf{\alpha_2\cdot x}\). Escogiendo la última de las ecuaciones y considerando que la restricción de normalización es necesaria otra vez, la cantidad a ser maximizada es

\[\mathbf{\alpha_2'\Sigma\alpha_2}-\lambda (\mathbf{\alpha_2'\lambda_2}-1)-\phi \mathbf{\alpha_2'\alpha_1},\]

donde \(\lambda\), \(\phi\) son multiplicadores de Lagrange. Al derivar respecto a \(\mathbf{\alpha_2}\), resulta

\[\mathbf{\Sigma \alpha_2}-\lambda \mathbf{\alpha_2}-\phi \mathbf{\alpha_1}= \mathbf{0}\]

multiplicando esta ecuación por la izquierda con \(\mathbf{\alpha_1'}\) se tiene

\[\mathbf{\alpha_1'\Sigma \alpha_2}-\lambda \mathbf{\alpha_1'\alpha_2}-\phi \mathbf{\alpha_1'\alpha_1}= \mathbf{0}\]

la cual, dado que los primeros dos términos son cero y \(\mathbf{\alpha_1'\alpha_1}=1\), se reduce a \(\phi=0\). Por lo tanto \(\mathbf{\Sigma\alpha_2}-\lambda \mathbf{\alpha_2}=\mathbf{0}\), o de manera equivalente \((\mathbf{\Sigma}-\lambda \mathbf{I}_p)\mathbf{\alpha_2}=\mathbf{0}\), luego, \(\lambda\) es un eigenvalor de \(\mathbf{\Sigma}\) y \(\mathbf{\alpha_2}\) el eigenvector correspondiente.

Puede ser demostrado que la tercera, cuarta, …, \(p\)-ésima componentes principales, los vectores de coeficientes \(\mathbf{\alpha_3}, \mathbf{\alpha_4}, ...,\mathbf{\alpha_p}\) son eigenvectores de \(\mathbf{\Sigma}\) correspondientes a \(\lambda_3, \lambda_3, ..., \lambda_p\) ordenados de mayor a menor, respectivamente. Además

\[var(\mathbf{\alpha_k'x})=\lambda_k\quad \text{ para } k=1,2,..., p.\]

En ocasiones los vectores \(\mathbf{\alpha_k}\) son referidos como componentes principales, pero es preferible reservar tal término para las variables \(\mathbf{\alpha_k'x}\) y referirnos a \(\mathbf{\alpha_k}\) como el vector de coeficientes (o loadings) para la \(k\)-ésima componente principal.

Ejemplo. Consideraremos el conjunto de datos iris, de las cuatro variables numéricas: Sepal.Length, Sepal.Width, Petal.Length y Petal.Width; se realizan combinaciones lineales que serán las componentes principales de tal manera que describan de mejor manera la variabilidad del conjunto de datos y que estás nuevas variables tengan correlación cero.

iris_num <- dplyr::select(iris, Sepal.Length:Petal.Width) #Se seleccionan las variables numéricas
iris_pca <- PCA(iris_num, ncp=4, graph=FALSE) #Se lleva a cabo el PCA
iris_eig <- get_eigenvalue(iris_pca) # Se obtiene la tabla de eigenvalores
iris_var <- get_pca_var(iris_pca)

En la siguiente tabla se muestran los eigenvalores, nótese que las dos primeras componentes capturan el 95.81% de la varianza acumulada, por lo tanto son suficientes para una descripción óptima de los datos. Esto también se ve reflejado en el screeplot.

kable(iris_eig, align = "c") %>% kable_styling("condensed", full_width = T, html_font = "helvetica") %>% scroll_box(width="100%", height="300px", fixed_thead = TRUE)
eigenvalue variance.percent cumulative.variance.percent
Dim.1 2.9184978 72.9624454 72.96245
Dim.2 0.9140305 22.8507618 95.81321
Dim.3 0.1467569 3.6689219 99.48213
Dim.4 0.0207148 0.5178709 100.00000 
fviz_screeplot(iris_pca, addlabels=TRUE)

La tabla y la gráfica que se muestran a continuación, contienen la contribución de cada variable para cada componente principal. En la primera contribuyen en mayor medida las variables Sepal.Length, Petal.Length y Petal.Width. En la segunda la que tiene una mayor contribución es Sepal.Width

kable(iris_var$contrib, align = "c") %>% kable_styling("condensed", full_width = T, html_font = "helvetica") %>% scroll_box(width="100%", height="300px", fixed_thead = TRUE)
Dim.1 Dim.2 Dim.3 Dim.4
Sepal.Length 27.150969 14.2444057 51.777574 6.827052
Sepal.Width 7.254804 85.2474875 5.972245 1.525463
Petal.Length 33.687936 0.0599839 2.019991 64.232089
Petal.Width 31.906291 0.4481230 40.230190 27.415396 
corrplot(iris_var$contrib, is.corr =FALSE, col=species_color)

En la siguiente Figura se muestra el biplot para las dos primeras componentes principales. se grafican los scores de cada una de las observaciones y los vectores para las 4 variables. Se observa una alta correlación entre las variables Petal.Length y Petal.Width, pues los respectivos vectores son casi paralelos, por otro lado, estas dos variables tienen una baja correlación con Sepal.Width (los vectores tienen un ángulo cercano a los 90 grados). El biplot también indica que la especie setosa tiene los valores más bajos para Sepal.Length, Petal.Length y Petal.Width; pero los más altos para Sepal.Width. Las flores de la especie virginica tienen las observaciones mayores para Sepal.Length, Petal.Length y Petal.Width; mientras que la especie versicolor tiene valores cercanos a la media para las cuatro variables. A pesar que el análisis de componentes principales no es un método para distinguir subconjuntos, en este caso son muy evidentes los grupos determinados por la variable Species.

fviz_pca_biplot(iris_pca, geom.ind = "point", fill.ind = iris$Species, pointshape=21,
                point_size=2.5, color.var="black", label="var", repel=TRUE, palette = species_color)

Bioconductor

Bioconductor es un proyecto de código abierto para el análisis de datos en Genómica. La misión del proyecto Bioconductor es desarrollar, respaldar y difundir software gratuito de código abierto que facilite el análisis riguroso y reproducible de datos de ensayos biológicos. La mayoría de los componentes de Bioconductor se distribuyen como paquetes del lenguaje R. El alcance funcional de los paquetes de bioconductores incluye el análisis de microarreglos de ADN, secuencias, flujo, SNP y otros datos.

Los objetivos generales del proyecto Bioconductor son: Proporcionar acceso generalizado a una amplia gama de potentes métodos estadísticos y gráficos para el análisis de datos genómicos. Para facilitar la inclusión de metadatos biológicos en el análisis de datos genómicos, por ejemplo, datos de literatura de PubMed, datos de anotación de genes Entrez. Proporcionar una plataforma de software común que permita el rápido desarrollo e implementación de software extensible, escalable e interoperable. Fomentar la comprensión científica mediante la producción de documentación de alta calidad e investigación reproducible. Formar investigadores en métodos computacionales y estadísticos para el análisis de datos genómicos.

A manera de ejemplo se considera un conjunto de datos para determinar cuáles organismos (a diferente nivel taxonómico) se encuentran en mayor abundancia en individuos (79 observaciones), de dos distintas nacionalidades (MX mexicanos, ESP españoles) y diferente condición médica (Saludables denotados por H y pacientes con colitis ulcerosa clínica inespecífica ICUC y enfermedad de Crohn CD); y extraer información relevante respecto a las diferencias encontradas, principalmente respecto a la estatus sano-paciente. Para ello se lleva a cabo un análisis de coordenadas principales.

Análisis de coordenadas principales (PCoA)

El PCA que está basado en la matriz de correlación o cavarianza, lo cual no siempre es una medida apropiada de asociación. El análisis de coordenadas principales (PCoA) que, como el PCA, se fundamenta en los eigenvalores de una matriz, pero puede cuantificar asociación entre observaciones.

El objetivo del PCoA es calcular una matriz de distancias y producir una configuración gráfica en un espacio euclídeo de dimensión baja (por lo regular 2 o 3), tal que las distancias entre los puntos (observaciones) en la gráfica reflejen las distancias originales tanto como sea posible (Zuur, 2007).

En las siguientes figuras se muestran las gráficas de los dos primeros ejes del análisis de coordenadas principales (PCoA), que explican el 28.7% de la variación de la matriz de distancias. En la primera figura se distingue la nacionalidad, no existe intersección entre las elipses a un nivel de confianza del 95% (dichas elipses determinan el área en la cual se encuentra la media de las coordenadas de cada grupo con una probabilidad de 0.95). De manera análoga en la segunda figura se distinguen los grupos respecto a la enfermedad donde tampoco existe intersección entre las elipses. Se llevó a cabo un análisis de similitud (ANOSIM) obteniendo valores de la estadística \(R=0.1535\) en el caso de nacionalidad y \(R=0.3041\) respecto a los grupos H-ICUC-CD (es decir que es mayor diferencia entre los grupos respecto a la enfermedad).

ggplotly(pcoa_data$graph_Country)
ggplotly(pcoa_data$graph_Disease)

Referencias.

Fisher, R. A. (1936) The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7, Part II, 179–188.

Lemon J.(2006). “Plotrix: a package in the red light district of R.” R-News, 6(4), 8-12.

R-Project (2022) https://www.r-project.org/about.html

Rincón, L. (2012a). Curso elemental de probabilidad y estadística. Serie Textos. Sociedad Matemática Mexicana.

Shahbaba, B. (2012). Biostatistics with R. An Introductión to Statistics Throuhg Biological Data. Springer, New York, NY. ISBN:978-1-4614-1301-1. DOI: https://doi.org/10.1007/978-1-4614-1302-8

Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.

Zuur, A. F.; Ieno, E. N. (2007). Analysing Ecological Data. Springer, New York, NY. ISBN: 978-0-387-45972-1. DOI: https://doi.org/10.1007/978-0-387-45972-1 .