En este trabajo vamos a estar explorando el conjunto de datos NCI-60, que es utilizado en un gran número de investigaciones de biología y farmacología. Fue creado por el Instituto Nacional del Cáncer de EE. UU. y tiene información sobre 64 líneas celulares humanas de varios tipos de cáncer, como el de mama, colon, pulmón, entre otros. Estas células han sido clave para estudiar cómo diferentes fármacos afectan a los distintos tipos de cáncer. Nos vamos a enfocar en específico en el cáncer de colon y vamos a hacer varias pruebas, como gráficas, análisis univariados y multivariados, pruebas analíticas y pruebas de hipótesis, para entender mejor estos datos y sus implicaciones.

# Insertamo data de COLON 

library(ISLR)
## Warning: package 'ISLR' was built under R version 4.3.3
data(NCI60)

NCI60 <- NCI60
data1 <- data.frame(t(NCI60$data))
names <- NCI60$labs
a <- which(names=="COLON")
colon <- data1[,a]
colon
# Visualización de los datos

head(colon)
##           V42    V43           V44         V45    V46    V47        V48
## 1 -0.32003900 -0.620 -4.900000e-01  0.07001953 -0.120 -0.290 -0.8100195
## 2  0.08996101  0.080  4.200000e-01 -0.82998050  0.000  0.030  0.0000000
## 3 -0.29003900  0.140 -3.400000e-01 -0.59998050 -0.010 -0.310  0.2199805
## 4 -1.03003900  0.740  7.000000e-02 -0.90998050  0.130  1.500  0.7399805
## 5  0.09496101  0.205 -2.050000e-01  0.24501950  0.555  0.005  0.1149805
## 6  0.05996101  0.000 -1.387779e-17 -0.43998050 -0.550 -0.540  0.1199805

1. Presentar por lo menos 2 gráficos que ayuden a explicar un poco su respectiva base de datos.

Estaremos utilizando un diagrama resumen y su correlación y una gráfica de boxplot. Explicaremos e interpretaremos las mismas.

# 1. Diagrama resumen y su correlación

library(GGally)
## Warning: package 'GGally' was built under R version 4.3.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(colon)

Lo más notable en esta gráfica son los picos extremadamente altos en los histogramas de la diagonal, lo que sugiere que, en nuestro análisis univariado, es probable encontrar valores de curtosis mayores a 3 (esto lo confirmaremos más adelante). Estos picos altos nos dice que hay una gran concentración de valores en un rango específico. Esta concentración también se observa en los gráficos de dispersión, donde los datos tienden a agruparse en un rango específico.

Sin embargo, los puntos en los gráficos de dispersión no presentan un alineamiento claro, lo que sugiere que no hay relaciones lineales fuertes entre las variables. Esto se confirma con los valores de correlación, donde la mayoría de las correlaciones son débiles o moderadas. Este patrón nos dice que, aunque existe cierta relación entre las variables, no es lo suficientemente fuerte como para considerar una dependencia directa entre ellas.

# 2. Gráfica de Boxplot

boxplot(colon)

Las medianas parecen estar centradas alrededor de 0 para todas las variables. En los boxplots que proporcionaste, parece que ninguna variable tiene una inclinación evidente hacia los valores altos de manera significativa. Sin embargo, la V46 y V47 tienen más outliers hacia los valores altos en comparación con los valores bajos, lo que podría indicar una ligera inclinación hacia la derecha, pero no es muy pronunciada ya que los bigotes son bastante simétricos.

Estos bigotes simétricos indican que la dispersión en los datos fuera del IQR es bastante uniforme.La presencia de numerosos outliers sugiere la existencia de valores extremos que podrían requerir una mayor investigación, como su origen o si son errores. Las distribuciones parecen en su mayoría simétricas, con algunas diferencias en las colas de las distribuciones.

2. Realizar el análisis univariado y multivariado de las variables correspondientes.

Comenzaremos nuestro analisis con el analisis univariado. Este tipo de analisis se enfoca en una única variable sin consideraar la influencia de otras. Es una herramienta util para comprender mejor nuestros datos antes de otros análisis.

Análisis Univariado:

# Calculamos la media, varianza, desviación estándar, asimetría y curtosis para la variable COLON

library(e1071) # Descargamos la paquetería e1071 para sacar la asimetría y curtosis
## Warning: package 'e1071' was built under R version 4.3.3
media <- apply(colon, 2, mean)
varianza <- apply(colon, 2, var)
desviacion_std <- apply(colon, 2, sd)
asimetria <- apply(colon, 2, skewness)
curtosis <- apply(colon, 2, kurtosis)

Descriptiva <- rbind(media,varianza,desviacion_std,asimetria,curtosis)
round(Descriptiva, 2)
##                 V42   V43   V44  V45   V46  V47  V48
## media          0.02 -0.02 -0.05 0.01  0.01 0.01 0.01
## varianza       0.32  0.62  0.45 0.64  0.52 0.90 0.66
## desviacion_std 0.56  0.79  0.67 0.80  0.72 0.95 0.81
## asimetria      0.88 -0.16 -0.64 0.40  1.11 0.42 0.39
## curtosis       8.32  4.43  7.61 7.36 12.64 6.45 7.67

Podemos apreciar cómo V43 (-0.02) y V44 (-0.05) tienen medias negativas, lo que significa que los valores de estas variables se encuentran ligeramente por debajo de cero. En contraste, podemos ver que las demás variables, como V42 (0.02), V45 (0.01), V46 (0.01), V47 (0.01) y V48 (0.01), tienen medias positivas, lo que indica que los valores de esas variables están relativamente por encima de cero.

