User Tools

Site Tools


r:r-tutorials:calculatekappa

Calculate Kappa and Disagreement Indices

Concept of the Kappa Index

Cohen's Kappa Index {[:aval:Cohen1960]} vergleicht den Anteil übereinstimmender Pixel zweier Raster mit der erwarteten Übereinstimmung bei zufälliger Pixelanordnung. Die Indexwerte liegen zwischen -1 und 1 wobei ein Indexwert von 0 einem Ergebnis entspricht, das bei zufälliger Pixelanordnung erwartet werden kann. Ein Indexwert von 1 bedeutet genaue Übereinstimmung. Auch wenn die Interpretation der Werte sehr von der Fragestellung abhängt bieten {[:aval:landis1977]} und {[:aval:monserud1992]} eine Orientierung zur Interpretation:

Kappa Übereinstimmung
<0.05keine
0.05-0.2sehr gering
0.2-0.4gering
0.4-0.55ausreichend
0.55-0.7gut
0.7-0.85sehr gut
0.85-0.99exzellent
0.99-1.0perfekt

Berechnung vom Kappa Index

Die Berechnung vom Kappa Index basiert auf einer Kreuztabelle von zwei Karten. Dabei ist die eine Karte in in Zeilen und die andere in Spalten angeordnet. Die Zeilensummen zeigen den Anteil der Pixel für jede Risikostufe in Karte A, die Spaltensummen für Karte B. Die Tabelle gibt an, wie viel % der Pixel in beiden Karten übereinstimmen (Diagonale) sowie wie viel % der Pixel in jeder Risikostufe nicht übereinstimmen. Die Kreuztabelle der Lawinenrisikostufen ist folgendermaßen aufgebaut:

Kategorie A/B geringes Risiko(1) erhöhtes Risiko(2) hohes Risiko(3) Summe
geringes Risiko(1) $$p_{1, 1}$$ $$p_{1, 2}$$ $$p_{1, 3}$$ $$ Tp_1= \sum_{i=1}^3 p_{1,i} $$
erhöhtes Risiko(2) $$p_{2, 1}$$ $$p_{2, 2}$$ $$p_{2, 3}$$ $$ Tp_2=\sum_{i=1}^3 p_{2,i} $$
hohes Risiko(3) $$p_{3, 1}$$ $$p_{3, 2}$$ $$p_{3, 3}$$ $$ Tp_3= \sum_{i=1}^3 p_{3,i} $$
Summe $$ pT_1= \sum_{i=1}^3 p_{i,1} $$ $$ pT_2= \sum_{i=1}^3 p_{i,2} $$ $$ pT_3= \sum_{i=1}^3 p_{i,3} $$ 1

Aus der Kreuztabelle werden PA (Anteil an Übereinstimmung), PE (erwarteter Anteil an Übereinstimmung) und Pmax (maximal erreichbarer Anteil an Übereinstimmung) berechnet: $$ PA= \sum_{i=1}^3 p_{i,i} $$ $$ PE= \sum_{i=1}^3 pT_i*Tp_i $$ $$ Pmax= \sum_{i=1}^3 min(pT_i,Tp_i) $$

der Kappaindex K lässt sich in Khisto und Klocation aufteilen {[:aval:pontius2000]},{[:aval:hagen2002]} und ist das Produkt aus beiden, ergibt sich also aus der quantitativen Ähnlichkeit und der Ähnlichkeit im Bezug auf die Lage der Kategorien. Khisto gibt die quantitative Ähnlichkeit von 2 Karten an und ergibt sich aus dem Histogramm der einzelnen Kategorien. Der Index sagt aus wie gut die Verteilung der Pixel auf die einzelnen Kategorien übereinstimmt. Klocation gibt die Ähnlichkeit der räumlichen Lage der Kategorien an. Die Ähnlichkeit der Karten wird skaliert auf die maximal erreichbare Ähnlichkeit entsprechend der gegebenen Pixelzahlen der Kategorien. D.h. Klocation ist der Kappaindex der erreicht werden könnte wenn Khisto 1 wäre, also die Pixelzahlen der Klassen übereinstimmen. $$K=\frac{PA-PE}{1-PE}=K_{histo}*K_{location}$$

