Introducción

A continuación explicaré a grandes rasgos lo que hice.

Primero hice análisis decriptivo sobre la base de datos, analicé la composición de la base, así como las correlaciones entre las variables numéricas y calculé los cuantiles de las variables numérocas para ver si las variables están en la misma escala o no; luego utilicé gráficos de caja y Componentes Principales para ver la distribución de los datos y ver si había datos atípicos. Luego apliqué un árbol de clasificación para ver si había alguna agrupación en los datos.

Luego, apliqué los modelos de clasificación en el siguiente orden: Bagging, Random Forest, Regresión Logística y LDA, para cada uno de ellos calculé las tasas de error por grupo y globales, tanto apartentes como las tasas sobre el conjunto de prueba. Además estimé la tasa de error con métodos de resmuestreo.

Por último hice la comparación de los tres métodos para ver cuál se equivoca menos y concluir.

Análisis Descriptivo.

La base de datos OJ contiene 1070 registros con 18 variables. Los registros son ventas de jugos de naranja, en los cuales el comprador eligió entre la marca Citrus Hill y la marca Minute Maid Orange. Las variables son mediciones de características del cliente y del producto. La base de datos está en la paquetería ISLR de R.

Cabe resaltar que no tomé en cuenta a todas las variables para el análiisis; en la base de datos hay una variable que indica el id de la tienda en que se hace la venta y otra variable que idica el número de semana en que se hizo la venta, éstas son las variables que eliminié para hacer el análisis. Así, quedan 16 variables, de las cuales 11 son numéricas y las demás son categóricas. Una de las cinco variables categóricas no es dummy mientras que las otras cuatro son dummies, es decir, que toman valores en el conjunto \(\{0,1\}\). A continuación, la descripción de las variables:

Variable Descripción
Purchase Marca del jugo(variable respuesta)MM/CH
PriceCH Precio del jugo para la marca CH
PriceMM Precio del jugo para la marca MM
DiscCH Descuento para la marca CH
DiscMM Descuento para la marca MM
SpecialCH Especial para la marca CH
SpecialMM Especial para la marca MM
LoyalCH Lealtad del cliente con la marca CH
SalePriceMM Precio final(descuento) para MM
SalePriceCH Precio final(descuento) para CH
PriceDiff Precio final MM menos precio final CH
Store7 Indica si la venta se hizo en Store 7
PctDiscMM Porcentaje de descuento para MM
PctDiscCH Porcentaje de descuento para CH
ListPriceDiff Precio de MM menos precio de CH
STORE Indica en cuál tienda se hizo la venta

A grandes rasgos, las variables numéricas son características del precio del jugo y las variables categóricas indican la tienda en que se hizo la venta y si existe especial. Las variable Purchase es la variable respuesta de la base, es una variable dummy. Las variables SpecialCH, SpecialMM y Store7 son variables dummies y la variable STORE tiene 5 categorías. Las demás variables son numéricas y tienen escalas distintas, por ejemplo, la media de la variable PctDiscCH es 0.027, mientras que la media de la variable PriceMM es de 2.085, casi 77 veces más que 0.027. También pude observar que hay variables que están correlacionadas.

Variable Variables Correlacionadas
PriceCH +PriceMM +SPCH
PriceMM +SPMM +PriceDiff
DiscCH -SPCH +PctDiscH
DiscMM -SPMM -PriceDiff +PctDiscMM
SPMM +PriceDiff -PctDiscMM
SPCH -PctDisH
PiceDiff -PtDiscMM

Donde: SPMM=SalePriceMM, SPCH=SalePriceCH, +: correlación positva, -: Correlación negativa y las variables correlacionadas si \(|\rho|\geq0.5\). En la tabla omití las repeticiones, por ejemplo, PriceMM está correlacionada con PriceCH por lo tanto no incluí a la variable PriceCH en las variables correlacionadas con PriceMM.

