III Modelización Predictiva

Trabajando con estudio de fumadores y no fumadores

diegotriana11@gmail.com Mi Perfil en Linkedin

El siguiente informe nace como parte de las actividades de la asignatura Estadistica Avanzada del Master en Data Science en la Universitat Oberta de Catalunya, con esto busco ofrecer un blog personal de informes que permitan evidenciar el crecimiento académico que se puede ir adquiriendo con el paso del Master y como R Studio es una herramienta de gran uso.

En estos informes hag uso de lenguaje R, R markdown, HTML y CSS

#libraries
library(ggplot2)
library(GGally)
library(dplyr)
library(plotly)
library(pROC)

Descripción de la Actividad

En esta actividad se usará el fichero Fumadores_clean_5Y_1.csv ya preparado, es decir, después del preproceso (II Preprocesamiento) que se ha realizado en la primera actividad. En nuestro archivo fuente encontraremos almacenados los _datos de una investigación médica sobre la capacidad pulmonar de varias personas, con el objetivo de estudiar si los hábitos de salud y los hábitos como fumadores influencian la capacidad pulmonar. Para realizar el estudio se recogió una muestra de 300 personas. A cada persona, se le preguntó a través de un cuestionario su género, hábitos de deporte, si era fumadora, y en caso de que lo fuera, cuántos cigarrillos al día de promedio fumaba y los años que hacía que fumaba. Además, se midió la capacidad pulmonar de cada persona a partir de un test de aire expulsado, desde donde se tomó como capacidad pulmonar la medida FEF (forced expiratory flow), que es la velocidad del aire saliendo del pulmón durante la porción central de una espiración forzada. Se mide en litros / segundo. Otros datos personales recogidos son: la altura, peso y ciudad donde vive. Se incluye en el archivo una columna adicional “PC5Y” que es la capacidad pulmonar de cada persona medida al cabo de 5 años de realizar el primer test. Se asume que la persona no ha cambiado sus condiciones personales significativamente en este tiempo.

Desrcargar Set

Carga deL Archivo

fumadores<- read.csv("Fumadores_clean_5Y_1.csv", sep = "," ) #Lectura del archivo
fumadores<- data.frame(fumadores) # Casteo a dataframe

Análisis de variables y registros

Inicialmente, el archivo cuenta con 300 observaciones y 10 variables (Sex, Sport, Years, Cig, PC, City, Weight, Age, Height, PC5Y). Ahora, se imprimen los primeros registros:

head(fumadores)

Funciones

#functions
regressions<- function(coeficientes, test){
  intercepto<- if(!is.null(coeficientes[1])){coeficientes[1]}
  beta_0<- if(!is.null(coeficientes[2])){coeficientes[2]}
  beta_1<- if(!is.null(coeficientes[3])){coeficientes[3]}
  beta_2<- if(!is.null(coeficientes[4])){coeficientes[4]}
  beta_3<- if(!is.null(coeficientes[5])){coeficientes[5]}
  probabilidad<- (intercepto + (beta_0 * as.double(test$Weight)) + (beta_1 * as.double(test$Cig)) + (beta_2 * as.double(test$Years)) + (beta_3*as.double(test$Sex)) )
  beta_0
  return(probabilidad) 
}

reg_logistica<-function(modelo, test){
  coeficientes<- summary(modelo)$coefficients
  intercepto<- if(!is.null(coeficientes[1])){coeficientes[1]}
  beta_0<- if(!is.null(coeficientes[2])){coeficientes[2]}
  beta_1<- if(!is.null(coeficientes[3])){coeficientes[3]}
  beta_2<- if(!is.null(coeficientes[4])){coeficientes[4]}
  p_logi<- exp(intercepto+(beta_0*test$PC)+ (beta_1*test$Weight) + (beta_2 * as.numeric(test$Sex) )) / (1+exp(intercepto + (beta_0*test$PC) + (beta_1*test$Weight) + (beta_2*as.numeric(test$Sex) )))
  return(p_logi)
}

Modelo de regresión lineal

Inicialmente, estudiaremos la asociación entre La Capacidad Pulmonar y algunas de las caracteristicas de cada atributo.

Modelo de regresiÓn lineal mÚltiple (regresores cuantitativos)

Diagrama de DispersiÓn PC-Weight-Cig-Years

Para empezar, se aplico la funciÓn pairs() para graficar y representar a manera de matriz, la dispersión de las variables de manera general. Se puede obervar a simple vista que PC-Years y PC-Cig parecen tener una relaciÓn lineal clara, caso que no es claro entre PC-Weight. Por otra parte, se aplico la funcion Para empezar, se aplico la funciÓn cor() para calcular la matriz de coeficientes de correlacion entre las variables.

