Introduccion

En este anĆ”lisis, trabajamos con datos de expresiĆ³n gĆ©nica de las lĆ­neas celulares del conjunto NCI-60 del paquete ISLR en R, que incluye informaciĆ³n de diversas lĆ­neas celulares de cĆ”ncer humano. Con un enfoque en las lĆ­neas celulares de ovario para llevar a cabo un estudio detallado. Para evaluar la normalidad de estas variables, se realizarĆ”n pruebas de normalidad y se explorarĆ”n transformaciones que podrĆ­an ayudar a alcanzar la normalidad, lo que facilitarĆ­a inferencias estadĆ­sticas mĆ”s sĆ³lidas en anĆ”lisis posteriores.

# Identificar lĆ­neas celulares de ovario
DATA <- NCI60$data
label <- NCI60$labs
ovarian_cells <- which(NCI60$labs == "OVARIAN")

# Extraer los datos de expresiĆ³n gĆ©nica para estas lĆ­neas celulares
ovarian_data <-NCI60$data[ovarian_cells, ]
ovarianma <- t(ovarian_data)
X22<-ovarianma[,1] 
X25<-ovarianma[,2] 
X26<-ovarianma[,3] 
X27<-ovarianma[,4] 
X28<-ovarianma[,5] 
X29<-ovarianma[,6]

Medidas descriptivas

Medias    <- apply(ovarianma,2,mean)
Des_st    <- apply(ovarianma,2,sd)
Asimetria <- apply(ovarianma,2,skewness)
Curtosis  <- apply(ovarianma,2,kurtosis)

Descriptivas <-rbind(Medias,Des_st,Asimetria,Curtosis)
round(Descriptivas,2)    
##            V22  V25  V26   V27  V28  V29
## Medias    0.02 0.06 0.08  0.04 0.03 0.05
## Des_st    0.77 0.81 0.77  0.66 0.72 0.64
## Asimetria 0.63 1.19 1.60  1.77 0.99 0.90
## Curtosis  5.50 8.90 8.22 15.13 7.39 5.37

Estos datos, nos demuestran que las variables tienen medias bastante similares, con desviaciones estandar bajas, sin embargo estas variables contienen niveles de Curtosis bastante altos, lo que nos dice que hay una prescencia notable de datos atĆ­picos.

Visualizacion de datos

# Histograma
par(mfrow=c(2,3))  
for(i in 1:ncol(ovarianma)) {
  hist(ovarianma[,i], main=paste("Histograma de V", i+21, sep=""), 
       xlab=paste("V", i+21, sep=""), col="lightblue", border="black")}

# QQplot
par(mfrow=c(2,3))
for(i in 1:ncol(ovarianma)) {
    qqnorm(ovarianma[,i], main=paste("Q-Q plot de V", i+21, sep=""))
    qqline(ovarianma[,i], col="red")}

# Boxplot
boxplot(ovarianma)

Matriz de correlaciĆ³n

ggpairs(ovarianma)

SegĆŗn estos visuales, estos a simple vista, pareciera que estos datos no siguen una distribuciĆ³n normal, por lo que, para comprobar esto, se deberĆ­a hacer una prueba de hipĆ³tesis que ponga esto a prueba.

1. Establecer nuestra Prueba de HipĆ³tesis:

\(H_0\):Los datos son normales.

\(H_a\):Los datos no son normales.

2. Establecer nuestro nivel de confianza. \(Ī±= 0.05\)

Pruebas de normalidad univariada

# Cramer-Von Mises
CVM <- mvn(ovarianma, univariateTest = "CVM",desc=T)
## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately

## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately

## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately

## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately

## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately

## Warning in FUN(newX[, i], ...): p-value is smaller than 7.37e-10, cannot be
## computed more accurately
CVM$univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Cramer-von Mises    V22      12.8095  <0.001      NO    
## 2 Cramer-von Mises    V25      31.1092  <0.001      NO    
## 3 Cramer-von Mises    V26      19.9689  <0.001      NO    
## 4 Cramer-von Mises    V27      37.5216  <0.001      NO    
## 5 Cramer-von Mises    V28      19.7132  <0.001      NO    
## 6 Cramer-von Mises    V29      13.8992  <0.001      NO
# Mardia
Mardia <- mvn(ovarianma, mvnTest = "mardia")
Mardia$multivariateNormality
##              Test        Statistic p value Result
## 1 Mardia Skewness 14256.3146996139       0     NO
## 2 Mardia Kurtosis 419.076947796215       0     NO
## 3             MVN             <NA>    <NA>     NO

