Exámen Parcial II

Estadística Avanzada IV

PRESENTADO POR:

ROSMERY CHARRIS

KAREN CORONELL

ASTRID DÍAZ

#load functions
source('https://www.dropbox.com/s/gnwm7vuvhojg9q4/logistic-functions.R?dl=1')

Esta es una colección de funciones creadas 
por el profesor Jorge Vélez <https://jivelez.github.io/> 
para el ajuste y validación de modelos de Regresión Logística 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a jvelezv@uninorte.edu.co 
 

source('https://www.dropbox.com/s/gnwm7vuvhojg9q4/logistic-functions.R?dl=1')

Esta es una colección de funciones creadas 
por el profesor Jorge Vélez <https://jivelez.github.io/> 
para el ajuste y validación de modelos de Regresión Logística 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a jvelezv@uninorte.edu.co 
 

## load functions
source('https://www.dropbox.com/s/gnwm7vuvhojg9q4/logistic-functions.R?dl=1')

Esta es una colección de funciones creadas 
por el profesor Jorge Vélez <https://jivelez.github.io/> 
para el ajuste y validación de modelos de Regresión Logística 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a jvelezv@uninorte.edu.co 
 

Contexto

Un grupo de investigadores recolecta información de 386 individuos; de estos, algunos son diagnosticados con la enfermedad (columna y = yes) y otros no (columna y = no). Para todos los individuos se registró la presencia o ausencia de 20 síntomas diferentes; cada síntoma \(x_j = 1\) si el síntoma está presente y \(x_j = 0\) si no lo está \((j=1,2,\ldots,20)\).

El Problema

Se requiere determinar, con una alta precisión, el diagnóstico de la próxima persona a partir de los síntomas reportados. Por la naturaleza de la variable respuesta, pueden utilizarse los modelos de Regresión Logística, CART, Random Forests y SVM.

El otro requerimiento es que, en la medida de lo posible, el número de variables predictoras utilizadas para determinar si el diagnóstico de la enfermedad es yes o no debe ser el menor.

Contexto Analítico

El conjunto de datos pueden leerse en R utilizando los siguientes comandos:

## Presentación de los Datos
url <- 'https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1'
d <- read.table(url, header = TRUE)
k<-head(d,6)
kable( k, caption = "Presentación de las primeras 10 observaciones, de la base de datos divorce", row.names = TRUE
       , digits = 1
       , format.args = list( decimal.mark = ","))%>%
   kable_paper() %>%
  scroll_box(width = "650px", height = "400px")%>%
  footnote(general = "https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1")

Table 1: Presentación de las primeras 10 observaciones, de la base de datos divorce

	id	y	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	x13	x14	x15	x16	x17	x18	x19	x20
1	1	no	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0
2	2	no	0	1	0	1	0	0	0	1	0	0	0	1	1	0	0	0	1	0	0	0
3	3	yes	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	0	0	1	1	0
4	4	no	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0
5	5	no	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
6	6	no	0	0	0	0	1	0	1	0	1	0	0	0	0	0	0	1	0	0	1	0

Note: https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1

## lectura de datos
url <- 'https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1'
d <- read.table(url, header = TRUE)
require(stringr)
reemplazos<-c("no"="0","yes"="1")
d$y<-factor(str_replace_all( d$y, reemplazos ))
levels(d$y)

[1] "0" "1"

Las variables relevantes en este caso son:

1. id: identificador del individuo
2. y: diagnóstico (yes, no)
3. xj: síntoma \(j\) (1: presente; 0: ausente)

Análisis previo escogencia de las variables del modelo

out <- sapply(1:20, function(j) chisq.test(with(d, table(d[,2+j], y)))$statistic)
names(out) <- paste0("X", 1:20)
sort(out, decreasing = TRUE)

       X2        X3        X6        X8        X1       X14       X10 
100.55892  96.00006  92.33030  90.34970  85.51428  80.79859  78.31227 
      X15       X13        X4       X12        X7       X19        X5 
 78.31112  65.44953  65.15814  62.36677  59.28506  56.90504  54.83594 
      X11        X9       X18       X17       X16       X20 
 50.66607  42.70207  41.93533  35.53293  33.74325  18.28584 

out1 <- sapply(1:20, function(j) chisq.test(with(d, table(d[,2+j], y)))$p.value)
names(out1) <- paste0("X", 1:20)
sort(out1, decreasing = TRUE)

         X20          X16          X17          X18           X9 
