• Fenotipo: SM300-Efficiency
  • Modelo: Random Forest
  • Significancia: No
  • Datos: 20 SNP’s (Más importantes)
  • Finalidad: Arma progresivamente modelos, agregando genes con mayor importancia para el fenotipo (20 más importantes)

Preprocesamiento de datos

Librerias

library(readxl)
library(randomForest)
library(caret)
library(dplyr)
library(ggplot2)

Lectura de los datos

fenotipo <- "SM300_Efficiency"
fenotipos_path <- "C:/Users/ca_re/OneDrive/Desktop/Memoria/Desarrollo/fenotipos.xlsx"
lectura_datos <- read_excel(fenotipos_path, 
                            col_types = c(rep("text", 55), rep("numeric", 12)),
                            na = "NA")
top_path <- "C:/Users/ca_re/OneDrive/Desktop/Memoria/Desarrollo/Top20.xlsx"
datos_top <- read_excel(top_path, 
                            col_types = c(rep("text", 1), rep("numeric", 1)),
                            na = "NA")

Preparación de los datos

Se obtiene el valor para fenotipo SM300 Efficiency y valor para cada cepa en los 54 SNP’s. Se eliminan valores que no poseen información.

datos <- setNames(data.frame(lectura_datos[[fenotipo]],lectura_datos[, 2:55]), 
                  c(fenotipo, names(lectura_datos[, 2:55]))) #Se obtiene valor fenotipo y bases de los SNP's (genes)

datos <- na.omit(datos) # Se eliminan filas que no contengan datos
datos <- datos %>% mutate_at(vars(2:55), as.factor) # Se indican las columnas que son factor
datos

Se obtienen los nombres de los genes con más relevancia para el fenotipo SM300 Efficiency.

datos_top
nombres_genes <- datos_top$Nombre #Se almacenan los nombres de los genes que tienen importancia

Funciones Modelo de regresión

Converción de datos

Debido a que los valores de cada uno de los SNP’s son bases nitrogenadas se deben transformar los datos a variables dummy para poder ser utilizadas en el modelo de regresión.

# Descripción: Convierte las columnas categoricas de un dataframe a valores a dummy
# Argumentos: datos (df a convertir), nombre del fenotipo 
# Valor retorno: dataFrame (columnas del dataFrame convertidas a variables dummy)
convertir_dummy <- function(datos, nombre_fenotipo, nombre_gen){
  columnas <- sapply(datos, is.factor)# Seleccia columnas catégoricas
  df_dummies <- as.data.frame(model.matrix(~.-1, data = datos))
  return(df_dummies)
}

División datos

Al ser un conjunto de datos reducidos se utiliza dividen los datos en 70% entrenamiento y un 30% de prueba.

# Descripción: Toma un conjunto de datos y los divide en datos de entrenamiento y prueba
# Argumentos: datos(conjunto de datos a separar)
# Valor retorno: lista(con el conjunto de entrenamieto y prueba)
separacion_datos <- function(datos){
  set.seed(123)
  train_index <- sample(1:nrow(datos), 0.7 * nrow(datos))
  train <- datos[train_index, ]
  test <- datos[-train_index, ]
  return(list(train = train, test = test))
}

Entrenamiento

Se realiza un entrenamiento de un modelo Random Forest, tomando los datos de entrenamiento (70% de los datos totales) con una validación cruzada simple debido a la cantidad reducida de datos, evitando el sobreajuste o subajuste de los datos, a través de la creación de varias particiones de los datos, entrenando y probando el modelo con cada una de ellas, para esta ocasión se utilizan 10 particiones.

# Descripción: Entrena un modelo Random Forest sobre un conjunto de datos de entrenamiento
# Argumentos: train(conjunto de datos de entrenamiento), fenotipo (nombre del fenotipo, variable independiente)
# Valor retorno: modelo entrenado
entrenamiento <- function(train, fenotipo){
  model <- train(as.formula(paste(fenotipo, "~ .")), 
                 data = train, 
                 method = "rf", 
                 trControl = trainControl(method = "cv", number = 10),
                 tuneLength = 2)
  return(model)
}

Predicciones

# Descripción: Se predicen valores con el modelo anteriormente entrenado
# Argumentos: model (modelo entrenado), test (datos de prueba)
# Valor retorno: predicciones
predicciones <- function(model, test){
  predictions <- predict(model, newdata = test)
  return(predictions)
}