Luego de observar la prueba de normalidad de estos datos, podemos concluir, a base de que los p-valores de esta prueba se acercan a 0, que fallamos en rechazar nuestra hipĆ³tesis nula, lo que nos dirĆ­a que esta base de datos sobre las celulas de ovarios no son normales a nivel univariado ni multivariado.

Al estos datos no ser normales, deberĆ­a aplicarse una transformaciĆ³n para ver si se puede alcanzar dicha normalidad, en esto utilizamos la formula de bestNormalize:

# TransformaciĆ³n de variables

bt22 <- bestNormalize(X22)
## 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.
bt22
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.8481
##  - Center+scale: 3.8182
##  - Double Reversed Log_b(x+a): 5.9288
##  - Exp(x): 166.7043
##  - Log_b(x+a): 4.8977
##  - orderNorm (ORQ): 1.1943
##  - sqrt(x + a): 3.8143
##  - Yeo-Johnson: 3.7631
## 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
##  - 1947 unique values 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## -4.34 -0.38  0.00  0.41  6.45
bt25 <- bestNormalize(X25)
bt25
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 5.2556
##  - Center+scale: 9.1464
##  - Double Reversed Log_b(x+a): 12.6872
##  - Exp(x): 311.7133
##  - Log_b(x+a): 8.8336
##  - orderNorm (ORQ): 1.9355
##  - sqrt(x + a): 8.6894
##  - Yeo-Johnson: 8.2715
## 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
##  - 1914 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -5.340 -0.330  0.000  0.360  7.915
bt26 <- bestNormalize(X26)
bt26
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.1543
##  - Center+scale: 4.6084
##  - Double Reversed Log_b(x+a): 9.2499
##  - Exp(x): 250.444
##  - Log_b(x+a): 3.4312
##  - orderNorm (ORQ): 1.1093
##  - sqrt(x + a): 3.5213
##  - Yeo-Johnson: 3.0447
## 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
##  - 1882 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -4.230 -0.360  0.000  0.405  6.140
bt27 <- bestNormalize(X27)
bt27
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 3.5264
##  - Center+scale: 8.0377
##  - Double Reversed Log_b(x+a): 12.7233
##  - Exp(x): 356.3692
##  - Log_b(x+a): 10.3314
##  - orderNorm (ORQ): 1.175
##  - sqrt(x + a): 7.5511
##  - Yeo-Johnson: 7.0446
## 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
##  - 1726 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -4.053 -0.260  0.000  0.270  7.180
bt28 <- bestNormalize(X28)
bt28
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 6.4132
##  - Center+scale: 8.7189
##  - Double Reversed Log_b(x+a): 11.968
##  - Exp(x): 171.9059
##  - Log_b(x+a): 9.4507
##  - orderNorm (ORQ): 3.749
##  - sqrt(x + a): 8.5315
##  - Yeo-Johnson: 8.4132
## 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
##  - 1820 unique values 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## -4.78 -0.34  0.00  0.35  5.69
bt29 <- bestNormalize(X29)
bt29
## Best Normalizing transformation with 6830 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 2.1748
##  - Center+scale: 3.7987
##  - Double Reversed Log_b(x+a): 6.71
##  - Exp(x): 84.9539
##  - Log_b(x+a): 4.1692
##  - orderNorm (ORQ): 1.2041
##  - sqrt(x + a): 3.3093
##  - Yeo-Johnson: 2.9483
## 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
##  - 1761 unique values 
##  - Original quantiles:
##     0%    25%    50%    75%   100% 
## -3.390 -0.315  0.000  0.355  4.790

VisualizaciĆ³n de datos multivariados

datatrans <- matrix(c(bt22$x.t,bt25$x.t, bt26$x.t, bt27$x.t,bt28$x.t, bt29$x.t),ncol=6)

# Histograma
par(mfrow=c(2,3))  
for(i in 1:ncol(datatrans)) {
  hist(datatrans[,i], main=paste("Histograma Transformado de V", i+21, sep=""), 
       xlab=paste("V", i+21, sep=""), col="lightblue", border="black")}

# QQplot
par(mfrow=c(2,3))
for(i in 1:ncol(datatrans)) {
    qqnorm(datatrans[,i], main=paste("Q-Q plot de V", i+21, sep=""))
    qqline(datatrans[,i], col="red")}

# Boxplot
boxplot(datatrans)

Encontramos que la mejor transformaciĆ³n para estas variables resultĆ³ ser la de OrderNorm, y al ver los graficos estas variables empezaron a tener formas normales, ahora ponemos a prueba si estas variables realmente son normales con una prueba de hipĆ³tesis:

1. Establecer nuestra Prueba de HipĆ³tesis:

\(H_0\):Los datos son normales.

\(H_a\):Los datos no son normales.

2. Establecer nuestro nivel de confianza.

\(Ī±= 0.05\)

