Ejercicio A.1

Apartado 1:

Simule o genere aleatoriamente datos de la expresión génica de un grupo de 10000 genes (en filas) para 12 momentos muestrales (en columnas):

Iniciamos la semilla y generamos los datos de genes y momentos muestrales.

set.seed(789)

num_genes <- 10000
mom_muestral_1 <- 6
mom_muestral_2 <- 6

mitad_1 <- matrix(runif(num_genes*mom_muestral_1,-2,2), nrow = num_genes, ncol = mom_muestral_1)
mitad_2 <- matrix(runif(num_genes*mom_muestral_2,-4,4), nrow = num_genes, ncol = mom_muestral_2)

datosaleaEG01 <- cbind(mitad_1,mitad_2)
dim(datosaleaEG01)
## [1] 10000    12
head(datosaleaEG01)
##             [,1]       [,2]        [,3]      [,4]       [,5]       [,6]
## [1,]  0.79957746  1.4943258  0.40687822  1.816425 -1.2727249 -1.1624271
## [2,] -1.62600451  0.5128004  0.87907447  1.077089  1.5792660 -0.9147813
## [3,] -1.95245273 -1.0258133 -0.06373976  1.447735 -0.5931999  0.3099358
## [4,]  0.36642541  1.6125608 -1.51397604  0.895905 -1.8915467 -1.1032403
## [5,] -0.03140225 -1.8225208  0.48047472 -1.562934 -0.3440984 -1.7929110
## [6,] -1.91934570 -0.5603091  0.47769946  1.185907  0.3340634 -0.2731466
##             [,7]       [,8]      [,9]     [,10]      [,11]      [,12]
## [1,] -0.14484139 -2.6215491 -1.864603  3.116052 -2.2350171 -2.2416845
## [2,]  3.59555582 -1.4881947  3.905833  3.002433  2.1280332 -3.7134398
## [3,]  0.05924968  1.8103155  1.259982 -3.680984 -1.2667359 -2.5608207
## [4,] -3.93094174 -0.3781865  1.681134  1.748022  1.7376943 -1.7946701
## [5,]  2.53134453 -1.9637646 -3.571918  2.462766 -0.1786467  1.4077508
## [6,] -0.44784125  1.9372120 -2.828150  3.601530  0.3162603 -0.8776999

La tabla generada contiene 10000 filas correspondientes a los 10000 genes y 12 columnas correspondientes a los 12 momentos muestrales.

Apartado 2.

Asignación de nombres a los genes y a los momentos muestrales.

nombre_filas <- sprintf("gen%.5d",1:num_genes)
rownames(datosaleaEG01) = nombre_filas

nombre_columnas <- c(paste0("notratado0",1:mom_muestral_1),paste0("tratado0",1:mom_muestral_2))
colnames(datosaleaEG01) = nombre_columnas

head(datosaleaEG01)
##          notratado01 notratado02 notratado03 notratado04 notratado05
## gen00001  0.79957746   1.4943258  0.40687822    1.816425  -1.2727249
## gen00002 -1.62600451   0.5128004  0.87907447    1.077089   1.5792660
## gen00003 -1.95245273  -1.0258133 -0.06373976    1.447735  -0.5931999
## gen00004  0.36642541   1.6125608 -1.51397604    0.895905  -1.8915467
## gen00005 -0.03140225  -1.8225208  0.48047472   -1.562934  -0.3440984
## gen00006 -1.91934570  -0.5603091  0.47769946    1.185907   0.3340634
##          notratado06   tratado01  tratado02 tratado03 tratado04  tratado05
## gen00001  -1.1624271 -0.14484139 -2.6215491 -1.864603  3.116052 -2.2350171
## gen00002  -0.9147813  3.59555582 -1.4881947  3.905833  3.002433  2.1280332
## gen00003   0.3099358  0.05924968  1.8103155  1.259982 -3.680984 -1.2667359
## gen00004  -1.1032403 -3.93094174 -0.3781865  1.681134  1.748022  1.7376943
## gen00005  -1.7929110  2.53134453 -1.9637646 -3.571918  2.462766 -0.1786467
## gen00006  -0.2731466 -0.44784125  1.9372120 -2.828150  3.601530  0.3162603
##           tratado06
## gen00001 -2.2416845
## gen00002 -3.7134398
## gen00003 -2.5608207
## gen00004 -1.7946701
## gen00005  1.4077508
## gen00006 -0.8776999

Apartado 3.

Conversión de “datosaleaEG01” en data.frame:

dfEG01 <- as.data.frame(datosaleaEG01)
head(dfEG01)
##          notratado01 notratado02 notratado03 notratado04 notratado05
## gen00001  0.79957746   1.4943258  0.40687822    1.816425  -1.2727249
## gen00002 -1.62600451   0.5128004  0.87907447    1.077089   1.5792660
## gen00003 -1.95245273  -1.0258133 -0.06373976    1.447735  -0.5931999
## gen00004  0.36642541   1.6125608 -1.51397604    0.895905  -1.8915467
## gen00005 -0.03140225  -1.8225208  0.48047472   -1.562934  -0.3440984
## gen00006 -1.91934570  -0.5603091  0.47769946    1.185907   0.3340634
##          notratado06   tratado01  tratado02 tratado03 tratado04  tratado05
## gen00001  -1.1624271 -0.14484139 -2.6215491 -1.864603  3.116052 -2.2350171
## gen00002  -0.9147813  3.59555582 -1.4881947  3.905833  3.002433  2.1280332
## gen00003   0.3099358  0.05924968  1.8103155  1.259982 -3.680984 -1.2667359
## gen00004  -1.1032403 -3.93094174 -0.3781865  1.681134  1.748022  1.7376943
## gen00005  -1.7929110  2.53134453 -1.9637646 -3.571918  2.462766 -0.1786467
## gen00006  -0.2731466 -0.44784125  1.9372120 -2.828150  3.601530  0.3162603
##           tratado06
## gen00001 -2.2416845
## gen00002 -3.7134398
## gen00003 -2.5608207
## gen00004 -1.7946701
## gen00005  1.4077508
## gen00006 -0.8776999
class(dfEG01)
## [1] "data.frame"

Adición 100 valores NA de forma aleatoria a momentos muestrales “notratados” y “tratados”.

set.seed(789)
num_NA = 100

muestra_fila_mitad_1 = sample(1:num_genes, num_NA, replace = TRUE)
muestra_columna_mitad_1 = sample(1:6, num_NA, replace = TRUE)
muestra_fila_mitad_2 = sample(1:num_genes, num_NA, replace = FALSE)
muestra_columna_mitad_2 = sample(7:12, num_NA, replace = TRUE)


dfEG01[[1]][muestra_fila_mitad_1[muestra_columna_mitad_1==1]] = NA
dfEG01[[2]][muestra_fila_mitad_1[muestra_columna_mitad_1==2]] = NA
dfEG01[[3]][muestra_fila_mitad_1[muestra_columna_mitad_1==3]] = NA
dfEG01[[4]][muestra_fila_mitad_1[muestra_columna_mitad_1==4]] = NA
dfEG01[[5]][muestra_fila_mitad_1[muestra_columna_mitad_1==5]] = NA
dfEG01[[6]][muestra_fila_mitad_1[muestra_columna_mitad_1==6]] = NA

paste("Número de NA totales en momento muestral No_tratados:")
## [1] "Número de NA totales en momento muestral No_tratados:"
sum( is.na(dfEG01[,1:6]) ) 
## [1] 100
dfEG01[[7]][muestra_fila_mitad_2[muestra_columna_mitad_2==7]] = NA
dfEG01[[8]][muestra_fila_mitad_2[muestra_columna_mitad_2==8]] = NA
dfEG01[[9]][muestra_fila_mitad_2[muestra_columna_mitad_2==9]] = NA
dfEG01[[10]][muestra_fila_mitad_2[muestra_columna_mitad_2==10]] = NA
dfEG01[[11]][muestra_fila_mitad_2[muestra_columna_mitad_2==11]] = NA
dfEG01[[12]][muestra_fila_mitad_2[muestra_columna_mitad_2==12]] = NA


paste("Número de NA totales en momento muestral tratados:")
## [1] "Número de NA totales en momento muestral tratados:"
sum( is.na(dfEG01[,7:12]) ) 
## [1] 100

Apartado 4.

Construya dos vectores de tipo factor de tamaño 12 asociados a los 12 momentos muestrales.

vtratamiento <- factor(c(rep("notratados",6),rep("tratados",6)))
vtratamiento
##  [1] notratados notratados notratados notratados notratados notratados
##  [7] tratados   tratados   tratados   tratados   tratados   tratados  
## Levels: notratados tratados
vsexo <- factor(rep(c("hombre","mujer"),6))
vsexo
##  [1] hombre mujer  hombre mujer  hombre mujer  hombre mujer  hombre mujer 
## [11] hombre mujer 
## Levels: hombre mujer

Creación del data.frame que tenga como columnas los dos vectores creados anteriormente:

dfCar01 <- data.frame(Tratamiento = vtratamiento, Sexo = vsexo)
dfCar01
##    Tratamiento   Sexo
## 1   notratados hombre
## 2   notratados  mujer
## 3   notratados hombre
## 4   notratados  mujer
## 5   notratados hombre
## 6   notratados  mujer
## 7     tratados hombre
## 8     tratados  mujer
## 9     tratados hombre
## 10    tratados  mujer
## 11    tratados hombre
## 12    tratados  mujer

Apartado 5.

Creación objeto RData con objetos dfEG01 y dfCar01:

datosEG01 <- c(dfEG01,dfCar01)
save(datosEG01, file =  "datosEG01.RData")

Apartado 6.

¿Qué momento muestral tiene mayor número de valores NA?