Las variables PriceCH y PriceMM tienen una correlación de 0.62, ésto es interesante ya que significa que si el precio de uno cambia entonces el del otro también.

Como parte del análisis descriptivo usé gráficos de caja y Componentes Principales para ver gráficamente si existía un patrón en la distribución de los datos o si había presencia de datos atípicos.

En los gráficos de caja de la Figura 1 es posible ver la distribución de las variables numéricas respecto a la marca que petenecen, las variables DiscCH, DiscMM, LoyalCH, PctDiscMM y PctDiscCH tienen una distribución distinta entre las dos marcas de jugo, en las demás variables la distribución es similar entre las marcas. En los gráficos de caja de las variables DiscCH, DiscMM, PctDiscMM y PctDiscCH hay presencia de datos atípidos, sin embargo no usé alguna tranformación para arreglarlo ya que no me parecen errores de registro.

Gráficos de caja de las variables numéricas respecto a su clase, en este caso, a la marca de jugo que pertenecen. Donde el color rojo son las centas de la marca CH y en azul las ventas de la marca MM

Gráficos de caja de las variables numéricas respecto a su clase, en este caso, a la marca de jugo que pertenecen. Donde el color rojo son las centas de la marca CH y en azul las ventas de la marca MM

Después calculé los componentes principales. Dado que las variables son mixtas tuve que transformar la variable STORE a dummy para poder calcular los componentes principales. Para hacer ésto usé la paquetería fastDummies que está en R. Tomando las variables dummies como numéricas apliqué el método Componentes Principales y grafiqué los primeros dos componentes y los últimos dos componentes los cuales están en la Figura 2.

Gráfica de los primeros dos componentes principales a la izquierda, a la derecha gráfica de los dos últimos componentes principales. Rojo: CH. Azul: MM

Gráfica de los primeros dos componentes principales a la izquierda, a la derecha gráfica de los dos últimos componentes principales. Rojo: CH. Azul: MM

Así, lo obtenido fue lo siguiente: en los gráficos de caja y en la gráfica de los primeros dos componentes principales se observa que hay datos atípicos, curiosamente los datos atípicos aparecen en los gráficos de caja de las variables correspondientes al descuento sobre el precio del jugo y aparecen por arriba de la mediana de los datos. Los decuentos de la marca CH llegan hasta 10 % ; mientras que los descuentos de la marca MM llegan hasta 40 %. En los gráficos de caja de éstas variables se observa que hay pocas ventas con descuentos mayores a 10 % y/o 40 %, es decir, que hubo pocas ventas con mucho descuento lo cual podría significar que había pocas ventas y eso los llevó a poner un descuento muy alto. Por eso mejor me parece que no son errores de registro.

También se observa que el descuento influye en la marca que se vendió, dicho de otra forma: si hay descuento en la marca CH entonces hay más ventas de la marca CH, lo mismo con la marca MM.

Ahora la distribución de los datos. La distribución de los datos es distinta, entre las dos marcas, en las variables que corresponden al descuento sobre el precio y en la variable de lealtad (LoyalCH). En la gráfica de los dos primeros componentes principales se observa que la distribución de los datos es mixta, es decir, que los datos están mezclados y que no se ve una diferencia marcada entre las dos categorías.

Además ajusté un árbol de clasificación, en la Figura 3 se puede apreciar. En principio se puede observar que la variable LoyalCH separa a las dos clases, ya que del lado izquierdo solo hay una rama que vale CH y las demás MM, del lado derecho aparece solo una rama que vale MM y las demás CH. Lo que se ve es que las observaciones tales que la variable LoyalCH tome valores menores que 0.5 serán las de la marca Minute Maid y las que son mayores o iguales a 0.5 serán Citrus Hill. En los gráficos de caja ya se había visto que la distribución de los datos, para la variable LoyalCH, es distinta respecto a las dos marcas.

Árbol de clasificación de tamaño 8, en las cuales cuatro ramas son MM y cuatro son CH.

