Fenotipo: SM300-Efficiency
Modelo: Random Forest
Significancia: No
Datos: 3 SNP’s
Finalidad: Realizar un grid search con la finalidad de mejorar los parámetros del modelo de regresión random forest entrenado anteriormente con 3 SNP’s.

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 <- head(datos_top$Nombre, 3) #Se almacenan los nombres de 3 primeros genes que tienen importancia para el fenotipo.
datos_modelo <- setNames(data.frame(dato = datos[[fenotipo]]),fenotipo)
# Se obtiene la información de los genes 
  for(i in 1:length(nombres_genes)){
    col_index <- match(nombres_genes[i], colnames(datos))
    dato_gen <- setNames(data.frame(dato = datos[,col_index]),nombres_genes[i])
    datos_modelo <- cbind(datos_modelo, dato_gen)
  }
datos_modelo

Modelo de regresión

Tansformación datos categoricos

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))
}

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 (ntree, mtry y nodesize)
# Valor retorno: dataframe (valores R2, MSE y valores para ntree, mtry y nodesize)
evaluacion_coeficientes <- function(predictions, test, fenotipo, n_tree, m_try, node_size){
  # 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(ntree = n_tree, mtry = m_try, nodesize= node_size, R2 = r2, MSE = mse)
  return(resultados)
}

Grid search

Con la finalidad de encontrar la mejor combinación de hiperparámetros para el algoritmo se utilzia grid search, si bien los hiperparámetros son valores que no se aprenden directamente del conjunto de datos, afectan el rendimiento y el comportamiento del modelo. En este caso de ven incolucrados los parametros ntree, mty y nodesize.

ntree: controla el número de árboles que se construirán en el bosque.
mtry: controla el número de variables predictoras (features) que se consideran al construir cada árbol de decisión en el bosque, indica el número máximo de variables que se seleccionan de forma aleatoria en cada partición de árbol.
nodesize: controla el tamaño mínimo de un nodo terminal en un árbol de decisión dentro del bosque.

# Descripción: entrena un modelo de regresión con random forest
# Argumentos: train (datos de entrenameinto), fenotipo (nombre del fenotipo),tune_control(contiene tipo de validación cruzada) y valores
# para ntree, mtry y nodesize
# Valor retorno: dataframe (valores R2, MSE y valores para ntree, mtry y nodesize)
entrenamiento_grid <- function(train, fenotipo, tune_control, n_tree, mtry, node_size){
  tune_grid <- expand.grid(.mtry=mtry)
  
  model <- train(as.formula(paste(fenotipo, "~ .")), 
                data = train, 
                method = "rf", 
                trControl = tune_control,
                tuneGrid = tune_grid,
                tuneLength = 2,
                ntree = n_tree,
                nodesize = node_size)
  return(model)
}

tune_control <- trainControl(method = "cv", number = 9) #Tipo de validación cruzada
resultados <- data.frame(matrix(ncol=5,nrow=0))

# ntree (Se realizan pruebas modificando el parámetro entre 1:20)
for (i in 1:20){
  # mtry (Se realizan pruebas modificando el parámetro entre 1:5)
  for (j in 1:5){
    # nodesize (Se realizan pruebas modificando el parámetro entre 1:10)
    for (k in 1:10){
      # Transformación variables categoricas 
      datos_transformados <- convertir_dummy(datos_modelo,fenotipo)
      datos_separados <- separacion_datos(datos_transformados)
      # Separación datos prueba | entrenameinto (30/70)
      train <- datos_separados$train
      test <- datos_separados$test
      # Se entrena el modelo
      model <- entrenamiento_grid(train,fenotipo,tune_control,i,j,k)
      # Se utlizan los datos de prueba
      predictions <- predicciones(model,test)
      # Se calcula el R2 y MSE
      evaluacion <- evaluacion_coeficientes(predictions,test,fenotipo,i,j,k)
      # Se alamcenan los resultados
      resultados <- rbind(resultados, evaluacion)
    }
  }
}

indice_mayorR2 <- which.max(resultados$R2)
cat("Mayor R2:", max(resultados$R2), ", ntree:", resultados$ntree[indice_mayorR2],
    ", mtry:", resultados$mtry[indice_mayorR2], ", nodesize:",resultados$nodesize[indice_mayorR2])

## Mayor R2: 0.31871 , ntree: 4 , mtry: 2 , nodesize: 9

indice_menorMSE <- which.min(resultados$MSE)
cat("\nMenor MSE:", max(resultados$MSE), ", ntree:", resultados$ntree[indice_menorMSE],
    ", mtry:", resultados$mtry[indice_menorMSE], ", nodesize:",resultados$nodesize[indice_menorMSE])

## 
## Menor MSE: 0.01584 , ntree: 4 , mtry: 2 , nodesize: 9

Resultados

	Modelo base(54 SNP’s)	Modelo con un gen con mayor R²	Modelo con un gen con menor MSE	Modelo con 3 SNP’s	Modelo con 3 SNP’s cpn parámetros ajustados
R²	0.23622	0.22245	0.22245	0.25972	0.31871
MSE	0.01403	0.01455	0.01455	0.01371	0.01285

Mejor modelo SM300-Efficiency Random Forest (Modificación parámetros)

Catalina Riquelme

2023-05-11