1.901152e-05 6.288695e-09 2.507754e-09 9.434266e-11 6.374474e-11 
         X11           X5          X19           X7          X12 
1.094961e-12 1.310230e-13 4.573713e-14 1.364082e-14 2.850938e-15 
          X4          X13          X15          X10          X14 
6.912229e-16 5.962143e-16 8.802516e-19 8.797393e-19 2.499361e-19 
          X1           X8           X6           X3           X2 
2.300341e-20 1.995757e-21 7.335269e-22 1.148803e-22 1.149268e-23 

library(vcd)
tabla1=table(d$y,x=d$x2)
assocstats(tabla1)

                    X^2 df P(> X^2)
Likelihood Ratio 108.04  1        0
Pearson          102.64  1        0

Phi-Coefficient   : 0.516 
Contingency Coeff.: 0.458 
Cramer's V        : 0.516 

tabla2=table(d$y,x=d$x3)
assocstats(tabla2)

                     X^2 df P(> X^2)
Likelihood Ratio 103.305  1        0
Pearson           98.032  1        0

Phi-Coefficient   : 0.504 
Contingency Coeff.: 0.45 
Cramer's V        : 0.504 

tabla3=table(d$y,x=d$x6)
assocstats(tabla3)

                     X^2 df P(> X^2)
Likelihood Ratio 102.900  1        0
Pearson           94.357  1        0

Phi-Coefficient   : 0.494 
Contingency Coeff.: 0.443 
Cramer's V        : 0.494 

Evaluación

Los ys deberán comprar los modelos de Regresión Logística, CART, Random Forests y SVM, y escoger aquel que proporcione el mejor Accuracy para el testing data set con el menor número de variables predictoras.

Inicialmente empleamos el software R studio, se determinó la correlación entre las variables independientes xj; síntoma \(j\) y la variable descenlace y; diagnóstico (yes, no). Con El fin de crear un sistema inteligente de clasificación, tomamos las variables \(x2\), \(x3\), \(x6\) pues son las que mayor peso tienen con respecto a las demás variables que se encuentran consignadas en nuestra base de datos.

La medida que se tiene en cuenta para comparar los diferentes Modelos es ACCURACY , pues esta nos indica que tan acertado es el modelo en la predicción del diagnostico en los pacientes; se refiere con exactitud, a cuán cerca del valor real se encuentra el valor medido. Otro criterio elegido, es explicar la mayor cantidad de información posible, pero con el menor número de variables.

MODELOS DE CLASIFICACIÓN

Regresión Logistica

#Convertimos las variables en factor, para un mejor manejo en las herramientas estadisticas.
d %<>% mutate_if(is.integer,as.factor)
head(glimpse(d),3)

Rows: 386
Columns: 22
$ id  <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16...
$ y   <fct> 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0,...
$ x1  <fct> 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0,...
$ x2  <fct> 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0,...
$ x3  <fct> 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0,...
$ x4  <fct> 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0,...
$ x5  <fct> 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1,...
$ x6  <fct> 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0,...
$ x7  <fct> 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,...
$ x8  <fct> 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0,...
$ x9  <fct> 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0,...
$ x10 <fct> 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0,...
$ x11 <fct> 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0,...
$ x12 <fct> 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0,...
$ x13 <fct> 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0,...
$ x14 <fct> 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0,...
$ x15 <fct> 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,...
$ x16 <fct> 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0,...
$ x17 <fct> 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1,...
$ x18 <fct> 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0,...
$ x19 <fct> 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0,...
$ x20 <fct> 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0,...

  id y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18
1  1 0  0  0  0  0  0  0  0  0  0   0   1   0   0   0   0   0   1   0
2  2 0  0  1  0  1  0  0  0  1  0   0   0   1   1   0   0   0   1   0
3  3 1  1  1  1  1  1  1  1  1  1   1   1   1   1   1   1   0   0   1
  x19 x20
1   0   0
2   0   0
3   1   0

#Subconjunto con las variables de mayor peso
require(dplyr)
d_sub <- d %>%dplyr::select(y, x2, x3, x6)
## Subconjuntos Training y Testing
intrain <- createDataPartition(y = d_sub$y, p= 0.7, list = FALSE)
training <- d_sub[intrain,]

testing <- d_sub[-intrain,]

Modelo 1