sub_fumadores<- subset(fumadores, select = c("PC","Weight","Cig","Years")) #Extración de variables de interes
pairs(sub_fumadores)

ggpairs(sub_fumadores, lower = list(continuous = "smooth"), diag = list(continuous = "bar"), axisLabels = "none")

Los coeficientes de Years y Cig muestran que ambas variables tienen coficientes mayor a |0.5|, reflejando, aparentemente, una correlacion fuerte.

De igual forma, en el siguiente apartado se procederá a detallar y explicar detenidamente si existe relacion o no, y en que grado. Para empezar, se realizara el calculo y posterior representacion de la estimacion por minimos cuadrados ordinarios en un modelo lineal entre PC y Weight.

Para cada una de las siguientes representaciones se hizo uso de la función lm() (Linear Model), la cual permite calcular y representar la recta de minimos cuadrados dada por la ecuación:
\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_nX_n \]

PC vs Weigh, Cig & Years

reg_mul1_PC <- lm(sub_fumadores$PC ~ sub_fumadores$Weight + sub_fumadores$Cig +sub_fumadores$Years, data = sub_fumadores)
summary(reg_mul1_PC) # Resumen del modelo lineal

## 
## Call:
## lm(formula = sub_fumadores$PC ~ sub_fumadores$Weight + sub_fumadores$Cig + 
##     sub_fumadores$Years, data = sub_fumadores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97891 -0.18424 -0.01939  0.19799  0.78591 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           3.677888   0.284701  12.918   <2e-16 ***
## sub_fumadores$Weight  0.001283   0.004178   0.307    0.759    
## sub_fumadores$Cig    -0.032711   0.001923 -17.008   <2e-16 ***
## sub_fumadores$Years  -0.023139   0.001583 -14.613   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2734 on 296 degrees of freedom
## Multiple R-squared:  0.812,  Adjusted R-squared:  0.8101 
## F-statistic: 426.3 on 3 and 296 DF,  p-value: < 2.2e-16

r2_wei<-summary(reg_mul1_PC)$r.squared * 100  #R cuadrado - bondad del ajuste
r2_wei_res<- 100-r2_wei # valor no explicado por el modelo

Como se puede observar, tanto el intercepto como la pendiente en Cig y Years tienen una influencia significativa a más del 99.9% (***), caso contratio, la pendiente (Weight) no presenta influencia significativa. La bondad del ajuste (R cuadrado) es de 0.8120466, la varaibilidad de la variable dependiente (PC) con respecto a la variable independiente (Weight) es de 81.2046647 % , revelando que el 18.7953353 % no esta siendo explicado por el modelo.

La ecuación de la recta de minimos cuadrados que presenta la dependencia/relación entre las variables PC y Weight es:

\[ Y = 3.67788 + 0.0012 Weight + (-0.03271 Cig) + (-0.02313 Years) \]

Como se pudó observar, tanto la relación PC-Cig y PC-Years, presentan valores atipicos, que deben pertenecer a observaciones exepcionales situadas en el punto 0 de la variable independiente.

La variable Height no debe ser tenida en cuenta porque puede presentar problemas de Multicolinealidad en la regresión, incrementando la varianza de los coeficientes de regresion, logrando que estos sean inestables.

Revisar colinealidad http://www.hrc.es/bioest/Reglin_15.html uuu

Modelo de regresión lineal múltiple (regresores cuantitativos y cualitativos)

Para empezar, se analizaron las variables cuantitativas mediante el diagrama de dispersion multiple y la relacion entre variables. Posteriormente, se procede a estudiar las variables categoricas con el fin de generar un boxplot con sus relaciones que permitan identificar la influencia en la variable dependiente. Se utilizó la función relevel() para reordenar las variables acorde al valor de la referencia.

sexR<- relevel(fumadores$Sex, ref = "F") # Reordenación
sportR<- relevel(fumadores$Sport, ref="N") # Reordenación
sub_fumadores2<- cbind(sub_fumadores, sexR, sportR) # Combinar set con variables
pairs(sub_fumadores2) # graficar

En los siguientes gráficos se obsevan la relación entre las Variable PC-Sex y PC-Sport(variable numerica y nominal). Para este apartado se ha hecho uso de lo que hoy por hoy es una de las mejores librerías de graficos en r y python ggplot2()

ggplot(data = sub_fumadores2, mapping=aes(x = sub_fumadores2$sexR, y = sub_fumadores2$PC,
                                          color=sub_fumadores2$sexR)) +geom_boxplot() + 
  geom_jitter(width = 0.1) + theme_bw() + theme(legend.position = "none")

