The expected Confusion Matrix

The model used in this simulation is a Bayesian network. Like many kinds of classification models its output is the probability that the subject is in each category. So the columns Reading.Novice, Reading.Intermediate, and Reading.Advanced represent the probability based on the observed evidence that the simulee’s reading ability is in those three categories. This is called the marginal distribution as it is one side of a big multi-way table with all of the varaibles. The mode.Reading, which we used as the score before, just tells which category has the highest probability. 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",
                     "mode.Reading")],
     caption="Reading data from first five simulees.",
     digits=3) |>
  kable_classic()

Look at the first row of the data. The first simulee was randomly assigned a reading ability of Intermediate and then the data were generated based on this estimate. The classification of the first simulee based on the 16 item form is a 56.5% chance of Novice and a 43.5% chance of Intermediate. So, the “modal” (highest probability) score is Novice.

When building the confusion matrix, each observation add a count of one to a cell based on the true and estimated value. So the first simulee contributes 1 towards the cell in the second row and first column in the table. In building the expected table, this value is split across the row according to the probabilities, so the first simulee adds .565 to the first cell in the second row, .435 to the second cell and nothing to the third cell.

The table below shows the calculations for the first five simulees. The first row is the sum of the Reading.Novice, Reading.Intermediate and Reading.Advanced variable for simulees whose true Reading value is Novice, the second for true Intermediate, and third row for true Advanced.

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 \(i\)th simulee, and let \(p_{im} = P(\hat{X_i}=m|e)\) be the estimated probability that \(X_i=m\) for the \(i\)th 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")

Based on this definition of the expected confusion matrix, here are the values for the 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")

Raw Agreement

The sum of the diagonal of the confusion matrix, \(\sum_k a_{kk}\), gives a count of how many cases are exact agreements (in this case between the simulation and estimation). Let \(N=\sum_k\sum_m a_{km}\); then the agreement rate is \(\sum_k a_{kk}/N\). For the reading data using the MAP scores, this is 831 out of 1000, so over 80% agreement. The function accuracy() calculates the agreement rate.

acc.tab <- data.frame(MAP=sapply(cm,accuracy),
                      EAP=sapply(em,accuracy))

Raw agreement can be easy to achieve if there is not much variability in the population. For example, if 80% of the target population was intermediate, a classifier that simply classified each respondent as intermediate would achieve 80% accuracy, and one that guessed intermediate randomly 80% of the time would achieve at least 64% accuracy. For that reason, two adjusted agreement rates, lambda and kappa, are often used.

Goodman-Kruskal Lambda

There is a problem with the agreement rate: depending on the proportions of cases in the population, it could be very easy to get a good agreement rate. Suppose the question is whether or not a randomly chosen person has traveled beyond the Earth’s atmosphere. As there are only a few dozen astronauts and cosmonauts, just guessing that everybody has not been to space is likely to get 100% accuracy.

@Goodman1954 suggested a correction. Pick the row which has the biggest proportion. For the Reading data, this is Intermediate, there are 480 intermediate readers in the sample, so simply classifying everybody as an intermediate would correctly classify 48% of the sample.

Let \(p_{ij}\) be the proportion of cases where the “true” category was \(i\) and the estimated category was \(j\). Let \(p_{i+} = \sum_{j} p_{ij}\) be the row sums and \(\sum_{i} p_{ii}\) is the agreement. Lambda discounts the agreement by the biggest of the row sums.

\[ \lambda = \frac{\sum p_{ii} - \max_{i} p_{i+}}{1 - \max_{i} p_{i+}}\] 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))

This is pretty good. There is a 67.5% improvement in classification from the one-size-fits-all model of classifying everybody at Intermediate.

Suppose the reading test was being used as a placement test.
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, Intermediate. On the other hand, by giving the reading test, the classification rate improves by 67%; this looks like it might be useful.

Fleis-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 \(\sum_{k} a_{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 \ge w_{km} \ge 0\) is the desirability of classifying a subject who is \(k\) as \(m\). In this case, weighted agreement is \(\sum_k\sum_m w_{km} a_{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)^2/(K-1)^2\)

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")
)