La regresión logística es una técnica que se utiliza cuando la variable dependiente es categórica (o nominal). Para la regresión logística binaria, el número de variables dependientes es dos, mientras que el número de variables dependientes para la regresión logística multinomial es más de dos.
Ejemplos: Los consumidores toman la decisión de comprar o no comprar, un producto puede pasar o no el control de calidad, existen buenos o malos riesgos crediticios y el empleado puede ser ascendido o no.
En la regresión logística, una transformación logística de las probabilidades (conocida como logit) sirve como la variable \[ \log(odds) = logit(P) = \ln \left(\frac{P}{1-P} \right) = a + b_{1} x_{1} + b_{2} x_{2} + b_{3} x_{3} + \ldots \] o \[ p = \frac{\exp \left(a + b_{1} X_{1} + b_{2} X_{2} + b_{3} X_{3} + \ldots \right)} {1+ \exp \left(a + b_{1} X_{1} + b_{2} X_{2} + b_{3} X_{3} + \ldots \right)} \]
Dónde:
\(p =\) la probabilidad de que un caso esté en una categoría particular,
\(\exp=\) el exponencial (aprox. \(2.72\) ).
\(\mathrm{a} =\) la constante de la ecuación y
\(b =\) el coeficiente del predictor o variables independientes.
Logits o Log Odds.
El valor de las probabilidades puede variar de 0 a infinito y le indica cuánto más probable es que una observación sea un miembro del grupo objetivo en lugar de un miembro del otro grupo.
Si la probabilidad es 0,80, las probabilidades son 4 a 1 o 0,80 / 0,20; si la probabilidad es 0.25, las probabilidades son .33 (.25 / .75).
La razón de probabilidades (OR), estima el cambio en las probabilidades de pertenencia al grupo objetivo para un aumento de una unidad en el predictor. Se calcula utilizando el coeficiente de regresión del predictor como exponente o exp.
En el ejemplo anterior, supongamos que estábamos prediciendo el éxito contable mediante un predictor de competencia matemática que b = 2,69. Por tanto, la razón de posibilidades es exp (2,69) o 14,73. Por lo tanto, las probabilidades de aprobar son 14,73 veces mayores para un estudiante, por ejemplo, que obtuvo una puntuación previa a la prueba de 5 que para un estudiante cuya puntuación previa a la prueba fue 4.
En regresión logística, las hipótesis son de interés:
La hipótesis nula, que es cuando todos los coeficientes de la ecuación de regresión toman el valor cero, y
La hipótesis alternativa de que el modelo actualmente en consideración es preciso y difiere significativamente del nulo de cero, es decir, da significativamente mejor que el nivel de predicción aleatoria o aleatoria de la hipótesis nula.
Evaluación de hipótesis
Luego calculamos la probabilidad de observar los datos que realmente observamos bajo cada una de estas hipótesis. El resultado suele ser un número muy pequeño y, para facilitar su manejo, se utiliza el logaritmo natural, lo que produce una probabilidad logarítmica (LL). Las probabilidades son siempre menores que uno, por lo que los LL son siempre negativos. La verosimilitud logarítmica es la base para las pruebas de un modelo logístico.
La prueba de razón de verosimilitud se basa en la razón -2LL. Es una prueba de la significancia de la diferencia entre la razón de verosimilitud (-2LL) para el modelo del investigador con predictores (llamado modelo chi cuadrado) menos la razón de verosimilitud para el modelo de línea de base con solo una constante en él.
La significancia en el nivel .05 o menor significa que el modelo del investigador con los predictores es significativamente diferente del modelo con la constante solamente (todos los coeficientes “b” son cero). Mide la mejora de ajuste que realizan las variables explicativas frente al modelo nulo.
El chi cuadrado se utiliza para evaluar la importancia de esta relación (consulte Información de ajuste del modelo en la salida de SPSS).
\(H_0\): No hay diferencia entre el modelo nulo y el modelo final.
\(H_1\): Hay diferencia entre el modelo nulo y el modelo final
Simplemente ejecute una “regresión lineal” después de asumir la variable dependiente categórica como variable continua.
Si el VIF (factor de inflación de la varianza) más grande es mayor que 10, entonces hay un motivo de preocupación (Bowerman y O’Connell, 1990).
La tolerancia por debajo de 0,1 indica un problema grave.
La tolerancia por debajo de 0,2 indica un problema potencial (Menard, 1995).
Si el índice de condición es mayor que 15, se asume la multicolinealidad.
Regresión logística multinomial para predecir la pertenencia a más de dos categorías. (Básicamente) funciona de la misma manera que la regresión logística binaria. El análisis divide la variable de resultado en una serie de comparaciones entre dos categorías.
Por ejemplo, si tiene tres categorías de resultados (A, B y C), el análisis consistirá en dos comparaciones que elija:
Compare todo con su primera categoría (por ejemplo, A frente a B y A frente a C),
O su última categoría (por ejemplo, A frente a C y B frente a C),
O una categoría personalizada (por ejemplo, B frente a A y B frente a C).
Las partes importantes del análisis y la salida son muy parecidas a las que acabamos de ver para la regresión logística binaria.
El conjunto de datos (hsbdemo.sav) contiene variables sobre 200 estudiantes. La variable de resultado es prog, tipo de programa (1 = general, 2 = académico y 3 = vocacional). Las variables predictoras son ses, estatus socioeconómico (1 = bajo, 2 = medio y 3 = alto), matemáticas, puntaje en matemáticas y ciencias, puntaje en ciencias: ambas son variables continuas.
(Pregunta de investigación): Cuando los estudiantes de secundaria eligen el programa (programas generales, vocacionales y académicos), ¿cómo afectan sus calificaciones en matemáticas y ciencias y su estatus económico social (SES) su decisión?
# Starting our example by import the data into R
library(haven)
hsbdemo <- read_sav("E:/hsbdemo.sav")
hsb<-hsbdemo # Get a new copy of data
summary(hsb)## id female ses schtyp
## Min. : 1.00 Min. :0.000 Min. :1.000 Min. :1.00
## 1st Qu.: 50.75 1st Qu.:0.000 1st Qu.:2.000 1st Qu.:1.00
## Median :100.50 Median :1.000 Median :2.000 Median :1.00
## Mean :100.50 Mean :0.545 Mean :2.055 Mean :1.16
## 3rd Qu.:150.25 3rd Qu.:1.000 3rd Qu.:3.000 3rd Qu.:1.00
## Max. :200.00 Max. :1.000 Max. :3.000 Max. :2.00
## prog read write math
## Min. :1.000 Min. :28.00 Min. :31.00 Min. :33.00
## 1st Qu.:2.000 1st Qu.:44.00 1st Qu.:45.75 1st Qu.:45.00
## Median :2.000 Median :50.00 Median :54.00 Median :52.00
## Mean :2.025 Mean :52.23 Mean :52.77 Mean :52.65
## 3rd Qu.:2.250 3rd Qu.:60.00 3rd Qu.:60.00 3rd Qu.:59.00
## Max. :3.000 Max. :76.00 Max. :67.00 Max. :75.00
## science socst honors awards cid
## Min. :26.00 Min. :26.00 Min. :0.000 Min. :0.00 Min. : 1.00
## 1st Qu.:44.00 1st Qu.:46.00 1st Qu.:0.000 1st Qu.:0.00 1st Qu.: 5.00
## Median :53.00 Median :52.00 Median :0.000 Median :1.00 Median :10.50
## Mean :51.85 Mean :52.41 Mean :0.265 Mean :1.67 Mean :10.43
## 3rd Qu.:58.00 3rd Qu.:61.00 3rd Qu.:1.000 3rd Qu.:2.00 3rd Qu.:15.00
## Max. :74.00 Max. :71.00 Max. :1.000 Max. :7.00 Max. :20.00# Load the jmv package for frequency table
library(jmv)Ahora hagamos el análisis descriptivo
# Use the descritptives function to get the descritptive data
descriptives(hsb, vars = vars(ses, prog, math, science), freq = TRUE)##
## DESCRIPTIVES
##
## Descriptives
## -----------------------------------------------------------
## ses prog math science
## -----------------------------------------------------------
## N 200 200 200 200
## Missing 0 0 0 0
## Mean 2.055000 2.025000 52.64500 51.85000
## Median 2.000000 2.000000 52.00000 53.00000
## Minimum 1.000000 1.000000 33.00000 26.00000
## Maximum 3.000000 3.000000 75.00000 74.00000
## -----------------------------------------------------------
##
##
## FREQUENCIES
##
## Frequencies of ses
## --------------------------------------------------
## Levels Counts % of Total Cumulative %
## --------------------------------------------------
## 1 47 23.50000 23.50000
## 2 95 47.50000 71.00000
## 3 58 29.00000 100.00000
## --------------------------------------------------
##
##
## Frequencies of prog
## --------------------------------------------------
## Levels Counts % of Total Cumulative %
## --------------------------------------------------
## 1 45 22.50000 22.50000
## 2 105 52.50000 75.00000
## 3 50 25.00000 100.00000
## --------------------------------------------------# To see the crosstable, we need CrossTable function from gmodels package
library(gmodels)
# Build a crosstable between admit and rank
CrossTable(hsb$ses, hsb$prog)##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 200
##
##
## | hsb$prog
## hsb$ses | 1 | 2 | 3 | Row Total |
## -------------|-----------|-----------|-----------|-----------|
## 1 | 16 | 19 | 12 | 47 |
## | 2.783 | 1.305 | 0.005 | |
## | 0.340 | 0.404 | 0.255 | 0.235 |
## | 0.356 | 0.181 | 0.240 | |
## | 0.080 | 0.095 | 0.060 | |
## -------------|-----------|-----------|-----------|-----------|
## 2 | 20 | 44 | 31 | 95 |
## | 0.088 | 0.692 | 2.213 | |
## | 0.211 | 0.463 | 0.326 | 0.475 |
## | 0.444 | 0.419 | 0.620 | |
## | 0.100 | 0.220 | 0.155 | |
## -------------|-----------|-----------|-----------|-----------|
## 3 | 9 | 42 | 7 | 58 |
## | 1.257 | 4.381 | 3.879 | |
## | 0.155 | 0.724 | 0.121 | 0.290 |
## | 0.200 | 0.400 | 0.140 | |
## | 0.045 | 0.210 | 0.035 | |
## -------------|-----------|-----------|-----------|-----------|
## Column Total | 45 | 105 | 50 | 200 |
## | 0.225 | 0.525 | 0.250 | |
## -------------|-----------|-----------|-----------|-----------|
##
## A continuación, usamos la función multinom del paquete nnet para estimar un modelo de regresión logística multinomial. Hay otras funciones en otros paquetes de R capaces de regresión multinomial. Elegimos la función multinom porque no requiere que los datos se modifiquen (como lo hace el paquete mlogit) y para reflejar el código de ejemplo que se encuentra en los modelos de regresión logística de Hilbe.
Primero, debemos elegir el nivel de nuestro resultado que deseamos usar como nuestra línea de base y especificarlo en la función de relevamiento. Luego, ejecutamos nuestro modelo usando multinom. El paquete multinom no incluye el cálculo del valor p para los coeficientes de regresión, por lo que calculamos los valores p utilizando las pruebas de Wald (aquí, las pruebas z).
# Load the multinom package
library(nnet)
# Since we are going to use Academic as the reference group, we need relevel the group.
hsb$prog2 <- relevel(as.factor(hsb$prog), ref = 2)
hsb$ses <- as.factor(hsb$ses)
levels(hsb$prog2)## [1] "2" "1" "3"# Give the names to each level
levels(hsb$prog2) <- c("academic","general","vocational")
# Run a "only intercept" model
OIM <- multinom(prog2 ~ 1, data = hsb)## # weights: 6 (2 variable)
## initial value 219.722458
## final value 204.096674
## convergedsummary(OIM)## Call:
## multinom(formula = prog2 ~ 1, data = hsb)
##
## Coefficients:
## (Intercept)
## general -0.8472980
## vocational -0.7419374
##
## Std. Errors:
## (Intercept)
## general 0.1781742
## vocational 0.1718249
##
## Residual Deviance: 408.1933
## AIC: 412.1933# Run a multinomial model
multi_mo <- multinom(prog2 ~ ses + math + science + math*science, data = hsb,model=TRUE)## # weights: 21 (12 variable)
## initial value 219.722458
## iter 10 value 173.831002
## iter 20 value 167.382760
## final value 166.951813
## convergedsummary(multi_mo)## Call:
## multinom(formula = prog2 ~ ses + math + science + math * science,
## data = hsb, model = TRUE)
##
## Coefficients:
## (Intercept) ses2 ses3 math science
## general 5.897618 -0.4081497 -1.1254491 -0.1852220 0.01323626
## vocational 22.728283 0.8402168 -0.5605656 -0.5036705 -0.28297703
## math:science
## general 0.001025283
## vocational 0.006185571
##
## Std. Errors:
## (Intercept) ses2 ses3 math science math:science
## general 0.002304064 0.2613732 0.2134308 0.02694593 0.02953364 0.0004761369
## vocational 0.003856861 0.2959741 0.1984775 0.02681947 0.03142872 0.0004760567
##
## Residual Deviance: 333.9036
## AIC: 357.9036# Check the Z-score for the model (wald Z)
z <- summary(multi_mo)$coefficients/summary(multi_mo)$standard.errors
z## (Intercept) ses2 ses3 math science math:science
## general 2559.659 -1.561559 -5.273132 -6.873837 0.4481759 2.153336
## vocational 5892.948 2.838819 -2.824328 -18.780028 -9.0037711 12.993348# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p## (Intercept) ses2 ses3 math science
## general 0 0.118391943 1.341147e-07 6.249667e-12 0.6540263
## vocational 0 0.004528089 4.737981e-03 0.000000e+00 0.0000000
## math:science
## general 0.03129229
## vocational 0.00000000Estos son los coeficientes logit relativos a la categoría de referencia. Por ejemplo, en “matemáticas”, el -0,185 sugiere que para un aumento de una unidad en la puntuación de “ciencia”, el coeficiente logit para “bajo” en relación con “medio” bajará en esa cantidad, -0,185.
# La función anova está conflictiva con la función anova de JMV, por lo que necesitamos desactivar de la función JMV antes de usar la función anova.
# detach ("paquete: jmv", unload = TRUE)
# Compare nuestro modelo de prueba con el modelo "Solo interceptar"
# anova (OIM, multi_mo)detach("package:jmv", unload=TRUE)
# Compare the our test model with the "Only intercept" model
anova(OIM,multi_mo)Interpretación de la información de ajuste del modelo
La probabilidad logarítmica es una medida de cuánta variabilidad inexplicable hay en los datos. Por lo tanto, la diferencia o cambio en la probabilidad logarítmica indica cuánta nueva varianza ha sido explicada por el modelo.
La prueba de chi-cuadrado prueba la disminución en la varianza inexplicable desde el modelo de línea de base (408.1933) al modelo final (333.9036), que es una diferencia de 408.1933 - 333.9036 = 74.29. Este cambio es significativo, lo que significa que nuestro modelo final explica una cantidad significativa de la variabilidad original.
La razón de verosimilitud chi-cuadrado de 74,29 con un valor de p <0,001 nos dice que nuestro modelo en su conjunto se ajusta significativamente mejor que un modelo vacío o nulo (es decir, un modelo sin predictores).
# Check the predicted probability for each program
head(multi_mo$fitted.values,30)## academic general vocational
## 1 0.18801940 0.17122451 0.6407561
## 2 0.12019189 0.10715542 0.7726527
## 3 0.52212681 0.08123771 0.3966355
## 4 0.23683979 0.23125435 0.5319059
## 5 0.10132130 0.12329032 0.7753884
## 6 0.38079544 0.10780793 0.5113966
## 7 0.32321815 0.16454057 0.5122413
## 8 0.09033932 0.08381233 0.8258484
## 9 0.02336687 0.09050704 0.8861261
## 10 0.32321815 0.16454057 0.5122413
## 11 0.16304678 0.29839918 0.5385540
## 12 0.22842326 0.27539161 0.4961851
## 13 0.32747927 0.28141483 0.3911059
## 14 0.05717483 0.12540921 0.8174160
## 15 0.15003741 0.52649953 0.3234631
## 16 0.12638004 0.20495962 0.6686603
## 17 0.25269654 0.39841475 0.3488887
## 18 0.05771613 0.08029142 0.8619924
## 19 0.27404420 0.16131436 0.5646414
## 20 0.25197679 0.20299587 0.5450273
## 21 0.15870561 0.17100945 0.6702849
## 22 0.27404420 0.16131436 0.5646414
## 23 0.16304678 0.29839918 0.5385540
## 24 0.16340250 0.31572103 0.5208765
## 25 0.26080538 0.36665885 0.3725358
## 26 0.74715288 0.12007209 0.1327750
## 27 0.33135572 0.26860688 0.4000374
## 28 0.32321815 0.16454057 0.5122413
## 29 0.40025162 0.32553566 0.2742127
## 30 0.19518342 0.30892507 0.4958915# We can get the predicted result by use predict function
head(predict(multi_mo),30)## [1] vocational vocational academic vocational vocational vocational
## [7] vocational vocational vocational vocational vocational vocational
## [13] vocational vocational general vocational general vocational
## [19] vocational vocational vocational vocational vocational vocational
## [25] vocational academic vocational vocational academic vocational
## Levels: academic general vocational# Test the goodness of fit
chisq.test(hsb$prog2,predict(multi_mo))## Warning in chisq.test(hsb$prog2, predict(multi_mo)): Chi-squared approximation
## may be incorrect##
## Pearson's Chi-squared test
##
## data: hsb$prog2 and predict(multi_mo)
## X-squared = 47.841, df = 4, p-value = 1.019e-09# Saque el signo "#" para ejecutar el código
# Cargue el paquete DescTools para calcular el cuadrado R
# library("DescTools")
# Calcule el cuadrado R
# PseudoR2 (multi_mo, which = c ("CoxSnell", "Nagelkerke", "McFadden"))library("DescTools")## Registered S3 method overwritten by 'DescTools':
## method from
## reorder.factor gdata# Calcule el cuadrado R
PseudoR2(multi_mo, which = c ("CoxSnell", "Nagelkerke", "McFadden"))## CoxSnell Nagelkerke McFadden
## 0.3102655 0.3565873 0.1819964Interpretación de la R-Square:
Estos son tres valores pseudo R cuadrado. La regresión logística no tiene un equivalente al R cuadrado que se encuentra en la regresión MCO; sin embargo, muchas personas han intentado encontrar uno. Estas estadísticas no significan exactamente lo que significa R cuadrado en la regresión MCO (la proporción de varianza de la variable de respuesta explicada por los predictores), sugerimos interpretarlos con mucha precaución.
El R-Square de Cox y Snell imita múltiples R-Square en función de la “probabilidad”, pero su máximo puede ser (y suele ser) inferior a 1,0, lo que dificulta su interpretación. Aquí está indicando que existe una relación de \(31\%\) entre la variable dependiente y las independientes. O indica que el \(31\%\) de la variación en la variable dependiente se explica por el modelo logístico.
La modificación de Nagelkerke que varía de 0 a 1 es una medida más confiable de la relación. El R2 de Nagelkerke normalmente será más alto que la medida de Cox y Snell. En nuestro caso es 0,357, lo que indica una relación del \(35,7\%\) entre los predictores y la predicción.
McFadden = {LL (nulo) - LL (completo)} / LL (nulo). En nuestro caso es 0,182, lo que indica una relación del \(18,2\%\) entre los predictores y la predicción.
# Use the lmtest package to run Likelihood Ratio Tests
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlrtest(multi_mo, "ses") # Chi-Square=12.922,p=0.01166*## # weights: 15 (8 variable)
## initial value 219.722458
## iter 10 value 176.645309
## iter 20 value 173.670728
## iter 30 value 173.430200
## final value 173.412997
## convergedlrtest(multi_mo, "math") # Chi-Square=10.613,p=0.004959*## # weights: 18 (10 variable)
## initial value 219.722458
## iter 10 value 172.387155
## final value 172.258318
## convergedlrtest(multi_mo, "science") # Chi-Square=5.3874,p=0.06763## # weights: 18 (10 variable)
## initial value 219.722458
## iter 10 value 169.972735
## final value 169.645522
## convergedlrtest(multi_mo, "math:science") # Chi-Square=5.249,p=0.072## # weights: 18 (10 variable)
## initial value 219.722458
## iter 10 value 170.088834
## final value 169.576379
## convergedInterpretación de las pruebas de razón de verosimilitud
\(X^2(2) = 10.613\), \(p = .005\) para SES.
# Let's check our model again
summary(multi_mo)## Call:
## multinom(formula = prog2 ~ ses + math + science + math * science,
## data = hsb, model = TRUE)
##
## Coefficients:
## (Intercept) ses2 ses3 math science
## general 5.897618 -0.4081497 -1.1254491 -0.1852220 0.01323626
## vocational 22.728283 0.8402168 -0.5605656 -0.5036705 -0.28297703
## math:science
## general 0.001025283
## vocational 0.006185571
##
## Std. Errors:
## (Intercept) ses2 ses3 math science math:science
## general 0.002304064 0.2613732 0.2134308 0.02694593 0.02953364 0.0004761369
## vocational 0.003856861 0.2959741 0.1984775 0.02681947 0.03142872 0.0004760567
##
## Residual Deviance: 333.9036
## AIC: 357.9036# Check the Wald Z again
z <- summary(multi_mo)$coefficients/summary(multi_mo)$standard.errors
z## (Intercept) ses2 ses3 math science math:science
## general 2559.659 -1.561559 -5.273132 -6.873837 0.4481759 2.153336
## vocational 5892.948 2.838819 -2.824328 -18.780028 -9.0037711 12.993348# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1))*2
p## (Intercept) ses2 ses3 math science
## general 0 0.118391943 1.341147e-07 6.249667e-12 0.6540263
## vocational 0 0.004528089 4.737981e-03 0.000000e+00 0.0000000
## math:science
## general 0.03129229
## vocational 0.00000000Tenga en cuenta que la tabla está dividida en dos filas. Esto se debe a que estos parámetros comparan pares de categorías de resultados.
Especificamos la segunda categoría (2 = académica) como nuestra categoría de referencia; por lo tanto, la primera fila de la tabla denominada General compara esta categoría con la categoría “Académica”. la segunda fila de la tabla denominada Vocacional también compara esta categoría con la categoría “Académica”.
Debido a que solo estamos comparando dos categorías, la interpretación es la misma que para la regresión logística binaria:
# extract the coefficients from the model and exponentiate
exp(coef(multi_mo))## (Intercept) ses2 ses3 math science math:science
## general 3.641690e+02 0.6648794 0.3245067 0.8309198 1.0133243 1.001026
## vocational 7.426219e+09 2.3168692 0.5708861 0.6043085 0.7535371 1.006205Las probabilidades logarítmicas relativas de estar en un programa general frente a un programa académico disminuirán en 1.125 si se pasa del nivel más alto de SES (SES = 3) al nivel más bajo de SES (SES = 1), b = -1.125, Wald \(\chi^2( 1) = -5.27\), p< 0.001.
\(Exp(-1.1254491) = 0.3245067\) significa que cuando los estudiantes pasan del nivel más alto de SES (SES = 3) al nivel más bajo de SES (1 = SES), la razón de probabilidades es 0.325 veces más alta y, por lo tanto, los estudiantes con el nivel más bajo de Los SES tienden a elegir el programa general contra el programa académico más que los estudiantes con el nivel más alto de SES.
Las probabilidades logarítmicas relativas de estar en un programa vocacional versus un programa académico disminuirán en 0.56 si se pasa del nivel más alto de SES (SES = 3) al nivel más bajo de SES (SES = 1), \(b = -0.56\), Wald \(\chi^2= -2.82,\quad p <0.01\).
\(Exp(-0.56) = 0.57\) significa que cuando los estudiantes pasan del nivel más alto de SES (SES = 3) al nivel más bajo de SES (SES = 1), la razón de probabilidades es 0.57 veces más alta y, por lo tanto, los estudiantes con el nivel más bajo de Los SES tienden a elegir un programa vocacional contra un programa académico más que los estudiantes con el nivel más alto de SES.
# Load the summarytools package to use the classification function
library(summarytools)## Registered S3 method overwritten by 'pryr':
## method from
## print.bytes Rcpp## For best results, restart R session and update pander using devtools:: or remotes::install_github('rapporter/pander')# Build a classification table by using the ctable function
ctable <- table(hsb$prog2,predict(multi_mo))
ctable##
## academic general vocational
## academic 88 4 13
## general 24 12 9
## vocational 21 4 25