summary(dfEG01)
##   notratado01         notratado02        notratado03         notratado04       
##  Min.   :-1.999864   Min.   :-1.99852   Min.   :-1.998455   Min.   :-1.999938  
##  1st Qu.:-0.975421   1st Qu.:-1.00692   1st Qu.:-0.988136   1st Qu.:-0.988395  
##  Median : 0.013646   Median :-0.02164   Median : 0.000659   Median :-0.020727  
##  Mean   : 0.009644   Mean   :-0.01131   Mean   : 0.005422   Mean   :-0.008069  
##  3rd Qu.: 1.010238   3rd Qu.: 0.98089   3rd Qu.: 0.983029   3rd Qu.: 0.990407  
##  Max.   : 1.999802   Max.   : 2.00000   Max.   : 1.999921   Max.   : 1.999953  
##  NA's   :9           NA's   :21         NA's   :20          NA's   :15         
##   notratado05        notratado06         tratado01           tratado02        
##  Min.   :-1.99981   Min.   :-1.99939   Min.   :-3.999126   Min.   :-3.999282  
##  1st Qu.:-0.96817   1st Qu.:-0.98651   1st Qu.:-1.953489   1st Qu.:-2.034113  
##  Median : 0.02784   Median : 0.02738   Median : 0.023218   Median : 0.005791  
##  Mean   : 0.02294   Mean   : 0.01374   Mean   : 0.009655   Mean   : 0.013956  
##  3rd Qu.: 1.02664   3rd Qu.: 1.02123   3rd Qu.: 1.976667   3rd Qu.: 2.010523  
##  Max.   : 1.99916   Max.   : 1.99977   Max.   : 3.999996   Max.   : 3.998916  
##  NA's   :17         NA's   :18         NA's   :25          NA's   :13         
##    tratado03           tratado04           tratado05          tratado06        
##  Min.   :-3.998357   Min.   :-3.999820   Min.   :-3.99987   Min.   :-3.999772  
##  1st Qu.:-2.011652   1st Qu.:-1.975206   1st Qu.:-2.02836   1st Qu.:-2.001994  
##  Median : 0.035676   Median : 0.005697   Median :-0.05363   Median : 0.019204  
##  Mean   :-0.001685   Mean   : 0.017464   Mean   :-0.02608   Mean   :-0.006787  
##  3rd Qu.: 1.982843   3rd Qu.: 2.057326   3rd Qu.: 1.96745   3rd Qu.: 1.963359  
##  Max.   : 3.999944   Max.   : 3.999746   Max.   : 3.99898   Max.   : 3.999814  
##  NA's   :13          NA's   :15          NA's   :18         NA's   :16
frec_NA_columnas <- colSums(is.na(dfEG01))
frec_NA_columnas
## notratado01 notratado02 notratado03 notratado04 notratado05 notratado06 
##           9          21          20          15          17          18 
##   tratado01   tratado02   tratado03   tratado04   tratado05   tratado06 
##          25          13          13          15          18          16
mom_muestral_max_NA <- frec_NA_columnas[which.max(frec_NA_columnas)]
mom_muestral_max_NA
## tratado01 
##        25
colmax =  which.max(frec_NA_columnas)
colmax
## tratado01 
##         7

El momento muestral con número máximo de NA es “tratado01” con 25 valores NA.

¿En qué posiciones se encuentran en ese momento muestral?

posiciones_NA_colmax <- which(is.na(dfEG01[[colmax]]))
posiciones_NA_colmax
##  [1]   52   58  892 1503 1902 2292 2357 3009 3176 3749 4781 4807 4813 4877 5034
## [16] 5207 5456 5876 6192 6322 6366 6376 7275 9280 9912

Los números indicados son las posiciones donde se encuentran los valores NA en la columna de momento muestral “tratado01”.

Apartado 7.

Obtenga en un vector la suma de los valores de las expresiones génicas para los genes “gen00048” y “gen00257”.

suma_valores_gen48 <- sum(dfEG01["gen00048",])
suma_valores_gen257 <- sum(dfEG01["gen00257",])

suma_genes = suma_valores_gen48+suma_valores_gen257

vector_suma_genes = c(suma_genes)
(vector_suma_genes)
## [1] -13.74306

Apartado 8.

Construya un data.frame que contenga únicamente los genes que tienen sus valores completos o dicho de otro modo todos sus valores son no NA. Guárdelo en un fichero: “dfEG01noNA.RData”:

Con la función complete.cases se forma un vector sin los “missing values” NA.

paste("Dimensiones data.frame sin procesar valores NA:")
## [1] "Dimensiones data.frame sin procesar valores NA:"
dim(dfEG01)
## [1] 10000    12
dfEG01noNA = dfEG01[complete.cases(dfEG01),]
class(dfEG01noNA)
## [1] "data.frame"
paste("Dimensiones data.frame tras procesar valores NA:")
## [1] "Dimensiones data.frame tras procesar valores NA:"
dim(dfEG01noNA)
## [1] 9801   12
head(dfEG01noNA)
##          notratado01 notratado02 notratado03 notratado04 notratado05
## gen00001  0.79957746   1.4943258  0.40687822    1.816425  -1.2727249
## gen00002 -1.62600451   0.5128004  0.87907447    1.077089   1.5792660
## gen00003 -1.95245273  -1.0258133 -0.06373976    1.447735  -0.5931999
## gen00004  0.36642541   1.6125608 -1.51397604    0.895905  -1.8915467
## gen00005 -0.03140225  -1.8225208  0.48047472   -1.562934  -0.3440984
## gen00006 -1.91934570  -0.5603091  0.47769946    1.185907   0.3340634
##          notratado06   tratado01  tratado02 tratado03 tratado04  tratado05
## gen00001  -1.1624271 -0.14484139 -2.6215491 -1.864603  3.116052 -2.2350171
## gen00002  -0.9147813  3.59555582 -1.4881947  3.905833  3.002433  2.1280332
## gen00003   0.3099358  0.05924968  1.8103155  1.259982 -3.680984 -1.2667359
## gen00004  -1.1032403 -3.93094174 -0.3781865  1.681134  1.748022  1.7376943
## gen00005  -1.7929110  2.53134453 -1.9637646 -3.571918  2.462766 -0.1786467
## gen00006  -0.2731466 -0.44784125  1.9372120 -2.828150  3.601530  0.3162603
##           tratado06
## gen00001 -2.2416845
## gen00002 -3.7134398
## gen00003 -2.5608207
## gen00004 -1.7946701
## gen00005  1.4077508
## gen00006 -0.8776999
save(dfEG01noNA, file = "dfEG01noNA.RData")

Apartado 9.

Utilizando el data.frame del apartado anterior, construya un nuevo data.frame que contenga únicamente las expresiones génicas de los momentos muestrales de: “tratados y mujeres”, utilizando el objeto: “dfCar01” (llámelo: “dfEG01_red01”). Muestre la dimensión del data.frame obtenido y sus primeras 6 filas.

Data.frame con expresiones génicas de momentos muestrales “tratados” y “mujeres”.

dfEG01_red01 <- dfEG01noNA[, dfCar01$Tratamiento == "tratados" & dfCar01$Sexo == "mujer"]
paste("Dimensiones dfEG01_red01:")
## [1] "Dimensiones dfEG01_red01:"
dim(dfEG01_red01)
## [1] 9801    3
head(dfEG01_red01)
##           tratado02 tratado04  tratado06
## gen00001 -2.6215491  3.116052 -2.2416845
## gen00002 -1.4881947  3.002433 -3.7134398
## gen00003  1.8103155 -3.680984 -2.5608207
## gen00004 -0.3781865  1.748022 -1.7946701
## gen00005 -1.9637646  2.462766  1.4077508
## gen00006  1.9372120  3.601530 -0.8776999

Apartado 10.

Orden decreciente de los 6 genes que suman mayor expresión génica en todos los momentos muestrales.

suma_todos_genes <- apply(dfEG01noNA, FUN = sum, MARGIN = 1)
sort(suma_todos_genes, decreasing = TRUE)[1:6]
## gen01983 gen03832 gen01602 gen07599 gen01814 gen09992 
## 21.68290 21.22535 20.64785 19.68284 19.19854 19.02482

Ejercicio A.2

Se consideran los datos obtenidos en un estudio nutrigenómico sobre el ratón, los cuales se encuentran en el paquete del proyecto Bioconductor: “mixOmics”, en el que el objeto tipo “list”: “nutrimouse” cuyos elementos contiene toda la información relacionada (ver: “Descripción de los datos nutrimouse”): dos primeros son data. frame con información de expresión génica (en filas observaciones muestrales de los 40 ratones estudiados y en las columnas se encuentran una representación de los genes más relevantes en células hepáticas y concentraciones (%) de ácidos grasos hepáticos): (1) “gene”, (2) “lipid”, y los dos últimos son factores: (3) “genotype”, (4) “diet”.

Con dicha información:

Apartado 1.

Cargue la información del objeto: “nutrimouse”. Obtenga las dimensiones de los dos data.frame que contiene. Y construya una tabla de frecuencias de los dos factores: “diet” y “genotype”, para ver cómo fueel diseño en la clasificación de los 40 ratones. Cuántas observaciones NA existen en cada uno de los data.frame:

  • Carga la información del objeto “nutrimouse” procedente del paquete de Bioconductor “mixOmics”.
library(mixOmics)
## Loading required package: MASS
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Loaded mixOmics 6.24.0
## Thank you for using mixOmics!
## Tutorials: http://mixomics.org
## Bookdown vignette: https://mixomicsteam.github.io/Bookdown
## Questions, issues: Follow the prompts at http://mixomics.org/contact-us
## Cite us:  citation('mixOmics')
data(nutrimouse)
class(nutrimouse)
## [1] "list"

Dentro de la base de datos “nutrimouse” existen 2 data.frames.

El primero de ellos llamado “gene” representa la expresión génica de 120 genes expresados en células hepáticas de los 40 ratones estudiados.