$$K_{histo}=\frac{Pmax-PE}{1-PE}$$ $$K_{location}=\frac{PA-PE}{Pmax-PE}$$

Kappa Indizes können sowohl für die Gesamtkarte als auch pro Kategorie angegeben werden. Für die Berechnung des Kappas pro Kategorie werden die Kreuztabellen reklassifiziert, sodass sie nur aus der jeweiligen Zielkategorie und der Kombination aller anderen Kategorien als zweite Kategorie bestehen, wie von {[:aval:monserud1992]} vorgeschlagen.

Disagreements of quantity and allocation

{[:aval:pontius2011]} zeigen Schwächen des Kappaindex auf und schlagen daher die Verwendung von “disagreement of quantity” und “disagreement of allocation” als Maß für den Vergleich kategorialer Karten vor.

Berechnung der Disagreements

Die disagreements basieren ebenfalls auf der Kreuztabelle.

$$allocation \ disagreement=100*(PMax-PA)$$

$$quantity \ disagreement=100*(1-PMax)$$

zusätzlich können noch chance agreement, allocation agreement und quantity agreement angegeben werden, die vollständigkeitshalber auch hier erwähnt werden sollen, auch wenn nach {[:aval:pontius2011]} die Angabe der disagreements zur Beurteilung der Ähnlichkeit zweier Karten ausreicht.

$$chance \ agreement= 100*min(1/N,PA,PE)$$

$$quantity \ agreement= ifelse(min(1/N,PE,PA)==1/N,100*min((PE-1/N),PA-1/N),0)$$ $$allocation \ agreement= 100*max(PA-PE,0)$$

mit $N$ der Anzahl der Kategorien, in diesem Fall also 3.

….wird vervollständigt…

Implementierung in R

Das folgende Skript kann verwendet werden um aus 2 Klassifikationen die Kappaindizes sowie disagreements zu berechnen. Als Input werden 2 klassifizierte Raster benötogt.

Beispiel:In den Beispieldaten sind 2 Raster enthalten die räumlich das Lawinenrisiko aus 2 Achtung Lawinen darstellen enthalten. Hierfür soll jetzt der Kappa-Index berechnet werden.

Zunächst ist die unten stehende Funktion “kstat” nach R zu kopieren. Die Raster weisen die Klassen 0,1,2,3,4 auf, wobei 0 und 4 nicht von Interesse sind, da sie für beide Raster identisch als Maske überlagert wurden. Sie sollen in der Kappaberechnung ausmaskiert werden.

Um den Kappaindex zu berechnen kann folgendermaßen vorgegangen werden:

setwd("PFAD ZU DEN DATEN")
library(raster,rgdal)
r1=raster("SC.tif")
r2=raster("GRM2.tif")
kstat(r1,r2,mask=c("0","4"))

Die kstat Funktion sieht wie folgt aus:

