Estadística Multivariada Examen 1
Mario Becerra 124362
Felipe Domínguez 127041
05/03/2015
Descripción de la base de datos
Los datos están relacionados con las campañas de marketing de una institución bancaria de Portugal. Abarca 17 campañas que ocurrieron entre mayo 2008 y noviembre 2010, las campañas consisten en llamadas a clientes sumando 40,690 contactos en esta base de datos. En las llamadas se ofrecía un contrato a largo plazo para invertir dinero, con tasas de interés atractivas (de acuerdo al banco). En cada contacto se registraban 16 atributos que incluyen la respuesta del cliente, esto es si el cliente contrató el producto o no.
Las variables son las siguientes:
La institución bancaria quisiera conocer como están relacionadas estas variables, en particular, quisieran entender la relación con la variable respuesta y.
library(dplyr)
library(ggplot2)
library(bnlearn)
library(gRain)
library(tidyr)
library(SDMTools)
load('./Datos/mkt.Rdata')Dividimos la muestra en un conjunto de prueba y uno de entrenamiento.
set.seed(18937)
mkt_train <- sample_n(mkt, 30000)
mkt_test <- anti_join(mkt, mkt_train)## Joining by: c("age", "job", "marital", "education", "default", "balance", "housing", "loan", "contact", "month", "duration", "campaign", "previous", "poutcome", "y", "id")Tenemos en la descripción de la base de datos que hay algunas variables numéricas. Vamos a categorizar estas variables, pues si hay pocas observaciones sobre algún número en particular, puede haber estimaciones ruidosas. Esto se puede ver gráficamente en los histogramas siguientes.
qplot(mkt_train$age) +
ggtitle("Distribución de la edad") +
theme(plot.title = element_text(lineheight=.8, face="bold")) +
xlab("Edad") +
ylab("Num. de personas")
qplot(mkt_train$balance) +
ggtitle("Distribución del balance") +
theme(plot.title = element_text(lineheight=.8, face="bold")) +
xlab("Balance en euros") +
ylab("Num. de personas")
qplot(mkt_train$duration) +
ggtitle("Distribución de la duración") +
theme(plot.title = element_text(lineheight=.8, face="bold")) +
xlab("Duración") +
ylab("Num. de personas")
qplot(mkt_train$campaign) +
ggtitle("Distribución de las veces que se \n contactó con el cliente durante la campaña") +
theme(plot.title = element_text(lineheight=.8, face="bold")) +
xlab("Veces que se contactó con el cliente") +
ylab("Num. de personas")
qplot(mkt_train$previous) +
ggtitle("Distribución de las veces que se había \n contactado con el cliente antes de la campaña") +
theme(plot.title = element_text(lineheight=.8, face="bold")) +
xlab("Veces que se contactó con el cliente") +
ylab("Num. de personas")
Como se puede apreciar en los histogramas, todas las variables están sesgadas hacia algún lado y tienen valores extremos con pocas observaciones.
Para evitar esto, vamos a categorizar cada una de las variables anteriores en intervalos tales que la cantidad de clientes en cada intervalo sea más o menos igual.
A age se le dividirá en 15 grupos de edad, de esta forma hay una distancia promedio de 5.13 años de diferencia entre en grupo y otro.
A balance se le dividirá en 15 grupos, de esta forma hay una distancia promedio de 7343.07 euros de diferencia entre en grupo y otro.
A duration se le dividirá en 5 grupos, de esta forma hay una distancia promedio de 983.6 segundos (4.1 minutos) de diferencia entre en grupo y otro.
A campaign se le dividirá en 3 grupos.
A previous se le dividirá en 2 grupos, ahora dirá si se había contactado con el cliente en la campaña anterior o no. Tendrá el valor TRUE si había sido contactado anteriormente y FALSE si no.
mkt_train$age_cat <- cut(mkt_train$age, quantile(mkt_train$age, probs=(0:15)/15))
mkt_train$balance_cat <- cut(mkt_train$balance, quantile(mkt_train$balance, probs=(0:15)/15))
mkt_train$duration_cat <- cut(mkt_train$duration, quantile(mkt_train$duration, probs=(0:5)/5))
mkt_train$campaign_cat <- cut(mkt_train$campaign, c(0,1,3,max(mkt_train$campaign)))
mkt_train$previous_cat <- factor(mkt_train$previous!=0)
train <- mkt_train[,sapply(mkt_train, class)=='factor']
names(train) <- c("job","marital","education","default","housing","loan","contact","month","poutcome","y", "age","balance", "duration","campaign","previous")
train <- train[complete.cases(train),]Vamos a aprender la estructura de una red bayesiana a partir de la muestra de entrenamiento. Debido a que la variable de interés es si contrató el plan o no (\(y\)), podemos considerarla como la variable respuesta, por lo que en la red bayesiana vamos a ‘obligar’ al algoritmo a que no le asigne hijos a la \(y\) en la red. Ajustaremos primero dos redes con distinto score, una con AIC y otra con BIC.
black <- expand.grid('y',names(train))
net1aic <- hc(train, score = 'aic', blacklist = black)
strength1aic <- arc.strength(net1aic, train)
strength.plot(net1aic, strength1aic, shape='ellipse', main='AIC')## Loading required namespace: Rgraphviz
net1bic <- hc(train, score = 'bic', blacklist = black)
strength1bic <- arc.strength(net1bic, train)
strength.plot(net1bic, strength1bic, shape='ellipse', main='BIC')
Tal vez lo primero que salta a la vista es que con BIC hay menos relaciones que con AIC, esto debido a que el primero es más ‘estricto’ con las relaciones y deja menos relaciones espurias en el modelo.
Vemos de la red resultante del BIC que las variables que afectan directamente a \(y\) son duration, poutcome y housing; es decir, según esta red, lo que afecta directamente a que acepten el contrato o no es el hecho de que tengan un préstamo de casa, cuánto tiempo duró la llamada de ofrecimiento del servicio y el resultado de la campaña de marketing anterior. Mientras que en AIC se agrega la variable la variable month.
Ahora, para tener redes más robustas, generaremos 1000 muestras bootstrap del conjunto de entrenamiento y creamos una red para cada muestra, después se crea una red promedio. Esto ayuda a que las relaciones que se mantienen sean más fuertes, pues se reduce el impacto de los máximos locales en el proceso de aprender la estructura. La red promedio usa arcos con fuerza mayor a \(0.8\), esto significa que se usan los arcos que están presentes al menos un 80% de las redes creadas.
net1aicboot <- boot.strength(train, R=100, algorithm='hc', algorithm.args=list(score='aic', blacklist=black), cpdag=FALSE)
net1aicavg<-averaged.network(net1aicboot, names(train), 0.8)
strength1aicavg <- arc.strength(net1aicavg, train)
strength.plot(net1aicavg, strength1aicavg, shape='ellipse', main='AIC promedio bootstrap')
net1bicboot <- boot.strength(train, R=100, algorithm='hc', algorithm.args=list(score='bic', blacklist=black), cpdag=FALSE)
net1bicavg<-averaged.network(net1bicboot, names(train), 0.8)
strength1bicavg <- arc.strength(net1bicavg, train)
strength.plot(net1bicavg, strength1bicavg, shape='ellipse', main='BIC promedio bootstrap')
Vemos ahora en esta nueva red que las flechas son todas más grusas, esto indica que la fuerza de los arcos es mayor ahora; y se puede ver que las variables que afectaban directamente a la \(y\) no han cambiado; tanto con BIC como con AIC.
A continuación se pueden ver los arcos con fuerza igual a 1 (el máximo) para el criterio con AIC.
arrange(net1aicboot[net1aicboot$strength==1,], from)## from to strength direction
## 1 age job 1 0.00
## 2 age marital 1 1.00
## 3 age education 1 1.00
## 4 age month 1 0.00
## 5 balance default 1 1.00
## 6 balance loan 1 0.46
## 7 balance month 1 0.00
## 8 campaign contact 1 0.17
## 9 campaign month 1 1.00
## 10 campaign poutcome 1 1.00
## 11 campaign duration 1 0.00
## 12 contact month 1 0.83
## 13 contact duration 1 0.83
## 14 contact campaign 1 0.83
## 15 default balance 1 0.00
## 16 duration contact 1 0.17
## 17 duration poutcome 1 1.00
## 18 duration y 1 1.00
## 19 duration campaign 1 1.00
## 20 education job 1 0.00
## 21 education age 1 0.00
## 22 housing marital 1 1.00
## 23 housing loan 1 1.00
## 24 housing month 1 0.00
## 25 housing y 1 1.00
## 26 job marital 1 1.00
## 27 job education 1 1.00
## 28 job month 1 0.00
## 29 job age 1 1.00
## 30 loan housing 1 0.00
## 31 loan month 1 0.00
## 32 loan balance 1 0.54
## 33 marital job 1 0.00
## 34 marital housing 1 0.00
## 35 marital age 1 0.00
## 36 month job 1 1.00
## 37 month housing 1 1.00
## 38 month loan 1 1.00
## 39 month contact 1 0.17
## 40 month poutcome 1 0.00
## 41 month y 1 1.00
## 42 month age 1 1.00
## 43 month balance 1 1.00
## 44 month campaign 1 0.00
## 45 poutcome month 1 1.00
## 46 poutcome y 1 1.00
## 47 poutcome duration 1 0.00
## 48 poutcome campaign 1 0.00
## 49 poutcome previous 1 1.00
## 50 previous poutcome 1 0.00
## 51 y housing 1 0.00
## 52 y month 1 0.00
## 53 y poutcome 1 0.00
## 54 y duration 1 0.00Para BIC.
arrange(net1bicboot[net1bicboot$strength==1,], from)## from to strength direction
## 1 age job 1 0.93
## 2 age marital 1 0.07
## 3 balance default 1 0.99
## 4 balance housing 1 0.00
## 5 balance loan 1 0.00
## 6 campaign month 1 0.02
## 7 campaign duration 1 0.98
## 8 contact month 1 0.90
## 9 default balance 1 0.01
## 10 duration y 1 1.00
## 11 duration campaign 1 0.02
## 12 education job 1 0.00
## 13 housing job 1 0.94
## 14 housing loan 1 1.00
## 15 housing month 1 0.00
## 16 housing balance 1 1.00
## 17 job education 1 1.00
## 18 job housing 1 0.06
## 19 job age 1 0.07
## 20 loan housing 1 0.00
## 21 loan month 1 0.00
## 22 loan balance 1 1.00
## 23 marital age 1 0.93
## 24 month housing 1 1.00
## 25 month loan 1 1.00
## 26 month contact 1 0.10
## 27 month campaign 1 0.98
## 28 poutcome previous 1 1.00
## 29 previous poutcome 1 0.00
## 30 y duration 1 0.00Como ya se había mencionado, el criterio BIC es más estricto, por lo que es natural que haya menos relaciones con fuerza igua a 1 usando BIC que usando AIC; y se puede ver que muchas de las relaciones de BIC están en AIC; de hecho, el número e elementos en la intersección de las dos tablas es 26.
Después de todo este análisis inicial de redes, decidiremos usar la red creada con el criterio del BIC, pues tiene menos relaciones, y como el objetivo final es describir el comportamiento de los clientes, no necesitamos de muchas relaciones que tal vez intuitivamente no tengan mucho sentido.
<<<<<<< HEAD Procederemos ahora a estimar los parámetros de la red. Podríamos primero probar con máxima verosimilitud, sin ninguna especie de suavizamiento ni regularización, pero si se observa la tabla siguiente (la cual muestra los conteos de la variable respuesta \(y\) condicionada a los nodos que apuntan a ella) se nota que muy seguramente habrá estimaciones ruidosas debido a que algunas probabilidades condicionales tienen pocas o ninguna obervación. ======= Procederemos ahora a estimar los parámetros de la red. Vamos primero a probar con máxima verosimilitud, sin ninguna especie de suavizamiento ni regularización, aunque si se observa la tabla siguiente (la cual muestra los conteos de la variable respuesta \(y\) condicionada a los nodos que apuntan a ella) se nota que muy seguramente habrá estimaciones ruidosas debido a que algunas probabilidades condicionales tienen pocas observaciones. >>>>>>> Mario
data.frame(table(train$y, train$duration, train$poutcome, train$housing))## Var1 Var2 Var3 Var4 Freq
## 1 no (0,90] failure no 147
## 2 yes (0,90] failure no 4
## 3 no (90,147] failure no 188
## 4 yes (90,147] failure no 22
## 5 no (147,223] failure no 146
## 6 yes (147,223] failure no 41
## 7 no (223,370] failure no 161
## 8 yes (223,370] failure no 68
## 9 no (370,4.92e+03] failure no 89
## 10 yes (370,4.92e+03] failure no 93
## 11 no (0,90] other no 55
## 12 yes (0,90] other no 1
## 13 no (90,147] other no 64
## 14 yes (90,147] other no 9
## 15 no (147,223] other no 56
## 16 yes (147,223] other no 18
## 17 no (223,370] other no 54
## 18 yes (223,370] other no 42
## 19 no (370,4.92e+03] other no 56
## 20 yes (370,4.92e+03] other no 57
## 21 no (0,90] success no 22
## 22 yes (0,90] success no 7
## 23 no (90,147] success no 54
## 24 yes (90,147] success no 38
## 25 no (147,223] success no 61
## 26 yes (147,223] success no 112
## 27 no (223,370] success no 38
## 28 yes (223,370] success no 167
## 29 no (370,4.92e+03] success no 30
## 30 yes (370,4.92e+03] success no 156
## 31 no (0,90] unknown no 2434
## 32 yes (0,90] unknown no 22
## 33 no (90,147] unknown no 2264
## 34 yes (90,147] unknown no 103
## 35 no (147,223] unknown no 1946
## 36 yes (147,223] unknown no 193
## 37 no (223,370] unknown no 1793
## 38 yes (223,370] unknown no 297
## 39 no (370,4.92e+03] unknown no 1419
## 40 yes (370,4.92e+03] unknown no 787
## 41 no (0,90] failure yes 488
## 42 yes (0,90] failure yes 2
## 43 no (90,147] failure yes 456
## 44 yes (90,147] failure yes 18
## 45 no (147,223] failure yes 485
## 46 yes (147,223] failure yes 20
## 47 no (223,370] failure yes 424
## 48 yes (223,370] failure yes 40
## 49 no (370,4.92e+03] failure yes 278
## 50 yes (370,4.92e+03] failure yes 104
## 51 no (0,90] other yes 207
## 52 yes (0,90] other yes 1
## 53 no (90,147] other yes 136
## 54 yes (90,147] other yes 6
## 55 no (147,223] other yes 146
## 56 yes (147,223] other yes 10
## 57 no (223,370] other yes 127
## 58 yes (223,370] other yes 12
## 59 no (370,4.92e+03] other yes 108
## 60 yes (370,4.92e+03] other yes 52
## 61 no (0,90] success yes 24
## 62 yes (0,90] success yes 2
## 63 no (90,147] success yes 30
## 64 yes (90,147] success yes 13
## 65 no (147,223] success yes 28
## 66 yes (147,223] success yes 30
## 67 no (223,370] success yes 38
## 68 yes (223,370] success yes 57
## 69 no (370,4.92e+03] success yes 25
## 70 yes (370,4.92e+03] success yes 59
## 71 no (0,90] unknown yes 2677
## 72 yes (0,90] unknown yes 10
## 73 no (90,147] unknown yes 2515
## 74 yes (90,147] unknown yes 24
## 75 no (147,223] unknown yes 2628
## 76 yes (147,223] unknown yes 37
## 77 no (223,370] unknown yes 2597
## 78 yes (223,370] unknown yes 85
## 79 no (370,4.92e+03] unknown yes 1970
## 80 yes (370,4.92e+03] unknown yes 706fit_mle <- bn.fit(net1bicavg, data = train, method = 'mle')Ahora las probabilidades estimadas por máxima verosimilitud.
probs_fit_mle <- data.frame(fit_mle[['y']]$prob)
probs_fit_mle## y housing poutcome duration Freq
## 1 no no failure (0,90] 0.973509934
## 2 yes no failure (0,90] 0.026490066
## 3 no yes failure (0,90] 0.995918367
## 4 yes yes failure (0,90] 0.004081633
## 5 no no other (0,90] 0.982142857
## 6 yes no other (0,90] 0.017857143
## 7 no yes other (0,90] 0.995192308
## 8 yes yes other (0,90] 0.004807692
## 9 no no success (0,90] 0.758620690
## 10 yes no success (0,90] 0.241379310
## 11 no yes success (0,90] 0.923076923
## 12 yes yes success (0,90] 0.076923077
## 13 no no unknown (0,90] 0.991042345
## 14 yes no unknown (0,90] 0.008957655
## 15 no yes unknown (0,90] 0.996278377
## 16 yes yes unknown (0,90] 0.003721623
## 17 no no failure (90,147] 0.895238095
## 18 yes no failure (90,147] 0.104761905
## 19 no yes failure (90,147] 0.962025316
## 20 yes yes failure (90,147] 0.037974684
## 21 no no other (90,147] 0.876712329
## 22 yes no other (90,147] 0.123287671
## 23 no yes other (90,147] 0.957746479
## 24 yes yes other (90,147] 0.042253521
## 25 no no success (90,147] 0.586956522
## 26 yes no success (90,147] 0.413043478
## 27 no yes success (90,147] 0.697674419
## 28 yes yes success (90,147] 0.302325581
## 29 no no unknown (90,147] 0.956485002
## 30 yes no unknown (90,147] 0.043514998
## 31 no yes unknown (90,147] 0.990547460
## 32 yes yes unknown (90,147] 0.009452540
## 33 no no failure (147,223] 0.780748663
## 34 yes no failure (147,223] 0.219251337
## 35 no yes failure (147,223] 0.960396040
## 36 yes yes failure (147,223] 0.039603960
## 37 no no other (147,223] 0.756756757
## 38 yes no other (147,223] 0.243243243
## 39 no yes other (147,223] 0.935897436
## 40 yes yes other (147,223] 0.064102564
## 41 no no success (147,223] 0.352601156
## 42 yes no success (147,223] 0.647398844
## 43 no yes success (147,223] 0.482758621
## 44 yes yes success (147,223] 0.517241379
## 45 no no unknown (147,223] 0.909770921
## 46 yes no unknown (147,223] 0.090229079
## 47 no yes unknown (147,223] 0.986116323
## 48 yes yes unknown (147,223] 0.013883677
## 49 no no failure (223,370] 0.703056769
## 50 yes no failure (223,370] 0.296943231
## 51 no yes failure (223,370] 0.913793103
## 52 yes yes failure (223,370] 0.086206897
## 53 no no other (223,370] 0.562500000
## 54 yes no other (223,370] 0.437500000
## 55 no yes other (223,370] 0.913669065
## 56 yes yes other (223,370] 0.086330935
## 57 no no success (223,370] 0.185365854
## 58 yes no success (223,370] 0.814634146
## 59 no yes success (223,370] 0.400000000
## 60 yes yes success (223,370] 0.600000000
## 61 no no unknown (223,370] 0.857894737
## 62 yes no unknown (223,370] 0.142105263
## 63 no yes unknown (223,370] 0.968307233
## 64 yes yes unknown (223,370] 0.031692767
## 65 no no failure (370,4.92e+03] 0.489010989
## 66 yes no failure (370,4.92e+03] 0.510989011
## 67 no yes failure (370,4.92e+03] 0.727748691
## 68 yes yes failure (370,4.92e+03] 0.272251309
## 69 no no other (370,4.92e+03] 0.495575221
## 70 yes no other (370,4.92e+03] 0.504424779
## 71 no yes other (370,4.92e+03] 0.675000000
## 72 yes yes other (370,4.92e+03] 0.325000000
## 73 no no success (370,4.92e+03] 0.161290323
## 74 yes no success (370,4.92e+03] 0.838709677
## 75 no yes success (370,4.92e+03] 0.297619048
## 76 yes yes success (370,4.92e+03] 0.702380952
## 77 no no unknown (370,4.92e+03] 0.643245694
## 78 yes no unknown (370,4.92e+03] 0.356754306
## 79 no yes unknown (370,4.92e+03] 0.736173393
## 80 yes yes unknown (370,4.92e+03] 0.263826607Para evitar estimaciones ruidosas, se utilizará un enfoque bayesiano con distintos tamaños de muestra imaginarios (\(5\), \(100\) y \(1000\)), los cuales, intuitivamente, son los pesos que le asigna a la distribución previa. Después comparamos el error cuadrático entre ellos para ver qué tanto cambian las estimaciones para los tamaños distintos de muestra (notar que esta comparación solo se hace para la variable de interés \(y\)).
fit1 <- bn.fit(net1bicavg, data = train, method = 'bayes', iss=5)
probs_fit1 <- data.frame(fit1[['y']]$prob)
fit2 <- bn.fit(net1bicavg, data = train, method = 'bayes', iss=100)
probs_fit2 <- data.frame(fit2[['y']]$prob)
fit3 <- bn.fit(net1bicavg, data = train, method = 'bayes', iss=1000)
probs_fit3 <- data.frame(fit3[['y']]$prob)
dif <- c(sum((probs_fit1$Freq - probs_fit2$Freq)^2),sum((probs_fit1$Freq - probs_fit3$Freq)^2),sum((probs_fit2$Freq - probs_fit3$Freq)^2))Las diferencias cuadráticas en las estimaciones son:
dif## [1] 0.005352825 0.256581331 0.188980325Comparamos ahora las estimaciones entre los distintos tamaños de muestra y el de máxima verosimilitud. Eso se puede ver en la gráfica siguiente.
probs <- full_join(probs_fit1, probs_fit2, by=c("y", "housing", "poutcome", "duration"))
probs <- full_join(probs, probs_fit3, by=c("y", "housing", "poutcome", "duration"))
probs <- full_join(probs, probs_fit_mle, by=c("y", "housing", "poutcome", "duration"))
names(probs) <- c("y","housing","poutcome", "duration", "iss5","iss100", "iss1000", "mle")
probs_l <- gather(probs, metodo, est,iss5,iss100, iss1000, mle)
probs2 <- full_join(probs_fit1, probs_fit_mle, by=c("y", "housing", "poutcome", "duration"))
probs2_l <- gather(probs, metodo, est,iss1000,mle)ggplot(probs_l, aes(x = y, y = est, colour = metodo)) +
geom_jitter(size=3,position=position_jitter(width=0.1, height=0))+
facet_wrap(~ duration+poutcome) +
geom_point(aes(shape=housing))
Se pueden ver disintas grupaciones, lo cual significa que la estimaciones no difieren mucho entre sí. Debido a esto y a que las diferencias entre los tamaños de muestra imaginario no son muy grandes, nos quedamos con las probabilidades estimadas por el tamaño de muestra más grande (\(1000\)).
Las comparaciones entre el modelo con muestra imaginaria de \(1000\) y con la estimación de máxima verosimilitud se pueden ver en la siguiente gráfica.
ggplot(probs2_l, aes(x = y, y = est, colour = metodo)) +
geom_jitter(size=3,position=position_jitter(width=0.1, height=0))+
facet_wrap(~ duration+poutcome) +
geom_point(aes(shape=housing))
Ahora vamos a predecir el valor de \(y\) en la muestra de prueba usando las probabilidades estimadas por el modelo elegido.
mkt_test$age_cat <- cut(mkt_test$age, quantile(mkt_train$age, probs=(0:15)/15))
mkt_test$balance_cat <- cut(mkt_test$balance, quantile(mkt_train$balance, probs=(0:15)/15))
mkt_test$duration_cat <- cut(mkt_test$duration, quantile(mkt_train$duration, probs=(0:5)/5))
mkt_test$campaign_cat <- cut(mkt_test$campaign, c(0,1,3,max(mkt_train$campaign)))
mkt_test$previous_cat <- factor(mkt_test$previous!=0)
test <- mkt_test[,sapply(mkt_test, class)=='factor']
names(test) <- c("job","marital","education","default","housing","loan","contact","month","poutcome","y", "age","balance", "duration","campaign","previous")
test <- train[complete.cases(test),]
net<-compile(as.grain(fit3))## Warning in as.grain.bn.fit(fit3): the gRain package does not support
## ordinal networks, disregarding the ordering of the levels.pred <- predict.grain(net, 'y', names(train)[-10],test)Vemos la matriz de confusión de la predicción, y el estadístico kappa de Cohen.
pred_1 <- factor(pred$pred$y)
levels(pred_1)<-c(0,1)
obs<-factor(test$y)
levels(obs) <- c(0,1)
cm <- confusion.matrix(obs, pred_1)## Warning in Ops.factor(pred, threshold): '>=' not meaningful for factors## Warning in Ops.factor(pred, threshold): '<' not meaningful for factorscm## obs
## pred 0 1
## 0 26090 2792
## 1 365 731
## attr(,"class")
## [1] "confusion.matrix"Kappa(cm)## [1] 0.2761494Se puede ver en la matriz de confusión que la predicción está muy cargada hacia los \(0\)’s, es decir, a rechazar la oferta. Esto tiene sentido pues la proporción de \(0\)’s en la muestra es mucho mayor comparada a la de \(1\)’s. Se puede ver la proporción en la siguiente tabla.
prop.table(table(obs))## obs
## 0 1
## 0.8824805 0.1175195Repetiremos la predicción anterior, pero ahora usando las estimaciones de máxima verosimilitud.
net_mle<-compile(as.grain(fit_mle))## Warning in as.grain.bn.fit(fit_mle): the gRain package does not support
## ordinal networks, disregarding the ordering of the levels.pred_mle <- predict.grain(net_mle, 'y', names(train)[-10],test)
pred_1_mle <- factor(pred_mle$pred$y)
levels(pred_1_mle)<-c(0,1)
cm2 <- confusion.matrix(obs, pred_1_mle)## Warning in Ops.factor(pred, threshold): '>=' not meaningful for factors## Warning in Ops.factor(pred, threshold): '<' not meaningful for factorscm2## obs
## pred 0 1
## 0 26090 2792
## 1 365 731
## attr(,"class")
## [1] "confusion.matrix"Kappa(cm2)## [1] 0.2761494write.net(file = './fit3.net', fit3)Vemos que las estimaciones son muy parecidas. Esto se debe a que las estimaciones de máxima verosimilitud y con iss son muy parecidas también, para términos prácticos son iguales. Con ambos métodos, del conjunto de personas con resultado positivo clasificamos correctamente el 20.75%.
Probamos ahora algunos queries de probabilidad condicional en SAMIAM.
Los queries de probabilidad condicional devueltos por SAMIAM son los siguientes:
P(Y=1)=0.1327 P(Y=1|Job=Unemployed,marital=Single,Education=Primary)=0.1478 P(Job=Admin|marital=married,Edu=primary,Balance=1000)=0.0310
El MAP(Maximum a Posteriori) de \(Y_i\) obtenido con la evidencia observada E={AGE=20, MARITAL=married, EDUCATION=secondary, DEFAULT=no, BALANCE=1000, MONTH=may, CAMPAIGN=1} es la siguiente:
MAP <- read.csv('MAP.csv')
MAP## Variable Resultado Probabilidad
## 1 contact unknown 0.5929
## 2 duration 223-370 0.2169
## 3 housing yes 0.9519
## 4 job blue-collar 0.2513
## 5 loan no 0.8626
## 6 poutcome unknown 0.8841
## 7 previous FALSE 0.8938
## 8 y no 0.9127Lo anterior se realizó con la funcion MAP de SAMIAM. Aunque en estricto sentido el MAP solo devuelve la variable más probable, no su probabilidad, en la tabla anterior viene la probabilidad solo para darnos una idea de cómo se comporta.
La probabilidad condicional nos devuelve una probabilidad, mientras que el MAP nos devuelve el nivel de la variable que es más probable con la evidencia conocida.
Así que la respuesta de cuándo utilizar una y cuándo utilizar otra varía, es decir, la probabilidad condicional la podemos utilizar para obtener la probabilidad de que un evento suceda. El MAP se refiere a el resultado que es más probable que suceda con la evidencia observada. Por ejemplo, si queremos encontrar el grupo de personas que es más probable que vayan a aceptar invertir, utilizaríamos la primera para condicionar las características de las personas que la compañía debería llamar (el target) para maximizar las inversiones, esto nos daría una probabilidad condicional que desearíamos fuera alta. La segunda la utilizaríamos para saber con evidencia dada, cuál es el resultado más probable que vamos a obtener para las variables que aún no conocemos.