Árbol de clasificación de tamaño 8, en las cuales cuatro ramas son MM y cuatro son CH.

Métodos de clasificación

Primero separé en un conjunto de entrenamiento y un conjunto de prueba, la base de datos tiene 1070 observaciones, utilicé el 75 % de la observacione en el conjunto de entrenamiento y el resto en el conjunto de prueba, quedando como resultado un conjnuot de entrenamiento con 802 observaciones.

Random Forest & Bagging

El primer modeló que utilicé fue el de Bagging, el cual apliqué con la paquetería randomForest que está en R, luego apliqué RandomForest con 3 variables por cada árbol (más adelante se verá porqué usé 3).

Para el conjunto de prueba ajusto el modelo Bagging 1000 veces usando a los 15 predictores, luego usando la función predict de R calculé la tasas de error aparentes y no aparentes por clase y global. Para el modelo Ranfom Forest utilicé 3 predictores para lo 1000 árboles ajustados, luego con la función predict calculé las tasas de error.

Método Error Marca CH Error Marca MM Error Global
Bagging Train 0.002 0.056 0.024
Bagging Test 0.138 0.255 0.179
Random Forest Train 0.052 0.173 0.101
Random Forest Test 0.109 0.223 0.149

Como se puede ver en la tabla, la tasa de error global de Random Forest sobre el conjunto de prueba es de 10 % el cual es menor que la tasa de error global de Bagging:14.9 %, sin embargo ambas tasas son muy bajas, ya se verá a continuación la comparación con otros modelos.

La ventaja de usar éstos métoso es su poder predictivo, además se puede estimar la tasa de error sobre el conjunto de prueba con el Out-of-Bag Error. El resultado se puede ver en la siguiente tabla:

Método Error Marca CH Error Marca MM Error Global
OOB-Bagging 0.17 0.28 0.21
OOB-Random Forest 0.14 0.29 0.20

A diferencia de la tabla pasada, aquí las tasas de error globales y sobre la marca MM solo difieren en 1 %, a excepción de la tasa de error sobre la marca CH, la cual es menjor para el método Random Forest

Otra ventaja de éstos métodos es que podemos ver la importancia que tienen las variables predictoras sobre la variable respuesta. En la figura 4 se puede ver que la variable con más importancia es la variable LoyalCH. La cual idica el porcentaje de lealtad del cliente a la marca CH, al ajustar un solo árbol de clasificación ya se había visto que ésta variable es la que separa a los datos en las dos marcas, a excepción de dos ramas.

Gráfica de importancia de las variables sobre la variable respuesta. Se escogieron 3 predictores. Hice Random Forest variando la cantidad de predictores usados en cada iteración y al usar 3 se minimizan las tasas de errores por grupo y global.

Regresión Logística

Para el modelo de Regresión Logística primero utilicé las 15 variables para predecir, sin embargo obtuve NA´s en summary, anova y drop1 del modelo, lo cual significa que hay multicolinealidad entre las variables, así que usando la función step que está en R calculé un modelo, tal que todas sus variables son significativas y no tiene NA´s, formado por las variables: PriceCH, PriceMM, DiscCH, DiscMM, LoyalCH, PctDiscMM, PctDiscCH y STORE. Éste es el que utilicé para predecir a las observaciones del conjunto de prueba. Con las probabilidades estimadas, si la proba es mayor que .5 entonces claifiqué en la marca MM. Los errores fueron los siguientes:

Regresión Logística Error Marca CH Error Marca MM Error Global
Train Error 0.13 0.24 0.17
Test Error 0.11 0.20 0.15
Bootstrap 0.12 0.25 0.17
Validation Approach 0.12 0.24 0.17

Para usar Bootstrap ocupé 1000 iteraciones, en cada una de ellas creé un conjunto de prueba a partir de una muestra con reemplazo de los números 1:1070, de tamaño 802, después ajusté el modelo con ese conjunto de prueba y predije con las observaciones restantes. Para Validation Approach hice lo mismo, solo que ocupé muestras sin reemplazo. La tasa de error tanto de Bootstrap como de Validation Approach sobre la Marca CH es menor en 1 % a la tasa de error sobre el conjunto de prueba, el remuestreo no disminuyó las tasas de error significativamente.