1 |h R Funktion zur Berechnung der Kappa Indizes und Disagreements |h
kstat=function(a,b,mask="",perCategory=TRUE) {
  require(raster)
	ctab=function(a,b,mask=""){
	require(raster)
	va=values(a)
	vb=values(b)
	name=sort(union(unique(a),unique(b)))
	result=matrix(0,ncol=length(name),nrow=length(name))
	colnames(result)=name
	row.names(result)=name
	for (i in 1:nrow(result)){
		for (k in 1:ncol(result)){
      			result[i,k]=sum(va==row.names(result)[i]&vb==colnames(result)[k])
    		}
  	}
	for (i in 1:length(mask)){
		result=result[row.names(result)!=mask[i],]
		if (ncol(as.matrix(result))==1||nrow(as.matrix(result))==1){
			stop ("Calculation of contingency table requires at least two categories")
		}
		result=result[,colnames(result)!=mask[i]] 
	}
	if (ncol(as.matrix(result))==1||nrow(as.matrix(result))==1){
		stop ("Calculation of contingency table requires at least two categories")
	}	
	return(result)
	}
  ##################################################################################################
  #Function: Calculate kappa indices and disagreements
	calculateKappa=function(ct){
		KappaResult=matrix()
		ct=ct/sum(ct)#percent of pixels 
		cmax=nrow(ct)#number of categories
		#Fraction of Agreement:
		PA=0
		for (i in 1:cmax) {
			PA=PA+ct[i,i]
		}
		#Expected Fraction of Agreement subject to the observed distribution:
		PE=0
		for (i in 1:cmax) {
			PE=PE+sum(ct[i,])*sum(ct[,i])
		}
		#Maximum  Fraction  of  Agreement  subject  to  the  observed  distribution:
		PMax=0
		for (i in 1:cmax) {
			PMax=PMax+min(sum(ct[i,]),sum(ct[,i]))
		}	
		#Kappa Index:
		K=(PA-PE)/(1-PE)
		#Kappa of location:
		Kloc=(PA-PE)/(PMax-PE)
		#Kappa of histogram:
		Khisto=(PMax-PE)/(1-PE)
		#chance agreement:
		CA=100*min((1/cmax),PA,PE)
		#quantity agreement:
		QA=ifelse(min((1/cmax),PE,PA)==(1/cmax),100*min((PE-1/cmax),PA-1/cmax),0)
		#allocation agreement:
		AA=100*max(PA-PE,0)
		#allocation disagreement:
		AD=100*(PMax-PA)
		#quantity disagreement:
		QD=100*(1-PMax)
		KappaResult=cbind(K,Kloc,Khisto,CA,QA,AA,AD,QD)	
		return (KappaResult)
	}
  ##################################################################################################
  #Initialise
  if (any(dim(a)!=dim(b))){
    stop ("Dimensions of a and b do not match!")
  }
  ct=ctab(a,b,mask=mask)
	result=list()
  #reclass to calculate kappa per category
	cttmp=ct
	if (ncol(cttmp)<=2||perCategory==FALSE){
		result[[1]]=calculateKappa(ct)
		if (perCategory==TRUE){
			print ("Warning: Kappa per category requires at least 3 categories")
		}
	}
	if (perCategory==TRUE&ncol(cttmp)>2) {
		for (ca in 1:(ncol(cttmp)+1)){
			if (ca==ncol(cttmp)+1) {
				ct=cttmp
			} 
			else {
				ct=cttmp
				ctNew=ct[1:2,1:2]
				ctNew[1,1]=ct[ca,ca]
				ctNew[1,2]=sum(ct[ca,])-ct[ca,ca]
				ctNew[2,1]=sum(ct[,ca])-ct[ca,ca]
				ctNew[2,2]=sum(ct)-sum(ct[ca,])-sum(ct[,ca])+ctNew[1,1]
				ct=ctNew
			}
			result[[ca]]=calculateKappa(ct)
		}
	}
  ##################################################################################################
  #arrange results	
	resultTab=result[[1]]
	if (length(result)>1){
		for (i in 2:length(result)){
			resultTab=rbind(resultTab,result[[i]])
		}
	}
	resultTab=as.data.frame(resultTab)
	colnames(resultTab)=c("K","Kloc","Khisto","CA","QA","AA","AD","QD")
	if (nrow(resultTab)>1){
		row.names(resultTab)=c(colnames(cttmp),"overall")	
	}
	else row.names(resultTab)="overall"
	return (resultTab)
}

Hilfe zur R Funktion

Die R Funktion ist im package Rsenal enthalten, die aus dem git repository der Umweltinformatik Marburg heruntergeladen werden. Im Folgenden ist die Hilfeseite der kstat Funktion aus dem package angegeben:

Description

Calculation of Kappa indices as well as quantity and location (dis)agreements to assess similarity of two categorical maps. Might be used mainly to compare a raster of ground truth sites with a classified map or to compare a model output with a reference map.

Usage

kstat(a, b, mask= “”, perCategory = TRUE)

Arguments

a: An object of type “raster” containing categorical data

b: An object of type “raster” containing categorical data. Must be of the same dimensions as a.

mask: string vector of category numbers which are are to be ignored in the calculations. The numbers must be identical to the values of the input rasters.

