Taller 7 Regresión Logística
Julian Moreno - Andres Arango
2022-06-24
R Markdown
#Librerias
library("easypackages")
lib_req<-c("MASS","visdat","hardhat","car","caret","pROC","class")
easypackages::packages(lib_req) ## Loading required package: MASS## Loading required package: visdat## Loading required package: hardhat## Loading required package: car## Loading required package: carData## Loading required package: caret## Loading required package: ggplot2## Loading required package: lattice## Loading required package: pROC## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var## Loading required package: class## All packages loaded successfullyEstructura inicial de los datos
De acuerdo a la estructura que se presenta a continuación, se observa que originalmente las variables “Acc” y “Sexo”, que son dicotómicas, estan declaradas como númericas por lo cual vamos a proceder a aplicar transformación para definirlas como factor.
#Cargue dataset
setwd("C:/Users/julia/OneDrive/Desktop/Maestria/Semestre 1/Metodos cuantitativos/6 Regresion Logistica/Taller 7")library(readxl)
accidentes <- read_excel("accidentes.xlsx")
str(accidentes)## tibble [35 x 5] (S3: tbl_df/tbl/data.frame)
## $ Acc : num [1:35] 0 0 0 1 0 1 0 1 0 1 ...
## $ Exp : num [1:35] 10 15 7 1 10 2 8 20 18 4 ...
## $ Edad: num [1:35] 30 40 25 21 29 20 40 25 43 23 ...
## $ Pot : num [1:35] 90 85 95 145 70 120 95 135 85 110 ...
## $ Sexo: num [1:35] 1 1 1 2 1 2 1 2 1 2 ...Luego de analizar los datos anteriores, se pudo detectar que los datos de las lineas 8, 13 y 16 son inconsistentes dado que la edad de los individuos sería menor de 10 años cuando aprendieron a manejar, razón por la cual los excluimos del análisis
accidentes[c(8,13,16),]## # A tibble: 3 x 5
## Acc Exp Edad Pot Sexo
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 20 25 135 2
## 2 1 15 24 110 2
## 3 0 15 25 90 2A continuación se muestra el dataset sin los datos inconsistentes detectados en el punto anterior
accidentes <- accidentes[c(-8,-13,-16),]Estructura de datos transformados
Se aplica transformación a las variables dicotómicas “Acc” y “Sexo” quedando definidas como factor.
accidentes = transform(accidentes,Sexo=factor(Sexo,levels=1:2,labels=c("Mujer","Hombre")))
accidentes = transform(accidentes,Acc=factor(Acc,levels=0:1,labels=c("No","Si")))
accidentes## Acc Exp Edad Pot Sexo
## 1 No 10 30 90 Mujer
## 2 No 15 40 85 Mujer
## 3 No 7 25 95 Mujer
## 4 Si 1 21 145 Hombre
## 5 No 10 29 70 Mujer
## 6 Si 2 20 120 Hombre
## 7 No 8 40 95 Mujer
## 8 No 18 43 85 Mujer
## 9 Si 4 23 110 Hombre
## 10 No 12 35 90 Mujer
## 11 No 14 32 95 Mujer
## 12 Si 3 23 90 Hombre
## 13 No 10 27 80 Hombre
## 14 No 11 29 95 Mujer
## 15 Si 5 24 150 Hombre
## 16 No 8 26 90 Mujer
## 17 No 6 25 85 Mujer
## 18 No 7 37 90 Mujer
## 19 No 15 38 85 Mujer
## 20 Si 5 24 95 Hombre
## 21 Si 8 23 130 Mujer
## 22 No 12 31 90 Mujer
## 23 No 10 27 105 Mujer
## 24 Si 4 45 120 Hombre
## 25 Si 9 56 110 Mujer
## 26 No 9 32 90 Mujer
## 27 No 10 29 90 Mujer
## 28 Si 8 34 150 Hombre
## 29 Si 5 45 110 Hombre
## 30 No 20 46 95 Hombre
## 31 Si 9 26 90 Mujer
## 32 Si 9 26 100 Mujerstr(accidentes)## 'data.frame': 32 obs. of 5 variables:
## $ Acc : Factor w/ 2 levels "No","Si": 1 1 1 2 1 2 1 1 2 1 ...
## $ Exp : num 10 15 7 1 10 2 8 18 4 12 ...
## $ Edad: num 30 40 25 21 29 20 40 43 23 35 ...
## $ Pot : num 90 85 95 145 70 120 95 85 110 90 ...
## $ Sexo: Factor w/ 2 levels "Mujer","Hombre": 1 1 1 2 1 2 1 1 2 1 ...#Acc=factor(Acc,ordered=(c("No","Si"))) en caso que noVisualización de Datos Faltantes
La prueba de datos faltantes arroja que no hay datos faltantes por lo cual se pude proseguir con los análisis requeridos del conjunto de datos.
#Prueba de datos faltantes
windows(height=10,width=15)
visdat::vis_miss(accidentes)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Convención de color
Hemos decidido aplicar color negro para representar los datos “No accidente” y color rojo para los “Accidente”
#Convención de color
color = ifelse(accidentes$Acc=="No","Black","Red")Análisis exploratorio de las variables
Como parte del proyecto de analítica, procedemos a analizar las variables del estudio de manera independiente aplicando los gráficos que se muestran a continuación.En este caso se comparan todas las variables contra la accidentalidad
Se establece como la variable a predecir “Acc”
#Grafico de dispersión
windows(width=10, height=7)
par(mfrow=c(2,3))
# Ajuste del modelo lineal
attach(accidentes)
plot(Exp,as.numeric(Acc)-1,col=color,pch=20,ylab="Acc",ylim=c(-0.5,1.5),xlim=c(0,20),
cex=2, main="Modelo de Regresión Lineal")
mod_lin_Exp=lm((as.numeric(Acc)-1)~Exp)
abline(h=c(0,1), lty=2,col="gray")
abline(mod_lin_Exp,lty=2,lwd=2,col="Blue")
summary(mod_lin_Exp)##
## Call:
## lm(formula = (as.numeric(Acc) - 1) ~ Exp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6164 -0.3240 0.0298 0.3105 0.6029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.05487 0.15767 6.690 2.06e-07 ***
## Exp -0.07308 0.01598 -4.572 7.77e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3894 on 30 degrees of freedom
## Multiple R-squared: 0.4107, Adjusted R-squared: 0.391
## F-statistic: 20.91 on 1 and 30 DF, p-value: 7.77e-05plot(Edad,as.numeric(Acc)-1,col=color,pch=20,ylab="Acc",ylim=c(-0.5,1.5),xlim=c(20,60),
cex=2, main="Modelo de Regresión Lineal")
mod_lin_Edad=lm((as.numeric(Acc)-1)~Edad)
abline(h=c(0,1), lty=2,col="gray")
abline(mod_lin_Edad,lty=2,lwd=2,col="Blue")
plot(Pot,as.numeric(Acc)-1,col=color,pch=20,ylab="Acc",ylim=c(-0.5,1.5),xlim=c(70,150),
cex=2, main="Modelo de Regresión Lineal")
mod_lin_Pot=lm((as.numeric(Acc)-1)~Pot)
abline(h=c(0,1), lty=2,col="gray")
abline(mod_lin_Pot,lty=2,lwd=2,col="Blue")
evento <- accidentes[Acc=="Si",]
tabla <- prop.table(table(evento$Sexo))
coord <- barplot(tabla,col=c("gray","blue"), ylim=c(0,1.1), main="Proporción de accidentados")
text(coord,tabla,labels=round(tabla,2), pos=3)
mujeres <- accidentes[Sexo=="Mujer",]
tabla <- prop.table(table(mujeres$Acc))
coord <- barplot(tabla,col=c("gray","blue"), ylim=c(0,1.1), main="Indice accidentalidad Mujeres")
text(coord,tabla,labels=round(tabla,2), pos=3)
hombres <- accidentes[Sexo=="Hombre",]
tabla <- prop.table(table(hombres$Acc))
coord <- barplot(tabla,col=c("gray","blue"), ylim=c(0,1.1), main="Indice accidentalidad Hombres")
text(coord,tabla,labels=round(tabla,2), pos=3)
detach(accidentes)Observaciones
- Se aplica regresión líneal para comparar las variables numéricas
- Se aplica gráfico de barras para visualizar la variable sexo por estar definida como factor.
- Se observa que a mayor experiencia se presentan menos accidentes.
- A mayor edad, se presentan menos accidentes.
- A mayor potencia del vehículo, se presentan más accidentes.
- El 69% de los accidentes son de hombres mientras que el 31% de mujeres
- Del total de hombres, el 82% ha sufrido un accidente en el conjunto de datos.
- Del total de mujeres, el 19% ha sufrido un accidente.
- Los modelos de regresión líneal aplicados a las variables numéricas presentan el inconveniente de que pueden tener valores mayores que 1 y menores que cero, lo cual no tiene sentido al momento de hacer predicción de la ocurrencia de un evento.
Modelo logístico
Dado que el modelo de regresión líneal presenta valores mayores que 1 y menores que 0, se procede a aplicar el modelo logístico que por su expresión garantiza valores de probabilidad entre 0 y 1.
par(mfrow=c(1,2))
#Modelo logístico Exp
attach(accidentes)
plot(Exp,as.numeric(Acc)-1,col=color,pch=20,ylab="Prob(Accidente)",ylim=c(-0.5,1.5),
xlim=c(0,20),cex=2,, main="Modelo de Regresión Logística")
mod_log=glm(Acc~Exp, family = "binomial") #glm es generalice el modelo lineal
points(Exp,predict(mod_log,type = "response"),pch=20,col="Blue")
abline(h=c(0,1), lty=2,col="gray")
abline(h=0.5, col="gray")
detach(accidentes)
summary(mod_log)##
## Call:
## glm(formula = Acc ~ Exp, family = "binomial")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6570 -0.5784 -0.1058 0.4417 1.6239
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.256 2.086 2.52 0.01175 *
## Exp -0.696 0.254 -2.74 0.00614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 24.014 on 30 degrees of freedom
## AIC: 28.014
##
## Number of Fisher Scoring iterations: 6#Modelo logístico Exp~Sexo
color2 = ifelse(accidentes$Sexo=="Mujer","Orange","Blue")
attach(accidentes)
plot(Exp,as.numeric(Acc)-1,col=color,pch=20,ylab="Prob(Accidente)",ylim=c(-0.5,1.5),
xlim=c(0,20),cex=2,xlab = "Exp~Sexo", main="Modelo de Regresión Logística")
mod_log10=glm(Acc~Exp+Sexo, family = "binomial") #glm es generalice el modelo lineal
points(Exp,predict(mod_log10,type = "response"),pch=20,col=color2)
abline(h=c(0,1), lty=2,col="gray")
abline(h=0.5, col="gray")
detach(accidentes)En el modelo 2 las mujeres se representan con puntos naranja y los hombres con puntos azules
Modelo ajustado
Modelo_RL1 <- glm(formula= Acc~Exp, data=accidentes,family = "binomial") # Ajusta un Modelo de regresión logistica simple
Modelo_RL10 <- glm(formula= Acc~Exp+Sexo, data=accidentes,family = "binomial")
summary(Modelo_RL1)##
## Call:
## glm(formula = Acc ~ Exp, family = "binomial", data = accidentes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6570 -0.5784 -0.1058 0.4417 1.6239
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.256 2.086 2.52 0.01175 *
## Exp -0.696 0.254 -2.74 0.00614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 24.014 on 30 degrees of freedom
## AIC: 28.014
##
## Number of Fisher Scoring iterations: 6summary(Modelo_RL10)##
## Call:
## glm(formula = Acc ~ Exp + Sexo, family = "binomial", data = accidentes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3013 -0.5392 -0.1463 0.2889 1.7643
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4988 2.4780 1.412 0.1580
## Exp -0.5354 0.2809 -1.906 0.0567 .
## SexoHombre 1.9451 1.3711 1.419 0.1560
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 21.893 on 29 degrees of freedom
## AIC: 27.893
##
## Number of Fisher Scoring iterations: 6Los resultados arrojados por los modelos 1 y 2, evidencian que en el modelo 1 tanto el intersepto como la variable Exp son significantes, mientras que en el modelo 2 a pesar de que su significancia es menor, presenta un valor menor de AIC.
Interpretación de los coeficientes - odds y odds ratio
Según los resultados del análisis de los coeficientes Odds y Odds ratio, se puede evidenciar que los bajos niveles de experiencia están asociados con altas probabilidades de accidentes
Adicionalmente, es importante destacar, que un incremento en la experiencia desfavorece la ocurrencia de accidentes
Tabla_coef=cbind(Coef=coef(Modelo_RL1),IC=confint(Modelo_RL1,level=0.95)) # Coeficientes ## Waiting for profiling to be done...Tabla_odds = exp(Tabla_coef)
colnames(Tabla_odds)[1]="e-beta"
Tabla_coef10=cbind(Coef=coef(Modelo_RL10),IC=confint(Modelo_RL10,level=0.95))## Waiting for profiling to be done...Tabla_odds10 = exp(Tabla_coef10)
colnames(Tabla_odds10)[1]="e-beta"
Tabla_coef;Tabla_odds## Coef 2.5 % 97.5 %
## (Intercept) 5.2563249 1.999345 10.5324024
## Exp -0.6959545 -1.334606 -0.3008925## e-beta 2.5 % 97.5 %
## (Intercept) 191.7753971 7.3842164 3.751148e+04
## Exp 0.4985983 0.2632619 7.401574e-01Tabla_coef10;Tabla_odds10## Coef 2.5 % 97.5 %
## (Intercept) 3.4988244 -0.3346055 9.3976594
## Exp -0.5353588 -1.2316418 -0.1357385
## SexoHombre 1.9450648 -0.6869887 5.1125988## e-beta 2.5 % 97.5 %
## (Intercept) 33.0765429 0.7156203 1.206012e+04
## Exp 0.5854592 0.2918131 8.730709e-01
## SexoHombre 6.9940847 0.5030887 1.661015e+02Ecuaciones del pronóstico
A continuación se presentan las 2 ecuaciones del pronóstico de acuerdo a los calculos de la tabla de coeficientes calculados en el punto anterior.
P(accidente)=1/(1+exp(-5.256+0.696Exp)) P(accidente)=1/(1+exp(-3.4988+0.5354Exp-1.9451Sexo))
Ajuste del Modelo - test razón verosimilitud
Según el análisis Anova, los resultados arrojan para el modelo 1, nivel de significancia menor que 0.05, por lo cual se determina que es considerablemente significante por lo cual se rechaza la hipotesis nula lo cual implica que ambos modelos son mejores que no tenerlos.
anova(Modelo_RL1,test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Acc
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 31 43.230
## Exp 1 19.216 30 24.014 1.167e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(Modelo_RL10,test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Acc
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 31 43.230
## Exp 1 19.2158 30 24.014 1.167e-05 ***
## Sexo 1 2.1209 29 21.893 0.1453
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análisis de bondad de ajuste
Los resultados arrojados en el análisis de bondad de ajuste presentado a continuación, muestran que los valores de los parámetros AIC, deviance y R cuadrado son mejores para el modelo 2, lo cual significa un mayor nivel de precisión del pronóstico del modelo para el conjunto de datos.
R2=function(Modelo){
R2=(Modelo$null.deviance - Modelo$deviance)/Modelo$null.deviance}
Bondad_Ajuste=data.frame(AIC=Modelo_RL1$aic,deviance=Modelo_RL1$deviance,R2=R2(Modelo_RL1))
row.names(Bondad_Ajuste)="Modelo_RL1"
Bondad_Ajuste## AIC deviance R2
## Modelo_RL1 28.0139 24.0139 0.444505R22=function(Modelo10){
R22=(Modelo10$null.deviance - Modelo10$deviance)/Modelo10$null.deviance}
Bondad_Ajuste10=data.frame(AIC=Modelo_RL10$aic,deviance=Modelo_RL10$deviance,R22=R22(Modelo_RL10))
row.names(Bondad_Ajuste10)="Modelo_RL2"
Bondad_Ajuste10## AIC deviance R22
## Modelo_RL2 27.89302 21.89302 0.4935657Predicción
A continuación se muestran las probabilidades de la ocurrencia del evento Accidente calculadas con ambos modelos para cada individuo según los valores de las variables aplicadas en cada modelo.
Estas probabilidades serán el insumo para una futura clasificación de acuerdo a un determinado punto de corte, lo cual se hará más adelante.
Score_RL1 = predict(Modelo_RL1) # Evalua los Score estimados para el modelo logistico
Prob_RL1 = predict(Modelo_RL1,type = "response")
Score_RL10 = predict(Modelo_RL10)
Prob_RL10 = predict(Modelo_RL10,type = "response")
Prob_RL1;Prob_RL10## 1 2 3 4 5 6
## 0.1540452138 0.0055798824 0.5949926008 0.9896500573 0.1540452138 0.9794557498
## 7 8 9 10 11 12
## 0.4227948123 0.0006950333 0.9221918446 0.0433086397 0.0111286914 0.9596300244
## 13 14 15 16 17 18
## 0.1540452138 0.0832356878 0.8552706548 0.4227948123 0.7466071510 0.5949926008
## 19 20 21 22 23 24
## 0.0055798824 0.8552706548 0.4227948123 0.0433086397 0.1540452138 0.9221918446
## 25 26 27 28 29 30
## 0.2675153779 0.2675153779 0.1540452138 0.4227948123 0.8552706548 0.0001728758
## 31 32
## 0.2675153779 0.2675153779## 1 2 3 4 5 6
## 0.135314544 0.010649267 0.438146616 0.992670782 0.135314544 0.987545880
## 7 8 9 10 11 12
## 0.313448503 0.002155372 0.964512896 0.050908205 0.018053468 0.978913550
## 13 14 15 16 17 18
## 0.522559915 0.083929006 0.940871584 0.313448503 0.571181060 0.438146616
## 19 20 21 22 23 24
## 0.010649267 0.940871584 0.313448503 0.050908205 0.135314544 0.964512896
## 25 26 27 28 29 30
## 0.210917308 0.210917308 0.135314544 0.761517750 0.940871584 0.005151579
## 31 32
## 0.210917308 0.210917308Clasificación inicial de las probabilidades
Para realizar una clasificación, inicialmente tomamos arbitrariamente el punto de corte 0.5, asumiendo que toda probabilidad mayor a este valor determina la ocurrencia del evento, y la no ocurrencia para valores menores.
par(mfrow=c(1,2))
hist(Prob_RL1, main="Histograma de probabilidades Exp")
abline(v=0.5, lwd=2, lty=2, col="red")
hist(Prob_RL10, main="Histograma de probabilidades Exp~Sexo")
abline(v=0.5, lwd=2, lty=2, col="red")
# Convertir probabilidades predichas en clasificaciones
pc=0.5
Class_RL1 = as.factor(ifelse(Prob_RL1>pc,"Si","No"))
Class_RL10 = as.factor(ifelse(Prob_RL10>pc,"Si","No"))
data.frame(accidentes$Acc,Class_RL1,Class_RL10)## accidentes.Acc Class_RL1 Class_RL10
## 1 No No No
## 2 No No No
## 3 No Si No
## 4 Si Si Si
## 5 No No No
## 6 Si Si Si
## 7 No No No
## 8 No No No
## 9 Si Si Si
## 10 No No No
## 11 No No No
## 12 Si Si Si
## 13 No No Si
## 14 No No No
## 15 Si Si Si
## 16 No No No
## 17 No Si Si
## 18 No Si No
## 19 No No No
## 20 Si Si Si
## 21 Si No No
## 22 No No No
## 23 No No No
## 24 Si Si Si
## 25 Si No No
## 26 No No No
## 27 No No No
## 28 Si No Si
## 29 Si Si Si
## 30 No No No
## 31 Si No No
## 32 Si No NoLa tabla anterior muestra el valor de la variable a predecir para cada individuo y el valor de la predicción arrojado por cada modelo, para tener una mejor idea de cual de los dos modelos estuvo más acertado, vamos a proceder a hacer la matriz de confusión.
Evaluación de la bondad de la clasificación
Modelo 1
Para el modelo 1, el análisis de bondad muestra que este es más específico que sensible, lo cual indica que es más acertado para pronósticar los no accidentados.
Respecto al indicador de precisión, a pesar de arrojar un valor de 0.75, sólo el 16% es atribuible al modelo, debido a que la tasa de no información es del 59%.
También se puede ver del indicador Kappa que al superar el 40%, es un valor considerado bueno.
Modelo 2
Para el modelo 2, el análisis de bondad al igual que el modelo anterior muestra que es más específico que sensible, pero es más sensible que el modelo 1, con un valor de 0.69 versus un 0.61 del modelo 1.
Respecto al indicador de precisión, arroja un valor de 0.81, con un 22% es atribuible al modelo, debido a que la tasa de no información es del 59%.
También se puede ver del indicador Kappa supera el 60%, considerado muy bueno.
# Evaluando la bondad de clasificación
table(accidentes$Acc,Class_RL1,dnn=c("Observado","Predicho")) # Matriz de confusión## Predicho
## Observado No Si
## No 16 3
## Si 5 8caret::confusionMatrix(Class_RL1,accidentes$Acc,positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 16 5
## Si 3 8
##
## Accuracy : 0.75
## 95% CI : (0.566, 0.8854)
## No Information Rate : 0.5938
## P-Value [Acc > NIR] : 0.04978
##
## Kappa : 0.4689
##
## Mcnemar's Test P-Value : 0.72367
##
## Sensitivity : 0.6154
## Specificity : 0.8421
## Pos Pred Value : 0.7273
## Neg Pred Value : 0.7619
## Prevalence : 0.4062
## Detection Rate : 0.2500
## Detection Prevalence : 0.3438
## Balanced Accuracy : 0.7287
##
## 'Positive' Class : Si
## ICC_RL1=caret::confusionMatrix(Class_RL1,accidentes$Acc,positive = "Si")
ICC_RL1$byClass # Almacena todos los indicadores de desempeño## Sensitivity Specificity Pos Pred Value
## 0.6153846 0.8421053 0.7272727
## Neg Pred Value Precision Recall
## 0.7619048 0.7272727 0.6153846
## F1 Prevalence Detection Rate
## 0.6666667 0.4062500 0.2500000
## Detection Prevalence Balanced Accuracy
## 0.3437500 0.7287449table(accidentes$Acc,Class_RL10,dnn=c("Observado","Predicho")) # Matriz de confusión## Predicho
## Observado No Si
## No 17 2
## Si 4 9caret::confusionMatrix(Class_RL10,accidentes$Acc,positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 17 4
## Si 2 9
##
## Accuracy : 0.8125
## 95% CI : (0.6356, 0.9279)
## No Information Rate : 0.5938
## P-Value [Acc > NIR] : 0.007565
##
## Kappa : 0.6017
##
## Mcnemar's Test P-Value : 0.683091
##
## Sensitivity : 0.6923
## Specificity : 0.8947
## Pos Pred Value : 0.8182
## Neg Pred Value : 0.8095
## Prevalence : 0.4062
## Detection Rate : 0.2812
## Detection Prevalence : 0.3438
## Balanced Accuracy : 0.7935
##
## 'Positive' Class : Si
## ICC_RL10=caret::confusionMatrix(Class_RL10,accidentes$Acc,positive = "Si")
ICC_RL10$byClass## Sensitivity Specificity Pos Pred Value
## 0.6923077 0.8947368 0.8181818
## Neg Pred Value Precision Recall
## 0.8095238 0.8181818 0.6923077
## F1 Prevalence Detection Rate
## 0.7500000 0.4062500 0.2812500
## Detection Prevalence Balanced Accuracy
## 0.3437500 0.7935223Graficos de densidad
A continuación se muestran los gráficos de densidad de ambos modelos, donde con una línea negra se muestran el comportamiento de los no accidentados y la roja de los accidentados.Recordemos que tenemos el punto de corte en 0.5 y que mientras más tienda el punto de corte hacia la derecha, tiende a atinarle más a los no accidentados haciendolo más específico y si tiende a la izquierda, es porque procura un modelo más sensible.
windows(height=10,width=15)
par(mfrow=c(2,1))
with(accidentes,{
dD=density(Prob_RL1[Acc=="Si"])
dS=density(Prob_RL1[Acc=="No"])
plot(NULL, xlim=range(Prob_RL1), ylim=range(c(dD$y,dS$y)),
type="n",ylab="Densidad",xlab="Probabilidad Exp")
lines(dD,col= color[4]);lines(dS,col= color[1])
text( 0.1,1.6, "No accidentados"); text(0.9,1.6, "Accidentados")
PC.option=c(0.5,0.25,0.75)
abline(v=PC.option,col="red",lty=2)
text(PC.option,2,paste("PC = ",PC.option),cex=0.7,adj=0)
})
with(accidentes,{
dD10=density(Prob_RL10[Acc=="Si"])
dS10=density(Prob_RL10[Acc=="No"])
plot(NULL, xlim=range(Prob_RL10), ylim=range(c(dD10$y,dS10$y)),
type="n",ylab="Densidad",xlab="Probabilidad Exp~Sexo")
lines(dD10,col= color[4]);lines(dS10,col= color[1])
text( 0.1,1.6, "No accidentados"); text(0.9,1.6, "Accidentados")
PC.option=c(0.5,0.25,0.75)
abline(v=PC.option,col="red",lty=2)
text(PC.option,2,paste("PC = ",PC.option),cex=0.7,adj=0)
})
#### Explorando un mejor punto de corte
Teniendo en cuenta que se debe buscar el mejor punto para equilibrar sensibilidad con especificidad, se realiza la curva ROC para cada modelo que nos indique el punto que busca armonía entre sensibilidad y especificidad.
roc <- pROC::roc(accidentes$Acc,Prob_RL1, auc = TRUE, ci = TRUE) #AUC áreabajo la curva## Setting levels: control = No, case = Si## Setting direction: controls < casesprint(roc)##
## Call:
## roc.default(response = accidentes$Acc, predictor = Prob_RL1, auc = TRUE, ci = TRUE)
##
## Data: Prob_RL1 in 19 controls (accidentes$Acc No) < 13 cases (accidentes$Acc Si).
## Area under the curve: 0.9008
## 95% CI: 0.7989-1 (DeLong)windows(height=10,width=15)
pROC::plot.roc(roc, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
col = 2, grid = TRUE)
roc10 <- pROC::roc(accidentes$Acc,Prob_RL10, auc = TRUE, ci = TRUE) #AUC áreabajo la curva## Setting levels: control = No, case = Si
## Setting direction: controls < casesprint(roc10)##
## Call:
## roc.default(response = accidentes$Acc, predictor = Prob_RL10, auc = TRUE, ci = TRUE)
##
## Data: Prob_RL10 in 19 controls (accidentes$Acc No) < 13 cases (accidentes$Acc Si).
## Area under the curve: 0.9008
## 95% CI: 0.7953-1 (DeLong)windows(height=10,width=15)
pROC::plot.roc(roc10, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
col = 2, grid = TRUE)
Al ver las gráficas anteriores se observa que cada modelo escogió un
lado, el primero queriendo ser más sensible ubicandose a un punto de
corte de 0.211 para procurar un 100% de aciertos en accidentes y el
segundo modelo queriendo ser más específico ubicandose a un punto de
corte de 0.666 para procurar un 100% de aciertos en los no
accidentados.
Actualización de Clasificaciones
A continuación se dejan como puntos de corte los puntos recomendados por las gráficas ROC y se generan de nuevo las matrices de confusión.
pc=0.211
Predict_RL1 = as.factor(ifelse(Prob_RL1>pc,"Si","No"))
data.frame(accidentes$Acc,Predict_RL1,Class_RL1)## accidentes.Acc Predict_RL1 Class_RL1
## 1 No No No
## 2 No No No
## 3 No Si Si
## 4 Si Si Si
## 5 No No No
## 6 Si Si Si
## 7 No Si No
## 8 No No No
## 9 Si Si Si
## 10 No No No
## 11 No No No
## 12 Si Si Si
## 13 No No No
## 14 No No No
## 15 Si Si Si
## 16 No Si No
## 17 No Si Si
## 18 No Si Si
## 19 No No No
## 20 Si Si Si
## 21 Si Si No
## 22 No No No
## 23 No No No
## 24 Si Si Si
## 25 Si Si No
## 26 No Si No
## 27 No No No
## 28 Si Si No
## 29 Si Si Si
## 30 No No No
## 31 Si Si No
## 32 Si Si Notable(accidentes$Acc,Predict_RL1,dnn=c("Observado","Predicho"))## Predicho
## Observado No Si
## No 13 6
## Si 0 13pc10=0.666
Predict_RL10 = as.factor(ifelse(Prob_RL10>pc10,"Si","No"))
data.frame(accidentes$Acc,Predict_RL10,Class_RL10)## accidentes.Acc Predict_RL10 Class_RL10
## 1 No No No
## 2 No No No
## 3 No No No
## 4 Si Si Si
## 5 No No No
## 6 Si Si Si
## 7 No No No
## 8 No No No
## 9 Si Si Si
## 10 No No No
## 11 No No No
## 12 Si Si Si
## 13 No No Si
## 14 No No No
## 15 Si Si Si
## 16 No No No
## 17 No No Si
## 18 No No No
## 19 No No No
## 20 Si Si Si
## 21 Si No No
## 22 No No No
## 23 No No No
## 24 Si Si Si
## 25 Si No No
## 26 No No No
## 27 No No No
## 28 Si Si Si
## 29 Si Si Si
## 30 No No No
## 31 Si No No
## 32 Si No Notable(accidentes$Acc,Predict_RL10,dnn=c("Observado","Predicho"))## Predicho
## Observado No Si
## No 19 0
## Si 4 9Como se mencionó anteriormente en la primer matriz de confusión se observa un ajuste perfecto de sensibilidad donde de los 13 accidentes observados, todos fueron aciertos, sin dejar falsos negativos y un total de 6 no aciertos. Por otra parte, en la segunda matriz de confusión se observa un ajuste perfecto pero de especificidad donde de los 19 casos que no se accidentaron, todos fueron aciertos, sin dejar falsos positivos y un total de 4 no aciertos Como puntos a favor del segundo modelo, se tiene que posee menor cantidad de errores totales; adicional a la bondad de ajuste donde todos los indicadores son, aunque por poco, mejores que los del primer modelo.Sin embargo, queda a criterio del analista debido que ambos modelos son altamente eficientes.
Actualización de Indicadores de desempeño.
Se actualizan los indicadores de desempeño acorde a los nuevos puntos de corte.
caret::confusionMatrix(Predict_RL1,accidentes$Acc,positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 13 0
## Si 6 13
##
## Accuracy : 0.8125
## 95% CI : (0.6356, 0.9279)
## No Information Rate : 0.5938
## P-Value [Acc > NIR] : 0.007565
##
## Kappa : 0.6377
##
## Mcnemar's Test P-Value : 0.041227
##
## Sensitivity : 1.0000
## Specificity : 0.6842
## Pos Pred Value : 0.6842
## Neg Pred Value : 1.0000
## Prevalence : 0.4062
## Detection Rate : 0.4062
## Detection Prevalence : 0.5938
## Balanced Accuracy : 0.8421
##
## 'Positive' Class : Si
## ICC_RL1=caret::confusionMatrix(Predict_RL1,accidentes$Acc,positive = "Si")
ICC_RL1$byClass # Se generan las predicciones ## Sensitivity Specificity Pos Pred Value
## 1.0000000 0.6842105 0.6842105
## Neg Pred Value Precision Recall
## 1.0000000 0.6842105 1.0000000
## F1 Prevalence Detection Rate
## 0.8125000 0.4062500 0.4062500
## Detection Prevalence Balanced Accuracy
## 0.5937500 0.8421053Ind_class=list()
Ind_class$RL1=ICC_RL1$byClass
caret::confusionMatrix(Predict_RL10,accidentes$Acc,positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 19 4
## Si 0 9
##
## Accuracy : 0.875
## 95% CI : (0.7101, 0.9649)
## No Information Rate : 0.5938
## P-Value [Acc > NIR] : 0.0005536
##
## Kappa : 0.7277
##
## Mcnemar's Test P-Value : 0.1336144
##
## Sensitivity : 0.6923
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.8261
## Prevalence : 0.4062
## Detection Rate : 0.2812
## Detection Prevalence : 0.2812
## Balanced Accuracy : 0.8462
##
## 'Positive' Class : Si
## ICC_RL10=caret::confusionMatrix(Predict_RL10,accidentes$Acc,positive = "Si")
ICC_RL10$byClass # Se generan las predicciones ## Sensitivity Specificity Pos Pred Value
## 0.6923077 1.0000000 1.0000000
## Neg Pred Value Precision Recall
## 0.8260870 1.0000000 0.6923077
## F1 Prevalence Detection Rate
## 0.8181818 0.4062500 0.2812500
## Detection Prevalence Balanced Accuracy
## 0.2812500 0.8461538Ind_class10=list()
Ind_class10$RL10=ICC_RL10$byClassBondad de ajuste
Con toda la información recogida hasta el momento se validará la bondad de ajuste de cada modelo para escoger el más adecuado.
Bondad_Ajuste=data.frame(AIC=Modelo_RL1$aic,deviance=Modelo_RL1$deviance,R2=R2(Modelo_RL1),AUC=roc$auc)
row.names(Bondad_Ajuste)="Modelo_RL1"
Bondad_Ajuste[2,]=data.frame(AIC=Modelo_RL10$aic,deviance=Modelo_RL10$deviance,R22=R22(Modelo_RL10),AUC=roc10$auc)
rownames(Bondad_Ajuste)[2]="Modelo_RL2"
Bondad_Ajuste## AIC deviance R2 AUC
## Modelo_RL1 28.01390 24.01390 0.4445050 0.9008097
## Modelo_RL2 27.89302 21.89302 0.4935657 0.9008097Teniendo en cuenta que en los indicadores mostrados en la tabla anterior no existe una diferencia significativa, que consideramos más importante pronosticar el que se va a accidentar y que en gran medida depende del criterio de acuerdo a lo requerido, escogerémos el primer modelo que tiene mayor acierto en sensibilidad.
Prediccion nuevo valor con modelo seleccionado
A continuación se busca predecir 5 valores con datos de prueba, para ver la reacción al nuevo modelo.
dato.nuevo=accidentes[25:30,c(2,3)]
Predict_new=predict(Modelo_RL1,newdata = dato.nuevo, type ="response")
as.factor(ifelse(Predict_new>pc,"Si","No")) ## 25 26 27 28 29 30
## Si Si No Si Si No
## Levels: No SiSe observa acorde a los datos reales, que falló solo en 1 dato que corresponde a un falso positivo.
Graficos individuales de las variables Cuantitativas
A continuación se observan las variables cuantitaticas contra la variable a predecir, accidentes (Acc). No es posible en graficar en este caso el sexo debido que corresponde a una categórica.
####Ajuste por regresión logística múltiple####
windows(height=10,width=15);par(mfrow=c(2,2)) # Partición de la ventana grafica 2x4
color = ifelse(accidentes$Acc=="No","Black","Red")
attach(accidentes)
tabla<-prop.table(table(Acc))
coord<-barplot(tabla, col=color[c(1,4)],ylim=c(0,1), main="Accidentes")
text(coord,tabla,labels=round(tabla,2), pos=3)
lapply(names(accidentes[,c(-1,-5)]),function(y){
boxplot(accidentes[,y]~accidentes[,"Acc"],ylab= y, xlab="Accidentes",boxwex = 0.5,col=NULL)
stripchart(accidentes[,y] ~ accidentes[,"Acc"], vertical = T,
method = "jitter", pch = 19,
col = color[c(1,4)], add = T) #stripchart es para que muestre los puntos sobrepuestos
})
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULLdetach(accidentes)Se puede observar en el gráfico de barras que los datos nos indican que aproximadamente 4 de cada 10 personas se accidentan, con respecto a los diagrama de caja se puede concluir que la edad es la que menos influencia tiene en los accidentes, mientras que la experiencia y potencia de motor si separan más los datos, lo cual se puede concluir viendo la diferencia en el eje Y.
Visualización bivariada
A continuación se muestran los gráficos de dispersión comparando las variables, diferenciando los accidentados (rojos) versus los no accidentados (negros).
windows(height=10,width=15)
par(mfrow=c(2,2),oma=c(1,1,1,1), mar=c(4,4,1,1))
attach(accidentes)
for (i in 2:3){
x=names(accidentes)[i]
rango.x=range(accidentes[,x])
for(j in (i+1):4){
y=names(accidentes)[j]
rango.y=range(accidentes[,y])
plot(accidentes[Acc=="No",x],accidentes[Acc=="No",y],xlab=x,ylab=y,xlim=rango.x,ylim=rango.y,col=color[1],pch=20)
points(accidentes[Acc=="Si",x],accidentes[Acc=="Si",y],col=color[4],pch=20)
}
}
detach(accidentes)
Estrategia de modelo Saturado
Se realiza un modelo saturado, pasando a esbelto eliminando las variables no significantes para obtener un mejor método de predicción.
Modelo_RLMsat = glm(formula= Acc~., data=accidentes,family = "binomial")
summary(Modelo_RLMsat)##
## Call:
## glm(formula = Acc ~ ., family = "binomial", data = accidentes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26813 -0.36434 -0.04645 0.02284 2.31957
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.82369 10.37123 -1.719 0.0857 .
## Exp -0.48020 0.38442 -1.249 0.2116
## Edad -0.01699 0.08921 -0.190 0.8489
## Pot 0.22186 0.11622 1.909 0.0563 .
## SexoHombre 2.78011 2.78173 0.999 0.3176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 12.327 on 27 degrees of freedom
## AIC: 22.327
##
## Number of Fisher Scoring iterations: 8## usando la función step de R
Modelo_RLM1=step(Modelo_RLMsat,direction = "backward")## Start: AIC=22.33
## Acc ~ Exp + Edad + Pot + Sexo
##
## Df Deviance AIC
## - Edad 1 12.363 20.363
## - Sexo 1 13.274 21.274
## <none> 12.327 22.327
## - Exp 1 15.834 23.834
## - Pot 1 21.285 29.285
##
## Step: AIC=20.36
## Acc ~ Exp + Pot + Sexo
##
## Df Deviance AIC
## - Sexo 1 13.392 19.392
## <none> 12.363 20.363
## - Exp 1 17.352 23.352
## - Pot 1 21.893 27.893
##
## Step: AIC=19.39
## Acc ~ Exp + Pot
##
## Df Deviance AIC
## <none> 13.392 19.392
## - Exp 1 21.297 25.297
## - Pot 1 24.014 28.014Esta opción step nos muestra que en el primer paso la variable elimienada es edad, algo que intuíamos después del análisis visual realizado a las cajas, donde se observaba como la menos significante. De manera posterior eliminó la variable sexo, motivo por el cual aumentar la potencia del motor si optiene una mejora significativa, pero agregar la edad no. Sin embargo, para la solución del punto vii se le adicionara la Edad a las variables seleccionadas experiencia y potencia (Exp + Pot + Edad).
summary(Modelo_RLM1)##
## Call:
## glm(formula = Acc ~ Exp + Pot, family = "binomial", data = accidentes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.15224 -0.27743 -0.02508 0.03644 2.31198
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -15.4824 9.9410 -1.557 0.1194
## Exp -0.7304 0.3891 -1.877 0.0605 .
## Pot 0.2162 0.1142 1.893 0.0583 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 13.392 on 29 degrees of freedom
## AIC: 19.392
##
## Number of Fisher Scoring iterations: 8Modelo_RLM1=update(Modelo_RLMsat, . ~ . - Sexo)
summary(Modelo_RLM1)##
## Call:
## glm(formula = Acc ~ Exp + Edad + Pot, family = "binomial", data = accidentes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.25677 -0.27816 -0.02113 0.03302 2.28821
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.14838 9.93984 -1.625 0.1042
## Exp -0.71989 0.39828 -1.807 0.0707 .
## Edad -0.02977 0.08590 -0.347 0.7289
## Pot 0.23177 0.12235 1.894 0.0582 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 13.274 on 28 degrees of freedom
## AIC: 21.274
##
## Number of Fisher Scoring iterations: 8Al ver este modelo anterior se observa que empeora la predicción del calculado con respecto al calculado con step debido que tiene un AIC mayor (21.27 contra 19.39) y Edad es la variable que daña la significancia es la edad con 0.72. Sin embargo, es mejor que los primeros 2 modelos examinados.
anova(Modelo_RLMsat,Modelo_RLM1,test="Chisq")## Analysis of Deviance Table
##
## Model 1: Acc ~ Exp + Edad + Pot + Sexo
## Model 2: Acc ~ Exp + Edad + Pot
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 27 12.327
## 2 28 13.274 -1 -0.94715 0.3304En el resultado anterior podemos ver que el modelo saturado es un poco mejor que el otro.
Bondad_Ajuste[2,]=c(Modelo_RLM1$aic,Modelo_RLM1$deviance,R2(Modelo_RLM1),NA);
rownames(Bondad_Ajuste)[2]="Modelo_Mult1"
Bondad_Ajuste## AIC deviance R2 AUC
## Modelo_RL1 28.01390 24.01390 0.4445050 0.9008097
## Modelo_Mult1 21.27402 13.27402 0.6929424 NACon este otro resultado es notorio el mejoramiento del modelo al que se le adicionó edad y potencia mucho mejor que el primer modelo escogido.
Interpretación de los parametros
Para completar la tabla anterior calculamos la curva ROC para agregar el AUC y realizar la comparación completa.
Tabla_coef=cbind(Coef=coef(Modelo_RLM1),IC=confint(Modelo_RLM1,level=0.95))## Waiting for profiling to be done...Tabla_odds = exp(Tabla_coef)
colnames(Tabla_odds)[1]="e-beta"
Tabla_coef;Tabla_odds## Coef 2.5 % 97.5 %
## (Intercept) -16.14838084 -40.66482616 -0.02501476
## Exp -0.71988538 -1.87506963 -0.13035832
## Edad -0.02976718 -0.20473233 0.15974753
## Pot 0.23176777 0.05201918 0.53377558## e-beta 2.5 % 97.5 %
## (Intercept) 9.701688e-08 2.185196e-18 0.9752955
## Exp 4.868081e-01 1.533443e-01 0.8777809
## Edad 9.706715e-01 8.148654e-01 1.1732146
## Pot 1.260827e+00 1.053396e+00 1.7053589# Predicciones del nuevo modelo
Prob_RLM1 = predict(Modelo_RLM1,type = "response")
# Buscando el punto de corte para el modelo múltiple
roc1 <- pROC::roc(accidentes$Acc,Prob_RLM1, auc = TRUE, ci = TRUE)## Setting levels: control = No, case = Si## Setting direction: controls < casesprint(roc1)##
## Call:
## roc.default(response = accidentes$Acc, predictor = Prob_RLM1, auc = TRUE, ci = TRUE)
##
## Data: Prob_RLM1 in 19 controls (accidentes$Acc No) < 13 cases (accidentes$Acc Si).
## Area under the curve: 0.9676
## 95% CI: 0.9131-1 (DeLong)windows(height=10,width=15)
pROC::plot.roc(roc1, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
col = 2, grid = TRUE)
Se obtiene el AUC y se adiciona a la tabla.
Bondad_Ajuste[2,4]=0.9676
Bondad_Ajuste## AIC deviance R2 AUC
## Modelo_RL1 28.01390 24.01390 0.4445050 0.9008097
## Modelo_Mult1 21.27402 13.27402 0.6929424 0.9676000Se modifica el punto de corte de acuerdo al nuevo ROC y se genera una nueva matriz de confusión.
# Calidad de la Clasificación en el modelo múltiple
pc=0.657
Predict_RLM1 = as.factor(ifelse(Prob_RLM1>pc,"Si","No"))
# Indicadores de correcta clasificación.
caret::confusionMatrix(Predict_RLM1,accidentes$Acc,positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction No Si
## No 19 2
## Si 0 11
##
## Accuracy : 0.9375
## 95% CI : (0.7919, 0.9923)
## No Information Rate : 0.5938
## P-Value [Acc > NIR] : 1.452e-05
##
## Kappa : 0.8672
##
## Mcnemar's Test P-Value : 0.4795
##
## Sensitivity : 0.8462
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9048
## Prevalence : 0.4062
## Detection Rate : 0.3438
## Detection Prevalence : 0.3438
## Balanced Accuracy : 0.9231
##
## 'Positive' Class : Si
## ICC_RLM1=caret::confusionMatrix(Predict_RLM1,accidentes$Acc,positive = "Si")
Ind_class$RLM1=ICC_RLM1$byClassEn los indicadores anteriores un excelente grado de predicción mayor al 93%, con un indice Kappa de 0.8672 una sensibilidad del 0.8462 y una especificidad de 1. Adicional en la matriz de confusión se ven solo 2 eventos no acertados. El modelo de mayor predicción total.
Evaluación de los residuos
Por último se muestran los residuos o diferencias de evento con respecto a lo predicho
windows()
residualPlots(Modelo_RLM1,type="deviance")
## Test stat Pr(>|Test stat|)
## Exp 0.0522 0.81923
## Edad 2.9773 0.08444 .
## Pot 0.0248 0.87477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Reporte de hallazgos
Se puede concluir que el último modelo mostrado es el de mayor acierto y que de las variables que afectan la siniestralidad solo la edad no tiene tanta importancia, porque la experiencia aporta toda la información que la edad podría dar a este conjunto de datos. Desde el inicio se observa cómo la potencia del motor influye significativamente, siendo los vehículos con mayor potencia son los más accidentados, por otra parte no es despreciable la persona que conduce, pues con el diagrama de barras se observa que el hombre se accidenta mucho más que las mujeres porcentualmente y la experiencia es inversamente proporcional a la accidentalidad, es decir a mayor experiencia, menor accidentalidad.En ese caso se podría decir que las personas que menor probabilidad de accidentes son las mujeres con un motor de baja potencia y con muchos años de experiencia y los potenciales accidentados son los hombres de poca experiencia con un motor de alta potencia.