Evalaución

  • Coeficiente de determinación (R2): medida estádistica que refelja la relación entre las variables y qué tan bien predice los valores de la variable independiente. Se aproxima el resultado con 4 dígitos. Su ecuación es: \[R^2 = \frac{\sum_{t=1}^{T}(y_t - \hat{y_t})^2}{\sum_{t=1}^{T}(y_t - \bar{y})^2}\]
  • Error cuadrático medio (MSE): utilizado para evaluar la calidad de un modelo de regresión. Indicando cuanto se desvia el modelo generado de los datos observados. Su ecuación es: \[MSE = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y_i})^2\]
# Descripción: Se calcula r2 y mse
# Argumentos: predictions (predicciones), test (df de prueba), fenotipo (nombre del fenotipo)
# Valor retorno: lista (con valores test y train)
evaluacion_coeficientes <- function(predictions, test, fenotipo, nombre_modelo){
  # Coeficiente de determinación
  n = length(predictions)
  x_mean = mean(predictions)
  y_mean = mean(test[[fenotipo]])
  ssx = sum((predictions - x_mean)^2)
  ssy = sum((test[[fenotipo]] - y_mean)^2)
  sxy = sum((predictions - x_mean)*(test[[fenotipo]] - y_mean))
  r2 = round((sxy^2) / (ssx*ssy),5)

  # Error cuadrático medio (MSE)
  mse <- round(mean((test[[fenotipo]] - predictions)^2),5)
  
  # Se organizan los resultados
  resultados <- data.frame(Modelo = nombre_modelo, R2 = r2, MSE = mse)
  return(resultados)
}

Ejecución

Se crea un total de 20 modelos de regresión con Random Forest para el fenotipo SM300 Efficiency, agregando progresivamente los 20 genes que tuvieron mayor importancia en el modelo de regresión Random Forest que involucraron el total de SNP’s (54).

  datos_top <- setNames(data.frame(dato = datos[[fenotipo]]),fenotipo)
  resultados <- data.frame(matrix(ncol=5,nrow=0)) # Dataframe donde se guardaran los resultados
  nombres <- "" # Se almacenan los nombres de los SNP's a trabajar en cada iteración
  for(i in 1:length(nombres_genes)){
    # Se obtiene la información del gen 
    col_index <- match(nombres_genes[i], colnames(datos))
    dato_gen <- setNames(data.frame(dato = datos[,col_index]),nombres_genes[i])
    datos_top <- cbind(datos_top, dato_gen)
    
    # Se transforma la columna del SNP a variable dummy para ser trabajada en el modelo
    datos_transformados <- convertir_dummy(datos_top,fenotipo, colnames(datos)[i+1])
    
    datos_separados <- separacion_datos(datos_transformados)
    train <- datos_separados$train
    test <- datos_separados$test
    
    # Se realiza un modelo Random forest (con validación cruzada simple)
    model <- entrenamiento(train,fenotipo)
    
     # Se utilizan los datos de prueba
    predictions <- predicciones(model,test)
    
    # Se realiza la evaluación de las predicciones del modelo
    evaluacion <- evaluacion_coeficientes(predictions, test, fenotipo, i)
    
    # Se agregan al df los nombres de los genes con los que se entreno el modelo
    nombres <- paste(nombres, nombres_genes[i], sep = "|")
    evaluacion <- cbind(evaluacion,nombres) 
    
    # Se almacenan los resultados obtenidos del gen
    resultados <- rbind(resultados, evaluacion)
  }
resultados

Resultados

Para el modelo base (entrenado con 54 SNP’s), se tiene:

  • R²: 0.23622
  • MSE: 0.01403

Mientras que para el entrenamiento de cada gen por separado para el fenotipo SM300 Efficiency se tiene:

  • Modelo con mayor R²: 0.22245
  • Modelo con menor MSE: 0.01455
resultados
ggplot(resultados, aes(x = Modelo)) + 
  geom_line(aes(y = R2, group = 1, color = "R2"), linetype = "solid") +
  geom_line(aes(y = MSE * 10, group = 1, color = "MSE"), linetype = "solid") +
  scale_y_continuous(name = "R2", sec.axis = sec_axis(~./10, name = "MSE")) +
  labs(title = "Modelos Top 20", x = "Número modelo", y = "Valor", color = "Leyenda")