I've generated a map of ecoregions based on some environmental data using a semi-supervised classification where I combined a statistical clustering with random forests to generate my map.

I created the map at a 30m resolution and to remove some of the small scale variability I've re-sampled (with the majority filter) to 2.5km. I was thinking of using McNemars test to investigate for any significant differences between my original (30m) and re-sampled (2.5km) data, but my validation dataset becomes significantly clustered at less than 2.5km (according to Ripley's K). Also, my validation dataset is quite large (~2000 samples) and from what I've noticed so far is that the difference is coming back as significant (p<<0.001) using McNemars test even though % agreement is ~95% with a kappa of .96. I think it's due to the large sample size of my validation set.

Would it be suitable to compare the differences in the map with something like a Mann-Whitney or Wilcoxon test on an independent random sample of pixels? i.e. create x amount of sampling points, extract class info from the 30m and 2.5km maps and then compare?

The other thing I've thought about doing, was to look at lumped class accuracy with McNemars test i.e. lump the class information from my validation sites on both maps and compare those counts across classes to the validation dataset. An issue I can see here though (apart from the clustering of my validation samples at less than 2.5km) is that there only needs to be a difference of 1 pixel between the validation set and my classification sets for it to counted as misclassified.

Any help would be greatly appreciated.

1 Answer 1


In case anybody else ever needs an answer to this question, I decided after a bit of research that a permutation based method comparing kappa values would be suitable. McNemars test would've been suitable if my classifications were independent of one another, but as they were obviously not the test wasn't. So I randomly sampled 300 of the 2000 validation points (making sure they were at a minimum 2.5 km away from the other points) and then assessed for differences using the R code below.

# load library and set random seed

## test statistics = difference between two kappa values

## This function uses a permutation test to assess if two kappa test results are significantly different.
perm_test <- function(input_data_frame,validation,classifier1,classifier2, iterations) {

## format input_data
scores <- input_data_frame[,c(validation, classifier1, classifier2)]

## Observer 1 kappa
classifier1.kappa <- cohen.kappa(scores[, c(validation, classifier1)], alpha=.05)

## Observer 2 kappa
classifier2.kappa <- cohen.kappa(scores[, c(validation, classifier2)], alpha=.05)

observed.kappa.difference <- round(classifier1.kappa$kappa,3) - round(classifier2.kappa$kappa,3)

## Permutation test
N.perm <- iterations

perm.kappa.differences <- c(rep(NA, length = N.perm), observed.kappa.difference)

for (i in 1:N.perm) {

    ## Under the null hypothesis that classifier1's and classifier2's scoring
    ## is identical, we can exchange classifier1's and classifier2's scores.
    perm.scores <- scores
    perm.scores[, c(classifier1, classifier2)] <-
        t(apply(perm.scores[, c(classifier1, classifier2)], 1,
                function(x) {x[sample(1:2)]} ))
    perm.kappa.differences[i] <-
            round(cohen.kappa(perm.scores[, c(validation, classifier1)], alpha=.05)$kappa,3) -
            round(cohen.kappa(perm.scores[, c(validation, classifier2)], alpha=.05)$kappa,3)

p.value <-
    mean(perm.kappa.differences <= (-1.0 * abs(observed.kappa.difference)) |
            perm.kappa.differences >= abs(observed.kappa.difference))

results <- p.value
attr(results, 'perm.kappa.differences') <- perm.kappa.differences
attr(results, 'observed.kappa.difference') <- observed.kappa.difference
return (results)

## read data file
## data file must be set-out: PointID, Validation, Classifier1, Classifier2...
data <- read.csv(file.choose(), header=TRUE)
names <- attr(data,'names')

## ALWAYS assume the following
## names[1] is the point ID column name
## names[2] is the validation column name
## names[3:length(names)] is the classifiers column names
validation <- names[2]
classifiers <- names[3:length(names)]

## set the number of iterations for permutation ##
N.perm <- 9999
print(paste("Number of iterations: ",N.perm+1), quote=FALSE)

## create empty matrices for test statistics ##
p.values <- matrix(rep(1,(length(classifiers)^2)), nrow=length(classifiers), ncol=length(classifiers))
observed.kappa.differences <- matrix(rep(0,(length(classifiers)^2)), nrow=length(classifiers), ncol=length(classifiers))
perm.kappa.differences <- array(rep(0,((length(classifiers)^2)*length(data[,1]))), c(length(classifiers), length(classifiers),(N.perm+1)))

## loop to do the perm test
row_count <- 1
for (classifier1 in classifiers[1:(length(classifiers)-1)]) {
col_count <- row_count+1
for (classifier2 in classifiers[(row_count+1):length(classifiers)]) {
    ## perform test
    results <- perm_test(data,validation,classifier1,classifier2,N.perm)
    ## record p-values
    p.values[row_count,col_count] <- results
    p.values[col_count,row_count] <- p.values[row_count,col_count]
    ## record observed.kappa.differences
    observed.kappa.differences[row_count,col_count] <- attr(results,'observed.kappa.difference')
    observed.kappa.differences[col_count,row_count] <- attr(results,'observed.kappa.difference')
    ## record perm kappa differences
    perm.kappa.differences[row_count,col_count,] <- attr(results, 'perm.kappa.differences')
    perm.kappa.differences[col_count,row_count,] <- attr(results, 'perm.kappa.differences')
    col_count <- col_count + 1  
row_count <- row_count+1
print(paste("Observed Kappa Difference =",round(observed.kappa.differences[1,2],3),"(after",N.perm+1,"permutations)" ), quote=FALSE)

if (p.values[1,2] <= 0.05)
    print("Observed Kappa difference is significant at .05", quote = FALSE)
    print(paste("p =", p.values[1,2]), quote=FALSE)
if (p.values[1,2] > 0.05)
    print("Observed Kappa difference is insignificant at .05", quote = FALSE)
    print(paste("p =", p.values[1,2]), quote=FALSE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.