LDA

Para ajustar un modelo LDA tuve que ocupar la paquetería InformationValue que se encuentra en R para poder aplicar Weights of Evidence y transformar las variables categóricas a numéricas. También intenté ajustar el modelo sin transformar a las variables, pero sí usándolas como dummies y declarándolas como numéricas, sin embargo no fue bueno ésto ya que hay colinealidad entre las variables, así que mejor opté por usar Weights of Evidence. Las tasas de error fueron las siguientes.

LDA Error Marca CH Error Marca MM Error Global
Train Error 0.24 0.13 0.17
Test Error 0.20 0.12 0.15
Bootstrap 0.24 0.13 0.17
Validation Approach 0.24 0.13 0.17

Conclusión

Para concluir hay que ver cuál modelo clasificó mejor, para ésto tengo que ver las tasas de error sobre el conjunto de prueba y las tasas de error estimadas.

Método Error Marca CH Error Marca MM Error Global
Bagging Test 0.14 0.26 0.18
Random Forest Test 0.11 0.22 0.15
Test Error Reg Log 0.11 0.20 0.15
Test Error LDA 0.20 0.12 0.15

Los modelos Random Forest y Regresión Logística tienen la misma tasa de error para la marca CH y ésta es más baja que en los demás métodos. El modelo LDA tiene la menor tasa de error sobre la marca MM y quitando el modelo Bagging los tres restantes tienen la misma tasa de error global.

Método Error Marca CH Error Marca MM Error Global
OOB-Bagging 0.17 0.28 0.21
OOB-Random Forest 0.14 0.29 0.20
Bootstrap RG 0.12 0.25 0.17
Validation Approach RG 0.12 0.24 0.17
Bootstrap LDA 0.24 0.13 0.17
Validation Approach LDA 0.24 0.13 0.17

Ahora me voy a fijar en las tasas de error con remuestreo. Las tasas de error, más chicas, sobre la marca CH son las del modelo de Regresión Logística al aplicar Bootstrap y Validation Approach. Sin embargo el modelo LDA tiene menor tasa de error sobre la marca MM al aplicar Bootstrap y Validation Approach. Tanto Regresión Logística como LDA obtuvieron las mismas tasas de error global que es 17 %. Por lo tanto, queda a punto de vista de la persona que quiera clasificar, ya que si se quiere un modelo que clasifique bien a la marca CH yo utilizaría el modelo de Regresión Logística o los modelos Bagging y Random Forest , ya que éstos tienen tasas de error sobre la marca CH menor que las del modelo LDA. Sin embargo si se quiere clasificar bien a la marca MM mejor utilizaría el modelo LDA. Cabe resaltar que Random Forest y Bagging tienen tasas de error global más altas que los demás métodos.

Código