lista_genes <- (nutrimouse$gene)
class(lista_genes)
## [1] "data.frame"
paste("Dimensiones data.frame gene:")
## [1] "Dimensiones data.frame gene:"
dim(lista_genes)
## [1]  40 120
head(lista_genes)
##   X36b4 ACAT1 ACAT2  ACBP  ACC1  ACC2 ACOTH ADISP ADSS1 ALDH3  AM2R   AOX  BACT
## 1 -0.42 -0.65 -0.84 -0.34 -1.29 -1.13 -0.93 -0.98 -1.19 -0.68 -0.59 -0.16 -0.22
## 2 -0.44 -0.68 -0.91 -0.32 -1.23 -1.06 -0.99 -0.97 -1.00 -0.62 -0.58 -0.12 -0.32
## 3 -0.48 -0.74 -1.10 -0.46 -1.30 -1.09 -1.06 -1.08 -1.18 -0.75 -0.66 -0.16 -0.32
## 4 -0.45 -0.69 -0.65 -0.41 -1.26 -1.09 -0.93 -1.02 -1.07 -0.71 -0.65 -0.17 -0.32
## 5 -0.42 -0.71 -0.54 -0.38 -1.21 -0.89 -1.00 -0.95 -1.08 -0.76 -0.59 -0.31 -0.31
## 6 -0.43 -0.69 -0.80 -0.32 -1.13 -0.79 -0.93 -0.97 -1.07 -0.75 -0.55 -0.23 -0.29
##    BIEN  BSEP Bcl.3 C16SR  CACP  CAR1   CBS CIDEA  COX1  COX2  CPT2 CYP24 CYP26
## 1 -0.89 -0.69 -1.18  1.66 -0.92 -0.97 -0.26 -1.21 -1.11 -1.18 -0.87 -1.37 -1.21
## 2 -0.88 -0.60 -1.07  1.65 -0.87 -0.92 -0.36 -1.17 -1.06 -1.06 -0.87 -1.14 -1.12
## 3 -0.89 -0.70 -1.17  1.57 -1.02 -0.98 -0.40 -1.29 -1.17 -1.14 -0.95 -1.30 -1.22
## 4 -0.77 -0.67 -1.12  1.61 -0.89 -0.97 -0.39 -1.18 -1.03 -1.13 -0.88 -1.27 -1.18
## 5 -0.97 -0.68 -0.93  1.66 -0.93 -1.06 -0.35 -1.15 -0.99 -1.10 -0.91 -1.20 -1.16
## 6 -0.84 -0.55 -1.08  1.70 -0.97 -1.03 -0.31 -1.14 -1.03 -1.16 -0.92 -1.11 -1.10
##   CYP27a1 CYP27b1 CYP2b10 CYP2b13 CYP2c29 CYP3A11 CYP4A10 CYP4A14 CYP7a CYP8b1
## 1   -0.71   -1.31   -1.23   -1.19   -0.06   -0.09   -0.81   -0.81 -0.77  -0.77
## 2   -0.62   -1.14   -1.20   -1.06   -0.20   -0.34   -0.88   -0.84 -0.71  -0.63
## 3   -0.78   -1.29   -1.32   -1.25   -0.30   -0.45   -0.71   -0.98 -0.93  -0.53
## 4   -0.71   -1.27   -1.23   -1.13   -0.07   -0.11   -0.65   -0.41 -0.80  -0.73
## 5   -0.69   -1.20   -1.22   -1.10   -0.29   -0.51   -1.16   -1.16 -0.71  -0.51
## 6   -0.60   -1.15   -1.10   -1.07   -0.28   -0.55   -0.99   -1.09 -0.74  -0.55
##     FAS   FAT  FDFT   FXR G6PDH G6Pase    GK    GS  GSTa GSTmu GSTpi2 HMGCoAred
## 1 -0.41 -1.03 -0.98 -0.93 -1.22  -0.46 -0.71 -1.24  0.00  0.02   0.45     -0.95
## 2 -0.37 -0.98 -0.92 -0.87 -1.09  -0.63 -0.67 -1.22 -0.05 -0.05   0.30     -0.86
## 3 -0.30 -1.03 -1.04 -1.00 -1.28  -1.06 -0.68 -1.36 -0.13 -0.19   0.18     -0.96
## 4 -0.59 -1.06 -1.00 -0.90 -1.19  -0.71 -0.75 -1.21 -0.09  0.03   0.36     -1.02
## 5 -0.06 -0.99 -0.99 -0.89 -1.16  -0.58 -0.62 -1.22 -0.02 -0.23   0.30     -0.70
## 6  0.18 -0.99 -1.00 -0.89 -0.96  -0.49 -0.59 -1.16 -0.11 -0.05   0.17     -0.76
##   HPNCL  IL.2 L.FABP   LCE  LDLr   LPK   LPL  LXRa  LXRb  Lpin Lpin1 Lpin2
## 1 -0.65 -0.94   0.24  0.09 -0.82 -0.32 -1.01 -0.82 -1.00 -0.87 -0.85 -0.85
## 2 -0.69 -0.94   0.27  0.06 -0.68 -0.39 -0.97 -0.82 -0.95 -0.97 -0.99 -0.87
## 3 -0.75 -1.16   0.17 -0.05 -0.82 -0.38 -1.11 -0.91 -1.16 -0.95 -0.94 -0.90
## 4 -0.61 -0.97   0.16  0.01 -0.94 -0.38 -0.99 -0.85 -1.01 -1.00 -1.02 -0.88
## 5 -0.66 -0.93   0.00 -0.07 -0.73 -0.17 -1.05 -0.83 -1.01 -0.57 -0.53 -0.72
## 6 -0.56 -0.96   0.23 -0.10 -0.74 -0.14 -0.99 -0.79 -0.99 -0.51 -0.51 -0.68
##   Lpin3 M.CPT1  MCAD  MDR1  MDR2  MRP6    MS MTHFR NGFiB NURR1  Ntcp OCTN2
## 1 -1.23  -1.15 -0.60 -1.15 -0.77 -0.99 -1.11 -0.96 -1.21 -1.21 -0.49 -1.15
## 2 -1.12  -1.06 -0.62 -1.10 -0.65 -0.85 -1.06 -0.99 -1.08 -1.10 -0.45 -1.15
## 3 -1.25  -1.26 -0.70 -1.26 -0.86 -0.90 -1.20 -1.10 -1.24 -1.32 -0.44 -1.20
## 4 -1.18  -1.10 -0.59 -1.13 -0.77 -0.95 -1.09 -0.95 -1.12 -1.11 -0.54 -1.17
## 5 -1.12  -1.11 -0.69 -1.11 -0.70 -0.91 -1.09 -0.93 -1.11 -1.14 -0.47 -1.19
## 6 -1.09  -1.14 -0.66 -1.09 -0.69 -0.84 -1.09 -0.96 -1.04 -1.18 -0.46 -1.11
##     PAL  PDK4  PECI  PLTP PMDCI   PON PPARa PPARd PPARg   PXR Pex11a  RARa
## 1 -1.32 -1.16 -0.68 -1.10 -0.52 -0.52 -0.93 -1.51 -1.06 -0.99  -1.00 -1.20
## 2 -1.25 -1.16 -0.69 -0.99 -0.52 -0.55 -0.86 -1.59 -1.02 -0.96  -1.02 -1.06
## 3 -1.16 -1.27 -0.92 -1.03 -0.60 -0.65 -0.95 -1.71 -1.14 -1.10  -1.20 -1.16
## 4 -1.25 -1.16 -0.71 -1.08 -0.52 -0.64 -0.97 -1.57 -1.05 -0.99  -1.00 -1.17
## 5 -1.24 -1.13 -0.83 -0.98 -0.71 -0.57 -0.94 -1.53 -1.09 -1.00  -0.95 -1.15
## 6 -1.02 -1.08 -0.81 -0.89 -0.69 -0.63 -0.95 -1.56 -1.01 -1.03  -1.07 -1.13
##   RARb2  RXRa RXRb2 RXRg1   S14  SHP1 SIAT4c SPI1.1 SR.BI   THB THIOL   TRa
## 1 -1.19 -0.67 -0.95 -1.16 -0.93 -1.10  -1.07   1.19 -0.84 -0.79 -0.18 -1.48
## 2 -1.11 -0.59 -0.95 -1.10 -0.86 -0.97  -0.97   1.15 -0.86 -0.85 -0.15 -1.46
## 3 -1.23 -0.68 -1.07 -1.21 -0.84 -1.09  -1.04   1.09 -0.95 -0.92 -0.24 -1.58
## 4 -1.16 -0.72 -0.95 -1.10 -1.05 -1.03  -0.99   1.07 -0.95 -0.79 -0.15 -1.54
## 5 -1.14 -0.78 -0.98 -1.11 -0.65 -1.13  -0.94   1.22 -1.06 -0.84 -0.35 -1.46
## 6 -1.07 -0.62 -0.94 -1.03 -0.40 -0.98  -0.93   1.05 -0.80 -0.86 -0.29 -1.44
##     TRb Tpalpha Tpbeta  UCP2  UCP3   VDR VLDLr  Waf1   ap2 apoA.I  apoB apoC3
## 1 -1.07   -0.69  -1.11 -1.06 -1.19 -1.17 -1.08 -1.18 -1.19   0.76 -0.12 -0.49
## 2 -1.00   -0.74  -1.09 -0.99 -1.10 -1.17 -0.99 -1.12 -1.16   0.86 -0.13 -0.35
## 3 -1.16   -0.81  -1.14 -1.05 -1.20 -1.30 -1.13 -1.30 -1.37   0.82 -0.19 -0.42
## 4 -1.11   -0.74  -1.04 -1.07 -1.12 -1.17 -1.01 -1.14 -1.19   0.71 -0.25 -0.39
## 5 -1.01   -0.82  -1.20 -0.98 -1.10 -1.18 -1.06 -1.15 -1.12   0.83 -0.17 -0.33
## 6 -1.00   -0.76  -1.05 -0.94 -1.07 -1.11 -1.04 -1.13 -1.17   0.90 -0.17 -0.30
##   apoE c.fos cHMGCoAS cMOAT eif2g hABC1 i.BABP i.BAT i.FABP i.NOS mABC1
## 1 1.08 -1.15    -1.15 -0.78 -1.23 -1.16  -0.78 -1.65  -1.14 -1.24 -0.85
## 2 1.12 -1.08    -0.95 -0.73 -1.02 -1.11  -0.73 -1.67  -1.03 -1.20 -0.84
## 3 1.04 -1.18    -0.93 -0.89 -1.14 -1.25  -0.89 -1.89  -1.16 -1.35 -0.96
## 4 1.05 -1.11    -1.10 -0.93 -1.10 -1.17  -0.82 -1.70  -1.17 -1.25 -0.87
## 5 1.08 -1.08    -0.94 -0.84 -1.10 -1.14  -0.71 -1.68  -1.15 -1.20 -0.82
## 6 1.12 -1.06    -0.83 -0.81 -1.05 -1.08  -0.76 -1.68  -1.08 -1.16 -0.81
##   mHMGCoAS
## 1    -0.03
## 2    -0.12
## 3    -0.12
## 4    -0.12
## 5    -0.29
## 6    -0.15
paste("Número de missing values:")
## [1] "Número de missing values:"
sum(is.na(lista_genes))
## [1] 0

No hay datos NA en el data.frame “gene”.

El segundo data.frame “lipid” representa la concentración de 21 ácidos grasos hepáticos medidos en los ratones.