#Modelo sin ninguna variable predictora
mod.ini <- glm(y ~ 1, data = training, family =binomial)
###Resumen modelo inicial
summary(mod.ini)

Call:
glm(formula = y ~ 1, family = binomial, data = training)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.316  -1.316   1.045   1.045   1.045  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   0.3200     0.1231   2.601   0.0093 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 368.83  on 270  degrees of freedom
Residual deviance: 368.83  on 270  degrees of freedom
AIC: 370.83

Number of Fisher Scoring iterations: 4

Modelo 2

#Modelo con todas las variables independientes escogidas como predictores 
logmodel<-glm(y~., data = training, family = binomial)
###Resumen modelo con todos los predictores
summary(logmodel)

Call:
glm(formula = y ~ ., family = binomial, data = training)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2959  -0.7775   0.3857   0.7677   1.6394  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.0416     0.2060  -5.056 4.29e-07 ***
x21           1.2062     0.3806   3.169  0.00153 ** 
x31           0.5449     0.3934   1.385  0.16598    
x61           1.8518     0.3991   4.640 3.48e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 368.83  on 270  degrees of freedom
Residual deviance: 264.43  on 267  degrees of freedom
AIC: 272.43

Number of Fisher Scoring iterations: 5

Modelo 3

###Modelo ajustado, usando el metódo stepwise
mod.final <- stepAIC(mod.ini, 
                     scope = list(upper = logmodel), 
                     direction = "both")

Start:  AIC=370.83
y ~ 1

       Df Deviance    AIC
+ x6    1   287.35 291.35
+ x2    1   295.74 299.74
+ x3    1   312.78 316.78
<none>      368.83 370.83

Step:  AIC=291.35
y ~ x6

       Df Deviance    AIC
+ x2    1   266.30 272.30
+ x3    1   274.43 280.43
<none>      287.35 291.35
- x6    1   368.83 370.83

Step:  AIC=272.3
y ~ x6 + x2

       Df Deviance    AIC
<none>      266.30 272.30
+ x3    1   264.43 272.43
- x2    1   287.35 291.35
- x6    1   295.74 299.74

El primer modelo se inició con todas las variables seleccionadas como predictores. sin embargo se utilizó un modelo alternativo con el método stepwise mixto (both) que selecciona los mejores predictores del modelo con el criterio Akaike (AIC), con el fin de determinar la calidad del modelo.

Algo realmente interesantes es que este criterio de selección, escogió basandose en ACCURACY, el modelo con todas las variables como predictores. Todas contrubuyen apropiadamente en el modelo

Supuestos del modelo para el training

#logverosimilitud
dev <- logmodel$deviance
nullDev <- logmodel$null.deviance
modelChi <- nullDev - dev
modelChi

[1] 104.4038

#Devianza
chigl <- logmodel$df.null - logmodel$df.residual
chisq.prob <- 1 - pchisq(modelChi, chigl)
chisq.prob

[1] 0

#R^2
R2.hl <- modelChi/logmodel$null.deviance
R2.hl

[1] 0.2830645

Del modelo completo, se concluye, que a un nivel de significancia del 5%, las variables x31 presente, x2 presente y x6 presente son significativas. Además se observó que la devianza del modelo (Residual deviance: 254.78 ) es mucho menor que la devianza solo con la constante (Null deviance: 368.83), lo cual significa que estas variables contribuyen a predecir mejor el diagnostico del individuo . Como medida de ajuste del modelo se obtuvo un AIC de 262.78

A través de la información que proporcionó el estadístico de Wald (z-statistic), vimos que las variables predictoras tiene una muy buena contribución en el modelo para predecir el diagnostico de un individuo.

Para saber la eficacia del modelo prediciendo la probabilidad de ser diagnosticado correctamente, se utilizó el estadístico chi-cuadrado \(LR=-2(l(\Theta_0)-l(\hat \Theta))\) dado que \(LR\sim \mathcal{X}_p^2\), donde \(l(\Theta_0)\) es la función log-verosimilitud del modelo final y \(l(\hat \Theta))\) es la función log-verosimilitud del modelo resultante sin variables independientes, sólo con la constante. Con \(p\) grados de libertad equivalente a la diferencia entre el número de parámetros del modelo final y el modelo nulo.

Se obtuvo un estadistico chi-cuadrado de 114.05 y la probabilidad asociada al estadístico fue de \(p-valor=P(\mathcal{X}_5^2>LR)=0.001\) es menor de 0.05, se concluye que el efecto general del modelo es estadisticamente significativo al 5%.