ggplot(data = sub_fumadores2, mapping=aes(x = sub_fumadores2$sportR, y = sub_fumadores2$PC, color=sub_fumadores2$sportR)) + geom_boxplot() + geom_jitter(width = 0.1) + theme_bw() + theme(legend.position = "none")

A continuación, se aplica el modelo de regresión lineal multiple, en este, se tienen en cuenta 5 variables (Weight, Cig, Years, Sex, Sport y la variable a explicar, PC).

reg_PC_multiple <- lm(sub_fumadores2$PC ~ fumadores$Weight+ fumadores$Cig + fumadores$Years + sexR + sportR, data = fumadores)# Resultado Regresion
summary(reg_PC_multiple) #Resumen de la regresión

## 
## Call:
## lm(formula = sub_fumadores2$PC ~ fumadores$Weight + fumadores$Cig + 
##     fumadores$Years + sexR + sportR, data = fumadores)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73530 -0.10739 -0.00663  0.10509  0.48296 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.561531   0.236800  15.040  < 2e-16 ***
## fumadores$Weight -0.000601   0.003599  -0.167 0.867511    
## fumadores$Cig    -0.034035   0.001293 -26.315  < 2e-16 ***
## fumadores$Years  -0.022633   0.001065 -21.249  < 2e-16 ***
## sexRM             0.102336   0.027373   3.739 0.000223 ***
## sportRE           0.565597   0.032939  17.171  < 2e-16 ***
## sportRR           0.370267   0.031566  11.730  < 2e-16 ***
## sportRS           0.199592   0.025917   7.701 2.11e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1827 on 292 degrees of freedom
## Multiple R-squared:  0.9172, Adjusted R-squared:  0.9152 
## F-statistic: 462.2 on 7 and 292 DF,  p-value: < 2.2e-16

Como se puede observar en el anterior resumen, tanto el coeficiente de determinacion 0.9172211, como el coeficiente ajustado 0.9152367 de este modelo dee regresion multiple es mayor a los coeficientes de las regresiones anteriores.

Variable	Coef. Deter.	Coef. Deter. Ajustado
Modelo 1	0.8120466	0.8101417
Modelo 2	0.9172211	0.9152367

A manera comparativa, la regresion multiple del modelo 2 permite identificar un mejor ajuste, ya que obtiene un valor mayor y mas cercano a 1, ajustandose mejor al modelo lineal. Es decir, el presente modelo es capaz de explicar un 91.722% de la variabilidad de PC. De igual forma, el coeficiente de determinancion ajustado, con un valor de 91.524% nos indica que este modelo cuenta con predictores de gran utilidad. Otro aspecto importante de analizar es el pvalor, el cual, para todas las variables excepto weight (que no contribuye de forma significativa), tienen una influencia significativa de mas del 99%, identificado con (***).

Predicción de la capacidad pulmonar con los dos modelos

test<-data.frame(Age=c("30"),Sex=c("M"),Years=c("15"),Cig=c("10"),Sport=c("R"), City=("Lleida"),Height=c("175"),Weight=("68"), PC=c("")) # Sujeto de prueba
PC_test1<- regressions(coef(summary(reg_mul1_PC)), test) ##llamado a la funcion para predecir
# la funcion creada se hizo manual con fines de comprender la ecuación lineal
PC_test2<- regressions(coef(summary(reg_PC_multiple)), test)

Para el modelo 1, el PC hallado es 3.9080221, mientras que para el modelo 2 el PC hallado es 3.6065991

Modelo de regresión logística

Para el modelo de regresion logistica, la funcion esta descrita por:

\[ p = \frac{ e^{\beta_0 +\beta_1X} } { 1+ e^{\beta_0 +\beta_1X} } \]

Estimación de un modelo de regresión logística

El primer paso será crear una variable binaria (smoker) que indique la condición de fumador (smoker = 1) o no fumador (smoker = 0).

fumadores$Smoke[fumadores$Years>0]<-1 # nueva variable Smoke con valores en 1 si years es mayor a 0
fumadores$Smoke[fumadores$Years<1]<-0 # nueva variable Smoke con valores en 0 si years es menor a 0

Despues, a traves del modelo de la funcion glm() se crea el modelo logistico:

reg_logi<- glm(formula = fumadores$Smoke ~ fumadores$PC + fumadores$Weight + sexR, data = fumadores, family = binomial())
summary(reg_logi)