En las varianzas, vemos que el valor más alto se encuentra en V47 (0.90), lo que quiere decir que los valores de esta variable tienen una mayor dispersión en comparación con las demás variables. La variable con menor varianza es V42 (0.32), lo que nos dice que sus valores son más consistentes y se agrupan más cerca de la media.

Igualmente, con la desviación estándar. La variable más volátil y con mayor variabilidad es V47 (0.95), mientras que la variable más consistente y homogénea es V42 (0.56).

Respecto a la asimetría, V42 (0.88), V45 (0.40), V46 (1.11), V47 (0.42) y V48 (0.39) presentan una asimetría positiva, lo que indica que hay más valores extremos altos en estas variables. Por otro lado, V43 (-0.16) muestra una asimetría casi simétrica, y V44 (-0.64) tiene una asimetría negativa, sugiriendo que tiene más valores extremos bajos.

Finalmente, en cuanto a la curtosis, V46 (12.64) tiene la curtosis más alta, lo que indica una distribución con colas muy pesadas y sugiere la presencia de muchos outliers. Las variables V42 (8.32), V44 (7.61), V45 (7.36), V47 (6.45) y V48 (7.67) también presentan curtosis alta, indicando todas las variables tienen una gran concentración de valores en lugares en específico, cosa que discutimos al comienzo de este trabajo.

Ahora, otro análisis a considerar es el de la media en comparación a la mediana. Lo haremos a continuación.

# Análisis Media y Mediana

Medianas <-  apply(colon,2,median)
Medias   <- media
Tendencia <- rbind(Medias,Medianas)
Tendencia
##                 V42         V43         V44        V45        V46         V47
## Medias   0.01684022 -0.02274423 -0.04529071 0.01348149 0.00897246 0.009429658
## Medianas 0.00000000  0.00000000  0.00000000 0.00000000 0.00000000 0.000000000
##                  V48
## Medias   0.005350211
## Medianas 0.000000000

Podemos observar que la media de V42 es 0.02, mientras que la mediana es 0.00. Esto sugiere que, aunque la media es positiva, lo que nos dice que hay algunos valores que se encuentran por encima de cero, la mediana se mantiene en cero, lo que implica que la mayoría de los valores están concentrados alrededor de cero. Esta diferencia entre la media y la mediana puede señalar una ligera asimetría en la distribución de los datos.

En el caso de V43, la media es -0.02 y la mediana también es 0.00. Esto indica que los valores de esta variable tienden a estar por debajo de cero, y la media negativa se ve influenciada por algunos valores extremos que afectan el promedio. Sin embargo, la mediana en cero sugiere que al menos la mitad de los valores son cero o positivos, lo que refuerza la idea de que la distribución puede ser más compleja de lo que parece.

Para V44, la media es -0.05 y la mediana es 0.00. Aquí, la media negativa refuerza la noción de que la mayoría de los valores son menores o iguales a cero. La similitud entre la media y la mediana sugiere que la distribución está igualmente afectada por valores extremos negativos.

En V45, V46, V47 y V48, todas presentan medias positivas (0.01, 0.01, 0.01 y 0.01, respectivamente) con medianas en 0.00. Esto significa que, aunque hay algunos valores que se sitúan por encima de cero, la mayoría de los datos no logran superar este umbral, como lo indican las medianas. La media positiva en estas variables sugiere la existencia de algunos valores altos que influyen en el promedio, pero no lo suficiente como para que la mediana se desplace hacia un valor positivo.

En resumen, la relación entre la media y la mediana en estas variables indica que existen diferencias significativas en la distribución de los datos. Las medias tienden a ser más influenciadas por los valores extremos, mientras que las medianas ofrecen una visión más robusta de la tendencia central de los datos. Este comportamiento es clave para comprender la naturaleza de las distribuciones y cómo los valores extremos pueden alterar nuestra interpretación de las medidas de tendencia central.

Visualicemos estas relaciones y distribuciones.

# Visualización variables con Histogramas

par(mfrow=c(2,2))

hist(colon$V42, main = "V42",xlab=" ")
hist(colon$V43, main = "V43",xlab=" ")
hist(colon$V44, main = "V44",xlab=" ")
hist(colon$V45, main = "V45",xlab=" ")

hist(colon$V46, main = "V46", xlab = "")
hist(colon$V47, main = "V47", xlab = "")
hist(colon$V48, main = "V48", xlab = "")

Nuevamente, y con los análisis ya realizados, se puede observar el gran pico de los histogramas. Ya que analizamos las curtosis, sabemos que esto ocurre porque son curtosis altas. Podemos visualizar con mayor claridad que un gran volumen de valores se concentran en lugares en específico. Algo más que podemos abundar es que los valores no tienen sesgos. Se ven en su mayoría en el centro del histograma, esto hace sentido ya que todas las medianas dieron cero.

Aprovechemos para ver si estos histogramas cumplen con la normalidad.

# Normalidad de Histogramas

par(mfrow=c(2, 2))