Determinamos medidas de clasificación del modelo sobre el traininig y el testing

source('https://www.dropbox.com/s/gnwm7vuvhojg9q4/logistic-functions.R?dl=1')

Esta es una colección de funciones creadas 
por el profesor Jorge Vélez <https://jivelez.github.io/> 
para el ajuste y validación de modelos de Regresión Logística 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a jvelezv@uninorte.edu.co 
 

Curva ROC para el training

pred <- prediction(fitted.values(logmodel),training$y)
# con performance se selecciona tpr (true positive rate) 
# y fpr (false positive rate)
perf2 <- performance(pred, "tpr", "fpr")
AUC <- performance(pred, "auc")
plot(perf2, colorize = TRUE,main = "Curva ROC para el training") # mostramos colores según el punto de corte
# Añadimos la recta y=x que sería la correspondiente al peor clasificador
abline(a = 0, b = 1)
# añadimos el valor del área bajo la curva
text(0.4, 0.6, paste(AUC@y.name, "\n", round(unlist(AUC@y.values), 3)), cex = 0.7)

cutoff <- seq(.1, .9, length = 1000)
res <- data.frame(cutoff = cutoff, bycutoff2(logmodel,real = training$y, cutoff))
with(res, matplot(cutoff, cbind(sensitivity, specificity, class_rate),type = 'l', las = 1, ylim = c(0, 1), lty = 1, ylab = 'Value (%)'))
legend('bottomleft', c('Se', 'Sp', 'CR'), lty = 1, col = 1:3, bty = 'n')

Puntos de corte

res[which.min(abs(res$sensitivity - res$specificity)), ]

       cutoff   a  b  c  d sensitivity specificity       ppv
349 0.3786787 128 26 29 88   0.8152866   0.7719298 0.8311688
          npv       fdr       fpr class_rate     lift     delta
349 0.7521368 0.1688312 0.2280702   0.797048 1.434693 0.2934877

res[which.max(res$class_rate), ]

       cutoff   a  b  c  d sensitivity specificity       ppv
349 0.3786787 128 26 29 88   0.8152866   0.7719298 0.8311688
          npv       fdr       fpr class_rate     lift     delta
349 0.7521368 0.1688312 0.2280702   0.797048 1.434693 0.2934877

res[which.max(res$lift), ]

       cutoff  a  b  c   d sensitivity specificity       ppv
741 0.6925926 89 10 68 104    0.566879   0.9122807 0.8989899
          npv       fdr       fpr class_rate    lift     delta
741 0.6046512 0.1010101 0.0877193  0.7121771 1.55176 0.4419146

Indices de validez

require(ROCR)
require(pROC)
rocA1CNP <- with(training,roc(y,fitted.values(logmodel)))
coords(rocA1CNP,"b",ret=c("threshold","specificity","sensitivity","npv","ppv"),best.method="closest.topleft")

          threshold specificity sensitivity       npv       ppv
threshold 0.4596902   0.7719298   0.8152866 0.7521368 0.8311688

Matriz de confusión

prediccion <- ifelse(fitted.values(logmodel) >= 0.5028028, 1, 0)

clasificacion<-table(prediccion);clasificacion

prediccion
  0   1 
117 154 

### Tabla de clasificación 
tabla.clasif <- table(training$y, prediccion)
tcc <- 100 * sum(diag(tabla.clasif))/sum(tabla.clasif)
tcc

[1] 79.7048

### Matriz de confusión reorganizada
matriz<-as.matrix(table(training$y, prediccion))
a<-matriz[2,2]
b<-matriz[1,2]
c<-matriz[2,1]
d<-matriz[1,1]

Mconf<-matrix(c(a,b,c,d),ncol=2,nrow=2,byrow=TRUE);Mconf

     [,1] [,2]
[1,]  128   26
[2,]   29   88

Para el testing

Predichos

predtion<-sapply(1:115, function(i) predict(logmodel, newdata =testing[i,], type = 'response'))

Curva ROC para el testing