############ Proyecto final #############
#########################################
####
library(ISLR)
library(fastDummies)  #dummies
library(tree)  #árbol de clasificación
library(randomForest)
library(InformationValue) #Para WOE
options(scipen=100, digits=2)
data(OJ)
OJ<-OJ[,-c(2, 3)]; attach(OJ)
col=c("red", "blue")
## Visualización de datos
# Base con 1070 registros de ventas de jugos de naranja, en las cuales el comprador escogió
#entre la marca _Citrus Hill_ o _Minute Maid_
# Además se registró un número de carácterísticas tanto del comprados como del producto
str(OJ)
?OJ
par(mfrow=c(3,4))
for(i in c(2:5, 8:11, 13:15)){
  boxplot(OJ[,i]~Purchase, col=col, ylab=names(OJ)[i])    #boxplot de variable respuesta vs numéricas
}
par(mfrow=c(1,3))
boxplot(OJ[,4]~Purchase, col=col, ylab=names(OJ)[4])
boxplot(OJ[,5]~Purchase, col=col, ylab=names(OJ)[5])
boxplot(OJ[,8]~Purchase, col=col, ylab=names(OJ)[8])
summary(OJ[,c(2:5, 8:11, 13:15)])
# DiscCH, DiscMM, LolaylCH, PctDiscMM, PctDiscCH vairables donde hay diferencia en la marca 
# Hay pocos outliers en la base, no haré transformación para elmininarlos ya que me parece 
#que así deben de ser.
table(Purchase) #ventas por tipo
table(Purchase, OJ[, 7])
table(OJ[, 6])
table(Purchase, Store7)
table(STORE)
hist(DiscCH)
hist(DiscMM)
# Correlaciones
cor(OJ[, c(2:5, 8:11, 13:15)])
pairs(OJ[, c(2:5, 8:11, 13:15)], col=col[Purchase])
# Algunas variables relacionadas
# Las variables están en diferentes escalas, e.g. SalePriceCH vs PctDiscMM, máx: 2.29 y .4
apply(OJ[, c(2:5, 8:11, 13:15)], 2, mean)
# Otra característica de los datos: Hay 16 variables, 11 numéricas y 5 categóricas. De las 5 
#categóricas, una es la variabl
# respuesta, cuatro de ellas son dummies y la otra tiene 5 niveles que indican la tienda.
# Aplico PCA para ver si hay algún patron en los datos. 
######### Con variables dummys
# La variable "STORE" tiene 5 categorías que indican las tiendas donde ocurrieron las ventas
# La convierto en dummy, las otras ya son dummies
OJ_D<-dummy_cols(OJ, select_columns = c("STORE", "Store7"), remove_first_dummy = TRUE)
OJ_D<-OJ_D[,-c(1, 12, 16)]
str(OJ_D)
OJ_D$STORE_0<-as.numeric(OJ_D$STORE_0)
OJ_D$STORE_2<-as.numeric(OJ_D$STORE_2)
OJ_D$STORE_3<-as.numeric(OJ_D$STORE_3)
OJ_D$STORE_4<-as.numeric(OJ_D$STORE_4)
OJ_D$Store7_Yes<-as.numeric(OJ_D$Store7_Yes)
# variables categóricas: "SpecialCH", "SpecialMM", "Store7", y "STORE"
pr_d<-prcomp(OJ_D, scale=TRUE)
par(mfrow=c(1, 2))
plot(pr_d$x, col=col[Purchase], main = "PC con dummies") #CH- red #MM- Orange
plot(pr_d$x[,17], pr_d$x[,18], col=col[Purchase], xaxt="n", yaxt="n", ylab="PC18", 
     xlab="PC17", main = "PC con dummies")
# Como en los boxplots, los datos están poco separados
######## Solo con variables numéricas
OJ_C<-OJ[, -c(1, 6, 7, 12, 16)]
proj<-prcomp(OJ_C, scale=TRUE)
par(mfrow=c(1, 2))
plot(proj$x, col=col[Purchase], main = "PC sin categóricas") #CH- red #MM- blue
plot(proj$x[,10], proj$x[,11], col=col[Purchase], xaxt="n", yaxt="n", ylab="PC12",
     xlab="PC11", main= "PC sin categóricas")
biplot(proj, scale=0, col=c("blue", "red"))   
############# Árbol de clasif ###########
#######
## Variables como factores 
oj<-OJ
oj$SpecialCH<-as.factor(oj$SpecialCH)
oj$SpecialMM<-as.factor(oj$SpecialMM)
oj$Store7<-as.factor(oj$Store7)
oj$STORE<-as.factor(oj$STORE)
str(oj)
ojt<-tree(Purchase~., data=oj)
plot(ojt); text(ojt, pretty = 0)

