• Fenotipo: SM300-Efficiency
  • Modelo: Random Forest
  • Significancia: No
  • Datos: 54 SNP’s
  • Finalidad: contruir un modelo base con el total de SNP’s (54)

Preprocesamiento de datos

Librerias

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

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

Preparación de los datos

Se obtiene el valor para fenotipo SM300 Efficiency y valor para cada cepa en los 54 SNP’s. Ademá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

Modelo de regresión

Transformació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: dataFrame (datos a convertir), nombre del fenotipo 
# Valor retorno: dataFrame (columnas del dataFrame convertidas a variables dummy)
convertir_dummy <- function(datos, nombre_fenotipo){
  columnas <- sapply(datos, is.factor)# Seleccia columnas catégoricas
  df_dummies <- as.data.frame(model.matrix(~.-1, data = datos[, columnas]))
  df_final <- setNames(data.frame(datos[[nombre_fenotipo]], df_dummies),
                       c(nombre_fenotipo, names(df_dummies)))
  return(df_final)
}
datos_nuevos <- convertir_dummy(datos,fenotipo)
datos_nuevos

División datos

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

set.seed(123)
train_index <- sample(1:nrow(datos_nuevos), 0.7 * nrow(datos_nuevos))
train <- datos_nuevos[train_index, ]
test <- datos_nuevos[-train_index, ]

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.

model <- train(as.formula(paste(fenotipo, "~ .")), 
               data = train, 
               method = "rf", 
               trControl = trainControl(method = "cv", number = 10))

Evaluación

Predicciones

predictions <- predict(model, newdata = test)

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}\]

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))
r_squared = round((sxy^2) / (ssx*ssy),5)

print(paste("R²: ",r_squared))
## [1] "R²:  0.23622"

Coeficiente de Pearson

Cuantificar la relación lineal entre dos variables cuantitativas continuas.

pearson <- round(cor(test[[fenotipo]],predictions)^2,5)
print(paste("Pearson²: ",pearson))
## [1] "Pearson²:  0.23622"

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\]

mse <- round(mean((test[[fenotipo]] - predictions)^2),5)
print(paste("MSE: ",mse))
## [1] "MSE:  0.01403"

Resultados

# Importancia de las variables

importancia <- varImp(model)
importancia
## rf variable importance
## 
##   only 20 most important variables shown (out of 109)
## 
##          Overall
## YPL200Wt  100.00
## YOR222Wg   88.46
## YLR425Wt   87.74
## YBL052Cg   81.87
## YNR032Wg   76.71
## YPR025Cg   72.86
## YDL176Wg   71.88
## YGL222Cg   67.80
## YBR243Ct   66.76
## YJL099Wt   61.36
## YNL248Ct   59.03
## YAL042Wt   59.01
## YJL099Wk   57.53
## YGL127Cg   56.18
## YGR270Wt   55.09
## YDR006Cm   54.75
## YIL128Wg   54.03
## YOL098Ct   53.90
## YOR227Wt   51.54
## YBL060Wt   51.36
plot(importancia, main = "Importancia de variables", xlab = "Variables", ylab = "Importancia")