# Histograma y línea de normalidad para V42
hist(colon$V42, probability=TRUE, main = "V42", xlab=" ")
x <- seq(min(colon$V42)-1, max(colon$V42)+1, length = 100)
y <- dnorm(x, mean(colon$V42), sd(colon$V42))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V43
hist(colon$V43, probability=TRUE, main = "V43", xlab=" ")
x <- seq(min(colon$V43)-1, max(colon$V43)+1, length = 100)
y <- dnorm(x, mean(colon$V43), sd(colon$V43))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V44
hist(colon$V44, probability=TRUE, main = "V44", xlab=" ")
x <- seq(min(colon$V44)-1, max(colon$V44)+1, length = 100)
y <- dnorm(x, mean(colon$V44), sd(colon$V44))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V45
hist(colon$V45, probability=TRUE, main = "V45", xlab=" ")
x <- seq(min(colon$V45)-1, max(colon$V45)+1, length = 100)
y <- dnorm(x, mean(colon$V45), sd(colon$V45))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V46
hist(colon$V46, probability=TRUE, main = "V46", xlab=" ")
x <- seq(min(colon$V46)-1, max(colon$V46)+1, length = 100)
y <- dnorm(x, mean(colon$V46), sd(colon$V46))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V47
hist(colon$V47, probability=TRUE, main = "V47", xlab=" ")
x <- seq(min(colon$V47)-1, max(colon$V47)+1, length = 100)
y <- dnorm(x, mean(colon$V47), sd(colon$V47))
lines(x, y, col = "red", lwd = 3)

# Histograma y línea de normalidad para V48
hist(colon$V48, probability=TRUE, main = "V48", xlab=" ")
x <- seq(min(colon$V48)-1, max(colon$V48)+1, length = 100)
y <- dnorm(x, mean(colon$V48), sd(colon$V48))
lines(x, y, col = "red", lwd = 3)

A simple vista, se podría decir que las distribuciones aparentan ser normales. No obstante, ya conocemos que hay unos picos extremádamente altos, esto queriendo decir que hay unos valores concentrados. Hay una posibilidad de que no exista una distribución normal del todo, aunque lo parezca. Esto lo confirmaremos más adelante con las pruebas analíticas.

Ahora, veámos otra gráfica.

# Visualización con QQ-plots

par(mfrow=c(2, 2))

# Q-Q plot para V42
qqnorm(colon$V42, main = "Q-Q Plot V42", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V42, col = "pink", lwd = 2)

# Q-Q plot para V43
qqnorm(colon$V43, main = "Q-Q Plot V43", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V43, col = "pink", lwd = 2)

# Q-Q plot para V44
qqnorm(colon$V44, main = "Q-Q Plot V44", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V44, col = "pink", lwd = 2)

# Q-Q plot para V45
qqnorm(colon$V45, main = "Q-Q Plot V45", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V45, col = "pink", lwd = 2)

# Q-Q plot para V46
qqnorm(colon$V46, main = "Q-Q Plot V46", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V46, col = "pink", lwd = 2)

# Q-Q plot para V47
qqnorm(colon$V47, main = "Q-Q Plot V47", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V47, col = "pink", lwd = 2)

