head(ml)
id female ses schtyp prog read write math science socst honors awards cid
45 female low public vocation 34 35 41 29 26 not enrolled 0 1
108 male middle public general 34 33 41 36 36 not enrolled 0 1
15 male high public vocation 39 39 44 26 42 not enrolled 0 1
67 male low public vocation 37 37 42 33 32 not enrolled 0 1
153 male middle public vocation 39 31 40 39 51 not enrolled 0 1
51 female high public general 42 36 42 31 39 not enrolled 0 1

table(ml$prog)
academic general vocation
105 45 50
table(ml$ses)
high low middle
58 47 95

# multinom() from nnet package
model <- multinom(prog ~ -1 + ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 180.503647
## final  value 179.981726 
## converged

summary(model)
## Call:
## multinom(formula = prog ~ -1 + ses + write, data = ml)
## 
## Coefficients:
##           seshigh   seslow sesmiddle       write
## general  1.689318 2.852149  2.318866 -0.05792787
## vocation 4.235464 5.218146  5.509562 -0.11360180
## 
## Std. Errors:
##           seshigh   seslow sesmiddle      write
## general  1.226936 1.166438  1.170285 0.02141093
## vocation 1.204684 1.163547  1.168358 0.02221987
## 
## Residual Deviance: 359.9635 
## AIC: 375.9635

Generating predictions

preds <- predict(model, ml, "probs")
head(preds)
academic general vocation
0.1482810 0.3382491 0.5134700
0.1202039 0.1806273 0.6991688
0.4186852 0.2368081 0.3445067
0.1726931 0.3508414 0.4764655
0.1001253 0.1689368 0.7309380
0.3533677 0.2377985 0.4088338

cat_preds_ndx <- apply(preds, 1, which.max)
head(cat_preds_ndx)
## 1 2 3 4 5 6 
## 3 3 1 3 3 3
cat_preds <- colnames(preds)[cat_preds_ndx]
head(cat_preds)
## [1] "vocation" "vocation" "academic" "vocation"
## [5] "vocation" "vocation"

accuracy <- mean(cat_preds == ml$prog)
accuracy
## [1] 0.61
table(ml$prog)
academic general vocation
105 45 50

Changing the pivot class

ml$prog <- relevel(ml$prog, ref = "general")
model2 <- multinom(prog ~ -1 + ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 180.001912
## final  value 179.981726 
## converged

summary(model2)
## Call:
## multinom(formula = prog ~ -1 + ses + write, data = ml)
## 
## Coefficients:
##            seshigh    seslow sesmiddle       write
## academic -1.689292 -2.852130 -2.318847  0.05792773
## vocation  2.546166  2.366045  3.190722 -0.05567427
## 
## Std. Errors:
##          seshigh   seslow sesmiddle      write
## academic 1.22694 1.166441  1.170288 0.02141099
## vocation 1.27830 1.174252  1.199186 0.02333142
## 
## Residual Deviance: 359.9635 
## AIC: 375.9635

preds <- predict(model2, ml, "probs")
head(preds)
general academic vocation
0.3382420 0.1482800 0.5134780
0.1806251 0.1202042 0.6991707
0.2368056 0.4186894 0.3445050
0.3508345 0.1726922 0.4764733
0.1689346 0.1001255 0.7309400
0.2377961 0.3533715 0.4088324