#set.seed(101010)
#cv_tree<-cv.tree(ojt)
#plot(cv_tree$size, cv_tree$dev, type="b",xlab="Tamaño del árbol",
#ylab="MSE", col="blue") #El error sigue disminuyendo, pero voy a podar a 7 y con 11
#prune_simple_tree<-prune.tree(ojt, best= 6)
#plot(prune_simple_tree); text(prune_simple_tree, pretty=0)  

##### División entre conjunto de prueba 
#.75*1070=802
set.seed(1010)
train<-sample(1:1070, 802)
oj.test<-oj[-train,]
purchase.test<-Purchase[-train]
oj.train<-oj[train,]
#RF & Bagging
errors<-matrix(NA, ncol=3, nrow=8)    #Errores aparentes y no aparentes
colnames(errors)<-c("Error Clase CH", "Error Clase MM", "Error Global")
rownames(errors)<-c("Bagging Train", "Bagging Test", "RF Train", "RF Test", 
                    "RegLog Train", "RegLog Test",
                    "LDA Train", "LDA Test")
#Errores de remuestreo
est_err<-matrix(NA, ncol=3, nrow=6)   
colnames(est_err)<-c("Error Clase CH", "Error Clase MM", "Error Global")
rownames(est_err)<-c("OOB_Bagging", "OOB_RF", "Bootstrap RegLog", 
                     "Validation Approach RegLog",
                     "Bootstrap LDA", "Validation Approach LDA")

#######################################
################### Random Forest & Bagging
#Bagging 
set.seed(101010)
bagoj<-randomForest(Purchase~.,  data=oj, ntree=1000, mtry=15,
                    subset=train, importance=TRUE)
# Test Error
pred<-predict(bagoj, oj.test, type="class")
table(pred, purchase.test)
errors[2,1]<-round(1-150/(150+24),3)
errors[2,2]<-round(1-70/(24+70),3)
errors[2,3]<-round(1-(150+70)/(268),3)
#Train error
pred1<-predict(bagoj, oj.train, type="class")
table(pred1, Purchase[train])
errors[1,1]<-round(1-478/(478+1),3)
errors[1,2]<-round(1-305/(305+18),3)
errors[1,3]<-round(1-(478+305)/(1+478+18+305),3)
# OOB Error
est_err[1,3]<-mean(bagoj$err.rate[,1])
est_err[1,1]<-mean(bagoj$err.rate[,2])
est_err[1,2]<-mean(bagoj$err.rate[,3])
# Importance

varImpPlot(bagoj)

#RandomForest
best_tree<-function(m, test, real){
  set.seed(1010100)
  rf<-randomForest(Purchase~ ., data=oj, subset=train, importance=TRUE, 
                   ntree=1000, mtry=m)
  yhrf<-predict(rf, test, class="type") #árbol, datos a predecir
  t<-table(yhrf, real)
  print(c(1-t[1,1]/sum(t[,1]), 1-t[2,2]/sum(t[,2]), 1-sum(diag(t))/sum(t)))
  # Error marca CH, Error Marca MM, Error global
}
####### Obteniendo el tamaño de predictores
er<-matrix(NA, ncol=3, nrow = 15)
colnames(er)<-c("Error Clase CH", "Error Clase MM", "Error Global")
for(i in 1:15){
  b<-best_tree(i, test = oj.test, real = purchase.test)
  er[i,1]<-b[1]
  er[i,2]<-b[2]
  er[i,3]<-b[3]
}
x<-seq(1:15)
plot(x, er[,3], type="b") #3
plot(x, er[,2], type="b") #3
plot(x, er[,1], type = "b") #3-6

# RandomForest con m=3
set.seed(101010)
rfoj<-randomForest(Purchase~.,  data=oj, ntree=1000, mtry=3, 
                   subset=train, importance=TRUE)