# Q-Q plot para V48
qqnorm(colon$V48, main = "Q-Q Plot V48", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(colon$V48, col = "pink", lwd = 2)

Aquí se puede apreciar con mayor facilidad la “concentración de valores en un area en específico”. Se puede apreciar como en la mayoría de variables hay una gran concentración de valores en el centro de la gráfica, aparentando tener relaciones fuerte entre los valores. Sin embargo, se puede ver como esa línea se “descarrila” al principio y al final de todas las variables. Por eso en el histograma se ve que los valores van subiendo pero derrepende salta hacia pico extremo.

Visualmente, podemos concluir que las variables no son normales.

Ahora, hagamos las pruebas analíticas.

# Pruebas de normalidad univariadas

library(MVN)
## Warning: package 'MVN' was built under R version 4.3.3
# 1. Lilliefors

L <- mvn(colon, univariateTest = "Lillie", desc = T)
L$univariateNormality
##                              Test  Variable Statistic   p value Normality
## 1 Lilliefors (Kolmogorov-Smirnov)    V42       0.0818  <0.001      NO    
## 2 Lilliefors (Kolmogorov-Smirnov)    V43       0.0698  <0.001      NO    
## 3 Lilliefors (Kolmogorov-Smirnov)    V44       0.1042  <0.001      NO    
## 4 Lilliefors (Kolmogorov-Smirnov)    V45       0.0785  <0.001      NO    
## 5 Lilliefors (Kolmogorov-Smirnov)    V46       0.0994  <0.001      NO    
## 6 Lilliefors (Kolmogorov-Smirnov)    V47       0.0984  <0.001      NO    
## 7 Lilliefors (Kolmogorov-Smirnov)    V48       0.0940  <0.001      NO

Esta prueba es la alternativa cuando el numero es mayor a 30. El p-value al ser menor de 0.05 (que es lo que asumimos que es alfa), nos dice que la distribución de todas las variables no son normal.

# 2. Anderson Darling

AD <- mvn(colon, univariateTest = "AD", desc = T)
AD$univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling    V42     112.4558  <0.001      NO    
## 2 Anderson-Darling    V43      85.2220  <0.001      NO    
## 3 Anderson-Darling    V44     159.6318  <0.001      NO    
## 4 Anderson-Darling    V45     118.1201  <0.001      NO    
## 5 Anderson-Darling    V46     160.7872  <0.001      NO    
## 6 Anderson-Darling    V47     156.2484  <0.001      NO    
## 7 Anderson-Darling    V48     163.1686  <0.001      NO

Nuevamente, la prueba de Anderson Darling nos muestra todas las variables con p-value de (0.001). Este valor al ser menor que el alfa (0.05), nos vuelve a decir que las variables no son distribuidas de manera normal. Estas son las únicas pruebas que pudimos utilizar, ya que las restantes no lo permite por tener más de 5000 variables.

Ya que vimos que las distribuciones no son normales, vamos a utilizar una serie de transformaciones con bestNormaliz() para realizar un conjunto de transformaciones potenciales.No utilizaremos boxcox ya que tenemos valores negativos y cero.

# Cargar las librerías necesarias
library(bestNormalize)
## Warning: package 'bestNormalize' was built under R version 4.3.3
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:bestNormalize':
## 
##     boxcox
# Generar datos de una distribución gamma
n <- 250
set.seed(100)

# Extraer la variable V42

x42 <- colon$V42

# Mejor transformacion automatica V42

best_trans_42 <- bestNormalize(x42)
## Warning: `progress_estimated()` was deprecated in dplyr 1.0.0.
## ℹ The deprecated feature was likely used in the bestNormalize package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
best_trans_42
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 3.468
##  - Center+scale: 5.7889
##  - Double Reversed Log_b(x+a): 8.3844
##  - Exp(x): 96.0228
##  - Log_b(x+a): 6.2582
##  - orderNorm (ORQ): 1.3058
##  - sqrt(x + a): 5.62
##  - Yeo-Johnson: 5.5699
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1543 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -3.750 -0.265  0.000  0.270  5.200
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V43

# Extraer la variable V43

x43 <- colon$V43


#  Mejor transformacion automatica

best_trans_43 <- bestNormalize(x43)
best_trans_43
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.7769
##  - Center+scale: 3.9679
##  - Double Reversed Log_b(x+a): 5.0557
##  - Exp(x): 84.8446
##  - Log_b(x+a): 6.2838
##  - orderNorm (ORQ): 1.1809
##  - sqrt(x + a): 4.7266
##  - Yeo-Johnson: 3.9421
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1961 unique values 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## -5.22 -0.40  0.00  0.38  5.11
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V44

# Extraer la variable V44

x44 <- colon$V44


# 5. Mejor transformacion automatica

best_trans_44 <- bestNormalize(x44)
best_trans_44
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 3.2495
##  - Center+scale: 6.4924
##  - Double Reversed Log_b(x+a): 6.8927
##  - Exp(x): 97.2692
##  - Log_b(x+a): 12.4753
##  - orderNorm (ORQ): 1.1268
##  - sqrt(x + a): 8.2909
##  - Yeo-Johnson: 6.0373
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1777 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -4.580 -0.300  0.000  0.287  5.640
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V45

# Extraer la variable V45

x45 <- colon$V45

# 5. Mejor transformacion automatica

best_trans_45 <- bestNormalize(x45)
best_trans_45
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.3711
##  - Center+scale: 5.2788
##  - Double Reversed Log_b(x+a): 7.6968
##  - Exp(x): 231.2893
##  - Log_b(x+a): 8.2327
##  - orderNorm (ORQ): 1.1702
##  - sqrt(x + a): 5.7961
##  - Yeo-Johnson: 5.2184
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1862 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -4.555 -0.370  0.000  0.380  6.640
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V46

# Extraer la variable V46

x46 <- colon$V46


# 5. Mejor transformacion automatica

best_trans_46 <- bestNormalize(x46)
best_trans_46
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 3.0137
##  - Center+scale: 6.6046
##  - Double Reversed Log_b(x+a): 10.1515
##  - Exp(x): 367.7109
##  - Log_b(x+a): 7.1056
##  - orderNorm (ORQ): 1.3277
##  - sqrt(x + a): 6.4453
##  - Yeo-Johnson: 6.4776
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1810 unique values 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## -5.37 -0.32  0.00  0.31  7.17
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V47

# Extraer la variable V47

x47 <- colon$V47


# 5. Mejor transformacion automatica

best_trans_47 <- bestNormalize(x47)
best_trans_47
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.8804
##  - Center+scale: 7.0654
##  - Double Reversed Log_b(x+a): 11.0959
##  - Exp(x): 310.1319
##  - Log_b(x+a): 10.3007
##  - orderNorm (ORQ): 1.2284
##  - sqrt(x + a): 7.5852
##  - Yeo-Johnson: 6.9386
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 2082 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -5.215 -0.400  0.000  0.395  8.040
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# V48

# Extraer la variable V48

x48 <- colon$V48


# 5. Mejor transformacion automatica

best_trans_48 <- bestNormalize(x48)
best_trans_48
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 8.6322
##  - Center+scale: 12.4383
##  - Double Reversed Log_b(x+a): 15.1153
##  - Exp(x): 216.26
##  - Log_b(x+a): 13.8479
##  - orderNorm (ORQ): 4.7578
##  - sqrt(x + a): 12.4718
##  - Yeo-Johnson: 12.1358
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 6830 nonmissing obs and ties
##  - 1845 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -5.715 -0.350  0.000  0.340  5.770
# Restablecer la configuración gráfica

par(mfrow=c(2, 2))
# Mejores transformaciones histogramas

par(mfrow = c(2, 2))

hist(best_trans_42$x.t, probability = TRUE, main = "Mejor transformación V42", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_42$x.t) - 1, max(best_trans_42$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_42$x.t), sd(best_trans_42$x.t))
lines(x, y, col = "blue", lwd = 3)


hist(best_trans_43$x.t, probability = TRUE, main = "Mejor transformación V43", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_43$x.t) - 1, max(best_trans_43$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_43$x.t), sd(best_trans_43$x.t))
lines(x, y, col = "blue", lwd = 3)


