Comprehensive study of factors affecting sleep quality using supervised learning models (GLM, Random Forest, GBM)
El propósito de este trabajo es aplicar técnicas de modelización predictiva para resolver un problema realista, desde la preparación de datos hasta la evaluación de modelos, con el fin de comparar su rendimiento y justificar la elección final.
En primer lugar, se debía elegir una base de datos que tuviese una variable objetivo de tipo continuo, discreto, binario, nominal u ordinal (en nuestro caso ordinal, pero se explicará más en detalle a continuación)
Debe contener también al menos 12 variables predictoras y un mínimo de 1000 observaciones.
Antes de proceder con cualquier proceso de depuración o análisis, se cargan las distintas librerías que se necesitarán a lo largo del proyecto y los datos.
library(MASS)
library(stargazer)
library(caret)
library(car)
library(vcd)
library(purrr)
library(dplyr)
library(nnet)
library(lmtest)
library(rpart)
library(rpart.plot)
library(pROC)
library(rattle)
library(vcd)
library(partykit)
library(rpartScore)
library(randomForest)
library(OOBCurve)
library(furrr)
library(iml)
library(glmnet)
library(lubridate)
library(adabag)
library(gbm)
# Para paralelizar
library(doParallel)
cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)Una vez están cargadas todas las librerías, se cargan los datos.
datos<-read.csv("stress_detection_data.csv")La siguiente base de datos consiste en un estudio de individuos reales donde se tienen en cuenta sus hábitos de vida, salud e indicadores de estrés.
Cada observación representa a un individuo (se cuenta con 774 sujetos) donde se reflejan, entre otras cosas, sus consumiciones de café, comportamiento, patrones de sueño y variables de salud.
Cuenta con las siguientes 22 variables:
Age: Edad del individuo
Gender: Sexo del individuo (binaria)
Occupation: Profesión del individuo (toma 169 valores)
Marital_Status: Soltero, Casado o Divorciado
Sleep_Duration: Número de horas que duerme el individuo
Sleep_Quality: Calidad de sueño entre 1 y 5 que siente que ha experimentado el individuo
Wake_Up_Time: Hora a la que se despierta normalmente
Bed_Time: Hora a la que se acuesta normalmente
Physical_Activity: Frecuencia con la que hace deporte a la semana
Screen_Time: Horas diarias de uso de pantallas en media del individuo
Caffeine_Intake: Cantidad de cafeína consumida (tazas de café)
Alcohol_Intake: Cantidad de bebidas alcohólicas consumidas (bebidas por semana)
Smoking_Habit: Si fuma o no (binaria)
Work_Hours: Número de horas trabajadas
Travel_Time: Número de horas al día invertidas en transporte
Social_Interactions: Frecuencia de interacciones sociales (eventos sociales)
Meditation_Practice: Si el individuo medita o no (binaria)
Exercise_Type: Tipo de ejercicio que hace (Cardio, Yoga, Pilates, etc. 6 categorías)
Blood_Pressure: Niveles de presión en sangre
Cholesterol_Level: Nivel de colesterol (mg/dL).
Blood_Sugar_Level: Nivel de azúcar en sangre (mg/dL).
Stress_Detection: Nivel de estrés (Low, Medium o High)
Para el estudio se va a utilizar como variable objetivo Sleep_Quality en función de los distintos factores que muestra el individuo. Además, como son valores reales, podemos extrapolar los resultados a la población.
Antes de comenzar a trabajar en los distintos modelos, se deben revisar las distintas variables y valores de las mismas que toma la base de datos, pasar a factor aquellas que lo necesiten y reagrupar aquellas variables que tengan muy pocas observaciones por categoría.
Se procede, con la ayuda de la función str(), a ver cómo están declaradas las variables y, en caso de que sea necesario, recodificarlas.
str(datos)'data.frame': 773 obs. of 22 variables:
$ Age : int 30 35 40 35 29 45 32 37 50 40 ...
$ Gender : chr "Male" "Female" "Male" "Male" ...
$ Occupation : chr "Software Engineer" "Marketing Manager" "Data Scientist" "Software Engineer" ...
$ Marital_Status : chr "Single" "Married" "Divorced" "Single" ...
$ Sleep_Duration : num 7 6 7 7 8 6 7 7.5 5.5 6 ...
$ Sleep_Quality : num 4 3 4 4 5 3 4 4 2 3 ...
$ Wake_Up_Time : chr "7:00 AM" "6:00 AM" "7:00 AM" "7:00 AM" ...
$ Bed_Time : chr "10:00 PM" "11:00 PM" "10:00 PM" "10:00 PM" ...
$ Physical_Activity : num 2 1 2 2 3 2 1 2.5 2 3 ...
$ Screen_Time : num 4 3 4 4 2 4 6 3.5 8 4.5 ...
$ Caffeine_Intake : int 1 0 1 1 1 2 3 2 3 2 ...
$ Alcohol_Intake : int 0 1 0 0 0 1 0 0 2 0 ...
$ Smoking_Habit : chr "No" "No" "No" "No" ...
$ Work_Hours : int 8 9 8 8 7 10 9 8 12 9 ...
$ Travel_Time : num 1 2 1 1 1 2 0.5 1.5 2 1 ...
$ Social_Interactions: int 5 3 5 5 4 3 2 5 2 4 ...
$ Meditation_Practice: chr "Yes" "No" "Yes" "Yes" ...
$ Exercise_Type : chr "Cardio" "Yoga" "Strength Training" "Cardio" ...
$ Blood_Pressure : int 120 110 130 120 110 130 120 140 150 125 ...
$ Cholesterol_Level : int 180 160 200 180 180 220 190 210 240 200 ...
$ Blood_Sugar_Level : int 90 80 100 90 90 110 85 100 130 105 ...
$ Stress_Detection : chr "Low" "Medium" "High" "Low" ...De base vemos como se cuenta con 22 variables (21 predictoras y la variable objetivo Sleep Quality).
En primer lugar se deberán pasar a factor las variables Gender, Occupation, Marital_Status, Smoking_Habit, Meditation_Practice, Exercise_Type y Stress_Detection. En el caso de Sleep_Quality, también se pasará a factores pero primero debemos ver y reagrupar los valores que toma (puesto que hay valores con decimales).
También se deberán recodificar las variables Wake_Up_Time y Bed_Time para que tomen valores de 0 a 23 respectivamente para tratar las horas como una variable numérica.
A tener en cuenta que en esta segunda variable, en ocasiones toma valores 1:00 AM, que parece “temprano” pero en realidad es tarde ya que es el día siguiente. Para arreglarlo, se va a convertir esa variable a un horario continuo donde la noche se va a extender más allá de las 24h para poder tener en cuenta estos valores.
datos$Gender <- factor(
ifelse(datos$Gender == "Male", 1, 2),
levels = c(1, 2),
labels = c("Male", "Female"))
datos$Occupation <- as.factor(datos$Occupation)
datos$Marital_Status <- factor(
ifelse(datos$Marital_Status == "Single", 0,
ifelse(datos$Marital_Status == "Married", 1, 2)),
levels = c(0, 1, 2),
labels = c("Single", "Married", "Divorced")
)
#datos$Sleep_Quality <- as.factor(datos$Sleep_Quality)
datos$Wake_Up_Time <- hour(parse_date_time(datos$Wake_Up_Time, orders = c("h:M p", "H:M")))
datos$Bed_Time <- hour(parse_date_time(datos$Bed_Time, orders = c("h:M p", "H:M")))
# Ajustamos para que horas tempranas (0-5 am) sean consideradas como tarde
datos$Bed_Time <- ifelse(datos$Bed_Time < 6, datos$Bed_Time + 24, datos$Bed_Time)
datos$Smoking_Habit <- factor(
ifelse(datos$Smoking_Habit == "No", 1, 2),
levels = c(1, 2),
labels = c("No", "Yes"))
datos$Meditation_Practice <- factor(
ifelse(datos$Meditation_Practice == "No", 1, 2),
levels = c(1, 2),
labels = c("No", "Yes"))
datos$Exercise_Type <- as.factor(datos$Exercise_Type)
datos$Stress_Detection <- factor(
ifelse(datos$Stress_Detection == "Low", 1,
ifelse(datos$Stress_Detection == "Medium", 2, 3)),
levels = c(1, 2, 3),
labels = c("Low", "Medium", "High")
)Se confirma que todos los cambios se han hecho correctamente
str(datos)'data.frame': 773 obs. of 22 variables:
$ Age : int 30 35 40 35 29 45 32 37 50 40 ...
$ Gender : Factor w/ 2 levels "Male","Female": 1 2 1 1 2 1 2 1 1 2 ...
$ Occupation : Factor w/ 169 levels "Account Manager",..: 151 94 39 151 158 46 70 23 19 103 ...
$ Marital_Status : Factor w/ 3 levels "Single","Married",..: 1 2 3 1 1 2 1 2 2 3 ...
$ Sleep_Duration : num 7 6 7 7 8 6 7 7.5 5.5 6 ...
$ Sleep_Quality : num 4 3 4 4 5 3 4 4 2 3 ...
$ Wake_Up_Time : int 7 6 7 7 6 5 7 6 5 6 ...
$ Bed_Time : num 22 23 22 22 22 23 24 22 24 23 ...
$ Physical_Activity : num 2 1 2 2 3 2 1 2.5 2 3 ...
$ Screen_Time : num 4 3 4 4 2 4 6 3.5 8 4.5 ...
$ Caffeine_Intake : int 1 0 1 1 1 2 3 2 3 2 ...
$ Alcohol_Intake : int 0 1 0 0 0 1 0 0 2 0 ...
$ Smoking_Habit : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 2 1 1 2 1 ...
$ Work_Hours : int 8 9 8 8 7 10 9 8 12 9 ...
$ Travel_Time : num 1 2 1 1 1 2 0.5 1.5 2 1 ...
$ Social_Interactions: int 5 3 5 5 4 3 2 5 2 4 ...
$ Meditation_Practice: Factor w/ 2 levels "No","Yes": 2 1 2 2 2 1 1 2 1 2 ...
$ Exercise_Type : Factor w/ 7 levels "Aerobics","Cardio",..: 2 7 5 2 7 2 1 5 6 4 ...
$ Blood_Pressure : int 120 110 130 120 110 130 120 140 150 125 ...
$ Cholesterol_Level : int 180 160 200 180 180 220 190 210 240 200 ...
$ Blood_Sugar_Level : int 90 80 100 90 90 110 85 100 130 105 ...
$ Stress_Detection : Factor w/ 3 levels "Low","Medium",..: 1 2 3 1 1 3 2 2 3 2 ...Y se revisan los valores que toma la variable Bed_Time
histogram(datos$Bed_Time)
Al haber realizado la modificación en como se tratan las horas, se observa que efectivamente hay individuos que se van a la cama pasado media noche.
Primero y más importante, se deben comprobar los valores que toma la variable objetivo, ya que si desde el comienzo se encuentra muy desbalanceada, podemos tener muchos problemas de cara a futuro en la predicción.
table(datos$Sleep_Quality)/nrow(datos)
2 2.5 2.9 3 3.1 3.2
0.003880983 0.002587322 0.002587322 0.122897801 0.012936611 0.028460543
3.3 3.4 3.5 3.6 3.7 3.8
0.016817594 0.021992238 0.029754204 0.028460543 0.073738680 0.102199224
3.9 4 4.1 4.2 4.3 4.5
0.073738680 0.266494179 0.051746442 0.037516171 0.014230272 0.011642950
5
0.098318241 Observamos que toma muchos valores intermedios. Esto se debe a que es un indicador personal de cada individuo en función de lo que el siente que ha experimentado entre un rango de 1 y 5. Como hay muchos valores que cuentan con un porcentaje de individuos muy pequeño, se van a hacer nuevos rangos.
El primero tomará todos los valores por debajo o igual a 3.4 (Pobre), el segundo los valores desde 3.5 hasta 3.8 (Aceptable), aquellos que hayan calificado su descanso como 3.9 o 4 (Bueno) y un úlitmo rango que toma valores superiores a 4 hasta 5 (Excelente).
La partición se hace de esta manera porque, al ver la distribución de individuos original se observa como un 27% de los mismos han valorado su descanso con un 4, lo que haría que, en cualquier categoría donde se encontrase este valor estuviera muy desbalanceada.
datos$Sleep_Quality <- cut(
datos$Sleep_Quality,
breaks = c(-Inf, 3.4, 3.8, 4, Inf),
labels = c("Pobre", "Aceptable", "Buena", "Excelente")
)
datos$Sleep_Quality <- factor(datos$Sleep_Quality,
levels =c("Pobre", "Aceptable", "Buena", "Excelente"))
table(datos$Sleep_Quality)/nrow(datos)
Pobre Aceptable Buena Excelente
0.2121604 0.2341527 0.3402329 0.2134541 Con esta nueva partición se observa que las categorías están ligeramente desbalanceadas, pero no parece muy grave, entorno al 21% de los datos respectivamente se encuentran en Pobre y Excelente mientras que hay un 23% en Aceptable y el 34% restante se encuentra en Bueno
En caso de que esta partición no funcionase correctamente, quizás se reagruparían de otra manera.
La variable Occupation toma muchos valores, por ello, se van a estudiar dichos valores y reagruparlos
table(datos$Occupation)/nrow(datos)
Account Manager Accountant
0.001293661 0.014230272
Actor Advertising Executive
0.005174644 0.001293661
Advertising Manager Architect
0.001293661 0.021992238
Artist Baker
0.009055627 0.001293661
Bakery Owner Bank Manager
0.002587322 0.003880983
Banker Bartender
0.003880983 0.002587322
Biologist Blacksmith
0.001293661 0.003880983
Brand Manager Bus Driver
0.001293661 0.002587322
Business Analyst Business Consultant
0.010349288 0.005174644
Business Owner Carpenter
0.001293661 0.006468305
CEO Chef
0.001293661 0.033635188
Civil Engineer Civil Servant
0.015523933 0.001293661
Cleaner Cobbler
0.001293661 0.002587322
Construction Engineer Construction Manager
0.001293661 0.007761966
Construction Worker Consultant
0.018111255 0.003880983
Content Creator Content Strategist
0.002587322 0.001293661
Content Writer Copywriter
0.007761966 0.003880983
Courier Customer Support
0.001293661 0.003880983
Data Analyst Data Engineer
0.012936611 0.001293661
Data Scientist Database Administrator
0.012936611 0.001293661
Delivery Driver Dentist
0.001293661 0.003880983
Designer Developer
0.002587322 0.001293661
Digital Marketer Doctor
0.002587322 0.010349288
Driver Editor
0.005174644 0.002587322
Electrical Engineer Electrical Technician
0.003880983 0.001293661
Electrician Engineer
0.021992238 0.014230272
Entrepreneur Event Coordinator
0.003880983 0.001293661
Event Manager Event Planner
0.001293661 0.007761966
Executive Director Factory Worker
0.001293661 0.001293661
Farmer Fashion Designer
0.010349288 0.010349288
Financial Advisor Financial Analyst
0.002587322 0.009055627
Financial Planner Firefighter
0.001293661 0.001293661
Fisherwoman Fitness Instructor
0.002587322 0.001293661
Fitness Trainer Flower Seller
0.005174644 0.002587322
Freelancer Graphic Designer
0.006468305 0.023285899
Hair Stylist Handicrafts Maker
0.001293661 0.002587322
Healthcare Assistant HR Executive
0.001293661 0.002587322
HR Manager HR Specialist
0.016817594 0.010349288
Human Resources Human Resources Manager
0.003880983 0.001293661
Insurance Agent Interior Designer
0.001293661 0.009055627
IT Consultant IT Manager
0.007761966 0.001293661
IT Specialist IT Support
0.002587322 0.001293661
IT Support Specialist Janitor
0.001293661 0.001293661
Journalist Laboratory Technician
0.020698577 0.001293661
Lawyer Librarian
0.009055627 0.009055627
Manager Marketing Director
0.001293661 0.001293661
Marketing Executive Marketing Manager
0.010349288 0.014230272
Marketing Specialist Mechanic
0.009055627 0.010349288
Mechanical Engineer Medical Assistant
0.005174644 0.001293661
Musician Nanny
0.003880983 0.002587322
Network Administrator Network Engineer
0.002587322 0.001293661
Nurse Nurse Practitioner
0.029754204 0.001293661
Nutritional Specialist Nutritionist
0.001293661 0.001293661
Operations Manager Painter
0.003880983 0.002587322
Personal Trainer Pharmacist
0.001293661 0.007761966
Photographer Physician
0.028460543 0.005174644
Physicist Physiotherapist
0.001293661 0.002587322
Pilot Plumber
0.001293661 0.011642950
Police Officer Potter
0.007761966 0.002587322
Primary School Teacher Product Designer
0.002587322 0.001293661
Product Manager Program Manager
0.003880983 0.001293661
Project Coordinator Project Manager
0.001293661 0.019404916
Psychologist Public Relations Specialist
0.007761966 0.005174644
Real Estate Agent Receptionist
0.009055627 0.001293661
Research Analyst Research Assistant
0.001293661 0.003880983
Research Scientist Researcher
0.006468305 0.002587322
Restaurant Manager Retail Manager
0.002587322 0.005174644
Retail Worker Retired
0.002587322 0.001293661
Sales Executive Sales Manager
0.001293661 0.011642950
Sales Representative Salesperson
0.007761966 0.001293661
Scientist Seamstress
0.021992238 0.002587322
Secretary Security Guard
0.002587322 0.003880983
Security Officer SEO Specialist
0.001293661 0.006468305
Shopkeeper Social Worker
0.002587322 0.006468305
Software Architect Software Developer
0.002587322 0.019404916
Software Engineer Software Tester
0.020698577 0.002587322
Street Vendor Student
0.002587322 0.003880983
Surgeon Tailor
0.001293661 0.006468305
Taxi Driver Teacher
0.002587322 0.047865459
Technician Therapist
0.001293661 0.001293661
Truck Driver UX Designer
0.011642950 0.001293661
Vegetable Vendor Veterinarian
0.002587322 0.010349288
Waitress Warehouse Worker
0.003880983 0.001293661
Weaver Web Developer
0.002587322 0.018111255
Writer
0.006468305 A la vista de los puestos de trabajo, vamos a reagruparlos en 4 nuevas categorías (Ciencia, Público, Creativo y Práctico y Económico)
datos$Occupation <- case_when(
datos$Occupation %in% c(
"Data Analyst","Data Engineer","Data Scientist","Database Administrator",
"Software Developer","Software Engineer","Software Architect","Software Tester",
"Web Developer","Developer",
"IT Consultant","IT Manager","IT Specialist","IT Support","IT Support Specialist",
"Network Engineer","Network Administrator",
"Scientist","Physicist","Biologist",
"Research Analyst","Research Assistant","Research Scientist","Researcher",
"Engineer","Civil Engineer","Electrical Engineer","Mechanical Engineer",
"Architect","Laboratory Technician"
) ~ "Ciencia",
datos$Occupation %in% c(
"Doctor","Physician","Surgeon","Dentist",
"Nurse","Nurse Practitioner",
"Therapist","Physiotherapist","Psychologist",
"Pharmacist",
"Healthcare Assistant","Medical Assistant",
"Nutritionist","Nutritional Specialist",
"Veterinarian",
"Fitness Trainer","Fitness Instructor","Personal Trainer",
"Waitress","Bartender","Chef","Baker","Bakery Owner",
"Shopkeeper","Retail Worker","Retail Manager","Restaurant Manager",
"Event Planner","Event Manager","Event Coordinator",
"Receptionist","Customer Support","Courier",
"Hair Stylist","Nanny","Social Worker","Freelancer",
"Digital Marketer","Marketing Executive",
"Pilot", "Civil Servant", "Police Officer", "Secretary",
"CEO","Manager","Executive Director",
"Operations Manager","Product Manager","Program Manager","Project Manager","Project Coordinator",
"Business Owner","Entrepreneur","Retired","Student", "Marketing Manager", "Advertising Manager", "Advertising Executive", "Brand Manager"
) ~ "Público",
datos$Occupation %in% c(
"Teacher","Primary School Teacher",
"Artist","Musician","Writer","Photographer","Journalist","Editor",
"Graphic Designer","UX Designer","Designer","Fashion Designer",
"Content Creator","Content Writer","Content Strategist","Copywriter",
"Interior Designer","Product Designer","Handicrafts Maker","Weaver","Potter", "Actor",
"Mechanic","Electrician","Plumber",
"Carpenter","Blacksmith","Technician","Electrical Technician",
"Factory Worker","Construction Worker","Construction Manager","Construction Engineer",
"Driver","Taxi Driver","Truck Driver","Delivery Driver","Bus Driver",
"Farmer","Fisherwoman","Flower Seller","Vegetable Vendor","Street Vendor",
"Security Guard","Security Officer","Firefighter",
"Cleaner","Janitor","Cobbler","Tailor","Seamstress",
"Painter","Warehouse Worker"
) ~ "Creativo y práctico",
datos$Occupation %in% c(
"Accountant","Account Manager",
"Financial Analyst","Financial Advisor","Financial Planner",
"Bank Manager","Banker",
"Business Analyst","Business Consultant","Consultant",
"HR Manager","HR Specialist","HR Executive","Human Resources","Human Resources Manager",
"Marketing Specialist","Marketing Executive","Marketing Director","SEO Specialist",
"Sales Executive","Sales Representative","Sales Manager","Salesperson",
"Insurance Agent","Real Estate Agent","Lawyer","Librarian","Public Relations Specialist"
) ~ "Económico",
TRUE ~ "Sin Clasificar"
)
datos$Occupation <- factor(datos$Occupation)
table(datos$Occupation)/nrow(datos)
Ciencia Creativo y práctico Económico Público
0.2121604 0.3518758 0.1630013 0.2729625 Ahora que la variable se encuentra recodificada, se puede observar que las categorías se encuentran más o menos equilibradas y será más cómodo trabajar con ellas.
Además, se van a revisar también los valores y el número de individuos en las distintas categorías de la variable Excersice_Type para que no haya categorías muy desbalanceadas
table(datos$Exercise_Type)
Aerobics Cardio Meditation Pilates
5 215 24 88
Strength Training Walking Yoga
245 5 191 Se observa que en Aerobics, Walking y Meditation hay muy pocos individuos. Vamos a recodificar la variable de manera que Walking y Aerobics estarán incluidas en la categoría Cardio y Meditation y Pilates se combinarán con la categoría Yoga
datos$Exercise_Type <- case_when(
datos$Exercise_Type %in% c("Cardio", "Walking", "Aerobics") ~ "Cardio",
datos$Exercise_Type == "Strength Training" ~ "Fuerza",
datos$Exercise_Type %in% c("Yoga", "Pilates", "Meditation") ~ "Cuerpo y Mente",
TRUE ~ "Otros"
)
datos$Exercise_Type <- factor(datos$Exercise_Type)
table(datos$Exercise_Type)
Cardio Cuerpo y Mente Fuerza
225 303 245 Finalmente antes de comenzar a realizar los modelos y con ayuda de la función summary() se va a confirmar que no hay ninguna variable con datos ausentes o atípicos
summary(datos) Age Gender Occupation Marital_Status
Min. :18.00 Male :389 Ciencia :164 Single :353
1st Qu.:33.00 Female:384 Creativo y práctico:272 Married :360
Median :39.00 Económico :126 Divorced: 60
Mean :38.89 Público :211
3rd Qu.:45.00
Max. :60.00
Sleep_Duration Sleep_Quality Wake_Up_Time Bed_Time
Min. :3.500 Pobre :164 Min. :4.000 Min. :16.00
1st Qu.:6.000 Aceptable:181 1st Qu.:6.000 1st Qu.:18.00
Median :6.300 Buena :263 Median :7.000 Median :22.00
Mean :6.338 Excelente:165 Mean :6.743 Mean :20.54
3rd Qu.:7.000 3rd Qu.:7.000 3rd Qu.:23.00
Max. :8.000 Max. :9.000 Max. :25.00
Physical_Activity Screen_Time Caffeine_Intake Alcohol_Intake
Min. :1.000 Min. :2.000 Min. :0.000 Min. :0.0000
1st Qu.:2.000 1st Qu.:4.000 1st Qu.:1.000 1st Qu.:0.0000
Median :3.000 Median :4.000 Median :2.000 Median :1.0000
Mean :2.979 Mean :4.105 Mean :1.819 Mean :0.8887
3rd Qu.:4.000 3rd Qu.:5.000 3rd Qu.:2.000 3rd Qu.:1.0000
Max. :5.000 Max. :8.000 Max. :4.000 Max. :2.0000
Smoking_Habit Work_Hours Travel_Time Social_Interactions
No :358 Min. : 6.000 Min. :0.500 Min. :1.000
Yes:415 1st Qu.: 8.000 1st Qu.:2.000 1st Qu.:3.000
Median : 8.000 Median :3.000 Median :3.000
Mean : 8.259 Mean :2.858 Mean :3.197
3rd Qu.: 9.000 3rd Qu.:4.000 3rd Qu.:4.000
Max. :14.000 Max. :5.000 Max. :5.000
Meditation_Practice Exercise_Type Blood_Pressure Cholesterol_Level
No :292 Cardio :225 Min. :110.0 Min. :150.0
Yes:481 Cuerpo y Mente:303 1st Qu.:130.0 1st Qu.:210.0
Fuerza :245 Median :140.0 Median :220.0
Mean :137.9 Mean :220.8
3rd Qu.:150.0 3rd Qu.:230.0
Max. :170.0 Max. :290.0
Blood_Sugar_Level Stress_Detection
Min. : 80.0 Low :162
1st Qu.:105.0 Medium:310
Median :115.0 High :301
Mean :111.8
3rd Qu.:120.0
Max. :150.0 Se observa que no hay ningún dato atípico, las variables están bien codificadas y toman valores dentro de rangos normales.
Se puede comenzar a realizar el estudio.
Antes de comenzar a construir los modelos, se va a realizar la división entrenamiento-prueba para poder así contar con un subconjunto de datos diferente para evaluar el modelo más adelante.
Para ello, se asignará el 80% de los datos a entrenamiento y el 20% restante a prueba manteniendo las proporciones de cada categoría de la variable objetivo Sleep_Quality.
set.seed(12345)
trainIndex <- createDataPartition(datos$Sleep_Quality, p=0.8, list=FALSE)
data_train <- datos[trainIndex,]
data_test <- datos[-trainIndex,]La variable objetivo Sleep_Quality es una variable ordinal, es decir, es una variable categórica que cuenta con un orden natural entre sus categorías pero la distancia entre estas no es necesariamente constante.
Debido a la naturaleza de la misma, para seleccionar el mejor modelo GLM vamos a hacer regresión ordinal y multinomial.
La regresión ordinal aprovecha el orden de las categorías, suele ser más eficiente que la multinomial y cuenta con menos parámetros ya que requiere menos datos. Entre sus limitaciones se encuentra que no hay probabilidades individuales por categoría directamente sin transformaciones.
Por otra parte, la regresión multinomial es categórica pero no asume orden, calcula la probabilidad de cada cada categoría suceda frente a una de referencia. Es flexible puesto que permite efectos diferentes para cada categoría con respecto a la de referencia pero pierde la información del orden, requiere de más parámetros y su interpretación será más compleja.
Antes de comenzar a realizar el modelo, los nombres de las categorías deben cumplir una serie de requisitos para evitar problemas posteriores con las funciones. Como la variable que se va a utilizar ya se encuentran bien codificada ya que no cuenta con tíldes, espacios, números, etc. este paso no es necesario, pero es recomendable mencionarlo.
table(datos$Sleep_Quality)/nrow(datos)
Pobre Aceptable Buena Excelente
0.2121604 0.2341527 0.3402329 0.2134541 La variable objetivo es de tipo ordinal (pues los niveles se pueden ordenar de manera natural), por lo que es importante verificar que sus niveles están ordenados según su significado para que los modelos funcionen correctamente. La salida anterior permite verificar que es así
A continuación, vamos a construir un primer modelo de regresión logística ordinal con todas las variables independientes disponibles.
Se va a construir un primer modelo de regresión logísitca ordinal utilizando todas las variables dependientes disponibles.
modeloInicial <- polr(Sleep_Quality ~ ., data=data_train)
modeloInicial$edf[1] 29summary(modeloInicial)
Re-fitting to get HessianCall:
polr(formula = Sleep_Quality ~ ., data = data_train)
Coefficients:
Value Std. Error t value
Age 0.003568 0.01235 0.288952
GenderFemale -0.395441 0.17933 -2.205095
OccupationCreativo y práctico -0.101863 0.21775 -0.467802
OccupationEconómico 0.157991 0.26886 0.587634
OccupationPúblico 0.002096 0.23141 0.009056
Marital_StatusMarried 0.043855 0.18151 0.241610
Marital_StatusDivorced 0.721947 0.35843 2.014208
Sleep_Duration 0.344056 0.13788 2.495399
Wake_Up_Time 0.306777 0.09523 3.221587
Bed_Time 0.116501 0.05069 2.298533
Physical_Activity -0.256757 0.19235 -1.334853
Screen_Time -0.204107 0.16668 -1.224551
Caffeine_Intake 0.576920 0.18834 3.063256
Alcohol_Intake 0.650838 0.18499 3.518206
Smoking_HabitYes -0.258879 0.17843 -1.450877
Work_Hours -0.569782 0.12525 -4.549142
Travel_Time -0.354117 0.13752 -2.575012
Social_Interactions 1.229568 0.14188 8.666382
Meditation_PracticeYes 0.241184 0.20320 1.186953
Exercise_TypeCuerpo y Mente -0.039823 0.20609 -0.193230
Exercise_TypeFuerza 0.260712 0.20086 1.297984
Blood_Pressure 0.025111 0.01596 1.573721
Cholesterol_Level -0.010956 0.01454 -0.753495
Blood_Sugar_Level 0.042395 0.02191 1.934530
Stress_DetectionMedium -0.311657 0.21630 -1.440867
Stress_DetectionHigh -0.297472 0.25975 -1.145226
Intercepts:
Value Std. Error t value
Pobre|Aceptable 8.8044 0.0085 1035.1462
Aceptable|Buena 10.2481 0.1105 92.7496
Buena|Excelente 12.1914 0.1584 76.9533
Residual Deviance: 1453.384
AIC: 1511.384 En esta salida se puede observar el valor de cada uno de los parámetros estimados junto con sus errores estándar y el AIC del modelo.
Este primer modelo cuenta con 29 parámetros.
A destacar que esta salida no muestra los p-valores pero incluye los valores t en los que se basa. Como el número de observaciones es grande, la distribución a aplicar se asemeja mucho a una distribución normal por lo que podemos tomar como valor de referencia |1.96| para determinar si los parámetros son significativos al 95%.
En este caso, se observa que hay 10 parámetros significativos.
Respecto a los parámetros alpha (intercepts) no se interpretan ni contrasta si son significativametne distintos de 0.
Este tipo de modelos suelen ser más rígidos que modelos multinomiales y pueden resultar de peor calidad en ocasiones.
Para estudiar la significatividad de las variables y obtener una ordenación de utilidad vamos a aplicar el análisis de tipo II
Anova(modeloInicial, type = "II")Analysis of Deviance Table (Type II tests)
Response: Sleep_Quality
LR Chisq Df Pr(>Chisq)
Age 0.083 1 0.7726341
Gender 4.723 1 0.0297593 *
Occupation 1.172 3 0.7596725
Marital_Status 4.190 2 0.1230756
Sleep_Duration 5.156 1 0.0231658 *
Wake_Up_Time 7.807 1 0.0052038 **
Bed_Time 2.805 1 0.0939794 .
Physical_Activity 1.786 1 0.1813554
Screen_Time 1.502 1 0.2203055
Caffeine_Intake 9.320 1 0.0022664 **
Alcohol_Intake 12.424 1 0.0004239 ***
Smoking_Habit 2.112 1 0.1461190
Work_Hours 22.597 1 1.998e-06 ***
Travel_Time 6.449 1 0.0111031 *
Social_Interactions 82.588 1 < 2.2e-16 ***
Meditation_Practice 1.404 1 0.2361139
Exercise_Type 2.645 2 0.2664611
Blood_Pressure 2.157 1 0.1418787
Cholesterol_Level 0.564 1 0.4526430
Blood_Sugar_Level 3.740 1 0.0531293 .
Stress_Detection 2.150 2 0.3412181
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se observan 8 variables significativas al 5% y se ve que las variables más importantes son Social_Interactions, Work_Hours, Exercise_Type, Alcohol_Intake, Wake_Up_Time y Travel_Time.
Vamos a generar un segundo modelo con las 6 variables más importantes
modelo2<-polr(Sleep_Quality ~ Social_Interactions + Work_Hours + Alcohol_Intake + Wake_Up_Time + Caffeine_Intake + Travel_Time, data=data_train)
summary(modelo2)
Re-fitting to get HessianCall:
polr(formula = Sleep_Quality ~ Social_Interactions + Work_Hours +
Alcohol_Intake + Wake_Up_Time + Caffeine_Intake + Travel_Time,
data = data_train)
Coefficients:
Value Std. Error t value
Social_Interactions 1.2056 0.11425 10.552
Work_Hours -0.3436 0.09232 -3.723
Alcohol_Intake 0.4579 0.15784 2.901
Wake_Up_Time 0.2699 0.08398 3.213
Caffeine_Intake 0.6612 0.15022 4.401
Travel_Time -0.3843 0.10483 -3.665
Intercepts:
Value Std. Error t value
Pobre|Aceptable 1.7204 1.0679 1.6110
Aceptable|Buena 3.0853 1.0723 2.8774
Buena|Excelente 4.9365 1.0857 4.5470
Residual Deviance: 1509.828
AIC: 1527.828 Observamos ahora que no todos los parámetros son significativos como tipo de ejercicio cuando toma la categoría Fuerza o Cuerpo y Mente o variables que si parecían importantes como Travel_Time. Además, el AIC es mayor al del modelo con todo (pasa de 1511 a 1538).
Para verificar si estamos prescindiendo de variables necesarias, vamos a usar la función add1()
add1(modelo2, scope=formula(modeloInicial),test = "Chisq")Single term additions
Model:
Sleep_Quality ~ Social_Interactions + Work_Hours + Alcohol_Intake +
Wake_Up_Time + Caffeine_Intake + Travel_Time
Df AIC LRT Pr(>Chi)
<none> 1527.8
Age 1 1524.8 5.0632 0.024439 *
Gender 1 1517.0 12.8031 0.000346 ***
Occupation 3 1531.9 1.9502 0.582814
Marital_Status 2 1527.3 4.5622 0.102174
Sleep_Duration 1 1527.3 2.4941 0.114272
Bed_Time 1 1529.7 0.1078 0.742639
Physical_Activity 1 1529.6 0.2520 0.615678
Screen_Time 1 1529.4 0.3906 0.532004
Smoking_Habit 1 1526.0 3.8217 0.050593 .
Meditation_Practice 1 1529.3 0.5725 0.449250
Exercise_Type 2 1524.9 6.9324 0.031236 *
Blood_Pressure 1 1512.3 17.5608 2.783e-05 ***
Cholesterol_Level 1 1512.6 17.2487 3.279e-05 ***
Blood_Sugar_Level 1 1507.3 22.5385 2.060e-06 ***
Stress_Detection 2 1530.4 1.4235 0.490793
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se verifica que, efectivamente, hay varias variables que aportarían información relevante al modelo (las más importantes serían Caffeine_Intake o Blood_Preassure pero otras como Alcohol_Intake o Blood_Sugar_Level también podrían influir).
Dado que el modelo construido manualmente seguramente sea demasiado sencillo, se van a obtener otra serie de modelos aplicando los métodos automáticos de selección de variables.
Se deben construir dos modelos iniciales, uno máximo y otro mínimo y se aplicaran los 3 métodos de selección de variables (Forward, Backward y Step-Wise) y los 2 criterios de selección (AIC y BIC).
null<-polr(Sleep_Quality ~ 1, data=data_train)
full<-polr(Sleep_Quality ~ ., data=data_train)
modeloStepBIC<-step(null, scope=list(lower=null, upper=full), direction="both",
k=log(nrow(data_train)),trace=F)
modeloStepAIC<-step(null, scope=list(lower=null, upper=full), direction="both",trace=F)
modeloForwBIC<-step(null, scope=list(lower=null, upper=full), direction="forward",
k=log(nrow(data_train)),trace=F)
modeloForwAIC<-step(null, scope=list(lower=null, upper=full), direction="forward",trace=F)
modeloBackBIC1<-step(full, scope=list(lower=null, upper=full), direction="backward",
k=log(nrow(data_train)),trace=F)
modeloBackAIC<-step(full, scope=list(lower=null, upper=full), direction="backward",trace=F)
modelos<-list(modeloStepBIC,modeloStepAIC,modeloForwBIC,modeloForwAIC,modeloBackBIC1,
modeloBackAIC)
map(modelos,function(x) x$call)[[1]]
polr(formula = Sleep_Quality ~ Social_Interactions + Alcohol_Intake +
Work_Hours + Blood_Sugar_Level + Travel_Time + Caffeine_Intake +
Wake_Up_Time, data = data_train)
[[2]]
polr(formula = Sleep_Quality ~ Social_Interactions + Alcohol_Intake +
Work_Hours + Blood_Sugar_Level + Sleep_Duration + Travel_Time +
Caffeine_Intake + Wake_Up_Time + Gender + Screen_Time + Bed_Time +
Smoking_Habit + Marital_Status, data = data_train)
[[3]]
polr(formula = Sleep_Quality ~ Social_Interactions + Alcohol_Intake +
Work_Hours + Blood_Sugar_Level + Sleep_Duration + Travel_Time +
Caffeine_Intake + Wake_Up_Time, data = data_train)
[[4]]
polr(formula = Sleep_Quality ~ Social_Interactions + Alcohol_Intake +
Work_Hours + Blood_Sugar_Level + Sleep_Duration + Travel_Time +
Caffeine_Intake + Wake_Up_Time + Gender + Screen_Time + Bed_Time +
Smoking_Habit + Marital_Status, data = data_train)
[[5]]
polr(formula = Sleep_Quality ~ Wake_Up_Time + Caffeine_Intake +
Alcohol_Intake + Work_Hours + Travel_Time + Social_Interactions +
Blood_Sugar_Level, data = data_train)
[[6]]
polr(formula = Sleep_Quality ~ Gender + Marital_Status + Sleep_Duration +
Wake_Up_Time + Bed_Time + Physical_Activity + Caffeine_Intake +
Alcohol_Intake + Smoking_Habit + Work_Hours + Travel_Time +
Social_Interactions + Blood_Sugar_Level, data = data_train)map_int(modelos,function(x) x$edf)[1] 10 17 11 17 10 17Los 6 procesos dan como resultado a 4 modelos distintos, con 10, 11 y 17 parámetros respectivamente.
El modelo de 10 parámetros y el de 11 son iguales, salvo que el último contiene además la variable sleep_duration y los de 17 parámetros añaden además Gender, Bed_Time, Smoking_Habit, Marital_Status y uno de ellos añade Screen_Time mientras que el otro hace uso de Physical_Activity.
Para poder determinar cuál de los modelos es mejor, se va a recurrir a la validación cruzada repetida, gracias a la cual se podran evaluar los modelos mediante las medidas ya estudiadas.
Se van a generar los valores de múltiples métricas, como la tasa de acierto (acc) o el índice Kappa (tanto ponderado como sin ponderar), y guardar los resultados de los cuatro modelos en la lista “vcrTodosModelos” (el modelo generado con las 6 variables más importantes y los 4 modelos distintos que se acaban de obtener con selección automática).
my_summary <- function(data, lev = NULL, model = NULL){
a1 <- multiClassSummary(data, lev, model)
b1 <- vcd::Kappa(table(data[, "pred"],data[, "obs"]))$Weighted[1]
out <- c(a1, Wkappa=b1)
out
}
modelos<-list(modeloForwBIC, modeloForwAIC, modeloBackBIC1, modeloBackAIC)
vcrTodosModelos<-list()
for (i in 1:length(modelos)){
set.seed(1712)
vcr <- caret::train (formula(modelos[[i]]),
data = data_train, method = "polr",
tuneGrid = expand.grid(method="logistic"),
trControl = trainControl(method = "repeatedcv",
number = 5, repeats=20,
classProbs = TRUE,
summaryFunction = my_summary)
)
vcrTodosModelos[[i]]<-vcr
}A continuación, se representan en diagramas de cajas y se obtienen algunos estadísticos de las repeticiones
names(vcrTodosModelos)<-paste0("Model",1:length(modelos), "_",sapply(modelos,function(x) x$edf))
bwplot(resamples(vcrTodosModelos),metric=c("AUC","Kappa","Accuracy","Wkappa.value"),
scales = list(x = list(relation = "free")))
summary(resamples(vcrTodosModelos),metric=c("AUC","Kappa","Accuracy","Wkappa.value"))
Call:
summary.resamples(object = resamples(vcrTodosModelos), metric =
c("AUC", "Kappa", "Accuracy", "Wkappa.value"))
Models: Model1_11, Model2_17, Model3_10, Model4_17
Number of resamples: 100
AUC
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_11 0.5852485 0.6702389 0.6858567 0.6868145 0.7061416 0.7641823 0
Model2_17 0.6002285 0.6684757 0.6893721 0.6882314 0.7054301 0.7579099 0
Model3_10 0.5914135 0.6674673 0.6816280 0.6843980 0.7048366 0.7615330 0
Model4_17 0.5976624 0.6708832 0.6912632 0.6908514 0.7082052 0.7625811 0
Kappa
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_11 0.07412129 0.1646446 0.2164209 0.2122370 0.2574785 0.3502919 0
Model2_17 0.09872836 0.1546579 0.1879026 0.1900882 0.2214593 0.3037493 0
Model3_10 0.06905181 0.1725073 0.2177560 0.2138844 0.2532711 0.3714080 0
Model4_17 0.11021159 0.1545181 0.1869209 0.1914060 0.2281012 0.3214926 0
Accuracy
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_11 0.3387097 0.4056845 0.4453514 0.4388968 0.4715447 0.5365854 0
Model2_17 0.3571429 0.3920000 0.4227642 0.4198510 0.4408871 0.5000000 0
Model3_10 0.3467742 0.4080000 0.4417742 0.4397820 0.4658862 0.5520000 0
Model4_17 0.3600000 0.3943710 0.4176774 0.4210530 0.4480000 0.5161290 0
Wkappa.value
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_11 0.1620553 0.2963662 0.3508436 0.3422856 0.3840005 0.5194446 0
Model2_17 0.2068607 0.2950057 0.3319486 0.3315417 0.3734170 0.4755100 0
Model3_10 0.1556116 0.3030629 0.3497293 0.3491425 0.3942766 0.5021070 0
Model4_17 0.2109375 0.2933794 0.3254369 0.3298699 0.3693378 0.4605263 0Para determinar el mejor modelo, se estudia la capacidad predictiva, variabilidad y número de parámetros de los modelos.
La capacidad predictiva se evalúa observando la mediana. Queremos la mediana más alta para un mejor rendimiento. Observamos como los modelos de 10 y 11 parámetros muestran medianas ligeramente superiores en Kappa y Wkappa y Accuracy
La variabilidad se mide por la dispersión de los resultados. Un modelo será más estable cuanto más estrechos sean su caja y bigotes. Todos los modelos muestran una variabilidad similar pero en general los modelos de más paráemtros parecen tener cajas más estrechas, lo que sugiere mayor estabilidad
Finalmente, el principio de parsimonia sugiere elegir el modelo más simple (con menor número de parámetros) que explique tan bien los datos como un modelo más complejo.
Realizado el estudio, y buscando el mejor equilibrio entre capacidad predictiva, variabilidad y número de parámetros, el modelo más prometedor es el Model3_10, es decir BackBIC, ya que es de los que menos parámetros tiene, su rendimiento es muy similar a los modelos de 17 parámetros aún reduciendo en un 30% su número de parámetros y todos los valores entre modelos son muy cercanos, por lo que nos quedamos con el más sencillo.
Seleccionado el mejor modelo, se obtienen los ODDSratio del mismo y se va a realizar una evaluación y análisis final mediante los datos de prueba
En primer lugar, se observa la exponencial de los coeficientes junto con su significatividad.
stargazer(modeloBackBIC1, type="text", report=('vc*p'), apply.coef=exp, p.auto=F)
===============================================
Dependent variable:
---------------------------
Sleep_Quality
-----------------------------------------------
Wake_Up_Time 1.367***
p = 0.0003
Caffeine_Intake 1.564***
p = 0.006
Alcohol_Intake 1.526***
p = 0.008
Work_Hours 0.638***
p = 0.00001
Travel_Time 0.586***
p = 0.00001
Social_Interactions 3.478***
p = 0.000
Blood_Sugar_Level 1.041***
p = 0.00001
-----------------------------------------------
Observations 620
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01Todos los parámetros son significativos.
La interpretación de los parámetros es la siguiente
Wake_Up_Time: Por cada hora adicional en que tarde en despertarse, las posibilidades de tener una percepción de calidad de sueño mayor aumentan en un 37%.
Caffeine_Intake: Por cada taza de café adicional ingerida en el día, las posibilidades de tener una percepción de calidad de sueño mayor aumentan en un 56%.
Alcohol_Intake: Por cada bebida alcohólica adicional que se tome durante la semana, las posibilidades de tener una percepción de calidad de sueño mayor aumentan en un 53%.
Work_Hours: Por cada hora adicional que se trabaje, las posibilidades de tener una percepción de calidad de sueño mayor disminuyen en un 36%.
Travel_Time: Por cada hora adicional que haya que emplear en transporte, las posibilidades de tener una percepción de calidad de sueño mayor disminuyen en un 41%.
Social_Interactions: Por cada interacción social adicional durante la semana, las posibilidades de tener una percepción de calidad de sueño mayor aumentan casi por 3,5.
Blood_Sugar_Level: Finalmente, por cada mg/dL adicional de nivel de azúcar en sangre, las posibilidades de tener una percepción de calidad de sueño mayor aumentan en un 4%.
Para completar el estudio de las variables que componen el modelo, se realiza nuevamente el análisis de tipo II, detectando las variables más importantes.
Anova(modeloBackBIC1,type = "II")Analysis of Deviance Table (Type II tests)
Response: Sleep_Quality
LR Chisq Df Pr(>Chisq)
Wake_Up_Time 13.470 1 0.0002424 ***
Caffeine_Intake 7.684 1 0.0055719 **
Alcohol_Intake 7.119 1 0.0076261 **
Work_Hours 22.784 1 1.813e-06 ***
Travel_Time 23.234 1 1.435e-06 ***
Social_Interactions 127.553 1 < 2.2e-16 ***
Blood_Sugar_Level 22.538 1 2.060e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se observa como la variable más importante sigue siendo Social_Interactions seguida de Work_Hours, Blood_Sugar_Level y Travel_Time
Finalmente, se obtienen las medidas de evaluación del modelo
cm_test<-confusionMatrix(table(predict(modeloBackBIC1,newdata=data_test),data_test$Sleep_Quality))
cm_test$table
Pobre Aceptable Buena Excelente
Pobre 14 1 2 0
Aceptable 3 8 5 0
Buena 13 26 39 26
Excelente 2 1 6 7cm_test$overall[1:2] Accuracy Kappa
0.4444444 0.2047331 cm_test$byClass[,1:2] Sensitivity Specificity
Class: Pobre 0.4375000 0.9752066
Class: Aceptable 0.2222222 0.9316239
Class: Buena 0.7500000 0.3564356
Class: Excelente 0.2121212 0.9250000Se observa que en términos generales no parece ser un muy buen modelo, pues cuenta con una tasa de acierto del 44% y un índice kappa del 20%.
A pesar de que es capaz de detectar de manera bastante sólida cuando una observacion no pertenece a las clases Pobre, Aceptable o Excelente, no se le da especialmente bien clasificar correctamente una observación en ninguna de las clases.
Obtenemos el kappa ponderado para ampliar el estudio
summary(Kappa(cm_test$table)) value ASE z Pr(>|z|)
Unweighted 0.2047 0.05332 3.84 1.233e-04
Weighted 0.3254 0.05679 5.73 1.003e-08
Weights:
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.6666667 0.3333333 0.0000000
[2,] 0.6666667 1.0000000 0.6666667 0.3333333
[3,] 0.3333333 0.6666667 1.0000000 0.6666667
[4,] 0.0000000 0.3333333 0.6666667 1.0000000Como se ha dicho al empezar el modelo de regresion ordinal, también vamos a evaluar el modelo multinomial. Para ello vamos a generar un primer modelo multinomial que contenga todas las variables.
modeloInicial<-multinom(Sleep_Quality ~ ., data=data_train)# weights: 112 (81 variable)
initial value 859.502504
iter 10 value 742.625609
iter 20 value 673.004120
iter 30 value 618.821732
iter 40 value 610.531223
iter 50 value 607.750078
iter 60 value 605.965876
iter 70 value 604.301282
iter 80 value 603.208883
iter 90 value 602.973731
final value 602.972999
convergedEste primer modelo cuenta con 81 parámetros ya que el valor que aparece entre paréntesis denominado como variable en realidad nos muestra el número de parámetros. Como el proceso converge, el resultado obtenido es óptimo
Para poder conocer más información del modelo vamos a recurrir a la función summary().
summary(modeloInicial)Call:
multinom(formula = Sleep_Quality ~ ., data = data_train)
Coefficients:
(Intercept) Age GenderFemale OccupationCreativo y práctico
Aceptable -9.639836 -0.006934332 -0.8495559 0.9666152
Buena -13.819237 0.015928162 -0.8194338 0.7164468
Excelente -20.330888 -0.003942244 -0.9808427 -0.0100921
OccupationEconómico OccupationPúblico Marital_StatusMarried
Aceptable 1.0376902 0.48750300 0.4833062
Buena 0.6517466 0.39469798 0.3731659
Excelente 0.5877008 -0.01837401 0.5812440
Marital_StatusDivorced Sleep_Duration Wake_Up_Time Bed_Time
Aceptable -10.3821956 0.5827337 0.2640550 0.040579964
Buena 0.8101448 0.6729822 0.6298243 0.007341599
Excelente 1.3826525 0.7074802 0.4774802 0.252203350
Physical_Activity Screen_Time Caffeine_Intake Alcohol_Intake
Aceptable 0.3624392 -0.3525726 1.075827 0.2651763
Buena -0.1470358 -0.2144884 1.153036 0.7576151
Excelente 0.2695785 -0.6775997 1.093747 1.3459492
Smoking_HabitYes Work_Hours Travel_Time Social_Interactions
Aceptable 0.9207181 -1.4927209 0.1848512 -0.2551996
Buena 0.1843457 -0.9899435 -1.0942389 0.8661008
Excelente -0.1627461 -1.1038984 -0.8072585 1.5396431
Meditation_PracticeYes Exercise_TypeCuerpo y Mente
Aceptable 0.5078746 -0.4924045
Buena 0.7083076 -0.5440943
Excelente 0.7698420 -0.1558168
Exercise_TypeFuerza Blood_Pressure Cholesterol_Level
Aceptable -0.3125424 0.03843048 0.05047919
Buena -0.5352393 0.04430337 0.03631243
Excelente 0.2622745 0.10722496 -0.03643194
Blood_Sugar_Level Stress_DetectionMedium Stress_DetectionHigh
Aceptable -0.02081721 -0.2928292 -1.2009229
Buena -0.01282188 -0.1993820 -0.3135391
Excelente 0.06093724 -0.6611625 -0.7370610
Std. Errors:
(Intercept) Age GenderFemale OccupationCreativo y práctico
Aceptable 0.02804016 0.02650572 0.3779018 0.4539714
Buena 0.02095464 0.02232101 0.3311888 0.3937103
Excelente 0.01818803 0.02489569 0.3804195 0.4525432
OccupationEconómico OccupationPúblico Marital_StatusMarried
Aceptable 0.5711518 0.4670821 0.3903810
Buena 0.4942071 0.4004470 0.3356230
Excelente 0.5442670 0.4515065 0.3761635
Marital_StatusDivorced Sleep_Duration Wake_Up_Time Bed_Time
Aceptable 1.596939e-05 0.3126537 0.2052379 0.10727459
Buena 5.693911e-01 0.2261870 0.1781259 0.09198937
Excelente 6.386569e-01 0.2594784 0.2040107 0.10404367
Physical_Activity Screen_Time Caffeine_Intake Alcohol_Intake
Aceptable 0.4773544 0.4593549 0.4594658 0.4143905
Buena 0.2857160 0.2548939 0.3047887 0.2899376
Excelente 0.3612844 0.3250287 0.3573751 0.3429740
Smoking_HabitYes Work_Hours Travel_Time Social_Interactions
Aceptable 0.3915810 0.3434452 0.3501883 0.3307150
Buena 0.3155611 0.2207800 0.2791428 0.2427394
Excelente 0.3626576 0.2667834 0.3097909 0.2775032
Meditation_PracticeYes Exercise_TypeCuerpo y Mente
Aceptable 0.4464856 0.4674660
Buena 0.3297740 0.4064315
Excelente 0.3883015 0.4658422
Exercise_TypeFuerza Blood_Pressure Cholesterol_Level
Aceptable 0.4515821 0.03434679 0.03201360
Buena 0.4107217 0.02834662 0.02539155
Excelente 0.4580134 0.03264315 0.02980697
Blood_Sugar_Level Stress_DetectionMedium Stress_DetectionHigh
Aceptable 0.04585855 0.4511677 0.5748418
Buena 0.03836782 0.3671794 0.4411922
Excelente 0.04525401 0.4385469 0.5082228
Residual Deviance: 1205.946
AIC: 1367.946 Como se mencionaba y debido a la naturaleza de este tipo de modelos, hablar sobre los valores de los parámetros y como afectan a cada categoría de la variable objetivo respecto a la categoría de referencia (Pobre) se hace un proceso costoso y largo.
Para poder comparar con el resultado obtenido a partir de la regresión ordinal se procede a crear los modelos automáticos y seleccionar el mejor de los mismos
null <- multinom(Sleep_Quality ~ 1, data = data_train, trace=F)
full <- multinom(Sleep_Quality ~ ., data = data_train, trace=F)
invisible(capture.output({ # Esto se incluye para evitar la aparición del texto "trying..."
modeloStepBIC<-step(null, scope=list(lower=null, upper=full), direction="both", k=log(nrow(data_train)), trace=F)
modeloStepAIC<-step(null, scope=list(lower=null, upper=full), direction="both", k=2, trace=F)
modeloForwBIC<-step(null, scope=list(lower=null, upper=full), direction="forward", k=log(nrow(data_train)),trace=F)
modeloForwAIC<-step(null, scope=list(lower=null, upper=full), direction="forward", k=2, trace=F)
modeloBackBIC<-step(full, scope=list(lower=null, upper=full), direction="backward", k=log(nrow(data_train)),trace=F)
modeloBackAIC<-step(full, scope=list(lower=null, upper=full), direction="backward", k=2, trace=F)
}))
modelos<-list(modeloStepBIC,modeloStepAIC,modeloForwBIC,modeloForwAIC,modeloBackBIC,modeloBackAIC)
map(modelos,function(x) x$call) [[1]]
multinom(formula = Sleep_Quality ~ Social_Interactions + Blood_Pressure +
Wake_Up_Time + Travel_Time + Work_Hours + Caffeine_Intake,
data = data_train, trace = F)
[[2]]
multinom(formula = Sleep_Quality ~ Social_Interactions + Blood_Pressure +
Wake_Up_Time + Travel_Time + Work_Hours + Caffeine_Intake +
Alcohol_Intake + Cholesterol_Level + Gender + Smoking_Habit +
Meditation_Practice + Sleep_Duration + Marital_Status, data = data_train,
trace = F)
[[3]]
multinom(formula = Sleep_Quality ~ Social_Interactions + Blood_Pressure +
Wake_Up_Time + Travel_Time + Work_Hours + Caffeine_Intake,
data = data_train, trace = F)
[[4]]
multinom(formula = Sleep_Quality ~ Social_Interactions + Blood_Pressure +
Wake_Up_Time + Travel_Time + Work_Hours + Caffeine_Intake +
Alcohol_Intake + Cholesterol_Level + Gender + Smoking_Habit +
Meditation_Practice + Sleep_Duration + Marital_Status, data = data_train,
trace = F)
[[5]]
multinom(formula = Sleep_Quality ~ Gender + Marital_Status +
Sleep_Duration + Wake_Up_Time + Caffeine_Intake + Alcohol_Intake +
Smoking_Habit + Work_Hours + Travel_Time + Social_Interactions +
Meditation_Practice + Blood_Pressure + Cholesterol_Level,
data = data_train, trace = F)
[[6]]
multinom(formula = Sleep_Quality ~ Gender + Marital_Status +
Sleep_Duration + Wake_Up_Time + Caffeine_Intake + Alcohol_Intake +
Smoking_Habit + Work_Hours + Travel_Time + Social_Interactions +
Meditation_Practice + Blood_Pressure + Cholesterol_Level,
data = data_train, trace = F)map_int(modelos,function(x) x$edf) [1] 21 45 21 45 45 45En esta ocasión se obtienen 2 modelos, con 21 y 45 parámetros respectivamente.
Para poder determinar cuál de los modelos es mejor, se va a recurrir a la validación cruzada repetida, gracias a la cual se podran evaluar los modelos mediante las medidas ya estudiadas.
Se van a generar los valores de múltiples métricas, como la tasa de acierto (acc) o el índice Kappa (tanto ponderado como sin ponderar), y guardar los resultados de los cuatro modelos en la lista “vcrTodosModelos” (el modelo generado con las 6 variables más importantes y los 4 modelos distintos que se acaban de obtener con selección automática).
my_summary <- function(data, lev = NULL, model = NULL){
a1 <- multiClassSummary(data, lev, model)
b1 <- vcd::Kappa(table(data[, "pred"],data[, "obs"]))$Weighted[1]
out <- c(a1, Wkappa=b1)
out
}
modelos<-list(modeloForwBIC, modeloForwAIC)
vcrTodosModelos<-list()
for (i in 1:length(modelos)){
set.seed(1712)
vcr <- caret::train (formula(modelos[[i]]),
data = data_train, method = "polr",
tuneGrid = expand.grid(method="logistic"),
trControl = trainControl(method = "repeatedcv",
number = 5, repeats=20,
classProbs = TRUE,
summaryFunction = my_summary)
)
vcrTodosModelos[[i]]<-vcr
}Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
: There were missing values in resampled performance measures.A continuación, se representan en diagramas de cajas y se obtienen algunos estadísticos de las repeticiones
names(vcrTodosModelos)<-paste0("Model",1:length(modelos), "_",sapply(modelos,function(x) x$edf))
bwplot(resamples(vcrTodosModelos),metric=c("AUC","Kappa","Accuracy","Wkappa.value"),
scales = list(x = list(relation = "free")))
summary(resamples(vcrTodosModelos),metric=c("AUC","Kappa","Accuracy","Wkappa.value"))
Call:
summary.resamples(object = resamples(vcrTodosModelos), metric =
c("AUC", "Kappa", "Accuracy", "Wkappa.value"))
Models: Model1_21, Model2_45
Number of resamples: 100
AUC
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_21 0.5898264 0.6669773 0.6879562 0.6869375 0.7102239 0.7604921 0
Model2_45 0.5851239 0.6690485 0.6868226 0.6854664 0.7059867 0.7508426 0
Kappa
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_21 0.06346136 0.1750472 0.2180600 0.2127896 0.2513802 0.3614447 0
Model2_45 0.07200876 0.1631936 0.2012339 0.2021320 0.2410507 0.3118881 0
Accuracy
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_21 0.3306452 0.4112903 0.4435484 0.4375301 0.4649355 0.5440000 0
Model2_45 0.3387097 0.4024194 0.4257097 0.4289588 0.4569194 0.5121951 0
Wkappa.value
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
Model1_21 0.1888423 0.2991006 0.3508623 0.3454874 0.3909556 0.4750284 0
Model2_45 0.2046462 0.3059767 0.3443147 0.3435121 0.3832081 0.5016734 0Por el principio de parsimonia y a la vista de la falta de diferencia tanto en capacidad predictiva como en variabilidad, optamos por el modelo con menor número de parámetros formulado a través de ForwBIC.
Seleccionado el mejor modelo, se obtienen los ODDSratio del mismo y se va a realizar una evaluación y análisis final mediante los datos de prueba
En primer lugar, se observa la exponencial de los coeficientes junto con su significatividad.
stargazer(modeloForwBIC, type="text", report=('vc*p'), apply.coef=exp, p.auto=F)
=======================================================
Dependent variable:
-----------------------------------
Aceptable Buena Excelente
(1) (2) (3)
-------------------------------------------------------
Social_Interactions 1.048 3.335*** 6.974***
p = 0.854 p = 0.000 p = 0.000
Blood_Pressure 1.090*** 1.093*** 1.118***
p = 0.00000 p = 0.000 p = 0.000
Wake_Up_Time 1.258* 1.971*** 1.790***
p = 0.093 p = 0.00000 p = 0.00005
Travel_Time 1.295 0.365*** 0.428***
p = 0.341 p = 0.00001 p = 0.0005
Work_Hours 0.395*** 0.613*** 0.538***
p = 0.00001 p = 0.0003 p = 0.0003
Caffeine_Intake 2.292*** 2.796*** 2.826***
p = 0.006 p = 0.00005 p = 0.0004
Constant 0.0003*** 0.00000*** 0.000***
p = 0.00001 p = 0.000 p = 0.000
-------------------------------------------------------
Akaike Inf. Crit. 1,364.713 1,364.713 1,364.713
=======================================================
Note: *p<0.1; **p<0.05; ***p<0.01Anova(modeloBackBIC, type = "II")Analysis of Deviance Table (Type II tests)
Response: Sleep_Quality
LR Chisq Df Pr(>Chisq)
Gender 9.165 3 0.0271791 *
Marital_Status 12.899 6 0.0446648 *
Sleep_Duration 8.625 3 0.0347140 *
Wake_Up_Time 18.581 3 0.0003337 ***
Caffeine_Intake 18.957 3 0.0002790 ***
Alcohol_Intake 14.253 3 0.0025807 **
Smoking_Habit 9.414 3 0.0242613 *
Work_Hours 48.947 3 1.339e-10 ***
Travel_Time 47.481 3 2.745e-10 ***
Social_Interactions 68.047 3 1.118e-14 ***
Meditation_Practice 7.452 3 0.0588145 .
Blood_Pressure 8.560 3 0.0357511 *
Cholesterol_Level 7.469 3 0.0583623 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En este caso, se obtiene para cada variable como influye en cada una de las categorías de la variable objetivo Sleep_Quality a partir de la categoría de referencia Pobre.
Respecto al Género, se observa como ser mujer disminuye las posibilidades de tener una percepción de calidad de sueño Aceptable respecto a Pobre en un 58%, de Buena respecto a Pobre en un 55% y de Excelente respecto a Pobre en un 62%.
Sin embargo, en variables como como la presión sanguínea, un incremento unitario en la presión no influye en las posibilidades de una percepcion de calidad de sueño Aceptable respecto a Pobre ni Buena respecto a Pobre pero si aumentan en un 9% de Excelente respecto a Pobre.
Finalmente, se obtienen las medidas de evaluación del modelo
cm_test<-confusionMatrix(table(predict(modeloForwBIC,newdata=data_test),data_test$Sleep_Quality))
cm_test$table
Pobre Aceptable Buena Excelente
Pobre 13 1 9 2
Aceptable 5 23 11 7
Buena 12 8 27 15
Excelente 2 4 5 9cm_test$overall[1:2] Accuracy Kappa
0.4705882 0.2739469 Kappa(cm_test$table) value ASE z Pr(>|z|)
Unweighted 0.2739 0.05487 4.993 5.959e-07
Weighted 0.2746 0.06184 4.440 8.989e-06cm_test$byClass[,1:2] Sensitivity Specificity
Class: Pobre 0.4062500 0.9008264
Class: Aceptable 0.6388889 0.8034188
Class: Buena 0.5192308 0.6534653
Class: Excelente 0.2727273 0.9083333Comparando los resultados de este modelo respecto a los obtenidos con el modelo de regresión ordinal, las métricas de tasa de acierto e indice kappa han aumentado en 7 puntos porcentuales respectivamente y aproximadamente pero, respecto a sensibilidad y especificidad todas las métricas se han acercado más a 0.5, por lo que el modelo no es especialmente bueno clasificando las observaciones.
Vamos a evaluar un modelo multinomial con las mismas variables que hemos obtenido en el mejor modelo ordinal
modelo_multinom<-multinom(formula(modeloBackBIC1), data=data_train)# weights: 36 (24 variable)
initial value 859.502504
iter 10 value 741.944346
iter 20 value 660.038664
iter 30 value 652.974084
final value 652.919269
convergedcm_test_multi<-confusionMatrix(table(predict(modelo_multinom,newdata=data_test),data_test$Sleep_Quality))
cm_test_multi$overall[1] Accuracy
0.5098039 Kappa(cm_test_multi$table) value ASE z Pr(>|z|)
Unweighted 0.3263 0.05407 6.035 1.592e-09
Weighted 0.3336 0.06007 5.554 2.789e-08lrtest(modeloBackBIC,modelo_multinom)Likelihood ratio test
Model 1: Sleep_Quality ~ Gender + Marital_Status + Sleep_Duration + Wake_Up_Time +
Caffeine_Intake + Alcohol_Intake + Smoking_Habit + Work_Hours +
Travel_Time + Social_Interactions + Meditation_Practice +
Blood_Pressure + Cholesterol_Level
Model 2: Sleep_Quality ~ Wake_Up_Time + Caffeine_Intake + Alcohol_Intake +
Work_Hours + Travel_Time + Social_Interactions + Blood_Sugar_Level
#Df LogLik Df Chisq Pr(>Chisq)
1 45 -624.94
2 24 -652.92 -21 55.965 5.106e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Visto como ninguno de los dos modelos funciona especialmente bien para nuestra base de datos, procedemos a trabajar modelos de predicción basados en árboles
En esta ocasión y de manera análoga a lo que sucedía con los modelos GLM podemos realizar árboles tanto ordinales como de clasificación.
Los árboles ordinales están especialmente diseñados para variables categóricas ordinales. Es importante destacar que aunque es posible construir modelos de clasificación para variables ordinales ignorando la componente ordinal, si el orden de los niveles es importante, como sucede en nuestro caso, es preferible optar por estos pues están diseñados para respetar dicho orden, por lo que sólo se realizarán árboles de este tipo
La función que genera los árboles ordinales precisa que las variables dependientes estén codificadas como numéricas asignando un número según el orden de menor a mayor.
Una vez este codificada correctamente, se realiza nuevamente la partición entrenamiento prueba
datos$Sleep_Quality<-car::recode(datos$Sleep_Quality, as.factor = FALSE,
"'Pobre'=1;'Aceptable'=2;'Buena'=3;'Excelente'=4")
set.seed(12345)
trainIndex <- createDataPartition(datos$Sleep_Quality, p=0.8, list=FALSE)
data_ord_train <- datos[trainIndex,]
data_ord_test <- datos[-trainIndex,]Para ver una primera modelización, vamos a partir de un árbol con tamaño de hoja mínimo del 5%, indicando que no realice la poda ni obtenga reglas sustitutas.
A partir de este árbol generaremos más adelante otros para poder determinar el tamaño mínimo de hoja, parámetro de la poda y número de reglas sustitutas que mejor nos puedan ayudar en un futuro a tomar decisiones
set.seed(12345)
arbolOrd5 <- rpartScore(Sleep_Quality~., data = data_ord_train,
minbucket=ceiling(0.05*nrow(data_ord_train)),
cp=0, maxsurrogate = 0, maxcompete=0)
sum(arbolOrd5$frame$var == "<leaf>") [1] 11rpart.plot(arbolOrd5, nn=TRUE, tweak = 1.2) 
Vemos que se ha construido un árbol con 11 hojas y una profundidad de 6, y que el resultado aparece con una gama de colores donde los claros corresponden con los niveles bajos de la variable, es decir calidad de sueño percibida Pobre o Aceptable y los oscuros son los altos, Buena o Excelente.
Cada nodo muestra la categoría predicha y observamos que al menos hay dos nodos para cada categoría.
Aquellos con mejor percepción de calidad de sueño son aquellos que han tenido más de 3 interacción sociales en la semana y han dormido más de 7.3 horas o, habiendo dormido menos de 7 horas, contaban con un nivel de presión en sangre mayor a 138 y habían empleado en transporte más de 4,5h.
Sin embargo, para asegurar el mejor árbol posible, es preferible dejar que crezca siendo más flexibles con el tamaño mínimo de hoja y asi poder podarlo. Por ello, construimos un arbol con el 1%.
set.seed(12345)
arbolOrd1 <- rpartScore(Sleep_Quality~., data = data_ord_train,
minbucket=ceiling(0.01*nrow(data_ord_train)),
cp=0, maxsurrogate = 0, maxcompete=0)
rpart.plot(arbolOrd1, nn=TRUE, tweak = 1.2) 
Ahora el árbol es prácticamente imposible de ver ya que es muy grande.
Obtenemos una serie de métricas para ver las diferencias entre este árbol y el anterior
modelos<-list(arbolOrd5,arbolOrd1)
sapply(modelos,
function(x) Kappa(table(factor(predict(x,newdata=data_ord_train),1:max(data_ord_train$Sleep_Quality)), data_ord_train$Sleep_Quality))$Weighted[1]) value value
0.5076934 0.6758188 sapply(modelos,
function(x) Kappa(table(factor(predict(x, newdata=data_ord_test),1:max(data_ord_test$Sleep_Quality)),data_ord_test$Sleep_Quality))$Weighted[1]) value value
0.3714540 0.4504643 sapply(modelos,function(x) sum(x$frame$var == "<leaf>"))[1] 11 31Hemos construido un árbol con 31 hojas pero que permite aumentar el índice kappa ponderado significativamente tanto en entrenamiento como en prueba.
De igual manera, ha aumentado el sobreajuste, por lo que vamos a podar el árbol para procurar conseguir un tamaño de árbol más manejable, con similar capacidad predictiva y, a poder ser, reduciendo el sobreajuste.
arbolOrd1$cptable CP nsplit rel error xerror xstd
1 0.144688645 0 1.0000000 1.0000000 0.02839714
2 0.049450549 1 0.8553114 0.8553114 0.02998857
3 0.045787546 4 0.7069597 0.8223443 0.03593199
4 0.031135531 5 0.6611722 0.7051282 0.02514566
5 0.024725275 6 0.6300366 0.6758242 0.01980225
6 0.019230769 8 0.5805861 0.6355311 0.02551644
7 0.016483516 10 0.5421245 0.6465201 0.02766114
8 0.008241758 12 0.5091575 0.6007326 0.03587126
9 0.007326007 14 0.4926740 0.5952381 0.03981147
10 0.006410256 15 0.4853480 0.5952381 0.03981147
11 0.005494505 17 0.4725275 0.5842491 0.03113014
12 0.004120879 20 0.4560440 0.5824176 0.03201725
13 0.003663004 24 0.4395604 0.5787546 0.03163784
14 0.001831502 28 0.4249084 0.5842491 0.03058663
15 0.000000000 30 0.4212454 0.6007326 0.03106002plotcp(arbolOrd1, las=2)
lines(arbolOrd1$cptable[,3], col="red", type = "b")
Se observa la mejora relativa. Se observa que a partir de las 13 hojas el error relativo no disminuye significativamente, por lo que podríamos podar el árbol a ese tamaño.
No obstante y teniendo en cuenta que en árboles se persigue el equilibrio entre la interpretabilidad y capacidad predictiva, un tamaño de 9 también podría ser razonable.
Respecto al sobreajuste, se observa que este es menor en 9 hojas respecto a 13.
Dado que la métrica graficada puede ser algo confusa, generamos una métrica para ver como varía el kappa ponderado en función del número de hojas
set.seed(12345)
CVresults<-xpred.rpart(arbolOrd1, return.all =T)
kappaPondCV<-apply(CVresults, 2,
function(x) Kappa(table(factor(x,1:max(data_ord_train$Sleep_Quality)),
data_ord_train$Sleep_Quality))$Weighted[1])Warning in sqrt((sum(p * sweep(sweep(W, 1, W %*% colSums(p) * (1 - kw)), : Se
han producido NaNsplot(arbolOrd1$cptable[,2]+1, kappaPondCV, type="b", xlab = "Número de hojas")
Como estamos trabajando con una variable ordinal, lo más conveniente es recurrir al kappa ponderado. Viendo el gráfico volvemos a los mismos resultados, a partir de las 13 hojas apenas hay mejoras pero con 9 hojas también se pueden obtener buenos resultados.
Por tanto, procedemos a llevar a cabo dicha poda escribiendo manualmente el valor de alpha para 9 hojas
arbolOrd1_podado<-prune(arbolOrd1,cp=0.022)
rpart.plot(arbolOrd1_podado, nn=TRUE, tweak=1.1, digits=3) 
Observamos un árbol de 9 hojas y profundidad de 5, con hojas de distinta tonalidades.
Se observa como 3 de ellas están asociadas a calidad de sueño percibida Excelente, y el resto de categorías tienen 2 hojas asociadas cada una de ellas.
Para conocer el grado de credibilidad de las predicciones, creamos una nueva representación gráfica del árbol donde se muestra un gráfico de barras para cada hoja
plot(as.party(arbolOrd1_podado), gp = gpar(fontsize = 4), terminal_panel = node_barplot)
Se observa que las hojas más creibles son las asociadas a la categoria 1, es decir, calidad de sueñño percibida Pobre y las dos más a la derecha asociadas a la categoría 4 Excelente. Esto tiene sentido puesto que, por lo general, aquellas poblaciones en los extremos suelen tener características muy diferenciativas respecto a categorías centrales. Aún así, a excepción de la hoja central, todas ellas tienen una categoría de predicción significativamente mayor en porcentaje frente al resto de categorías en la misma hoja bastante notoria.
Finalemnte, realizamos una evaluación del modelo generado tanto para el conjunto de entrenamiento como para el de prueba
cm<-confusionMatrix(table(factor(predict(arbolOrd1_podado, newdata=data_ord_train),1:max(data_ord_train$Sleep_Quality)),data_ord_train$Sleep_Quality))
cm$overall[1] Accuracy
0.5848142 Kappa(cm$table) value ASE z Pr(>|z|)
Unweighted 0.4312 0.02717 15.87 9.965e-57
Weighted 0.5502 0.02480 22.18 5.011e-109cm$byClass[,1:2] Sensitivity Specificity
Class: 1 0.7681159 0.8939709
Class: 2 0.5507246 0.8669439
Class: 3 0.6161137 0.6985294
Class: 4 0.3787879 0.9609856cm_test<-confusionMatrix(table(factor(predict(arbolOrd1_podado, newdata=data_ord_test),1:max(data_ord_test$Sleep_Quality)),data_ord_test$Sleep_Quality))
cm_test$overall[1]Accuracy
0.474026 Kappa(cm_test$table) value ASE z Pr(>|z|)
Unweighted 0.2826 0.05580 5.066 4.067e-07
Weighted 0.4109 0.05436 7.560 4.033e-14cm_test$byClass[,1:2] Sensitivity Specificity
Class: 1 0.7692308 0.8593750
Class: 2 0.4186047 0.8108108
Class: 3 0.4615385 0.6568627
Class: 4 0.3333333 0.9421488Al generar un árbol más grande y posteriormente podarlo, hemos obtenido unas métricas terminar de explicar
Por último, vamos a obtener la importancia de las variables del modelo donde, haciendo uso de la información sobre las divisiones que forman parte del árbol multiplicamos la mejora del índice y el número de observaciones del nodo y las agregamos por variable
auxImportancia<-data.frame(arbolOrd1_podado$splits,
variable=rownames(arbolOrd1_podado$splits))
impVar<-aggregate(count*improve~variable, auxImportancia, sum)
barplot(impVar$`count * improve`/sum(impVar$`count * improve`),
names.arg=impVar$var, las=2)
De las 6 variables que utiliza nuestro modelo, se observa que la que más influencia tiene es Sleep_Duration con casi el 35%, seguida de Blood_Preassure y Social_Interactions cercanas al 25%
Cabe destacar que se pueden generar reglas sustitutas y, en caso de tener datos faltantes para alguna de las “preguntas” puede ser útil en vez de directamente asociar la observación con la categoría que mayor porcentaje de observaciones contenga
set.seed(12345)
arbolOrd1_pod_miss <- rpartScore(Sleep_Quality~., data = data_ord_train,
minbucket=ceiling(0.01*nrow(data_ord_train)),
cp=0.022, maxsurrogate = 2, maxcompete=0)
summary(arbolOrd1_pod_miss)Call:
rpartScore(formula = Sleep_Quality ~ ., data = data_ord_train,
minbucket = ceiling(0.01 * nrow(data_ord_train)), cp = 0.022,
maxsurrogate = 2, maxcompete = 0)
n= 619
CP nsplit rel error xerror xstd
1 0.14468864 0 1.0000000 1.0000000 0.03817265
2 0.04945055 1 0.8553114 0.8553114 0.03162193
3 0.04578755 4 0.7069597 0.8553114 0.03183338
4 0.03113553 5 0.6611722 0.7472527 0.03549524
5 0.02472527 6 0.6300366 0.6721612 0.03397424
6 0.02200000 8 0.5805861 0.6318681 0.03438642
Variable importance
Sleep_Duration Travel_Time Work_Hours Alcohol_Intake
19 14 13 11
Blood_Pressure Blood_Sugar_Level Screen_Time Caffeine_Intake
11 8 8 8
Social_Interactions Cholesterol_Level Stress_Detection
5 4 1
Node number 1: 619 observations, complexity param=0.1446886
predicted score= 3 expected loss= 0.88206785
left son=2 (145 obs) right son=3 (474 obs)
Primary splits:
Social_Interactions < 2.5 to the left, improve=0.09347201, (0 missing)
Surrogate splits:
Work_Hours < 7.5 to the left, agree=0.887, adj=0.517, (0 split)
Cholesterol_Level < 172.5 to the left, agree=0.772, adj=0.028, (0 split)
Node number 2: 145 observations, complexity param=0.04945055
predicted score= 2 expected loss= 0.58620690
left son=4 (47 obs) right son=5 (98 obs)
Primary splits:
Alcohol_Intake < 0.5 to the left, improve=0.2331931, (0 missing)
Surrogate splits:
Screen_Time < 3.5 to the left, agree=0.903, adj=0.702, (0 split)
Caffeine_Intake < 1.5 to the left, agree=0.897, adj=0.681, (0 split)
Node number 3: 474 observations, complexity param=0.04945055
predicted score= 3 expected loss= 0.80590717
left son=6 (324 obs) right son=7 (150 obs)
Primary splits:
Sleep_Duration < 6.95 to the left, improve=0.09408623, (0 missing)
Surrogate splits:
Travel_Time < 2.75 to the right, agree=0.785, adj=0.320, (0 split)
Stress_Detection splits as RLL, agree=0.724, adj=0.127, (0 split)
Node number 4: 47 observations
predicted score= 1 expected loss= 0.25531915
Node number 5: 98 observations
predicted score= 2 expected loss= 0.46938776
Node number 6: 324 observations, complexity param=0.04945055
predicted score= 3 expected loss= 0.93827160
left son=12 (125 obs) right son=13 (199 obs)
Primary splits:
Blood_Pressure < 137.5 to the left, improve=0.1958097, (0 missing)
Surrogate splits:
Travel_Time < 3.5 to the left, agree=0.920, adj=0.792, (0 split)
Blood_Sugar_Level < 112.5 to the left, agree=0.877, adj=0.680, (0 split)
Node number 7: 150 observations, complexity param=0.03113553
predicted score= 3 expected loss= 0.52000000
left son=14 (131 obs) right son=15 (19 obs)
Primary splits:
Sleep_Duration < 7.75 to the left, improve=0.09622678, (0 missing)
Surrogate splits:
Blood_Pressure < 117.5 to the right, agree=0.907, adj=0.263, (0 split)
Cholesterol_Level < 187.5 to the right, agree=0.880, adj=0.053, (0 split)
Node number 12: 125 observations, complexity param=0.04578755
predicted score= 1 expected loss= 0.84800000
left son=24 (15 obs) right son=25 (110 obs)
Primary splits:
Sleep_Duration < 5.25 to the left, improve=0.2009542, (0 missing)
Surrogate splits:
Work_Hours < 6.5 to the left, agree=0.976, adj=0.8, (0 split)
Cholesterol_Level < 230 to the right, agree=0.904, adj=0.2, (0 split)
Node number 13: 199 observations, complexity param=0.02472527
predicted score= 3 expected loss= 0.72361809
left son=26 (164 obs) right son=27 (35 obs)
Primary splits:
Travel_Time < 4.5 to the left, improve=0.09755948, (0 missing)
Surrogate splits:
Social_Interactions < 3.5 to the left, agree=0.844, adj=0.114, (0 split)
Node number 14: 131 observations
predicted score= 3 expected loss= 0.45801527
Node number 15: 19 observations
predicted score= 4 expected loss= 0.05263158
Node number 24: 15 observations
predicted score= 4 expected loss= 0.66666667
Node number 25: 110 observations
predicted score= 1 expected loss= 0.64545455
Node number 26: 164 observations, complexity param=0.02472527
predicted score= 3 expected loss= 0.71951220
left son=52 (122 obs) right son=53 (42 obs)
Primary splits:
Work_Hours < 8.5 to the left, improve=0.05667594, (0 missing)
Surrogate splits:
Cholesterol_Level < 242.5 to the left, agree=0.933, adj=0.738, (0 split)
Blood_Sugar_Level < 127.5 to the left, agree=0.890, adj=0.571, (0 split)
Node number 27: 35 observations
predicted score= 4 expected loss= 0.37142857
Node number 52: 122 observations
predicted score= 3 expected loss= 0.61475410
Node number 53: 42 observations
predicted score= 2 expected loss= 0.69047619 A continuación, vamos a generar un modelo RF, donde se introduce cierta aleatoriedad a la hora de generar los árboles para que entre ellos haya menor correlación y, por tanto, una mejora de la capacidad predictiva.
Debido a la naturaleza de nuestra variable objetivo, haremos uso de bosques aleatorios para clasificación donde nos centraremos en predecir si los individuos pertenecen a una categoría u otra.
En este caso, nos tenemos que asegurar de que la variable es de tipo factor. Como en el apartado anterior la hemos codificado como númerica que toma valores 1,2,3 y 4 basta con pasarla nuevamente a factor
datos$Sleep_Quality <- as.factor(datos$Sleep_Quality)Nuevamente y antes de comenzar con la modelización, llevamos a cabo la partición entrenamiento-prueba
set.seed(12345)
trainIndex <- createDataPartition(datos$Sleep_Quality, p=0.8, list=FALSE)
data_clasif_train <- datos[trainIndex,]
data_clasif_test <- datos[-trainIndex,]Comenzamos creando un RF de 500 árboles con un tamaño de hoja equivalente al 5% con ntree y mtry por defecto, que es equivalente a la raíz cuadrada del número de variables explicativas disponibles
set.seed(12345)
RF5<-randomForest(Sleep_Quality~., data=data_clasif_train, keep.inbag = TRUE,
nodesize=ceiling(0.05*nrow(data_clasif_train)))
print(RF5)
Call:
randomForest(formula = Sleep_Quality ~ ., data = data_clasif_train, keep.inbag = TRUE, nodesize = ceiling(0.05 * nrow(data_clasif_train)))
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 4
OOB estimate of error rate: 41.77%
Confusion matrix:
1 2 3 4 class.error
1 99 12 21 0 0.2500000
2 27 95 19 4 0.3448276
3 26 53 117 15 0.4454976
4 4 18 60 50 0.6212121Como podemos observar, el modelo da lugar a una tasa de error de 41.77% en las observaciones OOB (out of the bag), es decir, con observaciones que no ha visto o seleccionado para la creación de árboles, una especie de “set de validación”.
Es una tasa bastante elevada puesto que es casi la mitad del dataset.
Vamos a ir modificando los distintos parámetros para que tenga mejor capacidad predictiva y elimine lo máximo posible el sobreajuste.
Para estudiar el número óptimo de árboles, vamos a recurrir a la librería OOBcurve que permite hacerlo a partir de medidas diferentes a la tasa de error.
auxtask <- makeClassifTask(data = data_clasif_train, target = "Sleep_Quality")
results_5 <- OOBCurve(RF5, measures = list(kappa),
task = auxtask, data = data_clasif_train)
plot(results_5$kappa, ylab = "oob-kappa", xlab = "ntrees")
results_5$kappa[c(1,50,100,200,300,400,500)][1] 0.09867621 0.42522310 0.41745584 0.41775200 0.43373405 0.42837765 0.43072865A partir de los resultados se observa que los mejores se obtienen para unos 300 árboles pero a partir de los 50 el cambio es bastante ligero.
Vamos a probar a continuación a reducir el tamaño de hoja para determinar si se puede dar lugar a mejores resultados
set.seed(12345)
RF1<-randomForest(Sleep_Quality~., data=data_clasif_train,
keep.inbag = TRUE,
nodesize=ceiling(0.01*nrow(data_clasif_train)))
results1 <- OOBCurve(RF1, measures = list(kappa),
task = auxtask, data = data_clasif_train)
set.seed(12345)
RF05<-randomForest(Sleep_Quality~., data=data_clasif_train,
keep.inbag = TRUE,
nodesize=ceiling(0.005*nrow(data_clasif_train)))
results05 <- OOBCurve(RF05, measures = list(kappa),
task = auxtask, data = data_clasif_train)
plot(RF1$err.rate[,1], ylab = "Tasa de fallo", ylim=c(0.35,0.57))
points(RF5$err.rate[,1], col="red")
points(RF05$err.rate[,1], col="green")
plot(results1$kappa, ylab = "oob-kappa", xlab = "ntrees", ylim=c(0.1,0.5))
points(results_5$kappa, col="red")
points(results05$kappa, col="green")
Disminuir el tamaño de hoja supone una gran mejora en las métricas, existiendo diferencias además entre hojas del 1% y 0.05%, siendo mejor en téerminos generales la hoja menor.
Así mismo, el mejor valor se alcanza cuando el número de árboles está alrededor de 150.
El último parámetro que queda popr determinar es mtry (número de variables entre las que puede elegir para hacer cada una de las divisiones del árbol)
future::plan(multisession, workers=detectCores() - 1)
mtryValues<-floor(seq(2,ncol(data_clasif_train)-1, length.out=10))
nt<-150
tunningRF<- future_map(mtryValues,
~ {set.seed(12345);
randomForest(Sleep_Quality~., data=data_clasif_train,
keep.inbag = TRUE,
nodesize=ceiling(0.005*nrow(data_clasif_train)),
ntree=nt,
mtry=.x)}
)Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-1>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-2>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-3>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-4>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-5>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-6>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-7>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-8>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-9>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: package 'randomForest' was built under R version 4.5.2Warning: UNRELIABLE VALUE: Future (<unnamed-10>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".plan(sequential)
error_rate<-sapply(tunningRF, function(x) x$err.rate[nt,1])
plot(mtryValues,error_rate,type="b")
Kappas<-sapply(tunningRF, function(x) OOBCurve(x, measures = list(kappa),
task = auxtask, data = data_clasif_train)$kappa[nt])
plot(mtryValues, Kappas, type="b")
Como podemos observar, ambos indicadores coninciden en el valor óptimo de mtry de 4. Además, se comprueba que esta fuente de aleatoriedad permite resultados mejores respecto al bagging (que sería utilizando las 21 variables predictoras restantes).
Seleccionados los parámetros de 150 árboles, de tamaño mínimo de hoja de 0.05% y valor mtry = 4, creamos el modelo ganador
RFdef<-tunningRF[[which(mtryValues==4)]]
print(RFdef)
Call:
randomForest(formula = Sleep_Quality ~ ., data = data_clasif_train, keep.inbag = TRUE, nodesize = ceiling(0.005 * nrow(data_clasif_train)), ntree = nt, mtry = .x)
Type of random forest: classification
Number of trees: 150
No. of variables tried at each split: 4
OOB estimate of error rate: 36.77%
Confusion matrix:
1 2 3 4 class.error
1 97 12 23 0 0.2651515
2 16 103 22 4 0.2896552
3 16 41 126 28 0.4028436
4 1 9 56 66 0.5000000Se observa que efectivameente la tasa de error ha disminuido de casi un 42% a un 37%, clasifica por lo general mejor las clases aunque en la que más falla es en la última Excelente, donde la mitad de los casos los clasifica bien y la otra mitad no.
Obtenemos una estimación de la calidad de los términos en otros indicadores para la partición prueba
multiclass.roc(data_clasif_train$Sleep_Quality,
predict(RFdef,data_clasif_train,type="prob"))$aucMulti-class area under the curve: 0.9981barplot(t(RFdef$importance)/sum(RFdef$importance), las=2)
Respecto a la importancia de variables se observa que la más importante sigue siendo Sleep_Duration, el resto de variables, aunque importantes, están muy repartidas en porcentaje entre ellas.
Esto sucede porque la aleatoriedad añadida con mtry en los RF permite que se aproveche más el potencial predictivo de las variables “no dominantes” o menos importantes en modelos anteriores.
Finalmente, obtenemos una estimación de calidad en la partición de prueba
multiclass.roc(data_clasif_test$Sleep_Quality,
predict(RFdef,data_clasif_test,type="prob"))$aucMulti-class area under the curve: 0.8658cm_test<-confusionMatrix(table((predict(RFdef, newdata=data_clasif_test)),data_clasif_test$Sleep_Quality))
cm_test$table
1 2 3 4
1 23 1 5 0
2 3 28 11 4
3 5 6 30 13
4 1 1 6 16cm_test$overall[1] Accuracy
0.6339869 Kappa(cm_test$table) value ASE z Pr(>|z|)
Unweighted 0.5026 0.05351 9.393 5.843e-21
Weighted 0.5752 0.05209 11.042 2.389e-28cm_test$byClass[,1:2] Sensitivity Specificity
Class: 1 0.7187500 0.9504132
Class: 2 0.7777778 0.8461538
Class: 3 0.5769231 0.7623762
Class: 4 0.4848485 0.9333333Como podemos obsevar, se ha creado un modelo con una calidad decente, el AUC estimado es de 0.87, un índice Kappa de 0.5 y tasa de acierto de 63.
Respecto a la sensibilidad y especificidad observamos que todas las categorías cuentan con una especificidad mayor que sensibilidad, es decir, son mejores indicando que observación no pertenece a su categoría que indicar cuales si pero, las categorías Pobre y Aceptable cuentan con una sensibilidad superior a 0.7, que no está mal, Buena llega a 0.57 y Excelente tiene una sensibilidad ligeramente inferior a 0.5. Como veníamos indicando, es la categoría a la que más le cuesta asociar observaciones
Aún así, vamos a graficar como las distintas variables importantes influyen a cada categoría
Partiendo de la variable más importate, Sleep_Duration
predictor <- Predictor$new(RFdef, data = data_clasif_train, y = data_clasif_train$Sleep_Quality)
pdp <- FeatureEffect$new(predictor, feature = "Sleep_Duration", method="pdp")
pdp$plot()
predictor <- Predictor$new(RFdef, data = data_clasif_train, y = data_clasif_train$Sleep_Quality)
pdp <- FeatureEffect$new(predictor, feature = "Travel_Time", method="pdp")
pdp$plot()
Se observa claramente como aquellos que califican su percepción de calidad de sueño como Buena o Excelente duermen entre 6 y 8 horas mientras que aquellos que sienten que han tenido una calidad de sueño Pobre o Aceptable se encuentran mayoritariamente entre las 4 y 6 horas de sueño.
Se puede hacer algun grafiquillo más, que siempre van bien
Finalmente el último tipo de modelo que vamos a realizar es Gradient Boosting. Sigue más o menos la temática de los Random Forest pero, en vez de crear un montón de árboles distintos entre ellos, van de manera secuencial capturando nueva información. Esto logra una aproximación iterativa al optimo de una determinada función, la cual permite evaliar la capacidad predicitiva del conjunto
Para poder realizar este tipo de modelos, es necesario que la variable objetivo esté codificada como 0, 1, 2 y 3 y, además, ser de tipo numérico, por lo que hay que hacer dicha modificación previa y luego se realiza la partición
datos$Sleep_Quality<-as.numeric(datos$Sleep_Quality)-1
set.seed(12345)
trainIndex <- createDataPartition(datos$Sleep_Quality, p=0.8, list=FALSE)
data_clasif_train <- datos[trainIndex,]
data_clasif_test <- datos[-trainIndex,]De manera análoga a los modelos anteriores, comenzamos creando un GB genérico.
Como utilizamos una variable que toma más de dos niveles, en la distribución de los datos indicamos “multinomial” en vez de “bernoulli”.
Vamos a realizar 5000 árboles con una tasa de aprendizaje de 0.1, tamaño mínimo de hoja del 5%, utilizando el 75% de los datos como entrenamiento y el 25% restante para validación y una profundidad del árbol 3
set.seed(12345)
data_clasif_train_reord <- data_clasif_train[sample(nrow(data_clasif_train)),]
set.seed(12345)
gbm_clasif_5 <- gbm(
formula = Sleep_Quality~., data=data_clasif_train_reord,
distribution = "multinomial",
train.fraction = 0.75,
n.trees = 5000,
shrinkage = 0.1,
bag.fraction = 1,
interaction.depth = 3,
n.minobsinnode = ceiling(0.05*nrow(data_clasif_train_reord))
)Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.print(gbm_clasif_5)gbm(formula = Sleep_Quality ~ ., distribution = "multinomial",
data = data_clasif_train_reord, n.trees = 5000, interaction.depth = 3,
n.minobsinnode = ceiling(0.05 * nrow(data_clasif_train_reord)),
shrinkage = 0.1, bag.fraction = 1, train.fraction = 0.75)
A gradient boosted model with multinomial loss function.
5000 iterations were performed.
The best test-set iteration was 69.
There were 21 predictors of which 21 had non-zero influence.Como podemos observar, se ha construido un GBM para clasificación con 5000 árboles, donde la partición de validación ha indicado que serían sólo necesarios los 69 primeros y se han utilizado todas las variables.
Vamos a representar gráficamente la evolución de la deviance tanto en entrenamiento (negro) como en validación (rojo)
opt_clasif_gbm5<-gbm.perf(gbm_clasif_5)
abline(h=gbm_clasif_5$valid.error[opt_clasif_gbm5], lty=2, col="darkgray")
La línea vertical indica donde se ubica el valor óptimo.
Se aprecia claramente como pasada dicha línea, se obtienen peores resultados y cada vez hay mayor sobreajuste.
A continuación, vamos a construir una rejilla de parámetros a evaluar para poder compararlos. Vamos a limitar las opciones de cada parámetro a 4 para no alargar mucho los procesos pero podemos y deberíamos ampliar este número
hyper_grid <- expand.grid(
shrinkage = c(0.1, 0.05, 0.01, 0.005),
interaction.depth = c(3, 5, 7, 9),
n.minobsinnode = c(ceiling(0.01*nrow(data_clasif_train_reord)), ceiling(0.05*nrow(data_clasif_train_reord)),
ceiling(0.02*nrow(data_clasif_train_reord)), ceiling(0.005*nrow(data_clasif_train_reord)))
)Creada la rejilla, creamos un modelo para cada una de las 64 combinaciones
future::plan(multisession, workers=detectCores() - 1)
set.seed(12345)
tunningGB<- future_map(1:nrow(hyper_grid),
~ {set.seed(12345);
gbm(formula = Sleep_Quality~., data=data_clasif_train_reord,
distribution = "multinomial",
n.trees = 5000,
train.fraction = 0.75,
shrinkage = hyper_grid$shrinkage[.x],
bag.fraction = 1,
interaction.depth = hyper_grid$interaction.depth[.x],
n.minobsinnode = hyper_grid$n.minobsinnode[.x]
)}
)Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-11>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-12>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-13>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-14>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-15>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-16>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-17>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-18>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-19>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-20>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-21>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-22>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-23>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-24>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.
Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.Warning: UNRELIABLE VALUE: Future (<unnamed-25>) unexpectedly generated random
numbers without specifying argument 'seed'. There is a risk that those random
numbers are not statistically sound and the overall results might be invalid.
To fix this, specify 'seed=TRUE'. This ensures that proper, parallel-safe
random numbers are produced. To disable this check, use 'seed=NULL', or set
option 'future.rng.onMisuse' to "ignore".plan(sequential)A continuación, lo representamos gráficamente
resultado <- map_df(1:nrow(hyper_grid),
function(x) {
data.frame(
shrinkage = as.factor(hyper_grid$shrinkage[x]),
interaction.depth = as.factor(hyper_grid$interaction.depth[x]),
n.minobsinnode = as.factor(hyper_grid$n.minobsinnode[x]),
ntrees = 1:5000,
dev = tunningGB[[x]]$valid.error
)
}
)
ggplot(data=resultado, aes(x=ntrees, y=dev, col=shrinkage))+
geom_line(lwd=1.1)+
geom_hline(yintercept = min(resultado$dev))+
facet_grid(interaction.depth ~ n.minobsinnode)
resultado %>%
group_by(interaction.depth, n.minobsinnode, shrinkage) %>%
slice_min(dev, n = 1)%>%
arrange(dev)# A tibble: 64 × 5
# Groups: interaction.depth, n.minobsinnode, shrinkage [64]
shrinkage interaction.depth n.minobsinnode ntrees dev
<fct> <fct> <fct> <int> <dbl>
1 0.05 3 13 125 0.885
2 0.005 3 13 1271 0.885
3 0.01 3 13 640 0.887
4 0.05 3 31 138 0.888
5 0.01 3 31 652 0.889
6 0.1 3 31 69 0.890
7 0.005 3 31 1263 0.890
8 0.1 3 13 64 0.893
9 0.01 3 7 676 0.901
10 0.005 3 7 1347 0.902
# ℹ 54 more rowsEn la representación vemos como cada uno de los parámetros impacta en el GBM.
Observamos como el mejor modelo es el que tiene una profundidad de 3, tasa de aprendizaje de 0.05, 31 observaciones mínimo por hoja y hace uso de 138. Otra opción también buena parece la de 69 árboles con profundidad y tamaño mínimo de hoja igual al anterior pero con una tasa de aprendizeje de 0.1.
En este caso, vamos a optar por el primero e intentar reducir su complejidad con el número de árboles manteniendo la tasa de aprendizaje
optComb<-which((hyper_grid$shrinkage==0.05)&(hyper_grid$interaction.depth==3)&(hyper_grid$n.minobsinnode==31))
print(tunningGB[[optComb]])gbm(formula = Sleep_Quality ~ ., distribution = "multinomial",
data = data_clasif_train_reord, n.trees = 5000, interaction.depth = hyper_grid$interaction.depth[.x],
n.minobsinnode = hyper_grid$n.minobsinnode[.x], shrinkage = hyper_grid$shrinkage[.x],
bag.fraction = 1, train.fraction = 0.75)
A gradient boosted model with multinomial loss function.
5000 iterations were performed.
The best test-set iteration was 138.
There were 21 predictors of which 21 had non-zero influence.opt_gbm<-gbm.perf(tunningGB[[optComb]])
abline(h=tunningGB[[optComb]]$valid.error[opt_gbm], lty=2, col="gray")
Parece que con 138 o incluso 100 para redondear vamos a obtener mñas o menos los mismos resultados y no hay ninguna otra solución que permita reducir significativamente el sobreajuste
Una vez descubierta la parametrización óptima, construimos el modelo final utilizando los datos de entrenamiento al completo para mejorar la capacidad predictiva y más adelante determinar su calidad final
set.seed(12345)
gbm_clasif_opt<- gbm(
formula = Sleep_Quality~., data=data_clasif_train,
distribution = "multinomial",
train.fraction = 1,
n.trees = 130,
shrinkage = 0.05,
bag.fraction = 1,
interaction.depth = 3,
n.minobsinnode = 31
)Warning: Setting `distribution = "multinomial"` is ill-advised as it is
currently broken. It exists only for backwards compatibility. Use at your own
risk.print(gbm_clasif_opt)gbm(formula = Sleep_Quality ~ ., distribution = "multinomial",
data = data_clasif_train, n.trees = 130, interaction.depth = 3,
n.minobsinnode = 31, shrinkage = 0.05, bag.fraction = 1,
train.fraction = 1)
A gradient boosted model with multinomial loss function.
130 iterations were performed.
There were 21 predictors of which 21 had non-zero influence.Obtenemos tanto los valores de AUC como la matriz de confusión
pred_array <- predict(
gbm_clasif_opt,
data_clasif_train,
n.trees = gbm_clasif_opt$n.trees, # Asegura usar el número de árboles entrenados (130)
type = "response"
)
# 2. Extracción de la matriz 2D de probabilidades (Dimensiones: 619 x 4)
# Seleccionamos la única "capa" (índice 1) de la tercera dimensión
pred_matrix <- pred_array[,,1]
# 3. Calcular el AUC Multiclase
AUC_multiclase <- pROC::multiclass.roc(data_clasif_train$Sleep_Quality, pred_matrix)$auc
# 4. Imprimir el resultado
print(AUC_multiclase)Multi-class area under the curve: 0.9417pred_array <- predict(
gbm_clasif_opt,
data_clasif_test,
n.trees = gbm_clasif_opt$n.trees, # Asegura usar el número de árboles entrenados (130)
type = "response"
)
# 2. Extracción de la matriz 2D de probabilidades (Dimensiones: 619 x 4)
# Seleccionamos la única "capa" (índice 1) de la tercera dimensión
pred_matrix <- pred_array[,,1]
# 3. Calcular el AUC Multiclase
AUC_multiclase <- pROC::multiclass.roc(data_clasif_test$Sleep_Quality, pred_matrix)$auc
# 4. Imprimir el resultado
print(AUC_multiclase)Multi-class area under the curve: 0.8408Accuracy, kappa y sensibilidades y especificidades.
n_trees_optimo <- gbm_clasif_opt$n.trees
pred_array_test <- predict(
gbm_clasif_opt,
newdata = data_test,
n.trees = n_trees_optimo,
type = "response"
)
pred_matrix_test <- pred_array_test[,,1]
niveles_clase <- levels(data_test$Sleep_Quality)
if (is.null(niveles_clase)) {
niveles_clase <- sort(unique(data_test$Sleep_Quality))
}
colnames(pred_matrix_test) <- niveles_clase
pred_clases_test <- colnames(pred_matrix_test)[apply(pred_matrix_test, 1, which.max)]
pred_clases_test <- factor(pred_clases_test, levels = niveles_clase)
cm_test <- confusionMatrix(
data = pred_clases_test, # Clases predichas (Factor)
reference = data_test$Sleep_Quality # Clases reales (Factor)
)
cm_test$table Reference
Prediction Pobre Aceptable Buena Excelente
Pobre 23 2 5 1
Aceptable 2 29 5 2
Buena 5 5 37 10
Excelente 2 0 5 20cm_test$overall[c("Accuracy", "Kappa")] Accuracy Kappa
0.7124183 0.6084681 Kappa(cm_test$table) value ASE z Pr(>|z|)
Unweighted 0.6085 0.05020 12.12 8.081e-34
Weighted 0.6438 0.05136 12.54 4.799e-36cm_test$byClass[, c("Sensitivity", "Specificity")] Sensitivity Specificity
Class: Pobre 0.7187500 0.9338843
Class: Aceptable 0.8055556 0.9230769
Class: Buena 0.7115385 0.8019802
Class: Excelente 0.6060606 0.9416667çY para train
pred_array <- predict(
gbm_clasif_opt,
data_clasif_train,
n.trees = gbm_clasif_opt$n.trees,
type = "response"
)
pred_matrix <- pred_array[,,1]
niveles_clase_ref <- levels(data_clasif_train$Sleep_Quality)
if (is.null(niveles_clase_ref)) {
niveles_clase_ref <- sort(unique(data_clasif_train$Sleep_Quality))
data_clasif_train$Sleep_Quality <- factor(data_clasif_train$Sleep_Quality,
levels = niveles_clase_ref)
}
colnames(pred_matrix) <- niveles_clase_ref
pred_clases_train_raw <- colnames(pred_matrix)[apply(pred_matrix, 1, which.max)]
pred_clases_train <- factor(pred_clases_train_raw, levels = niveles_clase_ref)
cm_train <- confusionMatrix(
data = pred_clases_train,
reference = data_clasif_train$Sleep_Quality
)
cm_train$table Reference
Prediction 0 1 2 3
0 111 9 10 0
1 8 109 22 8
2 17 18 159 34
3 2 2 20 90cm_train$overall[c("Accuracy", "Kappa")] Accuracy Kappa
0.7576737 0.6704864 cm_train$byClass[, c("Sensitivity", "Specificity")] Sensitivity Specificity
Class: 0 0.8043478 0.9604990
Class: 1 0.7898551 0.9209979
Class: 2 0.7535545 0.8308824
Class: 3 0.6818182 0.9507187Finalemnte, observamos la importancia de las variables
summary(gbm_clasif_opt, las=2, cex.names=.5)
var rel.inf
Sleep_Duration Sleep_Duration 26.8293713
Travel_Time Travel_Time 20.7319431
Cholesterol_Level Cholesterol_Level 10.5056008
Work_Hours Work_Hours 6.0417555
Social_Interactions Social_Interactions 4.2673385
Bed_Time Bed_Time 4.1042407
Blood_Pressure Blood_Pressure 3.9788158
Wake_Up_Time Wake_Up_Time 3.8465322
Age Age 3.8336118
Blood_Sugar_Level Blood_Sugar_Level 2.9821662
Caffeine_Intake Caffeine_Intake 2.9355375
Alcohol_Intake Alcohol_Intake 2.8540974
Physical_Activity Physical_Activity 1.7370493
Exercise_Type Exercise_Type 1.2634732
Stress_Detection Stress_Detection 1.1670641
Screen_Time Screen_Time 0.7916989
Occupation Occupation 0.5984244
Meditation_Practice Meditation_Practice 0.5802062
Marital_Status Marital_Status 0.4016970
Smoking_Habit Smoking_Habit 0.2966212
Gender Gender 0.2527551En este caso, vemos como ha cambiado mucho respecto a lo que obteníamos anteriormente, siendo el nivel de colesterol la variable más importante, sequido de la hora a la que se va a la cama y la edad.