## 
## Call:
## glm(formula = fumadores$Smoke ~ fumadores$PC + fumadores$Weight + 
##     sexR, family = binomial(), data = fumadores)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.47787  -0.31210  -0.04843   0.04271   3.13797  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      38.51385    7.64106   5.040 4.65e-07 ***
## fumadores$PC     -9.50804    1.36436  -6.969 3.20e-12 ***
## fumadores$Weight -0.09876    0.08155  -1.211    0.226    
## sexRM             0.76825    0.64182   1.197    0.231    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 411.06  on 299  degrees of freedom
## Residual deviance: 107.87  on 296  degrees of freedom
## AIC: 115.87
## 
## Number of Fisher Scoring iterations: 8

Como se puede oservar, el regresor PC tiene una influencia significativa en el modelo, a mas del 99.9% (***). De igual manera, se evidencia que el regresor es negativo, lo que implica que un aumento en la capacidad pulmonar provocara una disminucion en la probabilidad de la respuesta. Por otra parte, se puede evidenciar que tanto weight como sex no influyen significativamente en el suceso.

Prediccion en el modelo lineal generalizado (modelo de regresión logística)

Ahora, creo mi set de prueba, para este caso es hombre que tiene una capacidad pulmonar de 3.75 l/s, un peso de 68 kg y altura de 175 cm; seguido envio este nuevo set a la funcion de funcion de probabilidad

test_log<-data.frame(PC=c(3.75),Sex=c("M"),Height=c(175),Weight=(68)) #sujeto de prueba
pro_test<- reg_logistica(reg_logi,test_log ) # resultado regresion

La probabilidad de que una persona hombre, con PC 3.75 l/s, un peso de 68 kg y altura de 175 cm es 4.3561991 %

Mejora del modelo

reg_logi_2<- glm(formula = fumadores$Smoke ~ fumadores$PC + fumadores$Weight + sexR + fumadores$Age, data = fumadores, family = binomial())
reg_logi_3<- glm(formula = fumadores$Smoke ~ fumadores$PC + fumadores$Weight + sexR + sportR, data = fumadores, family = binomial())
reg_logi_4<- glm(formula = fumadores$Smoke ~ fumadores$PC + fumadores$Weight + sexR + fumadores$Age + sportR, data = fumadores, family = binomial())

Modelo	AIC
Age	66.6790289
SportR	54.2610429
Age + SportR	41.5586883

ANALISIS AIC

Calidad del ajuste Falsos Negativos y Falsos Positivos

Se calculará lamatriz de confusión del mejor modelo del mejor modelo anteriormente desarrollado, se supondrá un umbral de discriminación del 70 %. Observaremos cuantos falsos negativos hay e interpretaremos qué es un __ falso negativo__ en este contexto. Hacer lo mismo con los falsos positivos.

Para iniciar, se aplico la funcion predict() que permite invocar el modelo de regresion y el dataset de prueba -para este caso es el mismo dataset, pero alguans fuentes consideran importante separar entrenamiento y prueba para evitar cesgo-, seguido se selecciono del modelo la prediccion en la cual el resultado fue mayor al umbral ( 70% ).

modelo_predictivo <- predict(reg_logi_4,fumadores,type = "response")
discriminacion <- ifelse(test = modelo_predictivo > 0.7,yes= 1,no = 0)
obs_pre <- data.frame(observaciones=fumadores$Smoke,predicciones= discriminacion)
obs_positivas<- sum(obs_pre$observaciones == 1)
obs_negativas<- sum(obs_pre$observaciones == 0)
pre_positivas<- sum(obs_pre$predicciones == 1)
pre_negativas<- sum(obs_pre$predicciones == 0)
#data.frame(obs_positivas, obs_negativas,pre_positivas,pre_negativas)

Como se puede observar, a un umbral del 70% nuestro modelo regresor aparentemente ha clasificado correctamente las personas de la muestra, con una diferencia de 2 individuos. Sin embargo, a continuacion se realiza el analisis de falsos positivos y falsos negativos:

verdaderos_pos<- sum(obs_pre$observaciones == 1 & obs_pre$predicciones == 1) # Verdaderos Positivos
verdaderos_neg<- sum(obs_pre$observaciones == 0 & obs_pre$predicciones == 0) # Verdaderos Negativos
falsos_pos<- sum(obs_pre$observaciones == 0  & obs_pre$predicciones == 1) # Falsos Positivos
falsos_neg<- sum(obs_pre$observaciones == 1 & obs_pre$predicciones == 0) # Falsos Negativos
sensibilidad<- round((verdaderos_pos/(verdaderos_pos+falsos_neg)*100),2) #Sensibilidad
especificidad<- round((verdaderos_neg/(verdaderos_neg+falsos_pos)*100),2) #Especificidad