hist(best_trans_44$x.t, probability = TRUE, main = "Mejor transformación V44", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_44$x.t) - 1, max(best_trans_44$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_44$x.t), sd(best_trans_44$x.t))
lines(x, y, col = "blue", lwd = 3)


hist(best_trans_45$x.t, probability = TRUE, main = "Mejor transformación V45", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_45$x.t) - 1, max(best_trans_45$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_45$x.t), sd(best_trans_45$x.t))
lines(x, y, col = "blue", lwd = 3)

hist(best_trans_46$x.t, probability = TRUE, main = "Mejor transformación V46", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_46$x.t) - 1, max(best_trans_46$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_46$x.t), sd(best_trans_46$x.t))
lines(x, y, col = "blue", lwd = 3)


hist(best_trans_47$x.t, probability = TRUE, main = "Mejor transformación V47", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_47$x.t) - 1, max(best_trans_47$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_47$x.t), sd(best_trans_47$x.t))
lines(x, y, col = "blue", lwd = 3)


hist(best_trans_48$x.t, probability = TRUE, main = "Mejor transformación V48", col = "cyan", xlab = "", ylab = "")
x <- seq(min(best_trans_48$x.t) - 1, max(best_trans_48$x.t) + 1, length = 100)
y <- dnorm(x, mean(best_trans_48$x.t), sd(best_trans_48$x.t))
lines(x, y, col = "blue", lwd = 3)

# Mejores transformaciones QQ-plot

par(mfrow = c(1, 2))


qqnorm(best_trans_42$x.t, main = "Mejor transformación V42", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_42$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_43$x.t, main = "Mejor transformación V43", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_43$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_44$x.t, main = "Mejor transformación V44", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_44$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_45$x.t, main = "Mejor transformación V45", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_45$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_46$x.t, main = "Mejor transformación V46", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_46$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_47$x.t, main = "Mejor transformación V47", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_47$x.t, col = "pink", lwd = 2)