pred2 <- prediction(predtion,testing$y)
# con performance se selecciona tpr (true positive rate) 
# y fpr (false positive rate)
perf3 <- performance(pred2, "tpr", "fpr")
AUC <- performance(pred2, "auc")
plot(perf3, colorize = TRUE,main = "Curva ROC para el testing") # mostramos colores según el punto de corte
# Añadimos la recta y=x que sería la correspondiente al peor clasificador
abline(a = 0, b = 1)
# añadimos el valor del área bajo la curva
text(0.4, 0.6, paste(AUC@y.name, "\n", round(unlist(AUC@y.values), 3)), cex = 0.7)

Indices de validez

rocA1CNP <- with(testing,roc(y,predtion))
coords(rocA1CNP,"b",ret=c("threshold","specificity","sensitivity","npv","ppv"),best.method="closest.topleft")

          threshold specificity sensitivity  npv       ppv
threshold 0.6056775      0.8125   0.8059701 0.75 0.8571429

Matriz de confusión

prediccion <- ifelse(predtion >= 0.646576  , 1, 0)
class<-table(prediccion); class

prediccion
 0  1 
52 63 

### Tabla de clasificación 
tabla.clasif <- table(testing$y,prediccion)
tcc <- 100 * sum(diag(tabla.clasif))/sum(tabla.clasif)
tcc

[1] 80.86957

### Matriz de confusión reorganizada
matriz2<-as.matrix(table(testing$y,prediccion ))
a<-matriz2[2,2]
b<-matriz2[1,2]
c<-matriz2[2,1]
d<-matriz2[1,1]

Mconf2<-matrix(c(a,b,c,d),ncol=2,nrow=2,byrow=TRUE);Mconf2

     [,1] [,2]
[1,]   54    9
[2,]   13   39

Caret

## lectura de datos
url <- 'https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1'
d <- read.table(url, header = TRUE)
d <- d[,-1] # eliminar la primera colmuna

if(!require(caret)) install.packages('caret')

d$y <- factor(d$y, levels = c('yes', 'no'))
d[,-1] <- apply(d[,-1], 2, as.factor)
d$y <- factor(d$y, levels = c('yes', 'no'))

Crear la partición

intrain <- createDataPartition(y = d$y, p= 0.7, list = FALSE)
training <- d[intrain,]
testing <- d[-intrain,]

caret model

set.seed(123) 
cartmodel <- train(y ~x2+x3+x6, data = training, 
  method = "rpart",trControl = trainControl("cv", number = 10), 
  tuneLength = 10)

grafica Cp

plot(cartmodel)

## mejor modelo
cartmodel$bestTune

          cp
3 0.08966862

# grafico del modelo final 
par(xpd = NA)  
plot(cartmodel$finalModel)
text(cartmodel$finalModel, digits = 3)

## optimize cp parameter using the training data set
ps <- prop.table(table(training$y))
cps <- seq(0.01, 0.1, by = 0.01)
out <- sapply(cps, function(cpi){
  fit <- rpart(y ~ x2+x3+x6, data = training, y = TRUE, method = 'class', parms = list(prior = ps, split = "information"), control = rpart.control(minsplit = 20, minbucket = round(20/3), cp = cpi, maxcompete = 4, maxsurrogate = 5, usesurrogate = 2, xval = 20))
  classhat <- apply(predict(fit), 1, which.max)
  m <- table(factor(classhat, levels = 1:2), factor(fit$y, levels = 1:2))
  c(measures(m), auc = rft(m, graph = FALSE)$auc)
})
cpopt <- cps[which.min(abs(out[1,] - out[2,]))]
cpopt

[1] 0.01

## fit the CART model using cpopt
set.seed(123)
fit <- rpart(y ~ x2+x3+x6, data = training, y = TRUE, method = 'class', parms = list(prior = ps, split = "information"), control = rpart.control(minsplit = 20, minbucket = round(20/3), cp = cpopt, maxcompete = 4, maxsurrogate = 5, usesurrogate = 2, xval = 20))
classhat <- apply(predict(fit), 1, which.max)
m <- table(predclass = factor(classhat, levels = 1:2), real = factor(fit$y, levels = 1:2))
(mtraining <- m)

         real
predclass   1   2
        1 129  27
        2  28  87

measures(mtraining)

sensitivity specificity         ppv         npv         fdr 
  0.8216561   0.7631579   0.8269231   0.7565217   0.1730769 
        fpr  class_rate        lift       delta 
  0.2368421   0.7970480   1.4273640   0.2964806 

graficas para el modelo Caret , Curva Roc y matriz de confusión