lista_hepatico <- (nutrimouse$lipid)
class(lista_hepatico)
## [1] "data.frame"
paste("Dimensiones data.frame lipid:")
## [1] "Dimensiones data.frame lipid:"
dim(lista_hepatico)
## [1] 40 21
head(lista_hepatico)
##   C14.0 C16.0 C18.0 C16.1n.9 C16.1n.7 C18.1n.9 C18.1n.7 C20.1n.9 C20.3n.9
## 1  0.34 26.45 10.22     0.35     3.10    16.98     2.41     0.26     0.00
## 2  0.38 24.04  9.93     0.55     2.54    20.07     3.92     0.23     0.00
## 3  0.36 23.70  8.96     0.55     2.65    22.89     3.96     0.26     0.19
## 4  0.22 25.48  8.14     0.49     2.82    21.92     2.52     0.00     0.00
## 5  0.37 24.80  9.63     0.46     2.85    21.38     2.96     0.30     0.27
## 6  1.70 26.04  6.59     0.66     7.26    28.23     8.99     0.36     2.89
##   C18.2n.6 C18.3n.6 C20.2n.6 C20.3n.6 C20.4n.6 C22.4n.6 C22.5n.6 C18.3n.3
## 1     8.93     0.00     0.00     0.78     3.07     0.00     0.00     5.97
## 2    14.98     0.30     0.30     1.64    15.34     0.58     2.10     0.00
## 3    16.06     0.27     0.33     1.51    13.27     0.54     1.77     0.00
## 4    13.89     0.00     0.00     1.10     3.92     0.00     0.00     0.49
## 5    14.55     0.27     0.23     1.58    11.85     0.32     0.44     0.42
## 6     3.47     2.66     0.00     0.81     5.09     0.00     0.56     0.00
##   C20.3n.3 C20.5n.3 C22.5n.3 C22.6n.3
## 1     0.37     8.62     1.75    10.39
## 2     0.00     0.00     0.48     2.61
## 3     0.00     0.00     0.22     2.51
## 4     0.00     2.99     1.04    14.99
## 5     0.00     0.30     0.35     6.69
## 6     0.00     0.00     2.13     2.56
paste("Número de missing values:")
## [1] "Número de missing values:"
sum(is.na(lista_hepatico))
## [1] 0

No hay datos NA en el data.frame “lipid”.

  • Tabla de frecuencias factores “diet” y “genotipe”
dieta <- nutrimouse$diet
tabla_dieta <- table(dieta)
genotipo <- nutrimouse$genotype
tabla_genotipo <- table(genotipo)
table(nutrimouse$diet,nutrimouse$genotype)
##       
##        wt ppar
##   coc   4    4
##   fish  4    4
##   lin   4    4
##   ref   4    4
##   sun   4    4

En el estudio se emplearon 2 genotipos, salvaje(wt) y PPAR.

Para cada uno de estos dos genotipos se experimentó con 5 dietas diferentes: ref(aceite maíz/colza 50/50), coc(aceite coco hidrogenado), sun(aceite de girasol), lin(aceite linaza) y fish(aceites pescado enriquecidos con maíz).

Siendo dos genotipos distintos y 5 dietas experimentales, para realizar 4 repeticiones en cada uno de los factores se hace necesario un total de 40 ratones para el estudio.

Apartado 2.

  • Selección de ratones con dieta “fish” y genotipo “ppar”:
ratones <- (nutrimouse$diet == "fish" & nutrimouse$genotype == "ppar")
(posicion <- which(ratones))
## [1] 24 29 35 38

Los ratones 24, 29, 35 y 38 son los de genotipo “ppar” que fueron alimentados con dieta “fish”.

  • Compare (gráficos y medidas estadísticas) las expresiones génicas en “lipid” considerando los 21 valores de los 4 ratones.:
lipid_raton <- nutrimouse$lipid[ratones,]
lipid_raton
##    C14.0 C16.0 C18.0 C16.1n.9 C16.1n.7 C18.1n.9 C18.1n.7 C20.1n.9 C20.3n.9
## 24  0.44 22.73  4.71     0.75     2.27    25.10     1.91     0.38        0
## 29  0.40 21.70  6.12     0.78     2.07    22.49     1.82     0.44        0
## 35  0.39 20.71  6.41     0.71     2.18    24.09     1.83     0.40        0
## 38  0.36 25.23 10.33     0.45     1.59    16.49     1.53     0.37        0
##    C18.2n.6 C18.3n.6 C20.2n.6 C20.3n.6 C20.4n.6 C22.4n.6 C22.5n.6 C18.3n.3
## 24    22.65        0     0.27     0.38     1.27        0     0.21     1.16
## 29    22.05        0     0.24     0.47     1.87        0     0.19     1.31
## 35    23.09        0     0.24     0.40     1.90        0     0.00     1.45
## 38    18.25        0     0.26     0.57     3.02        0     0.00     0.71
##    C20.3n.3 C20.5n.3 C22.5n.3 C22.6n.3
## 24        0     2.17     2.04    11.56
## 29        0     2.65     2.16    13.26
## 35        0     2.82     1.80    11.57
## 38        0     2.91     1.64    16.28
lipid_num <- as.numeric(unlist(lipid_raton))
hist(lipid_num,main = "Histograma porcentaje hepático", breaks=30,xlab = "% hepático")

caja_bigote <- boxplot(lipid_num, main="Boxplot porcentaje hepático", horizontal = TRUE)

Apartado 3.

Estudiar ratón 38 en los datos “lipid”:

(nutrimouse$lipid[38,])
##    C14.0 C16.0 C18.0 C16.1n.9 C16.1n.7 C18.1n.9 C18.1n.7 C20.1n.9 C20.3n.9
## 38  0.36 25.23 10.33     0.45     1.59    16.49     1.53     0.37        0
##    C18.2n.6 C18.3n.6 C20.2n.6 C20.3n.6 C20.4n.6 C22.4n.6 C22.5n.6 C18.3n.3
## 38    18.25        0     0.26     0.57     3.02        0        0     0.71
##    C20.3n.3 C20.5n.3 C22.5n.3 C22.6n.3
## 38        0     2.91     1.64    16.28
raton_38 <- as.numeric(nutrimouse$lipid[38,])
raton_38
##  [1]  0.36 25.23 10.33  0.45  1.59 16.49  1.53  0.37  0.00 18.25  0.00  0.26
## [13]  0.57  3.02  0.00  0.00  0.71  0.00  2.91  1.64 16.28
  • a: Crear una tabla de frecuencias en 7 intervalos de todos los 21 valoes:
intervalos <- cut(raton_38,breaks = 7)
absoluta=table(intervalos)
absacum=cumsum(absoluta)
relativa= prop.table(absoluta) 
relacum= cumsum(relativa)
tablaAg= cbind(absoluta,absacum,relativa,relacum)
tablaAg
##               absoluta absacum   relativa   relacum
## (-0.0252,3.6]       16      16 0.76190476 0.7619048
## (3.6,7.21]           0      16 0.00000000 0.7619048
## (7.21,10.8]          1      17 0.04761905 0.8095238
## (10.8,14.4]          0      17 0.00000000 0.8095238
## (14.4,18]            2      19 0.09523810 0.9047619
## (18,21.6]            1      20 0.04761905 0.9523810
## (21.6,25.3]          1      21 0.04761905 1.0000000

El intervalo que mayor número de observaciones presenta es el intervalo (-0.0252,3.6].

El valor 4.18 es el que deja un porcentaje menor o igual al 90.48%. Es decir, el 90.48% de los valores son menores o iguales a 4.18.

  • b: Gráfico simultáneo histograma y diagrama de caja y bigote.
par(mfrow = c(1,2))
caja_big_raton_38 <- boxplot(raton_38, horizontal = TRUE, main = "Boxplot datos lipid ratón 38")
hist_raton_38 = hist(raton_38, breaks = 20, main = "Histograma datos lipid ratón 38", ylab = "Frecuencia", xlab = "%hepático", col = "lightblue", border = "black")

  • c: Detección de outliers por método rango intercuartílico:
vX = raton_38
gbp01<-boxplot(vX, horizontal = TRUE)

Por este método se han obtenido 5 valores outliers como son:

head(gbp01$out, 50)
## [1] 25.23 10.33 16.49 18.25 16.28

Detección outliers método desviación típica:

vX = raton_38
vX
##  [1]  0.36 25.23 10.33  0.45  1.59 16.49  1.53  0.37  0.00 18.25  0.00  0.26
## [13]  0.57  3.02  0.00  0.00  0.71  0.00  2.91  1.64 16.28
k = 3
int_inlier <- c(mean(vX)-k*sd(vX), mean(vX)+k*sd(vX))
int_inlier
## [1] -18.12230  27.64516
ind_out = which( (vX < int_inlier[1]) | (vX > int_inlier[2]) )
vX_outliers = vX[ind_out]
vX_outliers
## numeric(0)
length(vX_outliers)
## [1] 0
range(vX_outliers)
## Warning in min(x): ningún argumento finito para min; retornando Inf
## Warning in max(x): ningun argumento finito para max; retornando -Inf
## [1]  Inf -Inf

No se ham podido detectar outliers con el método de la desviación típica.

Ejercicio B

El conjunto de datos “breast.TCGA” del paquete mixOmics contiene la expresión o abundancia de tres conjuntos de datos omicos coincidentes: mRNA y miRNA y proteómica para 150 muestras de cáncer de mama (Basal, Her2, Luminal A) en el conjunto de entrenamiento, y 70 muestras en el conjunto de Test.

En este ejercicio vamos a trabajar con los datos del conjunto Test pertenecientes al mrna:

Trabajaremos con los datos del conjunto Test pertenecientes al mRNA:

library(mixOmics)
data("breast.TCGA")

Transformación datos matriciales en dataset:

mrna_test <- breast.TCGA[["data.test"]][["mrna"]]
dim(mrna_test)
## [1]  70 200
class(mrna_test)
## [1] "matrix" "array"
datos <- as.data.frame(mrna_test)
head(datos)
##          RTN2    NDRG2  CCDC113   FAM63A    ACADS     GMDS    HLA-H   SEMA4A
## A54N 1.187128 7.519887 3.819383 4.827405 5.255628 5.867395 5.851815 3.422738
## A2NL 2.729705 7.064522 2.435517 4.329876 4.766857 6.775768 7.209585 4.621781
## A6VY 3.054289 5.978128 3.546771 4.478048 4.205232 5.043232 9.538554 6.182447
## A3XT 2.700442 7.061261 1.968639 2.671635 5.019725 5.063346 7.829317 7.024889
## A1F0 3.136745 7.660714 3.501533 4.269159 3.716986 3.847745 7.073876 5.957018
## A26I 1.555233 4.206059 3.436369 6.068822 3.573562 4.912124 5.794861 4.825844
##          ETS2    LIMD2     NME3      ZEB1    CDCP1    GIYD2    RTKN2   MANSC1
## A54N 5.168055 5.062627 6.934305 0.5514035 5.223230 5.537991 3.527975 3.976379
## A2NL 6.279341 6.983768 4.541756 2.6032285 3.389409 4.588594 4.450123 5.461346
## A6VY 6.418312 4.994062 6.162271 4.1038359 7.131570 5.529871 4.344942 4.385813
## A3XT 6.506575 5.131091 5.981112 3.5087997 6.311185 5.103600 4.157039 1.758701
## A1F0 6.701526 4.874302 3.706619 5.8543827 7.177210 4.638989 3.889588 5.813442
## A26I 6.282206 3.658165 3.778679 5.5782758 6.420026 3.207599 5.409406 4.724025
##          TAGLN    IFIT3    ARL4C    HTRA1   KIF13B   CPPED1    SKAP2     ASPM
## A54N  9.736904 4.630935 8.175230 4.528594 5.845351 1.187128 1.609557 3.781187
## A2NL  5.722427 5.631750 5.569918 5.241667 6.323421 3.052038 4.115293 6.533775
## A6VY  9.063117 6.462693 7.323727 7.865587 4.544098 4.859348 5.403746 5.505593
## A3XT  7.497878 5.755211 6.039350 6.261910 4.452044 5.213291 4.525320 6.509121
## A1F0  8.441572 8.876170 6.576875 7.464052 5.414650 5.918797 5.885867 5.129014
## A26I 10.723074 6.865341 5.455454 7.649439 5.052793 4.478539 4.351738 6.800977
##         KDM4B   TBXAS1     MT1X   MED13L   SNORA8     RGS1     CBX6       WWC2
## A54N 5.283221 2.646303 5.735780 5.967736 3.743502 6.570520 1.879170 -0.7303445
## A2NL 6.165904 2.440171 3.893984 6.701385 5.039805 3.684389 3.091034  4.5865281
## A6VY 6.746111 3.854454 5.988619 6.112593 4.017428 6.456311 3.865005  4.5052819
## A3XT 5.974125 4.669789 4.250094 6.649151 5.201179 6.812270 5.616148  3.3290889
## A1F0 4.846374 4.845783 4.579358 7.247279 4.471865 7.815029 4.493019  4.0416683
## A26I 5.159893 3.154965 4.262694 6.960758 4.554034 6.392959 3.939645  5.2499236
##      TNFRSF12A   ZNF552   MAPRE2    SEMA5A   STAT5A     FLI1  COL15A1  C7orf55
## A54N  6.830123 3.319912 5.277745 -2.241912 6.680443 1.403727 1.383473 4.207087
## A2NL  3.302618 2.208820 6.470343  1.959324 4.503354 3.576488 5.800971 5.309507
## A6VY  7.243114 4.122643 5.226903  3.409365 5.641183 3.770896 6.704371 3.637771
## A3XT  5.545963 2.020964 5.241129  2.862152 5.816116 4.224290 6.910697 4.528336
## A1F0  4.797755 4.097287 6.062403  5.490255 7.879610 6.139090 8.017701 2.749321
## A26I  5.546436 4.325393 4.949096  5.553437 5.833843 3.875443 6.598602 2.787239
##         ASF1B     FUT8    LASS4     SQLE     GPC4     AKAP12      AGL  ADAMTS4
## A54N 4.045684 3.064035 6.226267 4.772001 0.986887 -0.9813895 4.838607 1.849944
## A2NL 6.714172 3.636160 4.853150 5.575746 5.484062  2.3391301 5.240334 4.765003
## A6VY 5.726585 4.934309 5.224174 5.702275 8.686988  2.8849190 3.816908 5.876202
## A3XT 5.659429 4.287787 4.957693 7.467741 4.120203  3.2630375 3.664934 7.139604
## A1F0 4.218983 5.212094 4.707708 5.988112 6.145100  6.3909053 5.631897 4.134499
## A26I 4.300639 4.841710 3.624318 5.656612 6.350627  5.2243770 5.844229 4.282844
##         EPHB3   MAP3K1     PRNP    PROM2  SLCO3A1    SNHG1  PRKCDBP     MXI1
## A54N 8.500221 5.966886 6.243226 6.842696 0.986887 6.972519 5.359084 7.147861
## A2NL 7.640699 4.371894 5.747053 2.170119 3.650650 8.055622 2.235840 6.331003
## A6VY 6.360145 4.641411 7.881696 7.438401 4.993614 6.117381 4.854061 5.009551
## A3XT 7.962835 5.131038 7.849049 6.497884 6.365593 6.191967 4.739108 4.706641
## A1F0 6.911577 6.405200 7.723599 5.777677 4.897970 5.091300 3.251729 4.881359
## A26I 7.418132 6.123656 7.626591 6.941373 4.359748 5.951431 3.320072 5.227116
##         CSF1R    TANC2  SLC19A2     RHOU  C4orf34    LRIG1    DOCK8      BOC
## A54N 5.010240 4.488856 4.785570 5.837927 3.541748 5.297999 2.704609 4.517019
## A2NL 5.062336 5.038747 4.770558 2.931602 3.355317 6.075713 4.153871 6.062914
## A6VY 6.455439 6.022001 4.828233 4.081014 3.479703 5.908754 5.826684 6.682915
## A3XT 7.013994 4.850931 4.729532 5.054308 5.043809 5.898822 5.561749 3.690287
## A1F0 7.697743 4.459455 5.508647 5.341159 4.990070 6.441589 7.447124 5.173448
## A26I 6.058400 7.230066 4.935761 7.745330 4.045242 6.082100 5.186484 4.864360
##      C11orf52  S100A16    NRARP    TTC23   TBC1D4   DEPDC6    ILDR1      SDC1
## A54N 3.101716 8.965490 5.576516 4.105643 5.595350 5.357787 3.372241 10.442226
## A2NL 2.598740 3.229218 3.027504 3.937428 5.317617 2.951300 2.378460  5.918622
## A6VY 2.810454 8.711676 4.414826 4.478048 4.800104 4.815597 3.385135  8.851152
## A3XT 2.157077 7.526234 5.585107 3.635403 4.951715 2.546328 2.020964  7.828502
## A1F0 3.144542 6.116264 3.905538 3.698795 6.402994 4.935074 2.650023  7.678229
## A26I 2.538364 9.076867 3.112930 3.175534 6.000657 3.148044 2.865361  9.028143
##          STC2    DTWD2     TCF4    ITPR2     DPYD     NME1    EGLN3    CD302
## A54N 5.090524 2.093831 2.258579 7.709530 2.851926 7.127238 2.247415 4.851567
## A2NL 3.835123 2.402509 6.163443 6.473752 4.428034 6.485592 2.329127 5.241001
## A6VY 5.445937 1.896778 5.202845 5.959788 3.561982 6.614160 7.548588 3.203079
## A3XT 6.115187 1.932672 4.650542 5.576878 3.970075 8.146036 2.489692 4.815868
## A1F0 4.464847 2.641886 7.292928 7.587744 5.353670 5.743227 4.639861 6.458236
## A26I 3.871246 2.714235 6.628540 5.179825 4.150577 6.617857 5.468940 4.684787
##           AHR  LAPTM4B     OCLN HIST1H2BK   HDAC11  C18orf1 C6orf192    AMPD3
## A54N 4.019673 8.049912 2.791720  5.954067 4.947815 4.683653 2.923813 2.895486
## A2NL 4.923811 7.746536 2.392937  4.525423 4.273616 5.951595 6.133223 3.861387
## A6VY 6.607690 8.417615 3.312478  6.275678 4.291192 4.992812 5.040908 3.561982
## A3XT 6.030869 8.407960 3.305845  4.498916 6.216512 4.473064 4.739108 4.464693
## A1F0 7.634721 6.714357 4.819566  6.735704 3.816752 6.439253 5.924120 4.921692
## A26I 8.068331 7.587261 4.767577  7.131893 3.062287 5.046603 5.029129 3.561458
##         COL6A1  RAB3IL1  APBB1IP    PSIP1  EIF2AK2    CSRP2 EIF4EBP3      LYN
## A54N  8.026430 4.227116 2.494047 6.812241 4.367785 5.461665 4.660021 5.068491
## A2NL  7.995400 3.770066 3.091034 7.122138 5.020992 7.722498 3.846668 4.687849
## A6VY 10.398104 4.796434 3.589811 5.364627 5.157790 5.538921 2.583874 6.027099
## A3XT 11.303807 6.043572 4.609365 7.067148 4.780168 2.746059 4.321395 7.170656
## A1F0  8.803016 4.614440 5.433238 6.796996 6.377081 4.036504 2.708679 7.999637
## A26I  9.912655 3.605914 3.994696 6.091547 6.910860 4.331515 2.375975 5.865131
##         WDR76   SAMD9L     ASPH     RBL1  SLC43A3      HN1   TTC39A     MTL5
## A54N 1.935903 4.486484 6.581808 1.013491 5.745731 8.117577 5.846276 3.550858
## A2NL 4.656194 5.081083 7.438594 3.634496 7.545353 7.451381 3.389409 5.758744
## A6VY 4.091535 5.672553 7.608025 3.200304 7.340047 8.425730 4.141570 3.506906
## A3XT 4.733892 5.623687 6.749278 4.274246 8.550059 7.801542 4.542326 6.344360
## A1F0 3.741313 7.727129 7.599448 4.565084 7.137421 6.728517 4.868629 3.897585
## A26I 3.545743 6.712845 8.164311 4.068645 5.862078 7.776145 5.133854 2.202080
##           NES      APOD     RIN3    ALCAM  C1orf38    PLCD3    BSPRY     NTN4
## A54N 8.850995 -3.968811 4.808540 5.162122 5.476083 6.276557 4.045684 4.445552
## A2NL 4.155471  3.350380 4.076955 6.306678 3.511633 3.302618 4.087377 2.908277
## A6VY 8.891543  4.520934 5.055567 6.245715 5.316111 3.531399 5.031572 6.124716
## A3XT 8.883476  5.742715 5.052913 5.658512 6.169961 2.800291 6.605793 3.317514
## A1F0 6.933829  6.806887 5.267253 5.848047 6.117976 3.815547 6.411397 5.857742
## A26I 8.248330  3.047485 4.281791 7.641658 4.294383 3.470061 5.947499 3.386272
##         IL1R1     EMP3  ZKSCAN1    FMNL2   OGFRL1     IRF5    IGSF3      DBP
## A54N 4.702169 2.844536 3.432628 5.735780 5.894484 4.664896 6.054362 5.474181
## A2NL 4.973621 3.849190 5.106193 5.473317 4.287888 2.986727 6.133941 3.604084
## A6VY 6.181392 5.350305 4.943887 5.607545 5.486930 4.845207 6.821853 3.576584
## A3XT 5.773118 5.920424 5.624157 6.157679 5.570065 4.595054 6.496858 3.770161
## A1F0 7.433241 5.495157 6.058077 7.321446 6.164404 5.046134 7.349827 3.424478
## A26I 5.927091 4.403997 6.821424 6.674022 4.694242 4.597634 7.341074 1.429654
##          CNN2   CAMK2D   SIGIRR    AKAP9     ICA1      FGD5     DSG2     E2F1
## A54N 8.896371 5.356488 5.994692 5.081122 5.379681 0.6222589 6.447874 5.478472
## A2NL 5.985026 5.869000 4.261834 6.513767 4.916314 3.1028235 8.355314 6.367501
## A6VY 8.050750 5.540567 4.744049 6.218890 4.143015 4.0915346 8.071866 5.679045
## A3XT 5.903877 4.942701 5.753493 5.538005 4.295643 4.0782913 6.623286 6.224901
## A1F0 7.071984 6.437687 3.931365 7.036986 5.276482 5.4773546 7.488079 4.521392
## A26I 6.738799 5.138518 2.934033 7.400646 4.409796 4.2138134 8.658100 4.058838
##         QSOX1     TOB1    CSF3R  SHROOM3     CCDC80    FRMD6   CXCL12    CCNA2
## A54N 5.290038 5.438537 7.787710 6.872521 -0.7852309 3.476316 7.060861 4.763464
## A2NL 6.688729 4.206057 1.592055 5.286294  4.4524285 5.320144 5.256257 5.669763
## A6VY 8.396634 4.756381 2.471546 6.492201  6.6497392 5.657965 5.908436 6.116734
## A3XT 5.764192 6.716736 4.338654 5.304069  6.4339243 4.352332 6.473057 5.358724
## A1F0 7.078657 6.415172 4.089329 6.983514  8.3058455 6.093320 8.687017 5.287414
## A26I 7.385890 6.664888 2.778391 6.777249  8.0851881 6.317594 6.303906 5.746757
##         TIGD5  ALDH6A1     POSTN     FZD4   NCAPG2     SDC4    SNED1  PLEKHA4
## A54N 3.203428 2.354618 -3.311116 3.012211 3.591155 7.594471 3.082999 7.839741
## A2NL 4.626894 2.790783  7.254870 5.349716 6.356771 6.666061 2.722081 2.203355
## A6VY 5.551495 2.813544  9.419911 3.215335 5.756170 7.691160 3.948396 5.541115
## A3XT 4.125523 3.375489  8.737852 4.387566 5.521957 7.431553 3.595702 3.105035
## A1F0 2.723845 3.081031  8.743015 6.770625 5.691922 7.507996 5.128529 4.675513
## A26I 2.241631 3.528703 10.417133 4.868176 6.854575 6.480281 4.013859 5.762749
##        KCNAB2  SH3KBP1    IGSF9     DNLZ    S1PR3    PTPRE FLJ23867   PLSCR1
## A54N 2.603166 5.487989 6.660592 4.893349 2.859278 4.495948 7.114072 5.188629
## A2NL 4.060423 6.555923 6.450295 2.908277 2.679296 2.619120 5.263497 6.762821
## A6VY 4.853178 5.127438 6.946641 4.503032 5.009377 4.374781 6.642603 6.721406
## A3XT 5.624627 7.623331 6.527073 3.265449 4.303844 4.670699 1.652142 5.858269
## A1F0 4.888256 6.651397 5.265489 1.850072 5.783219 5.257965 5.974405 7.386008
## A26I 3.284747 5.953797 5.865834 2.011856 6.747125 5.076682 6.291397 6.384149
##          LMO4   IFITM2   LRRC25      TST     NCF4    NCOA7     IL4R   CCDC64B
## A54N 6.830591 4.802831 2.866592 4.919960 2.611897 4.528594 4.493588 6.5496541
## A2NL 8.957832 4.925471 2.035116 4.529796 3.289782 7.033673 3.866583 0.9096316
## A6VY 5.905457 6.685397 3.260148 4.895819 4.196904 5.950736 7.091513 5.2599236
## A3XT 8.490427 6.041111 4.208083 5.106361 3.683089 7.401442 5.713424 4.1802323
## A1F0 7.574977 7.402661 4.107173 3.786310 4.423509 7.933876 5.902142 4.0833318
## A26I 6.385376 5.609940 3.124729 3.515587 2.630634 5.503245 5.701427 1.9495285
##          SGPP1    RUNX3   SLC5A6    IFIH1    PREX1    PLAUR    CDK18  SLC43A2
## A54N 0.5514035 1.893564 9.111073 4.685722 5.825782 5.129050 5.182780 6.310476
## A2NL 3.3773258 5.690065 6.650079 5.229619 4.526517 3.559477 5.475586 3.587019
## A6VY 3.4728220 4.567972 7.771439 5.733319 5.332676 7.023484 4.567972 5.522909
## A3XT 4.2920990 4.517247 8.363555 6.759564 7.636201 6.543810 6.704858 6.085807
## A1F0 5.3685420 5.221568 7.951124 8.289767 6.815355 5.813140 4.915517 5.359058
## A26I 4.8629514 2.175097 9.251806 6.664150 4.989667 5.766896 2.138316 4.052674
##            GK    ICAM2    YPEL2     CBR1    MEX3A    ZNRF3    PTPRM C1orf162
## A54N 3.838108 3.849929 4.662796 5.949769 7.990682 1.864631 2.736897 3.495315
## A2NL 3.307721 5.097520 3.740183 6.384924 7.466574 3.478079 9.653613 2.832414
## A6VY 3.063480 3.737198 3.540203 5.830157 6.929857 3.249450 4.807416 3.484273
## A3XT 3.467499 3.916880 3.578347 6.176499 5.308169 3.652088 5.807867 4.843686
## A1F0 3.982746 4.821417 4.997531 5.767175 6.196013 3.477345 5.899585 5.375921
## A26I 4.116707 2.537165 4.328457 5.926043 7.231704 4.023976 5.643138 3.098640
##          GAS6     C1QB    PVRL4     CTSK     MRVI1      LEF1    PLCD4   ZNF37B
## A54N 2.484571 6.782579 8.007943 3.282119 4.9852572 -1.810061 2.540521 4.415944
## A2NL 5.299802 7.049180 1.414232 5.382532 0.5622033  4.396051 1.549611 4.554690
## A6VY 6.504344 7.922148 7.294104 8.126331 3.6129469  3.304748 1.930654 3.639820
## A3XT 7.307726 8.387484 6.456792 6.234166 1.8579275  2.987238 1.376781 3.610955
## A1F0 6.695641 8.342243 6.866889 8.028861 4.8973999  3.763987 1.840621 3.896445
## A26I 7.818992 6.177899 7.261003 7.730926 6.1239500  3.186834 1.210112 4.063750
##         MEGF9    GINS2   FAM13A    CPT1A    SNX10   TRIM45     ELP2    ALOX5
## A54N 3.729970 4.285585 5.341106 5.809744 3.892868 4.229954 6.168993 2.637779
## A2NL 4.529796 5.361217 4.025388 5.667418 3.678492 5.471044 5.598826 2.884569
## A6VY 5.051724 5.335208 5.367227 4.826435 4.716155 3.472822 5.875767 5.288291
## A3XT 4.335214 3.236238 3.932930 7.148125 5.488282 3.102339 5.154650 5.387812
## A1F0 5.382855 2.531686 5.430922 6.939373 6.022736 2.068202 5.457399 5.199096
## A26I 6.612497 3.366530 4.093261 5.257982 4.985143 3.488444 6.105906 5.443726
##          AMN1   CERCAM   SEMA3C      KRT8 TP53INP2      JAM3   ZNF680     PBX1
## A54N 3.330531 6.084473 3.392652 10.514946 5.567567 0.8160598 2.688190 7.674013
## A2NL 4.155471 4.357203 4.816030  7.536919 5.079592 7.7647397 3.614046 7.767636
## A6VY 2.002495 7.240247 6.953737  9.746922 5.076901 3.5554824 3.194739 6.279956
## A3XT 2.136144 5.419019 4.475149  7.839072 3.684892 3.4463966 2.865336 6.242147
## A1F0 2.784246 5.490255 5.805874  7.588995 4.899108 4.9149537 4.175481 7.422295
## A26I 1.446603 6.378000 8.351354  9.417362 5.118002 4.7862714 4.606918 8.173260
class(datos)
## [1] "data.frame"