qqnorm(best_trans_48$x.t, main = "Mejor transformación V48", xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
qqline(best_trans_48$x.t, col = "pink", lwd = 2)

Hemos transformado todas las variables. Visual y analíticamente, podemos observar como ahora si se cumple con la normalidad. En su totalidad, la mejor transformación fue la OQR.

Proseguiremos con el anáisis multivariado.

Análisis Multivariado:

# Matriz de Correlación- Buscamos correlacion entre variables 

correlacion <- cor(colon)

round(correlacion,2)
##      V42  V43  V44  V45  V46  V47  V48
## V42 1.00 0.25 0.31 0.28 0.23 0.18 0.17
## V43 0.25 1.00 0.37 0.33 0.30 0.41 0.26
## V44 0.31 0.37 1.00 0.43 0.43 0.34 0.28
## V45 0.28 0.33 0.43 1.00 0.50 0.46 0.33
## V46 0.23 0.30 0.43 0.50 1.00 0.51 0.41
## V47 0.18 0.41 0.34 0.46 0.51 1.00 0.40
## V48 0.17 0.26 0.28 0.33 0.41 0.40 1.00

Podemos observar cómo V42 tiene una relación moderada con V44 (0.31) y V45 (0.28). Por otro lado, V42 tiene una relación bastante débil con V43 (0.25), V46 (0.23), V47 (0.18) y V48 (0.17), en ese orden.

Respecto a V43, vemos una cantidad de relaciones moderadas mayor que en V42, siendo estas relaciones con V47 (0.41), V44 (0.37) y V45 (0.33), en ese respectivo orden. Con V46 (0.30) todavía podemos considerarla como moderada, pero con V48 (0.26) ya observamos una relación débil.

Con V44 vemos relaciones un poco más fuertes, por ejemplo, con V45 (0.43) y V46 (0.43). Vemos relaciones moderadas con V43 (0.37), V47 (0.34) y V42 (0.31). Finalmente, se observa una relación débil con V48 (0.28).

En V45, las relaciones relativamente fuertes son con V46 (0.50), V47 (0.46) y V44 (0.43). Vemos relaciones más moderadas con V43 (0.33) y V48 (0.33). También notamos una relación débil con V42 (0.28).

Respecto a V46, observamos una relación fuerte con V47 (0.51), una correlación moderada con V45 (0.50), V44 (0.43) y V43 (0.30), y una relación más débil con V42 (0.23) y V48 (0.41).

Finalmente, con V47 se observan relaciones fuertes con V46 (0.51) y moderadas con V45 (0.46), V44 (0.34) y V43 (0.41). La relación con V48 (0.40) es también moderada, mientras que con V42 (0.18) es bastante débil.

En V48, vemos correlaciones moderadas con V46 (0.41), V47 (0.40) y V45 (0.33), mientras que las relaciones con V44 (0.28), V43 (0.26) y V42 (0.17) son más débiles, siendo esta última la más débil de todas.

Culminando el analisis de correlación con todas las variables, podemos apreciar como las variables V45, V46 y V47 están más correlacionadas entre sí, lo que podría indicar que tienen patrones de comportamiento más similares. Por otro lado, las variables V42 y V48 muestran correlaciones más bajas con las otras variables, lo que indica que podrían tener patrones de comportamiento más independientes.

También podemos graficar esta matriz de correlación con un Heat map. Es una manera visual de ver estas relaciones entre estas variables.

library(ggplot2)
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.3.3
# Matriz de correlación

cor_matrix <- cor(colon[, c("V42", "V43", "V44", "V45", "V46", "V47", "V48")])

# Transformar 

cor_melt <- melt(cor_matrix)

# Graficar Heatmap

ggplot(cor_melt, aes(Var1, Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1, 1), space = "Lab", 
                       name="Correlación") +
  theme_minimal() +
  labs(title = "Heatmap de la Matriz de Correlación")

Nuevamente podemos ver que la diagonal de la matriz muestra correlaciones perfectas (valor 1, representado en tonalidad roja). Todas las correlaciones son positivas, como lo indican los colores cálidos. No se observan tonalidades azules, lo que nso dice que no hay que no hay correlaciones negativas. Esto implica que, a medida que una variable aumenta, la otra también lo hace. Además, el grado de la relación entre variables puede observarse a través de la intensidad del color: cuanto más intensa sea la tonalidad roja, mayor será la correlación positiva entre las variables.

Otra manera de analizar la relación entre estas variables es con el matriz de covarianza. Lo haremos a continuación.

# Matriz de Covarianza- Buscamos covarianzas entre cada variable

library(MASS)

# Calcular la matriz de covarianza

S <- round(cov(colon),2)
S
##      V42  V43  V44  V45  V46  V47  V48
## V42 0.32 0.11 0.12 0.13 0.10 0.09 0.08
## V43 0.11 0.62 0.20 0.21 0.17 0.31 0.17
## V44 0.12 0.20 0.45 0.23 0.21 0.22 0.15
## V45 0.13 0.21 0.23 0.64 0.29 0.35 0.21
## V46 0.10 0.17 0.21 0.29 0.52 0.35 0.24
## V47 0.09 0.31 0.22 0.35 0.35 0.90 0.31
## V48 0.08 0.17 0.15 0.21 0.24 0.31 0.66

Los valores en la diagonal principal de la matriz representan las varianzas de cada variable. Dicho esto, podemo ver como la variable V47 es la que presenta mayor varianza (0.90), lo que sugiere que es la que tiene la mayor dispersion es sus valores. Cabe destacar que es una variable muy volatil. Por otro lado, V42 tiene la menor varianza (0.32), lo que indica que sus valores están más concentrados alrededor de la media.

En general, podemos observar que la variable V42 tiene relaciones relativamente débiles con las demás. Sus covarianzas más altas son con V45 (0.13) y V44 (0.12), lo que sugiere una relación leve. Por otro lado, V42 muestra covarianzas más bajas con V48 (0.08) y V47 (0.09), lo que indica una relación débil entre estas variables.

En el caso de V43, observamos relaciones un tanto moderadas, como con V47 (0.31) y V45 (0.21). También tiene una covarianza moderada con V44 (0.20), mientras que sus relaciones más débiles son con V42 (0.11) y V48 (0.17).

La variable V44 presenta relaciones moderadas con V45 (0.23), V46 (0.21), y V43 (0.20), mientras que su relación con V42 (0.12) y V48 (0.15) es bastante baja. Esto sugiere que V44 comparte cierta similitud con las variables adyacentes, pero no una fuerte correlación.

Con V45, las covarianzas son más destacadas, especialmente con V47 (0.35) y V46 (0.29), lo que indica que estas variables pueden estar más relacionadas entre sí. Las relaciones más débiles de V45 son con V42 (0.13) y V48 (0.21).

La variable V46 también tiene una relación fuerte con V47 (0.35), mientras que las relaciones con V48 (0.24) y V42 (0.10) son más débiles. V46 parece estar más correlacionada con otras variables intermedias como V45 (0.29) y V44 (0.21).

En el caso de V47, se destaca su relación fuerte con V45 (0.35) y V46 (0.35), mientras que muestra correlaciones moderadas con V43 (0.31) y V44 (0.22). Las relaciones más débiles se observan con V42 (0.09) y V48 (0.31).

Finalmente, la variable V48 muestra relaciones moderadas con V47 (0.31), V46 (0.24), y V45 (0.21), mientras que sus relaciones con V44 (0.15), V43 (0.17), y V42 (0.08) son débiles.

Al culminar el análisis de covarianza entre todas las variables, se aprecia que V45, V46, y V47 están más correlacionadas entre sí, lo que podría indicar patrones de comportamiento similares. Por otro lado, V42 y V48 muestran covarianzas más bajas con el resto de las variables, sugiriendo que pueden ser más independientes en su comportamiento.

Ahora haremos una evaluación de la normal multivariada de esta distribución que analizamos. Solo utilizaremos Mardia y Henze-Zirkler debido a que las demás pruebas requiere que las variables sean menos de 5000. Ambas de estas pruebas analíticas se pueden utilizar para evaluar la normalidad multivariada de un conjunto de datos.

library(MVN)

# Evaluación de la normal multivariada

# 1. Mardia 

Mardia <- mvn(colon, mvnTest = "mardia")
Mardia$multivariateNormality
##              Test        Statistic p value Result
## 1 Mardia Skewness 11323.0870922812       0     NO
## 2 Mardia Kurtosis 463.704886476806       0     NO
## 3             MVN             <NA>    <NA>     NO

La prueba analítica de Mardia nos muestra la asimetría y la curtosis en un contexto multivariado. En este caso y con estos datos, podemos ver como Mardia nos muestra un sesgo positivo muy elevado (11323), lo que sugiere que la distribución tiene una asimetría hacia la derecha, con valores extremos que tienden a ser positivos. Además, observamos una curtosis extremadamente alta (463.70), lo que indica que la distribución tiene una gran concentración de valores cerca de la media, con colas más gruesas, lo que implica una mayor probabilidad de valores extremos en comparación con una distribución normal. Respecto a la normalidad de la distribución, la prueba analítica de Mardia nos dice que la distribución no es normal (MVN= NO).

Prosigamos con la prueba analítica de Henze-Zirkler.

# 2. Henze-Zirkler

HZ <- mvn(colon, mvnTest = "hz")
HZ$multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 34.51333       0  NO

La prueba Henze-Zirkler (HZ) tambien la podemos utilizar para evaluar la normalidad multivariada de un conjunto de datos. En el análisis observamos un estadístico HZ= 34.51 y un P-value= 0. Con el P-value ya podemos asumir que se rechazaría una posible hipotesis nula. Además, MVN siendo NO nos quiere decir que los datos no cumplen con la normalidad multivariada.

Además de estas pruebas analíticas de normalidad, también utilizaremos pruebas de hipótesis para investigar más a fondo las relaciones y la distribución de las variables en el conjunto de datos. Seguiremos considerando estas herramientas analíticas (como Mardia y Henze-Zirkler) para complementar el análisis de las pruebas de hipótesis.

3. Proponer por lo menos 3 pruebas de hipótesis de las trabajadas en clase.

1. Primera prueba de hipótesis

En nuestra primera prueba de hipótesis, utilizando un nivel de significancia de alfa= 0.05, probaremos si es posible suponer que las variables son independientes.

H0: Variables son independientes, con varianzas iguales. Ha: Las variables no son independientes.

# Matriz

head(colon)
##           V42    V43           V44         V45    V46    V47        V48
## 1 -0.32003900 -0.620 -4.900000e-01  0.07001953 -0.120 -0.290 -0.8100195
## 2  0.08996101  0.080  4.200000e-01 -0.82998050  0.000  0.030  0.0000000
## 3 -0.29003900  0.140 -3.400000e-01 -0.59998050 -0.010 -0.310  0.2199805
## 4 -1.03003900  0.740  7.000000e-02 -0.90998050  0.130  1.500  0.7399805
## 5  0.09496101  0.205 -2.050000e-01  0.24501950  0.555  0.005  0.1149805
## 6  0.05996101  0.000 -1.387779e-17 -0.43998050 -0.550 -0.540  0.1199805
# Calcular la matriz de covarianza

S <- round(cov(colon),2)
S
##      V42  V43  V44  V45  V46  V47  V48
## V42 0.32 0.11 0.12 0.13 0.10 0.09 0.08
## V43 0.11 0.62 0.20 0.21 0.17 0.31 0.17
## V44 0.12 0.20 0.45 0.23 0.21 0.22 0.15
## V45 0.13 0.21 0.23 0.64 0.29 0.35 0.21
## V46 0.10 0.17 0.21 0.29 0.52 0.35 0.24
## V47 0.09 0.31 0.22 0.35 0.35 0.90 0.31
## V48 0.08 0.17 0.15 0.21 0.24 0.31 0.66
# Tamaño de la muestra

n <- nrow(colon)
p <- ncol(colon)
# Matriz de covarianza fija

Sigma_0 <- matrix(c(0.32,0,0,0,0,0,0,0,0.62,0,0,0,0,0,0,0,0.45,0,0,0,0,0,0,0,0.64,0,0,0,0,0,0,0,0.52,0,0,0,0,0,0,0,0.90,0,0,0,0,0,0,0,0.66),ncol=7)
Sigma_0
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.32 0.00 0.00 0.00 0.00  0.0 0.00
## [2,] 0.00 0.62 0.00 0.00 0.00  0.0 0.00
## [3,] 0.00 0.00 0.45 0.00 0.00  0.0 0.00
## [4,] 0.00 0.00 0.00 0.64 0.00  0.0 0.00
## [5,] 0.00 0.00 0.00 0.00 0.52  0.0 0.00
## [6,] 0.00 0.00 0.00 0.00 0.00  0.9 0.00
## [7,] 0.00 0.00 0.00 0.00 0.00  0.0 0.66
# Estadístico de prueba 

Propios <- eigen(S%*%solve(Sigma_0))
val_prop <- Propios$values
val_prop
## [1] 3.1088230 0.9316472 0.7488588 0.6932365 0.5934416 0.4936209 0.4303720
lambda_0 <- -(n-1)*sum(log(val_prop))
lambda_0
## [1] 11357
a <- 0.05 
gl <- (1/2)*p*(p-1)
lambda_c <- qchisq(1-a,gl)
lambda_c
## [1] 32.67057

Con esta información, encontramos que el valor observado (11357) está por encima del valor crítico (32.67057), lo que indica que hay suficiente evidencia para rechazar la hipótesis nula. Como lambda0 > lambdaC, rechazamos la hipótesis nula. En otras palabras, se concluye que las variables no son independientes, lo que sugiere una relación de dependencia entre ellas.

2. Segunda prueba de hipótesis

En esta prueba de hipótesis analizaremos si es posible suponer que las variables son independientes con igual matrix de covarianza. Esto es lo que aprendimos en la clase que se le conoce como Esfericidad.

H0:Las variables son independientes con igual matrix de covarianza. Ha: Las variables no son indipendiente con igua matriz de covarianza.

# Matriz

head(colon)
##           V42    V43           V44         V45    V46    V47        V48
## 1 -0.32003900 -0.620 -4.900000e-01  0.07001953 -0.120 -0.290 -0.8100195
## 2  0.08996101  0.080  4.200000e-01 -0.82998050  0.000  0.030  0.0000000
## 3 -0.29003900  0.140 -3.400000e-01 -0.59998050 -0.010 -0.310  0.2199805
## 4 -1.03003900  0.740  7.000000e-02 -0.90998050  0.130  1.500  0.7399805
## 5  0.09496101  0.205 -2.050000e-01  0.24501950  0.555  0.005  0.1149805
## 6  0.05996101  0.000 -1.387779e-17 -0.43998050 -0.550 -0.540  0.1199805
# Calcular la matriz de covarianza

S <- round(cov(colon),2)
S
##      V42  V43  V44  V45  V46  V47  V48
## V42 0.32 0.11 0.12 0.13 0.10 0.09 0.08
## V43 0.11 0.62 0.20 0.21 0.17 0.31 0.17
## V44 0.12 0.20 0.45 0.23 0.21 0.22 0.15
## V45 0.13 0.21 0.23 0.64 0.29 0.35 0.21
## V46 0.10 0.17 0.21 0.29 0.52 0.35 0.24
## V47 0.09 0.31 0.22 0.35 0.35 0.90 0.31
## V48 0.08 0.17 0.15 0.21 0.24 0.31 0.66
# Tamaño de la muestra

n <- nrow(colon)
p <- ncol(colon)
# Calculamos el estadístico de prueba

sigma_est  <- sum(diag(S))/p
lambda_0   <- (n-1)*(log(sigma_est)-log(det(S)))
lambda_0 
## [1] 35289.78
# Calculamos el valor crítico

a          <-  0.05
gl         <- (1/2)*(p+2)*(p-1)
lambda_c   <- qchisq(1-a,gl)
lambda_c
## [1] 40.11327

Nuevamente, el estadístico de prueba (lambda0= 35289.78)es mayor que el valor crítico (lamdaC= 40.11327), por lo cual se rechaza la hipótesis nula. Es decir, no hay evidenia suficiente para suponer que las variables son independientes con igual matriz de covarianza.

3. Tercera prueba de hipótesis

Ahora haremos una prueba de hipótesis sobre covarianzas. Tenemos una covarianza estimada. Queremos ver si la covarianza de nuestros datos es igual a esta covarianza estimada.

H0:Σ=Σ0 Ha:Σ≠Σ0

# Tamaño de la muestra

n <- nrow(colon)
p <- ncol(colon)

# Matriz de covarianza fija

Cov0 <- matrix(c(
  0.40, 0.08, 0.10, 0.15, 0.09, 0.07, 0.06,
  0.08, 0.55, 0.22, 0.20, 0.14, 0.30, 0.12,
  0.10, 0.22, 0.48, 0.25, 0.23, 0.19, 0.14,
  0.15, 0.20, 0.25, 0.60, 0.32, 0.34, 0.26,
  0.09, 0.14, 0.23, 0.32, 0.50, 0.28, 0.22,
  0.07, 0.30, 0.19, 0.34, 0.28, 0.85, 0.33,
  0.06, 0.12, 0.14, 0.26, 0.22, 0.33, 0.72
), nrow = 7, byrow = TRUE)
Cov0
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.40 0.08 0.10 0.15 0.09 0.07 0.06
## [2,] 0.08 0.55 0.22 0.20 0.14 0.30 0.12
## [3,] 0.10 0.22 0.48 0.25 0.23 0.19 0.14
## [4,] 0.15 0.20 0.25 0.60 0.32 0.34 0.26
## [5,] 0.09 0.14 0.23 0.32 0.50 0.28 0.22
## [6,] 0.07 0.30 0.19 0.34 0.28 0.85 0.33
## [7,] 0.06 0.12 0.14 0.26 0.22 0.33 0.72
# Matriz de covarianza muestral

S <- round(cov(colon),2)
S
##      V42  V43  V44  V45  V46  V47  V48
## V42 0.32 0.11 0.12 0.13 0.10 0.09 0.08
## V43 0.11 0.62 0.20 0.21 0.17 0.31 0.17
## V44 0.12 0.20 0.45 0.23 0.21 0.22 0.15
## V45 0.13 0.21 0.23 0.64 0.29 0.35 0.21
## V46 0.10 0.17 0.21 0.29 0.52 0.35 0.24
## V47 0.09 0.31 0.22 0.35 0.35 0.90 0.31
## V48 0.08 0.17 0.15 0.21 0.24 0.31 0.66
# Calculamos el estadístico de prueba

Propios <- eigen(S%*%solve(Cov0))
val_prop <- Propios$values
lambda_0 <- (n-1)*(sum(val_prop-log(val_prop))-p)
lambda_0
## [1] 1206.19
# Calculamos el valor crítico

a <- 0.05
gl <- (1/2)*p*(p+1)
lambda_c <- qchisq(1-a, gl)
lambda_c
## [1] 41.33714

Como lambda_0 resulta ser mucho mayor a lambda_c, se rechaza la hipotesis nula. Se puede concluir con un 95% de confianza que la matriz de covarianza poblacional no es similar a la matriz de covarianza propuesta en H0.

Conclusión

Luego de realizar varias pruebas, como las gráficas, análisis univariado y multivado, pruebas analíticas y pruebas de hipotesis, nos encontramos con ciertos patrones o tendencias que explican las relaciones entre las variables mostrada en la báse de datos NCI60 del cáncer de colon.