## plotting the tree
par(mfrow = c(1,1), mar = c(2, 2, 2, 2))
fancyRpartPlot(fit, sub = "")

## nice colors
cols <- c("#0080ff", "#ff00ff", "darkgreen", "#ff0000", 
          "orange", "#00ff00", "brown")


## variable importance
par(mar = c(4, 6, 2, 1), mfrow = c(1, 2))
vimp <- fit$variable.importance
cart <- vimp/max(vimp)
names(cart) <- names(fit$variable.importance)
barplot(rev(cart), horiz = TRUE, las = 1, width = .2, space = .1, col = cols[1], border = cols[1], xlab = "", main = "", ylab = "", cex.axis = 1.1, cex.names = 1.1)
mtext("Score", 1, line = 2.3)
rft(mtraining, las = 1)$auc

[1] 0.792407

## testing data set
predclass <- predict(fit, newdata = testing, type = 'class')
m <- table(predclass = predclass, real = testing$y)
(mtesting <- m)

         real
predclass yes no
      yes  55 10
      no   12 38

measures(mtesting)

sensitivity specificity         ppv         npv         fdr 
  0.8208955   0.7916667   0.8461538   0.7600000   0.1538462 
        fpr  class_rate        lift       delta 
  0.2083333   0.8086957   1.4523536   0.2747384 

rft(mtesting, las = 1, add = TRUE, line.col = cols[2])$auc

[1] 0.8062811

## confusion matrices
mtraining 

         real
predclass   1   2
        1 129  27
        2  28  87

mtesting

         real
predclass yes no
      yes  55 10
      no   12 38

SVM

## lectura de datos
d <- read.table("https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1", 
                header = TRUE)

## verificar paquete caret
if(!require(caret)) install.packages('caret')

## random seed
set.seed(123)

## training & testing data sets
intrain <- createDataPartition(y = d$y, p= 0.7, list = FALSE)
training <- d[intrain,]
testing <- d[-intrain,]

## control parameters in caret
trctrl <- trainControl(method = "repeatedcv", 
                       number = 10, 
                       repeats = 3)

## parameters
grid <- expand.grid(C = c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2,5))
set.seed(123)
svm_Linear_Grid <- train(y ~x2+x3+x6, data = training, 
                         method = "svmLinear",
                         trControl = trctrl,
                         tuneGrid = grid,
                         tuneLength = 10)
plot(svm_Linear_Grid)

svm_Linear_Grid

Support Vector Machines with Linear Kernel 

271 samples
  3 predictor
  2 classes: 'no', 'yes' 

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 3 times) 
Summary of sample sizes: 243, 243, 244, 245, 244, 244, ... 
Resampling results across tuning parameters:

  C     Accuracy   Kappa    
  0.01  0.7675587  0.5347604
  0.05  0.7650047  0.5205467
  0.10  0.7662393  0.5248014
  0.25  0.7700346  0.5285598
  0.50  0.7650522  0.5231597
  0.75  0.7649166  0.5211492
  1.00  0.7724529  0.5331103
  1.25  0.7698955  0.5312293
  1.50  0.7663241  0.5229310
  1.75  0.7676536  0.5258181
  2.00  0.7635836  0.5165059
  5.00  0.7626272  0.5166582

Accuracy was used to select the optimal model using the
 largest value.
The final value used for the model was C = 1.

## training
test_pred_grid <- predict(svm_Linear_Grid, newdata = training)
m <- confusionMatrix(test_pred_grid, factor(training$y, levels = c('no','yes')))$table
(mtraining <- m[2:1, 2:1])

          Reference
Prediction yes  no
       yes 132  33
       no   25  81

measures(mtraining)

sensitivity specificity         ppv         npv         fdr 
  0.8407643   0.7105263   0.8000000   0.7641509   0.2000000 
        fpr  class_rate        lift       delta 
  0.2894737   0.7859779   1.3808917   0.3303801 

#Accuracy
round(rft(mtraining, graph = FALSE)$auc*100,2)

[1] 77.56

## testing
test_pred_grid <- predict(svm_Linear_Grid, newdata = testing)
m <- confusionMatrix(test_pred_grid, factor(testing$y, levels = c('no','yes')))$table
(mtesting <- m[2:1, 2:1])

          Reference
Prediction yes no
       yes  53 12
       no   14 36

measures(mtesting)