B.1.

Calcule la probabilidad de que una observación tenga NDRG2 mayor a 5.

paste("Cabecera de datos:") 
## [1] "Cabecera de datos:"
datos_NDRG2 <- datos$NDRG2
head(datos_NDRG2)
## [1] 7.519887 7.064522 5.978128 7.061261 7.660714 4.206059
paste("Tamaño muestal:") 
## [1] "Tamaño muestal:"
(n <- length(datos_NDRG2))
## [1] 70
intervalos=cut(datos_NDRG2,breaks=c(0,1,2,3,4,5,6,7,8,9))
frec_absoluta <- table(intervalos)
frec_relativa <- prop.table(frec_absoluta)
frec_acumu <- cumsum(frec_relativa)
(tablaAgrup <- cbind(frec_absoluta,frec_relativa,frec_acumu))
##       frec_absoluta frec_relativa frec_acumu
## (0,1]             0    0.00000000 0.00000000
## (1,2]             0    0.00000000 0.00000000
## (2,3]             2    0.02857143 0.02857143
## (3,4]             3    0.04285714 0.07142857
## (4,5]            12    0.17142857 0.24285714
## (5,6]            18    0.25714286 0.50000000
## (6,7]            18    0.25714286 0.75714286
## (7,8]            14    0.20000000 0.95714286
## (8,9]             3    0.04285714 1.00000000
data_frame_agrup <- as.data.frame(tablaAgrup)
(probabilidad <- data_frame_agrup[5,3])
## [1] 0.2428571
(prob_mayor_5 = 1-probabilidad)
## [1] 0.7571429

La probabilidad de que una observación tenga NDRG2 mayor a 5 es 0.7571429.

B.2.

Calcule la probabilidad de que una observación tenga NDRG2 y FAM63A mayor a 5.

Primero hay que calcular la probabilidad de que una observación tenga FAM63A mayor a 5:

datos_FAM63 <- datos$FAM63A
head(datos_FAM63)
## [1] 4.827405 4.329876 4.478048 2.671635 4.269159 6.068822
n <- length(datos_FAM63)

intervalos_FAM63 <- cut(datos_FAM63,breaks=c(0,1,2,3,4,5,6,7,8,9))
frec_absoluta_FAM63 <- table(intervalos_FAM63)
frec_relativa_FAM63 <- prop.table(frec_absoluta_FAM63)
frec_acumu_FAM63 <- cumsum(frec_relativa_FAM63)
(tablaAgrup_FAM <- cbind(frec_absoluta_FAM63,frec_relativa_FAM63,frec_acumu_FAM63))
##       frec_absoluta_FAM63 frec_relativa_FAM63 frec_acumu_FAM63
## (0,1]                   0          0.00000000       0.00000000
## (1,2]                   0          0.00000000       0.00000000
## (2,3]                   3          0.04285714       0.04285714
## (3,4]                   6          0.08571429       0.12857143
## (4,5]                  21          0.30000000       0.42857143
## (5,6]                  15          0.21428571       0.64285714
## (6,7]                  22          0.31428571       0.95714286
## (7,8]                   3          0.04285714       1.00000000
## (8,9]                   0          0.00000000       1.00000000
data_frame_agrup_FAM <- as.data.frame(tablaAgrup_FAM)
(probabilidad_FAM <- data_frame_agrup_FAM[5,3])
## [1] 0.4285714
(prob_mayor_5_FAM = 1-probabilidad_FAM )
## [1] 0.5714286