perCategory: If TRUE indices are calculated per category in addition to overall indices.

Details

Calculation of kappa index, kappa of quantity, Kappa of location (see Cohen 1960, Pontius 2000, Hagen 2002) for the overall map and for each category (Monserud and Leemans 1992) as well as quantity and location (dis)agreements according to Pontius (2011).

The first step to calculate the indices is to calculate a contingency table:

From the contingency table the following statistics are calculated:

From theses statistics, the kappa indices are calculated as:

Disagreements are calculated as:

with cmax is the number of categories.

The calculations of Kappa per category require a reclassification for each category:

The Kappa per category i is calculate in the same way as the normal Kappa but from the reclassified contingency table. Thus, the formula for the Kappa per category is

as introduced by Fleiss (1981) and applied in e. g. Monserud 1992 and Hagen (2002) The Disagreements per category are as well calculated as explained above but from the reclassified contingency table for each category.

Value

A data frame containing the values for overall kappa (K), Kappa of location (Kloc), Kappa of histogram (Khisto), chance agreement (CA), quantity agreement (QA), allocation agreement (AA), allocation disagreement (AD), quantity disagreement (QD).

Note

The calculation of kappa per category bases on the approach of Fleiss (1981). Dont get confused with the Kappa per category calculated in e.g. Idrisi which follows the approach of Light (1971) and is an asymmetric Kappa index depending on which map is treated as reference map.

Author(s)

Hanna Meyer

References

Cohen, J., (1960): A coefficient of agreement for nominal scales, Educational and Psychological Measurements, 20: 37-46.

Fleiss, J.L., (1981): Statistical Methods for Rates and Proportions (2nd edition). John Wiley and Sons, New York, NY, 321 pp.

Hagen, A. (2002): Multi-method assessement of map similarity. Proceeding of 5 th AGILE Conference on Geographic Information Science, April. 25-27 Spain, 171-182.

Light RJ (1971) Measure of response agreement for qualitative data: some generalizations and alternatives. Psychol Bull 76:365-377.

Monserud, R. A., and Leemans, R.(1992): Comparing global vegetation maps with the kappa statistic. Ecological Modeling, 62: 275-293.

Pontius, R.G. (2000): Quantifcation Error Versus Allocation Error in Comparison of Categorical Maps. Photogrammatric Engineering and Remote Sensing, 66: 1011-1016.

Pontius, R.G., Jr., and Millones, M. (2011): Death to Kappa: Birth of Quantity Disagreement and Allocation Disagreement for Accuracy Assessment. International Journal of Remote Sensing, 32: 4407-4429.

Examples

  library(Rsenal)
 
  #### Example 1: Calculate the indices to compare two Land cover maps
  #load data
  library(raster)
  LUC1990=raster(system.file("LUC1990.rst", package="Rsenal"))
  LUC2006=raster(system.file("LUC2006.rst", package="Rsenal"))
 
  #plot to get an impression:
  par(mfrow=c(1,2))
  plot(LUC1990,main="LUC Marburg 1990")
  plot(LUC2006,main="LUC Marburg 2006")
 
  #calculate indices:
  kstat(LUC1990,LUC2006)
 
  #### Example 2: Calculate the indices to compare a land cover map 
  ####with a map of ground truth sites
  #load data:
  LUC1990=raster(system.file("LUC1990.rst", package="Rsenal"))
  trainingSites=raster(system.file("training.rst", package="Rsenal"))
 
  #plot to get an impression:
  par(mfrow=c(1,2))
  plot(LUC1990,main="LUC Marburg 1990")
  plot(trainingSites,main="Training sites (Ground truth)")
 
  #calculate indices
  #(consider to mask the background since ground truth sites don't 
  #cover the whole area)
  kstat(LUC1990,trainingSites,mask="0")
 
  ###
  #For an example of an other field of application of these indices see 
  #http://giswerk.org/doku.php?id=people:c:creudenbach:alpenwerk:achtunglawine
 

Hanna Meyer

r/r-tutorials/calculatekappa.txt · Last modified: 2018/12/23 19:46 (external edit)