CVM22 <- cvm.test(bt22$x.t)     # Esta variable logrĆ³ la normalidad
CVM22
## 
##  Cramer-von Mises normality test
## 
## data:  bt22$x.t
## W = 0.0093468, p-value = 0.9996
CVM25 <- cvm.test(bt25$x.t)     # Esta variable logrĆ³ la normalidad
CVM25
## 
##  Cramer-von Mises normality test
## 
## data:  bt25$x.t
## W = 0.09933, p-value = 0.1149
CVM26 <- cvm.test(bt26$x.t)     # Esta variable logrĆ³ la normalidad
CVM26
## 
##  Cramer-von Mises normality test
## 
## data:  bt26$x.t
## W = 0.014544, p-value = 0.9951
CVM27 <- cvm.test(bt27$x.t)     # Esta variable logrĆ³ la normalidad
CVM27
## 
##  Cramer-von Mises normality test
## 
## data:  bt27$x.t
## W = 0.023887, p-value = 0.9246
CVM28 <- cvm.test(bt28$x.t)     # Esta variable no logrĆ³ la normalidad
CVM28
## 
##  Cramer-von Mises normality test
## 
## data:  bt28$x.t
## W = 0.31294, p-value = 0.0002367
CVM29 <- cvm.test(bt29$x.t)     # Esta variable logrĆ³ la normalidad
CVM29 
## 
##  Cramer-von Mises normality test
## 
## data:  bt29$x.t
## W = 0.010001, p-value = 0.9994

Al observar las pruebas a nivel univariado, se puede ver que 5 de las 6 variables de la celula de ovario si lograrĆ³n a obtener una distribuciĆ³n Normal, ahora se comprobarĆ” si las variables que alcanzaron esta normalidad, pudieran crear una matriz con Normalidad, para ver si pudieramos hacer inferencias sobre esta:

# Matriz tranformada con variables con distribucion normal.

datatri <- matrix(c(bt22$x.t,bt25$x.t, bt26$x.t, bt27$x.t, bt29$x.t),ncol=5)

1. Establecer nuestra Prueba de HipĆ³tesis:

\(H_0\):Los datos son normales.

\(H_a\):Los datos no son normales.

2. Establecer nuestro nivel de confianza.

\(Ī±= 0.05\)

# Mardia
Mardiatra <- mvn(datatri, mvnTest = "mardia")
Mardiatra$multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 502.560217991754 3.71330755841732e-84     NO
## 2 Mardia Kurtosis 35.0122996431687                    0     NO
## 3             MVN             <NA>                 <NA>     NO
# Henze-Zirkler
HZ <- mvn(datatri, mvnTest = "hz")
HZ$multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 2.834468       0  NO

Al observar las pruebas a nivel multivariado, se puede ver que la matriz creada con las variables transformadas no pudieron lograr una distribuciĆ³n Normal, debido a que el pvalor obtenido en estas es uno cercano a 0, por lo que se puede decir que no existe prueba significativa para rechazar H0 por lo que no se le pueden aplicar algunas propiedades para poder hacer inferencias sobre esta matriz de datos.

NOTA: Como la matriz no resulto normal, no se pudieran hacer supuestos sobre su matriz de covarianza, pero para propĆ³sitos del examen asumimos que los datos transformados sĆ­ lo fueron, a continuaciĆ³n se presentaran algunas pruebas las cuĆ”l se querĆ­an trabajar:

Supongamos que ahora queremos comprobar si las variables que conforman la matriz de datos son independientes,

1. Establecer nuestra Prueba de HipĆ³tesis:

\(H_0:Ī£=Ī£0\)

Las variables son independientes.

\(H_0:Ī£ā‰ Ī£0\)

Las variables no son independientes.

2. Establecer nuestro nivel de confianza.

\(Ī±= 0.05\)

# InformaciĆ³n para el ejercicio
n        <- 6830
p        <- 6
# Matriz de covarianza muestral
S <-round(cov(datatrans),2)
S
##       [,1]  [,2] [,3]  [,4] [,5] [,6]
## [1,]  1.00 -0.03 0.00 -0.03 0.15 0.01
## [2,] -0.03  1.00 0.28  0.11 0.07 0.01
## [3,]  0.00  0.28 1.00  0.16 0.14 0.04
## [4,] -0.03  0.11 0.16  1.00 0.12 0.08
## [5,]  0.15  0.07 0.14  0.12 1.00 0.15
## [6,]  0.01  0.01 0.04  0.08 0.15 1.00
diag(S)
## [1] 1 1 1 1 1 1
# Matriz de covarianza fija
Sigma0     <-matrix(c(1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1),ncol=6)
Sigma0
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    0    0    0    0
## [2,]    0    1    0    0    0    0
## [3,]    0    0    1    0    0    0
## [4,]    0    0    0    1    0    0
## [5,]    0    0    0    0    1    0
## [6,]    0    0    0    0    0    1
# Calculamos el valor crĆ­tico
a          <-  0.05
gl         <- (1/2)*p*(p-1)
lambda_c   <- qchisq(1-a,gl)
lambda_c
## [1] 24.99579
# Calculamos el estadĆ­stico de prueba
Propios    <- eigen(S%*%solve(Sigma0))
val_prop   <- Propios$values
val_prop 
## [1] 1.4916603 1.1527888 1.0042922 0.8761990 0.7675256 0.7075342
lambda_0   <- -(n-1)*sum(log(val_prop))
lambda_0 
## [1] 1340.929

