En esta investigación, nos enfocamos en analizar las líneas celulares derivadas de cáncer de mama utilizando la base de datos NCI-60, proporcionada por el Instituto Nacional del Cáncer de EE. UU. El cáncer de mama es una de las formas más prevalentes de cáncer y su estudio es crucial para entender su desarrollo, progresión y posibles tratamientos. A través de técnicas estadísticas, exploraremos datos relacionados con la expresión génica, sensibilidad a fármacos y otras características relevantes de estas líneas celulares.
Nuestro objetivo es aplicar métodos univariados y multivariados para identificar patrones y asociaciones dentro de las variables relacionadas con el cáncer de mama. Además, realizaremos pruebas de hipótesis para validar algunos de los hallazgos y generar conocimiento que contribuya a la investigación oncológica, particularmente en la personalización de terapias y el desarrollo de nuevos fármacos.
Los datos con los cuales estaremos trabajando son los siguientes:
library(ISLR)
NCI60 <- NCI60
data1 <- data.frame(t(NCI60$data))
names <- NCI60$labs
a <- which(names=="BREAST")
breast <- data1[,a]
A traves de estas graficas podremos observar el comportamiento de las lineas celulares humanas derivadas del cancer de mama (breast) con las que estaremos trabajando:
library(GGally)
ggpairs(breast)
En esta grafica podemos observar que algunas variables tienen correlaciones significativas, como V57 y V58 (0.832), Mientras que otras, como v5 Y v57, no estan correlacionadas (0.005) lo que indican independencia.
par(mfrow=c(1,2))
hist(breast$V5,probability=T,xlab="V5",ylab=" ",main=" ")
x <- seq(min(breast$V5)-1, max(breast$V5)+1, length = 100)
y <- dnorm(x, mean(breast$V5),sd(breast$V5))
lines(x, y, col = "red",lwd=3)
hist(breast$V8,probability=T,xlab="V8",ylab=" ",main=" ")
x <- seq(min(breast$V8)-1, max(breast$V8)+1, length = 100)
y <- dnorm(x, mean(breast$V8),sd(breast$V8))
lines(x, y, col = "red",lwd=3)
En esta grafica podemos observar como claramente las variables no son normales, pues su forma no muestra una perfecta campana de Gauss, si no que la parte superior en vez de ser suave se presenta en forma de pico.
Ahora realizaremos pruebas para estudiar la normalidad de estas variables.
Comenzamos con Pruebas Univariadas:
library(MVN)
AD <- mvn(breast[,-1], univariateTest = "AD",desc=T)
AD$univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling V8 94.9941 <0.001 NO
## 2 Anderson-Darling V18 64.9205 <0.001 NO
## 3 Anderson-Darling V50 117.7221 <0.001 NO
## 4 Anderson-Darling V52 87.5734 <0.001 NO
## 5 Anderson-Darling V57 101.7055 <0.001 NO
## 6 Anderson-Darling V58 132.4390 <0.001 NO
AD
## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 35.81118 0 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling V8 94.9941 <0.001 NO
## 2 Anderson-Darling V18 64.9205 <0.001 NO
## 3 Anderson-Darling V50 117.7221 <0.001 NO
## 4 Anderson-Darling V52 87.5734 <0.001 NO
## 5 Anderson-Darling V57 101.7055 <0.001 NO
## 6 Anderson-Darling V58 132.4390 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## V8 6830 0.034898165 0.7719684 0 -5.520000 5.65000 -0.3750000 0.3850000
## V18 6830 0.005414269 0.8192777 0 -5.380000 5.88000 -0.4300000 0.4150073
## V50 6830 -0.045211715 0.8783862 0 -6.939981 5.97998 -0.4300195 0.3899805
## V52 6830 0.014101316 0.8491402 0 -4.839961 6.61000 -0.4200000 0.4200195
## V57 6830 -0.004496054 0.7949380 0 -6.100020 6.50000 -0.3850000 0.3643676
## V58 6830 -0.019257425 0.8310602 0 -6.120039 5.99998 -0.3800195 0.3843481
## Skew Kurtosis
## V8 0.5434773 5.335457
## V18 0.2177262 4.431392
## V50 -0.2088429 5.256755
## V52 0.4517649 5.970943
## V57 0.2166695 6.973168
## V58 -0.3888450 7.296726
Pruebas Multivariadas:
library(MVN)
Mardia <- mvn(breast [,], mvnTest = "mardia")
Mardia$multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 22972.4630960486 0 NO
## 2 Mardia Kurtosis 404.134115759277 0 NO
## 3 MVN <NA> <NA> NO
Mardia
## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 22972.4630960486 0 NO
## 2 Mardia Kurtosis 404.134115759277 0 NO
## 3 MVN <NA> <NA> NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling V5 193.9133 <0.001 NO
## 2 Anderson-Darling V8 94.9941 <0.001 NO
## 3 Anderson-Darling V18 64.9205 <0.001 NO
## 4 Anderson-Darling V50 117.7221 <0.001 NO
## 5 Anderson-Darling V52 87.5734 <0.001 NO
## 6 Anderson-Darling V57 101.7055 <0.001 NO
## 7 Anderson-Darling V58 132.4390 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## V5 6830 0.148587415 0.9929087 0 -4.680000 7.42500 -0.4000000 0.5149878
## V8 6830 0.034898165 0.7719684 0 -5.520000 5.65000 -0.3750000 0.3850000
## V18 6830 0.005414269 0.8192777 0 -5.380000 5.88000 -0.4300000 0.4150073
## V50 6830 -0.045211715 0.8783862 0 -6.939981 5.97998 -0.4300195 0.3899805
## V52 6830 0.014101316 0.8491402 0 -4.839961 6.61000 -0.4200000 0.4200195
## V57 6830 -0.004496054 0.7949380 0 -6.100020 6.50000 -0.3850000 0.3643676
## V58 6830 -0.019257425 0.8310602 0 -6.120039 5.99998 -0.3800195 0.3843481
## Skew Kurtosis
## V5 1.8190590 7.972899
## V8 0.5434773 5.335457
## V18 0.2177262 4.431392
## V50 -0.2088429 5.256755
## V52 0.4517649 5.970943
## V57 0.2166695 6.973168
## V58 -0.3888450 7.296726
DH <- mvn(breast[,], mvnTest = "dh")
DH$multivariateNormality
## Test E df p value MVN
## 1 Doornik-Hansen 15791.54 14 0 NO
DH
## $multivariateNormality
## Test E df p value MVN
## 1 Doornik-Hansen 15791.54 14 0 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling V5 193.9133 <0.001 NO
## 2 Anderson-Darling V8 94.9941 <0.001 NO
## 3 Anderson-Darling V18 64.9205 <0.001 NO
## 4 Anderson-Darling V50 117.7221 <0.001 NO
## 5 Anderson-Darling V52 87.5734 <0.001 NO
## 6 Anderson-Darling V57 101.7055 <0.001 NO
## 7 Anderson-Darling V58 132.4390 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## V5 6830 0.148587415 0.9929087 0 -4.680000 7.42500 -0.4000000 0.5149878
## V8 6830 0.034898165 0.7719684 0 -5.520000 5.65000 -0.3750000 0.3850000
## V18 6830 0.005414269 0.8192777 0 -5.380000 5.88000 -0.4300000 0.4150073
## V50 6830 -0.045211715 0.8783862 0 -6.939981 5.97998 -0.4300195 0.3899805
## V52 6830 0.014101316 0.8491402 0 -4.839961 6.61000 -0.4200000 0.4200195
## V57 6830 -0.004496054 0.7949380 0 -6.100020 6.50000 -0.3850000 0.3643676
## V58 6830 -0.019257425 0.8310602 0 -6.120039 5.99998 -0.3800195 0.3843481
## Skew Kurtosis
## V5 1.8190590 7.972899
## V8 0.5434773 5.335457
## V18 0.2177262 4.431392
## V50 -0.2088429 5.256755
## V52 0.4517649 5.970943
## V57 0.2166695 6.973168
## V58 -0.3888450 7.296726
HZ <- mvn(breast[,], mvnTest = "hz")
HZ$multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 31.39747 0 NO
HZ
## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 31.39747 0 NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling V5 193.9133 <0.001 NO
## 2 Anderson-Darling V8 94.9941 <0.001 NO
## 3 Anderson-Darling V18 64.9205 <0.001 NO
## 4 Anderson-Darling V50 117.7221 <0.001 NO
## 5 Anderson-Darling V52 87.5734 <0.001 NO
## 6 Anderson-Darling V57 101.7055 <0.001 NO
## 7 Anderson-Darling V58 132.4390 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## V5 6830 0.148587415 0.9929087 0 -4.680000 7.42500 -0.4000000 0.5149878
## V8 6830 0.034898165 0.7719684 0 -5.520000 5.65000 -0.3750000 0.3850000
## V18 6830 0.005414269 0.8192777 0 -5.380000 5.88000 -0.4300000 0.4150073
## V50 6830 -0.045211715 0.8783862 0 -6.939981 5.97998 -0.4300195 0.3899805
## V52 6830 0.014101316 0.8491402 0 -4.839961 6.61000 -0.4200000 0.4200195
## V57 6830 -0.004496054 0.7949380 0 -6.100020 6.50000 -0.3850000 0.3643676
## V58 6830 -0.019257425 0.8310602 0 -6.120039 5.99998 -0.3800195 0.3843481
## Skew Kurtosis
## V5 1.8190590 7.972899
## V8 0.5434773 5.335457
## V18 0.2177262 4.431392
## V50 -0.2088429 5.256755
## V52 0.4517649 5.970943
## V57 0.2166695 6.973168
## V58 -0.3888450 7.296726
Segun las pruebas uni y multivariadas podemos confirmar que las variables no son normales. Por lo que procedemos a buscar el mejor metodo de transformacion. Utilizamos la funcion bestnormalize:
library(bestNormalize)
besttrans_V5 <- bestNormalize(breast$V5)
besttrans_V5
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 3.6123
## - Center+scale: 7.7261
## - Double Reversed Log_b(x+a): 13.8852
## - Exp(x): 445.1969
## - Log_b(x+a): 5.1397
## - orderNorm (ORQ): 1.3967
## - sqrt(x + a): 5.827
## - Yeo-Johnson: 4.3258
## 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
## - 2026 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -4.680 -0.400 0.000 0.515 7.425
besttrans_V8 <- bestNormalize(breast$V8)
besttrans_V8
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 2.153
## - Center+scale: 4.4864
## - Double Reversed Log_b(x+a): 6.7506
## - Exp(x): 124.6481
## - Log_b(x+a): 4.8154
## - orderNorm (ORQ): 1.253
## - sqrt(x + a): 4.2064
## - Yeo-Johnson: 4.1007
## 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
## - 1934 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -5.520 -0.375 0.000 0.385 5.650
besttrans_V18 <- bestNormalize(breast$V18)
besttrans_V18
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 4.4867
## - Center+scale: 6.1308
## - Double Reversed Log_b(x+a): 7.5286
## - Exp(x): 137.0376
## - Log_b(x+a): 7.7261
## - orderNorm (ORQ): 2.2616
## - sqrt(x + a): 6.2239
## - Yeo-Johnson: 6.1406
## 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
## - 1965 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -5.380 -0.430 0.000 0.415 5.880
besttrans_V50 <- bestNormalize(breast$V50)
besttrans_V50
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 3.8256
## - Center+scale: 6.8034
## - Double Reversed Log_b(x+a): 7.4842
## - Exp(x): 144.2676
## - Log_b(x+a): 8.6707
## - orderNorm (ORQ): 1.5648
## - sqrt(x + a): 7.3757
## - Yeo-Johnson: 6.5542
## 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
## - 1966 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -6.94 -0.43 0.00 0.39 5.98
besttrans_V52 <- bestNormalize(breast$V52)
besttrans_V52
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.9076
## - Center+scale: 4.1432
## - Double Reversed Log_b(x+a): 6.2762
## - Exp(x): 221.8626
## - Log_b(x+a): 6.065
## - orderNorm (ORQ): 1.2943
## - sqrt(x + a): 4.4237
## - Yeo-Johnson: 4.0827
## 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
## - 2003 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -4.84 -0.42 0.00 0.42 6.61
besttrans_V57 <- bestNormalize(breast$V57)
besttrans_V57
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 3.8586
## - Center+scale: 6.3021
## - Double Reversed Log_b(x+a): 7.9406
## - Exp(x): 213.9697
## - Log_b(x+a): 7.7176
## - orderNorm (ORQ): 1.7364
## - sqrt(x + a): 6.6036
## - Yeo-Johnson: 6.3693
## 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
## - 1941 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -6.100 -0.385 0.000 0.364 6.500
besttrans_V58 <- bestNormalize(breast$V58)
besttrans_V58
## Best Normalizing transformation with 6830 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 2.4397
## - Center+scale: 5.6775
## - Double Reversed Log_b(x+a): 6.7972
## - Exp(x): 148.5484
## - Log_b(x+a): 9.3145
## - orderNorm (ORQ): 1.2028
## - sqrt(x + a): 6.7399
## - Yeo-Johnson: 5.6222
## 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
## - 1895 unique values
## - Original quantiles:
## 0% 25% 50% 75% 100%
## -6.120 -0.380 0.000 0.384 6.000
La prueba bestNormalize nos indica que la mejor tranformacion para todas las variables es orderNorm. Ahora procedemos a transformarlas para convertirlas en normales.
library(MVN)
library(nortest)
V5_trans <- ad.test(besttrans_V5$x.t)
V5_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V5$x.t
## A = 0.13751, p-value = 0.9767
V8_trans <- ad.test(besttrans_V8$x.t)
V8_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V8$x.t
## A = 0.075227, p-value = 0.9992
V18_trans <- ad.test(besttrans_V18$x.t)
V18_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V18$x.t
## A = 0.6488, p-value = 0.09046
V50_trans <- ad.test(besttrans_V50$x.t)
V50_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V50$x.t
## A = 0.24629, p-value = 0.7566
V52_trans <- ad.test(besttrans_V52$x.t)
V52_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V52$x.t
## A = 0.051368, p-value = 0.9999
V57_trans <- ad.test(besttrans_V57$x.t)
V57_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V57$x.t
## A = 0.36322, p-value = 0.441
V58_trans <- ad.test(besttrans_V58$x.t)
V58_trans
##
## Anderson-Darling normality test
##
## data: besttrans_V58$x.t
## A = 0.086574, p-value = 0.9983
Podemos observar que todos los p-value son mayores que alfa y por ende concluimos que ya las variables que componen el cancer de mama (breast) son normales.
Nuestro siguiente paso sera formular y comprobar hipotesis sobre las repectivas variables.
Hipótesis 1: Prueba de hipótesis para un vector de medias con covarianza desconocida
Tenemos 6830 observaciones y 7 variables
n <- 6830
p <-7
mu0 <- as.vector(c(0.15,0.04,0.01,-0.05,0.02,-0.01,-0.02))
mu0
## [1] 0.15 0.04 0.01 -0.05 0.02 -0.01 -0.02
xbar <- colMeans(breast)
xbar
## V5 V8 V18 V50 V52 V57
## 0.148587415 0.034898165 0.005414269 -0.045211715 0.014101316 -0.004496054
## V58
## -0.019257425
S <- cov(breast)
S
## V5 V8 V18 V50 V52 V57
## V5 0.985867643 0.29188627 0.09625408 -0.20456245 -0.11453978 0.003670185
## V8 0.291886275 0.59593528 0.04657967 -0.10891561 -0.05654456 -0.033067190
## V18 0.096254080 0.04657967 0.67121601 -0.07503046 -0.07490351 -0.019812111
## V50 -0.204562452 -0.10891561 -0.07503046 0.77156235 0.35125689 -0.035140388
## V52 -0.114539781 -0.05654456 -0.07490351 0.35125689 0.72103914 -0.046439342
## V57 0.003670185 -0.03306719 -0.01981211 -0.03514039 -0.04643934 0.631926416
## V58 -0.025590365 -0.05403290 -0.02952381 -0.02150555 -0.02665193 0.549435404
## V58
## V5 -0.02559036
## V8 -0.05403290
## V18 -0.02952381
## V50 -0.02150555
## V52 -0.02665193
## V57 0.54943540
## V58 0.69066114
Asumimos un alpha de 0.05 y calculamos valor critico
a <- 0.05
Fc <- qf(1-a,p,n-p)
round(Fc,2)
## [1] 2.01
Calculamos el estadistico de prueba
T0 <- n*t((xbar-mu0)) %*% solve(S) %*% (xbar-mu0)
F0 <- ((n-p)/(p*(n-1)))*T0
round(F0,2)
## [,1]
## [1,] 0.33
Como f0 cae dentro de la zona de no rechazo, no rechazamos H0. Es decir, no existe evidencia suficiente para de decir que mu0 ≠ 0.15,0.04,0.01,-0.05,0.02,-0.01,-0.02.
Hipótesis 2: Prueba de independencia
n <- 6830
p <- 2
Creamos la matriz de covarianza muestral con las dos variables a estudiar
breast_nueva <- breast[,c(6,7)]
S <- cov(breast_nueva)
S
## V57 V58
## V57 0.6319264 0.5494354
## V58 0.5494354 0.6906611
Realizamos la matriz de covarianza fija
Sigma0 <- matrix(c(0.631926416,0,0,0.69066114),ncol=2)
Sigma0
## [,1] [,2]
## [1,] 0.6319264 0.0000000
## [2,] 0.0000000 0.6906611
Asumimos un alpha de 0.05 y procedemos a calcular el valor critico
a <- 0.05
gl <- (1/2)*p*(p-1)
lambda_c <- qchisq(1-a,gl)
lambda_c
## [1] 3.841459
Por ultimo, calculamos el estadistico de prueba
Propios <- eigen(S%*%solve(Sigma0))
val_prop <- Propios$values
val_prop
## [1] 1.8316697 0.1683303
lambda_0 <- -(n-1)*sum(log(val_prop))
lambda_0
## [1] 8034.996
Lambda0 cae dentro de la zona de rechazo, por lo que rechazamos la hipotesis. Confirmamos con un 95% de confianza que las lineas celulares de cancer de mama V57 y V58 son dependientes.
HIPóTESIS 3: Prueba de esferecidad, las Variables V8 y V18 son independientes y tienen la misma varianza.
Presentamos los datos:
n <- 6830
p <- 2
Formamos nuestra matriz de covarianza muestral
breast_nueva <- breast[,c(2,3)]
S <- cov(breast_nueva)
S
## V8 V18
## V8 0.59593528 0.04657967
## V18 0.04657967 0.67121601
Calculamos la matriz de covarianza fija
sigma_est <- sum(diag(S))/p
sigma_est
## [1] 0.6335756
Asumimos un alfa de 0.05 y calculamos valor critico
a <- 0.05
gl <- (1/2)*(p+2)*(p-1)
lambda_c <- qchisq(1-a,gl)
lambda_c
## [1] 5.991465
Por ultimo computamos el estadistico de prueba
lambda_0 <- (n-1)*(log(sigma_est)-log(det(S)))
lambda_0
## [1] 3177.879
La hipótesis nula se rechaza ya que el lambda0 cae dentro de la zona de rechazo. Es decir, con un 95% de confianza afirmamos que las variables V8 y V18 no tienen forma esferica.
En esta investigación, logramos transformar las variables, inicialmente no normales, a distribuciones normales. Posteriormente, realizamos pruebas de independencia, esfericidad y de vectores de medias para obtener resultados relevantes en el análisis de las líneas celulares de cáncer de mama.