La probabilidad de que una observación tenga FAM63 mayor a 5 es 0.5714286.

La probabilidad de que en una observación NDRG2 y FAM63A sea mayor a 5 es el producto de ambas probabilidades:

(prob_NDRG2_FAM63A_mayor_5 <- prob_mayor_5 * prob_mayor_5_FAM)
## [1] 0.4326531

Dicha probabilidad es 0.4326531

B.3.

Para una observación con FAM63A mayor a 5, calcule la probabilidad de que presente GMDS entre 3 y 5.

Para realizar este apartado es necesario dividir las observaciones FAM63A y GMDS en los intervalos convenientes.

Posteriormente, con la función prop.table podemos hallar la probabilidad de que GMDS esté entre 3 y 5 para una observación con FAM63A mayor a 5:

datos_GMDS <- datos$GMDS
head(datos_GMDS)
## [1] 5.867395 6.775768 5.043232 5.063346 3.847745 4.912124
intervalos_GMDS <- cut(datos_GMDS,breaks=c(0,3,5,10))
frec_absoluta_GMDS <- table(intervalos_GMDS)
frec_absoluta_GMDS
## intervalos_GMDS
##  (0,3]  (3,5] (5,10] 
##      4     42     24
datos_FAM63
##  [1] 4.827405 4.329876 4.478048 2.671635 4.269159 6.068822 4.846648 4.429770
##  [9] 4.215061 4.875190 4.590524 4.042204 6.216048 3.839578 4.226852 4.612914
## [17] 5.013307 4.880497 3.975930 5.249345 3.387204 4.745601 4.715738 5.158879
## [25] 4.823377 5.957478 4.448172 2.820693 3.159309 3.854134 3.505493 4.200345
## [33] 5.705105 2.520773 4.539731 6.322894 6.350165 6.742725 6.441197 6.025760
## [41] 6.536516 6.496451 5.993788 6.757154 6.096949 7.045815 6.124112 5.801687
## [49] 6.212728 5.983151 5.603054 4.791595 5.907861 6.053649 6.271782 7.437277
## [57] 5.418058 4.279430 6.722430 7.195024 6.225689 6.660829 5.823610 5.690087
## [65] 6.694670 5.636068 6.943967 6.578923 6.260935 5.124423
intervalos_FAM63 <- cut(datos_FAM63,breaks=c(0,5,10))
frec_absoluta_FAM63 <- table(intervalos_FAM63)
frec_absoluta_FAM63
## intervalos_FAM63
##  (0,5] (5,10] 
##     30     40
(tabla_condicionada = prop.table(table(intervalos_FAM63,intervalos_GMDS),1))
##                 intervalos_GMDS
## intervalos_FAM63 (0,3] (3,5] (5,10]
##           (0,5]   0.00  0.40   0.60
##           (5,10]  0.10  0.75   0.15

La condición inicial es FAM63A sea mayor que 5, fijamos la segunda fila. Luego, observamos la columna de interés que en nuestro caso es GMDS entre 3 y 5, que es la segunda columna.

La probabilidad es:

tabla_condicionada[2,2]
## [1] 0.75

B.4.

Sin utilizar ninguna función de R, a partir de los resultados de apartados anteriores, calcula la probabilidad de que una observación cualquiera presente GMDS entre 3 y 5.

(tabla_no_condicionada = prop.table(table(intervalos_FAM63,intervalos_GMDS)))
##                 intervalos_GMDS
## intervalos_FAM63      (0,3]      (3,5]     (5,10]
##           (0,5]  0.00000000 0.17142857 0.25714286
##           (5,10] 0.05714286 0.42857143 0.08571429

La probabilidad de que cualquier observación presente GMDS entre 3 y 5 puede calcularse realizando el sumatorio de las frecuencias de GMDS (3,5] obtenido en las tabla de probabilidades no condicionadas entre dos observaciones.

paste("Probabilidad de que cualquier observación presente GMDS entre 3 y 5:")
## [1] "Probabilidad de que cualquier observación presente GMDS entre 3 y 5:"
(tabla_no_condicionada[1,2]) + (tabla_no_condicionada[2,2])
## [1] 0.6

B.5.

Calcule y razone si son independientes los sucesos “tener NDRG2 mayor a 5” y “GMDS entre 3 y 5”.

Si dos sucesos son independientes, el producto de ambas probabilidades es igual a la probabilidad de que ocurran ambos sucesos a la vez, es decir la probabilidad de que ocurra uno de ellos no depende de que haya ocurrido el otro.

Vamos a calcular el producto de ambas probabilidades por separado:

paste("Probabilidad de observación NDRG2 mayor a 5")
## [1] "Probabilidad de observación NDRG2 mayor a 5"
(prob_NDRG2_mayor_5 = prob_mayor_5)
## [1] 0.7571429
frec_absoluta_GMDS <- table(intervalos_GMDS)
(tabla_relativa <- prop.table(frec_absoluta_GMDS))
## intervalos_GMDS
##      (0,3]      (3,5]     (5,10] 
## 0.05714286 0.60000000 0.34285714
paste("Probabilidad de observación GMDS entre 3 y 5")
## [1] "Probabilidad de observación GMDS entre 3 y 5"
(prob_GMDS_3_y_5 = tabla_relativa[2])
## (3,5] 
##   0.6
paste("Probabilidad de que se den ambas observaciones:")
## [1] "Probabilidad de que se den ambas observaciones:"
(prob_ambos = prob_NDRG2_mayor_5*prob_GMDS_3_y_5)
##     (3,5] 
## 0.4542857

Ahora calculamos la probabilidad NDRG2 mayor a 5 unión GMDS entre 3 y 5.

intervalos_NDRG2=cut(datos_NDRG2,breaks=c(0,5,10))

(tabla_condicionada = prop.table(table(intervalos_NDRG2,intervalos_GMDS),1))
##                 intervalos_GMDS
## intervalos_NDRG2      (0,3]      (3,5]     (5,10]
##           (0,5]  0.05882353 0.76470588 0.17647059
##           (5,10] 0.05660377 0.54716981 0.39622642
(tabla_condicionada_2 = prop.table(table(intervalos_NDRG2,intervalos_GMDS),2))
##                 intervalos_GMDS
## intervalos_NDRG2     (0,3]     (3,5]    (5,10]
##           (0,5]  0.2500000 0.3095238 0.1250000
##           (5,10] 0.7500000 0.6904762 0.8750000
valor_1 <- tabla_condicionada[1,2]
valor_2 <- tabla_condicionada_2[2,2]

Si ambas probabilidades son iguales los sucesos son independientes:

prob_ambos == valor_1
## (3,5] 
## FALSE
prob_ambos == valor_2
## (3,5] 
## FALSE

Los sucesos no son independientes.

B.6.

Calcula el nivel de expresión medio de RBL1 y su varianza.