Como \(Ī»_0ā‰„Ī»_c\), entonces se rechaza \(H_0\). Se puede concluir con un 95% de confianza que las variables transformadas de la matriz de covarianza de la celula de ovario no se pueden suponer que son independientes.

Ahora se quisiera estudiar si esta matriz de covarianza tiene forma esferica, que son independientes y tienen varianza en comĆŗn.

1. Establecer nuestra Prueba de HipĆ³tesis:

\(H_0:Ī£=Ļƒ^2I\)

Las variables son independientes y tienen la misma varianza. \(H_0:Ī£ā‰ Ļƒ^2I\)

Las variables no son independientes y no tienen la misma varianza.

2. Establecer nuestro nivel de confianza.

\(Ī±= 0.05\)

# InformaciĆ³n del problema
n        
## [1] 6830
p 
## [1] 6
# Matriz de covarianza muestral
S 
##       [,1]  [,2] [,3]  [,4] [,5] [,6]
## [1,]  1.00 -0.03 0.00 -0.03 0.15 0.01
## [2,] -0.03  1.00 0.28  0.11 0.07 0.01
## [3,]  0.00  0.28 1.00  0.16 0.14 0.04
## [4,] -0.03  0.11 0.16  1.00 0.12 0.08
## [5,]  0.15  0.07 0.14  0.12 1.00 0.15
## [6,]  0.01  0.01 0.04  0.08 0.15 1.00
# Calculamos el valor crĆ­tico
a          <-  0.05
gl         <- (1/2)*(p+2)*(p-1)
lambda_c   <- qchisq(1-a,gl)
lambda_c
## [1] 31.41043
# 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] 1340.929

Como \(Ī»_0ā‰„Ī»_c\), entonces se rechaza \(H_0\). Se puede concluir con un 95% de confianza que las variables transformadas de la matriz de covarianza de la celula de ovario no se pueden suponer que son independientes ni tienen varianza en comun.

Para conocer mas sobre la relaciĆ³n entre las celulas de ovarios quisieramos estudiar la densidad que tienen las variables transformadas X25 Y X26, debido a que estas tienen la mayor correlaciĆ³n, serĆ­a interesante observar como estas se vieran en un plano.

ma2526 <-matrix(c(bt25$x.t,bt26$x.t),ncol=2) 
mu <- colMeans(ma2526)     
Sigma2 <- cov(ma2526)           

densidad <- dmvnorm(ma2526, mean = mu, sigma = Sigma2)
graph <- plot_ly(x = ~ma2526[,1], y = ~ma2526[,2], z = ~densidad, type = "scatter3d", color=densidad, xlab="",ylab="", zlab="")

graphlab <- graph %>% layout(scene = list(
  xaxis = list(title = "X25 Transformada"),
  yaxis = list(title = "X26 Transformada"),
  zaxis = list(title = "Densidad")
))
graphlab
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## Warning: 'scatter3d' objects don't have these attributes: 'xlab', 'ylab', 'zlab'
## Valid attributes include:
## 'connectgaps', 'customdata', 'customdatasrc', 'error_x', 'error_y', 'error_z', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'projection', 'scene', 'showlegend', 'stream', 'surfaceaxis', 'surfacecolor', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'visible', 'x', 'xcalendar', 'xhoverformat', 'xsrc', 'y', 'ycalendar', 'yhoverformat', 'ysrc', 'z', 'zcalendar', 'zhoverformat', 'zsrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

Otra forma para observar esto serĆ­a un Diagrama de dispersiĆ³n

scatterplot3d(x = ma2526[,1], y = ma2526[,2], z = densidad, pch = 19, cex.lab = 1, highlight.3d = TRUE, xlab = "X25 Transformada ", ylab = "X26 Transformada ", zlab = "Densidad")

Al observar las graficas de densidad entre la variable X25 y X26 se puede notar que a simple vista, estos parecieran tener una distribuciĆ³n normal.