sensitivity specificity         ppv         npv         fdr 
  0.7910448   0.7500000   0.8153846   0.7200000   0.1846154 
        fpr  class_rate        lift       delta 
  0.2500000   0.7739130   1.3995408   0.3258255 

#Accuracy
round(rft(mtesting, graph = FALSE)$auc*100,2)

[1] 77.05

Curva ROC

## nice colors
cols <- c("#0080ff", "#ff00ff", "darkgreen", "#ff0000", 
          "orange", "#00ff00", "brown")

par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
rft(mtraining, las = 1)$auc

[1] 0.7756453

rft(mtesting, las = 1, add = TRUE, line.col =cols[2])$auc

[1] 0.7705224

Random Forest

source('https://www.dropbox.com/s/gnwm7vuvhojg9q4/logistic-functions.R?dl=1')

Esta es una colección de funciones creadas 
por el profesor Jorge Vélez <https://jivelez.github.io/> 
para el ajuste y validación de modelos de Regresión Logística 
 
En caso de presentar algún inconveniente, por favor 
repórtelo a jvelezv@uninorte.edu.co 
 

## lectura de datos
d <- read.table("https://www.dropbox.com/s/bxea58hj951vi7k/datos-p1-caso2.txt?dl=1", 
                header = TRUE)

require(stringr)
d$y<-factor(d$y)
levels(d$y)

[1] "no"  "yes"

x<-d[,2:20]

## Subconjuntos Training y Testing
intrain <- createDataPartition(y = d$y, p= 0.7, list = FALSE)
training <- d[intrain,]
testing <- d[-intrain,]

## Control de parámetros en caret
control <- trainControl(method = "repeatedcv", 
                       number = 10, 
                       repeats = 3)

## Método de comparación  elegido es el AUC
metric <- "Accuracy"
set.seed(123)
mtry <- sqrt(ncol(x))
tunegrid <- expand.grid(.mtry=mtry)
rf_default <- train(y~x2+x3+x6, 
                    data=training, 
                    method='rf', 
                    metric='Accuracy', 
                    tuneGrid=tunegrid, 
                    trControl=control)
print(rf_default)

Random Forest 

271 samples
  3 predictor
  2 classes: 'no', 'yes' 

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 3 times) 
Summary of sample sizes: 243, 243, 244, 245, 244, 244, ... 
Resampling results:

  Accuracy   Kappa    
  0.8011396  0.5928968

Tuning parameter 'mtry' was held constant at a value of 4.358899

Determinamos medidas de clasificación del modelo sobre el traininig y el testing

## training
test_pred_RF <- predict(rf_default, newdata = training)
m <- confusionMatrix(test_pred_RF, factor(training$y, levels = c('no','yes')))$table
(mtraining <- m[2:1, 2:1])

          Reference
Prediction yes  no
       yes 129  26
       no   28  88

measures(mtraining)

sensitivity specificity         ppv         npv         fdr 
  0.8216561   0.7719298   0.8322581   0.7586207   0.1677419 
        fpr  class_rate        lift       delta 
  0.2280702   0.8007380   1.4365728   0.2895213 

rft(mtraining, graph = FALSE)$auc

[1] 0.7967929

## testing
test_pred_RFT <- predict(rf_default, 
                                 newdata = testing)
m <- confusionMatrix(test_pred_RFT, factor(testing$y, levels = c('no','yes')))$table
(mtesting <- m[2:1, 2:1])

          Reference
Prediction yes no
       yes  55 11
       no   12 37

measures(mtesting)

sensitivity specificity         ppv         npv         fdr 
  0.8208955   0.7708333   0.8333333   0.7551020   0.1666667 
        fpr  class_rate        lift       delta 
  0.2291667   0.8000000   1.4303483   0.2908535 

Modelo Final

rf_default$finalModel

Call:
 randomForest(x = x, y = y, mtry = param$mtry) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 3

        OOB estimate of  error rate: 19.93%
Confusion matrix:
    no yes class.error
no  88  26   0.2280702
yes 28 129   0.1783439

Realizamos la grafica ROC

## nice colors
cols <- c("#0080ff", "#ff00ff", "darkgreen", "#ff0000", 
          "orange", "#00ff00", "brown")

par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
rft(mtraining, las = 1)$auc

[1] 0.7967929

rft(mtesting, las = 1, add = TRUE, line.col =cols[2])$auc

[1] 0.7958644

De los modelos analizados el que mejor clasifica con un 86.8% es el modelo de regresión logística para el testing.