The expected Confusion Matrix

The Bayes net scoring engine (like many scoring models) expresses its uncertainty about the abilities of the students by estimating the probability that the student is in each category. These are often called the marginal probabilities because they are one margin of the joint distribution over all four variables. The table below shows the “true” (simulated) reading ability and estimated marginal probabilities for the first five simulees.

kbl(language16[1:5,c("Reading","Reading.Novice",
                     "Reading.Intermediate","Reading.Advanced")],
     caption="Reading data from first five simulees.",
     digits=3) |>
  kable_classic()

The expected value of the confusion matrix, $\overline{\mathbf{A}}$ is calculated as follows: let X_i be the value of the first (simulated truth) variable for the ith simulee, and let $p_{im} = P(\hat{X_i}=m|e)$ be the estimated probability that X_i = m for the ith simulee is m. Then $\overline{a_{km}} = \sum_{i: X_i=k} p_{im}$.

In the running example, the first row (novice) is the sum of all rows for which the true scores is novice. The second row is the sum of the intermediate rows and the third row advanced.

The function expTable() does this work. Note that it expects the marginal probability to be in a number of columns marked var.state, where var is the name of the variable and state is the name of the state. If the data uses a different naming convention, this can be expressed with the argument pvecregex which is a regular expression with special symbols <var> to be substituted with the variable name and <state> to be substituted with the state name.

The table below shows the expected matrix from the first five rows of the reading data.

reading5 <- expTable(language16[1:5,],"Reading","Reading")

What follows are the expected confusion matrixes for all four proficiency variables.

em <- list()
em$Reading <- expTable(language16,"Reading","Reading")

em$Writing <- expTable(language16,"Writing","Writing")

em$Speaking <- expTable(language16,"Speaking","Speaking")

em$Listening <- expTable(language16,"Listening","Listening")

Goodman and Kruskal Lambda

If the test was not available, the only strategy would be a one-size-fits-all one, assuming that all subjects are at the same level of the variable. The best strategy is to treat all subjects as if they are in the modal (most likely) category. The marginal distribution for the variable is found from the row sums of A (or $\overline{\mathbf{A}}$), a_k+ = ∑_ma_km, and the best that can be done with the one-size-fits-all strategy is max_ka_k+.

Goodman and Kruskal (1952) suggest that the agreement rate be adjusted by subtracting out the best that can be done with mapping everybody to a single value. So they propose

$$\lambda = \frac{\sum_{k} a_{kk} - \max_{k} a_{k+}}{N - \max_{k} a_{k+}} \ .$$

This ranges from -1 to 1 with 0 representing doing no better than just guessing the most probable character and 1 indicating perfect agreement. (Negative numbers are doing worse than just guessing the mode.)

The function gkLambda() will do this calculation. Here are the values for the language test using both the MAP and expected agreements.

lambda.tab <- data.frame(MAP=sapply(cm,gkLambda),
                         EAP=sapply(em,gkLambda))

Flies-Cohen Kappa

Jacob Cohen (Flies, Levin & Paek, 2003) took a different approach which treats the two assignments of cases to categories more symmetrically. The idea is that these are two raters and the goal is to judge the extent of the agreement. Baseline here is to imagine two raters one of which assign categories randomly with probabilities a_k+/N and a_+k/N. Then the expected number of random agreements is ∑_ka_k+a_+k/N. So the agreement measure adjusted for random agreement is:

$$ \kappa = \frac{\sum_{k} a_{kk} - \sum_{k}a_{k+}a_{+k}/N}{N-\sum_{k}a_{k+}a_{+k}/N} \ .$$

Again this runs from -1 to 1. The function fcKappa calculates Cohen’s kappa.

kappa.tab <- data.frame(MAP=sapply(cm,fcKappa),
                         EAP=sapply(em,fcKappa))

Weighted versions

All three of kappa, lambda and raw agreement all assume that any missclassification is equally bad. Fleis suggest adding weights, were 1 ≥ w_km ≥ 0 is the desirability of classifying a subject who is k as m. In this case, weighted agreement is ∑_k∑_mw_kma_km/N. The weighted versions of kappa and lambda are given by:

$$ \lambda = \frac{\sum\sum w_{km} a_{km} - \max_k \sum_m w_{km} a_{km}}{N - \max_k \sum_m w_{km} a_{km}} \ ;$$

$$ \kappa = \frac{\sum\sum w_{km} a_{km} - \sum_k \sum_m w_{km} a_{k+}a_{+m}/N}{N - \sum_k \sum_m w_{km} a_{k+}a_{+k}/N} \ .$$

There are three commonly uses cases:

None

w_km = 1 if k = m, 0 otherwise.

Linear

w_km = 1 − |k − m|/(K − 1)

Quadratic

w_km = 1 − (k − m)²/(K − 1)²

Both linear and quadratic weight have increasing penalties for the number of categories of difference. So off-by-one has a lower penalty than off-by-two.

The accuracy, gkLambda and fcKappa have both a weights argument where no weights (“None”, default), “Linear” or “Quadratic” weights can be selected, and a w argument where a custom weight matrix can be entered.

wacc.tab <- data.frame(
  None = sapply(cm,accuracy,weights="None"),
  Linear = sapply(cm,accuracy,weights="Linear"),
  Quadratic = sapply(cm,accuracy,weights="Quadratic"))
wlambda.tab <- data.frame(
  None = sapply(cm,gkLambda,weights="None"),
  Linear = sapply(cm,gkLambda,weights="Linear"),
  Quadratic = sapply(cm,gkLambda,weights="Quadratic"))
wkappa.tab <- data.frame(
  None = sapply(cm,fcKappa,weights="None"),
  Linear = sapply(cm,fcKappa,weights="Linear"),
  Quadratic = sapply(cm,fcKappa,weights="Quadratic")
)