In this section we will focus on creating an confusion matrix in R. Additionally we will perform a significance test, and calculate confidence intervals as well as the kappa coefficient.<\/p>\n
<\/br><\/a><\/p>\n A confusion matrix, also known as error or accuracy matrix, is a specific type of table showing the performance of an algorithm. The name confusion matrix<\/em> comes from the fact that you can quickly see where the algorithm confuses two classes, which would indicate a misclassification.<\/p>\n Several metrics can be derived from the table:<\/p>\n Implementation in R:<\/strong>: Fine, then import your raster<\/span> library in R as usual:<\/p>\n Then import your classification image classification_RF.tif<\/span>, your training shapefile training_data.shp<\/span> and your validation shapefile validation_RF.shp<\/span>. In this example we use the output of the RF classification. However, the workflow is applicable to every classifier!<\/p>\n Our goal is first to generate two factor vectors, which we then compare in our confusion matrix:<\/p>\n For the reference vector, we need to address the validclass <\/em> column of our shapefile (we created this column in this section in QGIS<\/a>). As already mentioned, we need to transform the vectors into factors using as.factor()<\/span> (it becomes clear why in a second):<\/p>\n You may notice some NA<\/span> entries in your reference values. These arise, for example, if you have skipped some validation points during the labeling in QGIS<\/a> because you could not make a clear classification based on the available image data. We simply ignore them later on in the statistics – no problem!<\/p>\n For the predicted vector, we need to extract the classification values at the validation coordinates. Here again, we want a factor object via as.factor()<\/span>:<\/p>\n Once we prepared both factor vectors, we can utilize the table<\/span> function in an elegant way to build a contingency table of the counts at each combination of factor levels. We can additionally name the vectors so that they are also displayed in the table accordingly:<\/p>\n The numbers thus reflect the number of validation pixels. All pixels that have a NA<\/span> value in either reference<\/span> or predicted<\/span> were ignored here. We can now calculate the user’s accuracies UA<\/span>, producer’s accuracies PA<\/span>, and the overall accuracy OA<\/span>:<\/p>\n Actually we already extracted all information needed for a confusion matrix, so let us form a nicely formatted matrix in R:<\/p>\n <\/br><\/a><\/p>\n Furthermore, we can check if the result is purely coincidental, i.e., whether a random classification of the classes could have led to an identical result. We can use a binomial test<\/a> for this. We only need two values for this test: The p-value is much lower than 0.05, so the classification result is highly significant. If the classification were repeated under the same conditions, it can be assumed that the OA is 95% in the range of 81.2% to 90.0%.<\/p>\n <\/br><\/a><\/p>\n The Kappa Coefficient can be used to evaluate the accuracy of a classification. It evaluates how well the classification performs compared to map, in which all values are just randomly assigned. When you have the accuracy matrix as a table \\(m_{i, j}\\) with \\(c\\) different classes, then Kappa is<\/p>\n \\begin{equation} \\label{1}\\tag{1} That looks pretty complicated. Using R, we can write our own function, which calculates \\(\\kappa\\) for us! The calculation also looks much friendlier:<\/p>\n Use the accuracy matrix accmat<\/span> as argument for our new function to calculate the Kappa coefficient:<\/p>\n However, Kappas use has been questioned by many articles and is therefore not recommended (see Pontius Jr and Millones 2011<\/a>).<\/p>\n <\/br><\/br><\/p>\nAccuracy Matrix<\/h1>\n
\n
\nYou need your final classification map, as well as both the training and validation shapefiles created in the course of the Classification in R section<\/a>.
\nGood to go?<\/p>\n\r\n#install.packages(\"raster\")\r\nlibrary(raster)\r\n<\/pre>\n
\r\nsetwd(\"\/media\/sf_exchange\/landsatdata\/\")\r\nimg.classified <- raster("classification_RF.tif")\r\nshp.train <- shapefile("training_data.shp")\r\nshp.valid <- shapefile("validation_RF.shp")\r\n<\/pre>\n\n
\r\nreference <- as.factor(shp.valid$validclass)\r\nreference\r\n## [1] 3 1 4 4 6 4 5 3 6 2 1 5 4 2 \r\n## [15] 2 4 4 5 4 5 1 3 <NA> 2 6 3 3 6 \r\n## [29] 1 4 6 6 5 6 2 3 1 1 3 2 4 2 \r\n## [43] 5 4 2 5 2 4 5 5 3 1 4 1 3 5 \r\n## [57] 2 3 5 4 5 2 5 6 2 1 5 2 5 2 \r\n## [71] 2 6 4 1 5 3 6 4 6 6 6 6 1 5 \r\n## [85] 1 2 5 5 6 4 5 3 5 5 2 2 5 6 \r\n## [99] 2 3 5 6 6 1 1 2 1 2 1 <NA> 3 4 \r\n## [113] 6 1 6 6 5 5 6 2 3 6 6 5 3 5 \r\n## [127] 2 3 6 4 3 6 6 5 4 6 4 2 6 3 \r\n## [141] 2 2 2 3 2 6 6 5 3 3 4 2 6 <NA>\r\n## [155] 2 5 <NA> 3 6 1 3 6 1 2 5 2 2 5 \r\n## [169] 5 <NA> 4 6 5 5 4 4 2 5 1 6 3 5 \r\n## [183] 5 <NA> 5 5 2 2 2 5 3 2 6 2 5 3 \r\n## [197] 6 4 1 3 5 3 3 3 5 4 4 4 4 2 \r\n## [211] 3 5 2 4 1 1 3 2 4 2 2 3 3 5 \r\n## [225] 4 4 6 2 1 1 2 6 3 4 5 3 6 1 \r\n## [239] 6 5 2 2 6 3 6 5 5 1 4 1 2 5 \r\n## [253] 3 4 5 5 3 5 2 3 2 2 6 3 6 5 \r\n## [267] 5 5 3 4 3 <NA> 3 6 5 6 5 5 5 6 \r\n## [281] 2 3 4 5 1 1 3 4 4 1 4 5 6 4 \r\n## [295] 2 3 4 3 5 4 \r\n## Levels: 1 2 3 4 5 6\r\n<\/pre>\n
\r\npredicted <- as.factor(extract(img.classified, shp.valid))\r\npredicted\r\n## [1] 3 1 4 4 6 4 5 3 6 2 1 5 4 2 2 4 4 3 4 5 1 3 4 2 6 3 3 6 1 4 6 6 5 6 2\r\n## [36] 3 1 1 3 3 1 2 5 4 3 4 2 4 1 5 1 1 4 1 3 5 2 1 5 4 5 2 3 6 2 1 5 3 5 2\r\n## [71] 2 6 4 1 5 1 6 4 6 6 6 6 1 1 1 2 4 5 6 4 5 3 5 4 2 2 1 2 2 3 5 6 6 1 1\r\n## [106] 2 1 2 1 1 3 4 6 1 6 6 5 4 6 2 1 6 6 5 1 5 2 3 6 4 2 6 6 3 4 6 4 2 6 3\r\n## [141] 2 2 2 3 2 6 6 5 3 3 4 2 6 1 2 5 1 3 6 1 3 6 1 3 4 2 2 5 4 1 1 6 5 4 4\r\n## [176] 4 2 5 1 6 1 5 5 6 4 5 2 2 2 5 3 2 6 2 5 3 6 4 1 5 5 3 1 3 5 4 4 4 4 2\r\n## [211] 3 5 3 4 1 1 3 3 4 2 2 3 3 1 4 4 6 2 1 1 2 6 3 4 5 3 6 1 6 5 2 3 6 3 6\r\n## [246] 5 1 1 4 1 2 4 3 4 5 5 3 5 2 3 2 2 6 3 6 5 5 3 2 4 3 1 3 6 5 6 5 5 5 6\r\n## [281] 2 3 5 5 1 1 3 4 4 1 4 5 6 4 2 3 5 3 5 4\r\n## Levels: 1 2 3 4 5 6\r\n<\/pre>\n
\r\naccmat <- table("pred" = predicted, "ref" = reference)\r\naccmat\r\n## ref\r\n## pred 1 2 3 4 5 6\r\n## 1 31 0 7 2 5 0\r\n## 2 0 47 2 0 0 1\r\n## 3 0 7 39 0 4 0\r\n## 4 0 0 0 40 9 0\r\n## 5 0 0 1 2 47 0\r\n## 6 0 0 0 0 0 49\r\n<\/pre>\n
\nExcellent! This output already visualize if and where there are misclassifications in our map: all pixels located on the diagonale are correctly classified, all pixels off the diagonal are not.<\/p>\n\r\nUA <- diag(accmat) \/ rowSums(accmat) * 100\r\nUA\r\n## 1 2 3 4 5 6 \r\n## 68.88889 94.00000 78.00000 81.63265 94.00000 100.00000\r\n\r\nPA <- diag(accmat) \/ colSums(accmat) * 100\r\nPA\r\n## 1 2 3 4 5 6 \r\n## 100.00000 87.03704 79.59184 90.90909 72.30769 98.00000\r\n\r\nOA <- sum(diag(accmat)) \/ sum(accmat) * 100\r\nOA\r\n## [1] 86.34812\r\n<\/pre>\n
\r\naccmat.ext <- addmargins(accmat)\r\naccmat.ext <- rbind(accmat.ext, "Users" = c(PA, NA))\r\naccmat.ext <- cbind(accmat.ext, "Producers" = c(UA, NA, OA))\r\ncolnames(accmat.ext) <- c(levels(as.factor(shp.train$classes)), "Sum", "PA")\r\nrownames(accmat.ext) <- c(levels(as.factor(shp.train$classes)), "Sum", "UA")\r\naccmat.ext <- round(accmat.ext, digits = 1)\r\ndimnames(accmat.ext) <- list("Prediction" = colnames(accmat.ext),\r\n "Reference" = rownames(accmat.ext))\r\nclass(accmat.ext) <- "table"\r\naccmat.ext\r\n## Reference\r\n## Prediction baresoil forest grassland urban_hd urban_ld water Sum UA\r\n## baresoil 31.0 0.0 7.0 2.0 5.0 0.0 45.0 68.9\r\n## forest 0.0 47.0 2.0 0.0 0.0 1.0 50.0 94.0\r\n## grassland 0.0 7.0 39.0 0.0 4.0 0.0 50.0 78.0\r\n## urban_hd 0.0 0.0 0.0 40.0 9.0 0.0 49.0 81.6\r\n## urban_ld 0.0 0.0 1.0 2.0 47.0 0.0 50.0 94.0\r\n## water 0.0 0.0 0.0 0.0 0.0 49.0 49.0 100.0\r\n## Sum 31.0 54.0 49.0 44.0 65.0 50.0 293.0 \r\n## PA 100.0 87.0 79.6 90.9 72.3 98.0 86.3\r\n<\/pre>\nSignificance Test<\/h1>\n
\nx = <\/span> total number of correctly classified validation points, and
\nn = <\/span> the total number of validation points in our confusion matrix:<\/p>\n\r\nsign <- binom.test(x = sum(diag(accmat)),\r\n n = sum(accmat),\r\n alternative = c("two.sided"),\r\n conf.level = 0.95\r\n )\r\n\r\npvalue <- sign$p.value\r\npvalue\r\n## [1] 5.30128e-39\r\n\r\nCI95 <- sign$conf.int[1:2]\r\nCI95\r\n## [1] 0.8187710 0.9006459\r\n<\/pre>\nKappa<\/h1>\n
\nThe Kappa coefficient can range from -1 to 1.
\nA value of 0 indicates that the classification is as good as random values.
\nA value below 0 indicates the classification is significantly worse than random.
\nA value greater than 0 indicates that the classification is significantly better than random.<\/p>\n
\n \\kappa = \\frac{N_{o} – N_{e}}{N – N_{e}} , \\text{with} \\\\
\n N = \\sum_{i,j = 1}^{c} m_{i, j}\\\\
\n N_{o} = \\sum_{i=1}^c{m_{i, i}}\\\\
\n N_{e} = \\frac{1}{N}\\cdot\\sum_{l=1}^c\\left(\\sum_{j=1}^c{m_{l, j}} \\cdot \\sum_{i=1}^c{m_{i, l}}\\right)
\n\\end{equation}<\/p>\n\r\nkappa <- function(m) {\r\n N <- sum(m)\r\n No <- sum(diag(m))\r\n Ne <- 1 \/ N * sum(colSums(m) * rowSums(m))\r\n return( (No - Ne) \/ (N - Ne) )\r\n}\r\n<\/pre>\n\r\nkappacoefficient <- kappa(accmat)\r\nkappacoefficient\r\n## [1] 0.8359646\r\n<\/pre>\n
\n
\n