# Test Error
pred<-predict(rfoj, oj.test, type="class")
table(pred, purchase.test)
errors[4,1]<-1-155/(155+19)
errors[4,2]<-1-73/(21+73)
errors[4,3]<-1-(155+73)/(268)
#Train error
pred1<-predict(rfoj, oj.train, type="class")
table(pred1, Purchase[train])
errors[3,1]<-1-454/(454+25)
errors[3,2]<-1-267/(267+56)
errors[3,3]<-1-(454+267)/(1+478+18+305)
# OOB Error
est_err[2,3]<-mean(rfoj$err.rate[,1])
est_err[2,1]<-mean(rfoj$err.rate[,2])
est_err[2,2]<-mean(rfoj$err.rate[,3])
varImpPlot(rfoj)

#################
### Regresión Logística
#############
m1 = glm(Purchase ~ PriceCH+PriceMM+DiscCH+DiscMM+SpecialCH+SpecialMM+
           LoyalCH+SalePriceCH+SalePriceMM+PriceDiff+
           Store7+PctDiscCH+PctDiscMM+ListPriceDiff+STORE, 
         data = oj, family = binomial(link = "logit"), subset = train)
summary(m1); anova(m1, test = "Chisq"); drop1(m1, test = "Chisq")
#Modelo bueno
m2=glm(Purchase ~ PriceCH + PriceMM + DiscCH + DiscMM + 
         LoyalCH + PctDiscMM + STORE, family = binomial(link = "logit"),  
       data = oj, subset = train)
summary(m2); anova(m2, test = "Chisq"); drop1(m2, test = "Chisq")

#Test error
logpred<-predict(m2, oj.test, type="response") #Estima las probabilidades
prob_test<-ifelse(logpred>0.5, "MM", "CH") 
t_<-table(prob_test, purchase.test)
errors[6,1]<-1-t_[1,1]/sum(t_[,1])
errors[6,2]<-1-t_[2,2]/sum(t_[,2])
errors[6,3]<-1-sum(diag(t_))/sum(t_)
#Train error
logpr<-predict(m2, oj.train, type="response") #Estima las probabilidades
prob_train<-ifelse(logpr>0.5, "MM", "CH")
t<-table(prob_train, Purchase[train])
errors[5,1]<-1-t[1,1]/sum(t[,1])
errors[5,2]<-1-t[2,2]/sum(t[,2])
errors[5,3]<-1-sum(diag(t))/sum(t)
# Estimar el error
# Bootstrap (conjuntos de entrenamiento con reemplazo)
tr<-matrix(NA, ncol=3, nrow=1000)
colnames(tr)<-c("Error Clase CH", "Error Clase MM", "Error Global")
for(i in 1:1000){
  set.seed(10+i)
  tra<-sample(1:1070, 802, replace = T)
  p<-glm(Purchase ~ PriceCH + PriceMM + DiscCH + DiscMM + 
        LoyalCH + PctDiscMM + STORE, family = binomial(link = "logit"),  
      data = oj, subset = tra)
  pre<-predict(p, oj[-tra,], type="response")
  prob<-ifelse(pre>.5, "MM", "CH")
  tt<-table(prob, Purchase[-tra])
  tr[i,1]<-1-tt[1,1]/sum(tt[,1])
  tr[i,2]<-1-tt[2,2]/sum(tt[,2])
  tr[i,3]<-1-sum(diag(tt))/sum(tt)
}
est_err[3,1]<-mean(tr[,1])
est_err[3,2]<-mean(tr[,2])
est_err[3,3]<-mean(tr[,3])
# Validation Approach (conjuntos de entrenamiento sin reemplazo)
tra<-matrix(NA, ncol=3, nrow=1000)
colnames(tra)<-c("Error Clase CH", "Error Clase MM", "Error Global")
for(i in 1:1000){
  set.seed(10+i)
  trai<-sample(1:1070, 802, replace = F)
  p<-glm(Purchase ~ PriceCH + PriceMM + DiscCH + DiscMM + 
           LoyalCH + PctDiscMM + STORE, family = binomial(link = "logit"),  
         data = oj, subset = trai)
  pre<-predict(p, oj[-trai,], type="response")
  prob<-ifelse(pre>.5, "MM", "CH")
  tt<-table(prob, Purchase[-trai])
  tra[i,1]<-1-tt[1,1]/sum(tt[,1])
  tra[i,2]<-1-tt[2,2]/sum(tt[,2])
  tra[i,3]<-1-sum(diag(tt))/sum(tt)
}
est_err[4,1]<-mean(tra[,1])
est_err[4,2]<-mean(tra[,2])
est_err[4,3]<-mean(tra[,3])
errors
###########################
############ LDA
# Wights of Evidence