La anterior matriz, permite comparar las observaciones de la muestra con las predicciones resultantes del modelo. Para empezar, tenemos la siguiente matriz de confusión:

Matriz	Obs Positiva	Obs Negativa
Prediccion Positiva	127	2
Prediccion Negativa	4	167

Se puede observar que actualmente el modelo presenta 127 verdaderos positivos, es decir, la cantidad de muestras que el modelo acerto cuando las observaciones eran positivas. De igual forma, se identifican los 167 verdaderos negativos, que son las muestras de nofumadores que el modelo acerto cuando las observaciones indicaban que no eran fumadores. Sin embargo, el modelo tamb?en a presentado valores contradictorios, como es el caso de los 2 falsos positivos que son las muestras que el modelo identifica como fumadores pero la observaci?n no. Lo mismo pasa con los 4 falsos negativos , que son las muestras que el modelo identifica como nofumadores pero la observaci?n si lo hace.

Por otra parte, la proporcion de elementos no relevantes correctamente identificados fue de 98.82% (especificidad), mientras que la proporcion de elementos relevantes correctamente identificados fueron 96.95% (sensibilidad).

La selección de los individuos fumadores

Ahora, se establece un nivel de probabilidad (umbral de discriminación) a partir del cual se considera que un individuo tiene muchas posibilidades de ser un fumador. Despues, se comparará el nivel de probabilidad que da el modelo con el valor de capacidad pulmonar (PC) del individuo. Finalmente se identificarán los individuos que no se comportan según lo esperado, es decir que tienen elevada capacidad pulmonar y el modelo los clasifica como fumadores y reportar los valores de probabilidad de ser fumador y de PC. Utilizar como umbral para declarar un individuo con PC elevado el cuartil tercero de la variable PC.

Nivel de probabilidad = 70%

fumadores$Prediccion <- ifelse(test = modelo_predictivo > 0.70,modelo_predictivo,0)
plot_ly(data = fumadores, x = ~fumadores$PC, y = ~fumadores$Prediccion, name= 'Fumadores', color = ~fumadores$Smoke, colors = "Set1") %>%
   layout(title = 'Capacidad Pulmonar vs Prediccion', xaxis = list(title = "Capacidad Pulmonar"), yaxis =list(title = "Prediccion"))  

quantile(fumadores$PC, 0.75)

##     75% 
## 3.79325

anormalidades <- subset(fumadores, subset = (fumadores$PC >= quantile(fumadores$PC, 0.75)) & round(fumadores$Prediccion) ==1)

En la anterior grafica, se ha hecho uso de la libreria plotly con el fin de permitir interactuar con cada una de las muestras, conocer su PC y su clasificaci?n. el umbral ha sido el tercer cuartil de PC con un valor de 3.79 En la grafica se alcanzan a observar, en valor 1 y color rojo, los 2 falsos positivos; mientras en la parte inferior (0) se identifican los 4 falsos negativos.

Con un umbral predictivo de 70% y un umbral de capacidad pulmonar superior al tercer cuantil 3.79325 se tienen 1 muestra que no se comportan segun lo esperado.

head(anormalidades)

Curva ROC

Para la curva ROC, se ha hecho uso de la libreria pROC, a traves de esta se ha usado la funcion roc.

roc_ob<- roc(fumadores$Smoke, round(fumadores$Prediccion))
plot(roc_ob)

Como se puede observar, el grafico presenta un un modelo discriminante, ya que cuenta con un ?rea de 0.9788157. Para terminos de comprender mejor lo que representa el roc y la auc, se procede a analizar los modelos anteriores del punto 2.3

modelo_predictivo_A <- predict(reg_logi_2,fumadores,type = "response")
fumadores$modelo_A <- ifelse(test = modelo_predictivo_A > 0.70,modelo_predictivo_A,0)
roc_ob_A<- roc(fumadores$Smoke, round(fumadores$modelo_A))

modelo_predictivo_B <- predict(reg_logi_3,fumadores,type = "response")
fumadores$modelo_B <- ifelse(test = modelo_predictivo_B > 0.70,modelo_predictivo_B,0)
roc_ob_B<- roc(fumadores$Smoke, round(fumadores$modelo_B))

plot(roc_ob_A, colorize = TRUE, col = c("pink"))
plot(roc_ob_B, add = TRUE, colorize = TRUE, col = c("red"))
plot(roc_ob, add = TRUE, colorize = TRUE, col = c("blue"))

Modelo	AUC	Color
Age	0.9559149	Rosado
SportR	0.9817742	Rojo
Age + SportR	0.9788157	Azul

Referencias

Regression

Relevel

Curso R

Matriz de Confusion

Graficos Plotly