datos_RBL1 <- datos$RBL1
summary(datos_RBL1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -2.347   2.290   3.173   2.941   3.804   5.245
media_RBL1 <- mean(datos_RBL1)
media_RBL1
## [1] 2.941069
(n=length(datos_RBL1)) 
## [1] 70
cuasivar = var(datos_RBL1)
cuasivar
## [1] 1.832505
varianza  = ((n-1)/n)*(cuasivar)
varianza
## [1] 1.806326

Ejercicio C.1

1.

Para estudiar la regulación hormonal de una línea metabólica, se inyectan ratas albinas con un fármaco que inhibe la síntesis de proteínas del organismo. En general, tres de cada veinte ratas mueren a causa del fármaco antes de que el experimento haya concluido. Si se trata a 25 animales con el fármaco, ¿cuál es la probabilidad de que mueran al menos 5 en el experimento?

Estamos ante una distribución binomial X~B(25,3/20) Calcular P[X>=5].

n = 25
p = 3/20
1-pbinom(4,n,p)
## [1] 0.317893

2.

El tiempo necesario para llegar a la máxima concentración plasmática tras la administración de un fármaco sigue una distribución Exponencial de media 40 minutos. Calcule la probabilidad de que el fármaco tarde más de 50 minutos en alcanzar la máxima concentración plasmática.

Estamos ante una distribución exponencial X~Exp(l). E[X] = 1/l = 40. Prob[X>50].

media = 40
(l = 1/media)
## [1] 0.025
paste("Probabilidad de que el fármaco tarde más de 50 minutos en alcanzar la máxima concentración plasmática:")
## [1] "Probabilidad de que el fármaco tarde más de 50 minutos en alcanzar la máxima concentración plasmática:"
1-pexp(50,l)
## [1] 0.2865048

3.

La glucemia basal en individuos sanos de una población se considera que sigue una distribución N(90, 100).

  • a.Calcule la probabilidad de que la glucemia basal se encuentre comprendida entre 95 y 110.
pnorm(110,90,100)-pnorm(95,90,100)
## [1] 0.0593209
    1. Si se seleccionan muestras de tamaño 50. Probabilidad de que al media de la muestra sea inferior a 92.
n = 50
media_b = 90
varianza_nueva = 100/sqrt(n)
X = (92-media_b)/varianza_nueva
X
## [1] 0.1414214

P(Z<0.1414214)

pnorm(X,90,100)
## [1] 0.1844367

Ejercicio C.2

El dataset “liver.toxicity” del paquete “mixOmics” de Bioconductor contiene los niveles de expresión de 3116 genes y 10 medidas clínicas de 64 ratas que fueron expuestas a dosis no tóxicas, moderadamente tóxicas o severamente tóxicas de acetaminofén en un experimento controlado. Se supone que los individuos seleccionados constituyen una muestra aleatoria simple de la población de la que proceden.

Carga de datos:

library(mixOmics)
data("liver.toxicity")
datos_gen_A43 <- liver.toxicity[["gene"]][["A_43_P22018"]]
datos_gen_A43
##  [1]  0.00545  0.39802 -0.08280  0.00584  0.06808 -0.02287 -0.05686  0.01422
##  [9]  0.06941 -0.01738  0.34651  0.05166 -0.10649  0.07916  0.15065  0.12148
## [17]  0.15606  0.17269 -0.31992 -0.17538 -0.11689 -0.08529 -0.07735 -0.06310
## [25]  0.14567  0.10312  0.08260  0.10300 -0.14515  0.03420  0.28165 -0.05300
## [33]  0.26445 -0.48427 -0.04649  0.06660 -0.05235  0.17775  0.23898  0.03891
## [41]  0.05064  0.08118  0.09198  0.27578 -0.19264  0.05871  0.05190  0.05732
## [49] -0.20335  0.05094 -0.06902 -0.34582 -0.04918  0.35258  0.00235 -0.28462
## [57] -0.01347  0.08641 -0.07130  0.22619  0.13432  0.39030 -0.03404 -0.06248
class(datos_gen_A43)
## [1] "numeric"

1.

Utilizando estimadores insesgados, proporcione una estimación puntual para la media y la varianza del nivel de expresión del gen “A_43_P22018”.

Escogemos una muestra de los datos y calculamos la media y cuasivarianza muestral:

set.seed(123)
muestra <- sample(datos_gen_A43,10,replace=FALSE)
muestra
##  [1]  0.28165  0.15065 -0.06902  0.07916 -0.08280  0.08118  0.05094  0.35258
##  [9]  0.09198 -0.05235
(media_muestra <- mean(muestra))
## [1] 0.088397
(cuasivarianza <- var(muestra))
## [1] 0.02071848

Como estimadores puntuales se han obtenido una media muestral de 0.088397 y una cuasivarianza de 0.02071848

2.

Construya un intervalo de confianza al 99% para el nivel medio de albúmina en sangre (variable “ALB.g.dL.” de la componente “clinic”).

En primer lugar realizamos la selección de datos:

datos_albumina <- liver.toxicity[["clinic"]][["ALB.g.dL."]]
head(datos_albumina)
## [1] 4.9 5.0 5.0 5.0 5.1 5.1

Comando para intervalo de confianza:

t.test(datos_albumina, conf.level = 0.99)
## 
##  One Sample t-test
## 
## data:  datos_albumina
## t = 153.94, df = 63, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
##  4.927550 5.100575
## sample estimates:
## mean of x 
##  5.014062
t.test(datos_albumina, conf.level = 0.99)$conf.int
## [1] 4.927550 5.100575
## attr(,"conf.level")
## [1] 0.99

3.

Realice un contraste de hipótesis para determinar si el nivel medio de expresión del gen “A_43_P12995” es mayor que 0.

Selección de datos:

datos_gen_P12995 <- liver.toxicity[["gene"]][["A_43_P12995"]]
head(datos_gen_P12995)
## [1]  0.00676  0.08869 -0.07799 -0.02446 -0.01130 -0.06510

Función para contraste de hipótesis:

t.test(datos_gen_P12995,alternative="less",mu=0)
## 
##  One Sample t-test
## 
## data:  datos_gen_P12995
## t = -3.8079, df = 63, p-value = 0.0001601
## alternative hypothesis: true mean is less than 0
## 95 percent confidence interval:
##         -Inf -0.01846784
## sample estimates:
##   mean of x 
## -0.03288469

El p-valor obtenido es 0.0001601, un valor menor al de alpha,que por defecto es 0.05, en este caso, rechazamos la hipótesis nula, con lo cual el nivel medio de la expresión del gen “A_43_P12995” no es mayor que 0.

4.

Realice un contraste de hipótesis para determinar si el nivel de expresión medio de los genes “A_43_P17952” y “A_43_P21075” es el mismo.

En primer lugar seleccionamos los datos:

datos_gen_P17952 <- liver.toxicity[["gene"]][["A_43_P17952"]]
head(datos_gen_P12995)
## [1]  0.00676  0.08869 -0.07799 -0.02446 -0.01130 -0.06510
datos_gen_P21075 <- liver.toxicity[["gene"]][["A_43_P21075"]]
head(datos_gen_P21075)
## [1] 0.05093 0.07729 0.04081 0.02629 0.00059 0.01642

Antes de aplicar contrastes sobre la media de dos poblaciones hay que averiguar si las varianzas son iguales o no.

var.test(datos_gen_P17952,datos_gen_P21075,ratio = 1,alternative="two.sided")
## 
##  F test to compare two variances
## 
## data:  datos_gen_P17952 and datos_gen_P21075
## F = 2.1246, num df = 63, denom df = 63, p-value = 0.003231
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.290763 3.497177
## sample estimates:
## ratio of variances 
##           2.124624
var.test(datos_gen_P17952,datos_gen_P21075,ratio = 1,alternative="two.sided")$conf.int
## [1] 1.290763 3.497177
## attr(,"conf.level")
## [1] 0.95

El valor 1 no pertenece al intervalo de confianza además de que el p-valor es menor que alpha, con lo cual las varianzas de ambos conjuntos de datos son distintas.

Ahora comprobamos el constraste sobre la media de las dos poblaciones:

t.test(datos_gen_P17952,datos_gen_P21075,var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  datos_gen_P17952 and datos_gen_P21075
## t = -1.2655, df = 111.55, p-value = 0.2083
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.034487691  0.007603941
## sample estimates:
##    mean of x    mean of y 
## -0.005810781  0.007631094
t.test(datos_gen_P17952,datos_gen_P21075,var.equal = FALSE)$conf.int
## [1] -0.034487691  0.007603941
## attr(,"conf.level")
## [1] 0.95

A un nivel de confianza del 95%, se obtiene un intervalo donde se encuentra el valor 0, además de que el p-valor=0.2083 es mayor al valor de alpha=0.05.

Podemos concluir que no existen evidencias suficientes para rechazar la hipótesis nula, de modo que el nivel de expresión medio de los genes “A_43_P17952” y “A_43_P21075” se considera el mismo.

Ejercicio C.3

1.

Se quiere probar la efectividad de un antitérmico (reductor de la fiebre) en unos animales. Con tal fin se tomó la temperatura a 10 mamíferos afectados de una cierta enfermedad, antes y después de la administración del antitérmico resultando una diferencia de medias de 1.58 (X - Y) y una cuasidesviación típica de la diferencia de 1.71. Suponiendo condiciones de normalidad, estudie si el antitérmico es efectivo, a un nivel de confianza del 90 %. ¿La conclusión sería la misma a un nivel de confianza del 99 %?

  • muD = muX-muY = 1.58
  • Cuasidesv(D) = 1.71
  • Condiciones de normalidad
  • Estudiar si el fármaco es efectivo, a nivel de confianza del 90%

Para determinar si el antitérmico es efectivo o no, realizamos un contraste de hipótesis sobre la diferencia de medias, la H0 formulada será si la diferencia de medias es 0.

\[ H0:\mu 1-\mu2=0 \\ H1: \mu 1-\mu2 \neq 0\]

alpha = 0.1
n = 10
difer_medias = 1.58
cuasidesv = 1.71
grados_libertad = n-1

Calculamos el estadístico:

estadistico = difer_medias/(cuasidesv/sqrt(n))
estadistico
## [1] 2.921871

Determinar región crítica:

valor_critico = qt(alpha/2, df = grados_libertad, lower.tail = FALSE)
valor_critico
## [1] 1.833113
if (estadistico>valor_critico){
  cat("Rechazamos la hipótesis nula, el antitérmico no es efectivo")
} else {
  "Aceptamos"
}
## Rechazamos la hipótesis nula, el antitérmico no es efectivo
  • Parte 2: nivel confianza 99%:
alpha = 0.01
n = 10
difer_medias = 1.58
cuasidesv = 1.71
grados_libertad = n-1

Calculamos el estadístico:

estadistico = difer_medias/(cuasidesv/sqrt(n))
estadistico
## [1] 2.921871

Determinar región crítica:

valor_critico = qt(alpha/2, df = grados_libertad, lower.tail = FALSE)
valor_critico
## [1] 3.249836
if (estadistico>valor_critico){
  cat("Rechazamos la hipótesis nula, el antitérmico no es efectivo")
} else {
  "Aceptamos la hipótesis nula, el antitérmico es efectivo"
}
## [1] "Aceptamos la hipótesis nula, el antitérmico es efectivo"

Las conclusiones son distintas al cambiar el nivel de confianza del 90 al 99%

2.

En un ensayo clínico se han probado cinco medicamentos diferentes A, B, C, D y E sobre un grupo de pacientes afectados de una cierta enfermedad; el tamaño de cada una de las muestras es 51,54,48,49 y 48 respectivamente. Observándose una mejoría en 12 pacientes para el medicamento A, 8 para el medicamento B, 10 para el C, 15 para del D y 5 para el E. Contrastar si los 5 medicamentos son igualmente efectivos con un nivel de significación del 5%. Indicar el contraste utilizado.

  • Ensayo clínico con 5 medicamentos diferentes:
  • Medicamento A: 51 pacientes y 12 mejorías.
  • Medicamento B: 54 pacientes y 8 mejorías.
  • Medicamento C: 48 pacientes y 10 mejorías.
  • Medicamento D: 49 pacientes y 15 mejorías.
  • Medicamento E: 48 pacientes y 5 mejorías.

Contrastar si los 5 medicamentos son igualmente efectivos con nivel de significación del 5%.

En primer lugar proponemos los datos y la tabla de frecuencias

datos_medicamentos <- data.frame(Pacientes = c(51,54,48,49,48),mejora = c(12,8,10,15,5), no_mejora = c(51-12,54-8,48-10,49-15,48-5))

rownames(datos_medicamentos) = c("Medicamento A","Medicamento B","Medicamento C","Medicamento D","Medicamento E")

datos_medicamentos
##               Pacientes mejora no_mejora
## Medicamento A        51     12        39
## Medicamento B        54      8        46
## Medicamento C        48     10        38
## Medicamento D        49     15        34
## Medicamento E        48      5        43
class(datos_medicamentos)
## [1] "data.frame"

Para comprobar si los medicamentos son efectivos realizamos un test CHI-cuadtado de homogeneidad:

TablaFrec = matrix(c(12,8,10,15,5,39,46,38,34,43),nrow=5,
dimnames=list(c("Medicamento A","Medicamento B", "Medicamento C", "Medicamento D", "Medicamento E"),c("Mejoran", "No mejoran")))
TablaFrec
##               Mejoran No mejoran
## Medicamento A      12         39
## Medicamento B       8         46
## Medicamento C      10         38
## Medicamento D      15         34
## Medicamento E       5         43
chisq.test(TablaFrec)
## 
##  Pearson's Chi-squared test
## 
## data:  TablaFrec
## X-squared = 7.5295, df = 4, p-value = 0.1104

Al ser el p-valor = 0.1104 mayor que alpha = 0.05, no existen evidencias suficientes para rechazar la hipótesis nula con lo cual podemos decir que los 5 medicamentos son igualmente efectivos.

Contraste utilizado: Contraste de igualdad de poblaciones en muestras independientes, Test Chi-cuadrado de independencia.