oj_WOE<-oj
str(oj_WOE)
oj_WOE$Purchase<-ifelse(Purchase=="CH", 1, 0)   #CH:1 MM:0
oj_WOE$SpecialCH<-WOE(X=oj_WOE$SpecialCH, Y=oj_WOE$Purchase)
oj_WOE$SpecialMM<-WOE(X=oj_WOE$SpecialMM, Y=oj_WOE$Purchase)
oj_WOE$Store7<-WOE(X=oj_WOE$Store7, Y=oj_WOE$Purchase)
oj_WOE$STORE <-WOE(X=oj_WOE$STORE, Y=oj_WOE$Purchase)

woe.test<-oj_WOE[-train,]
pr.test<-oj_WOE$Purchase[-train]
woe.train<-oj_WOE[train,]

ldaoj<-lda(Purchase ~ ., data = oj_WOE, subset = train)
# Test Errors
ldap<-predict(ldaoj, woe.test, type="response")
tl<-table(ldap$class, pr.test)
errors[8,1]<-1-tl[1,1]/sum(tl[,1])
errors[8,2]<-1-tl[2,2]/sum(tl[,2])
errors[8,3]<-1-sum(diag(tl))/sum(tl)
#Train Errors
ldap1<-predict(ldaoj, woe.train, type="response")
tl1<-table(ldap1$class, oj_WOE$Purchase[train])
errors[7,1]<-1-tl1[1,1]/sum(tl1[,1])
errors[7,2]<-1-tl1[2,2]/sum(tl1[,2])
errors[7,3]<-1-sum(diag(tl1))/sum(tl1)
# Bootstrap (conjuntos de entrenamiento con reemplazo)
tr<-matrix(NA, ncol=3, nrow=1000)
colnames(tr)<-c("Error Clase CH", "Error Clase MM", "Error Global")
for(i in 1:1000){
  set.seed(10+i)
  tra<-sample(1:1070, 802, replace = T)
  p<-lda(Purchase ~ ., data = oj_WOE, subset = tra)
  pre<-predict(p, oj_WOE[-tra,], type="response")
  tt<-table(pre$class, oj_WOE$Purchase[-tra])
  tr[i,1]<-1-tt[1,1]/sum(tt[,1])
  tr[i,2]<-1-tt[2,2]/sum(tt[,2])
  tr[i,3]<-1-sum(diag(tt))/sum(tt)
}
est_err[5,1]<-mean(tr[,1])
est_err[5,2]<-mean(tr[,2])
est_err[5,3]<-mean(tr[,3])
# Validation Approach
tra<-matrix(NA, ncol=3, nrow=1000)
colnames(tra)<-c("Error Clase CH", "Error Clase MM", "Error Global")
for(i in 1:1000){
  set.seed(10+i)
  trai<-sample(1:1070, 802, replace = F)
  p<-lda(Purchase ~ ., data = oj_WOE, subset = trai)
  pre<-predict(p, oj_WOE[-trai,], type="response")
  tt<-table(pre$class, oj_WOE$Purchase[-trai])
  tra[i,1]<-1-tt[1,1]/sum(tt[,1])
  tra[i,2]<-1-tt[2,2]/sum(tt[,2])
  tra[i,3]<-1-sum(diag(tt))/sum(tt)
}
est_err[6,1]<-mean(tra[,1])
est_err[6,2]<-mean(tra[,2])
est_err[6,3]<-mean(tra[,3])