Nothing
################ MEVA.TABLE: Function originally from A. Guisan (Unil-ECOSPAT, Switzerland)
################ Modified from L. Maiorano, O. Broennimann and F. Collart (Unil-ECOSPAT, Switzerland)
ecospat.meva.table <- function(Pred, Sp.occ, th) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
# th: threshold used to cut the probability to binary values
{
# eva <- list()
tab <- table(Pred >= th, Sp.occ, dnn = list("Predicted values", "Observed values"))
a <- tab[4]
b <- tab[2]
c <- tab[3]
d <- tab[1]
N <- a + b + c + d
prev <- (a + c)/N
ccr <- (a + d)/N
mcr <- (b + c)/N
se <- a/(a + c)
sp <- d/(b + d)
fpr <- b/(b + d)
fnr <- c/(a + c)
ppp <- a/(a + b)
npp <- d/(c + d)
or <- (a * d)/(c * b)
kappa <- ecospat.cohen.kappa(tab)$kap
nmi <- (-(a * log(a)) - (b * log(b)) - (c * log(c)) - (d * log(d)) + ((a + b) *
log(a + b)) + ((c + d) * log(c + d)))/((N * log(N)) - ((a + c) * log(a +
c)) - ((b + d) * log(b + d)))
tss <- se + sp - 1
orss <- (a*d - b*c)/(a*d + b*c)
sed <- ( log(fpr) - log(se) - log(1-fpr) + log(1-se) ) /( log(fpr) + log(se) + log(1-fpr) + log(1-se))
mtrcNAME <- c("Prevalence", "Correct classification rate", "Mis-classification rate",
"Sensitivity", "Specificity", "False positive rate", "False negative rate",
"Positive predictive power", "Negative predictive power", "Odds Ratio", "Kappa", "Normalized mutual information",
"True skill statistic","Odds Ratio Skill Score","Symmetric Extremal Dependence Index")
mtrcVALUE <- c(round(prev, digits = 4), round(ccr, digits = 4), round(mcr, digits = 4),
round(se, digits = 4), round(sp, digits = 4), round(fpr, digits = 4), round(fnr,
digits = 4), round(ppp, digits = 4), round(npp, digits = 4), round(or,
digits = 4), round(kappa, digits = 4), round(nmi, digits = 4), round(tss,
digits = 4),round(orss, digits = 4),round(sed, digits = 4))
mtrcMATR <- cbind.data.frame(Metric=mtrcNAME, Value =mtrcVALUE)
return(list(CONTINGENCY_TABLE = tab, EVALUATION_METRICS = mtrcMATR))
}
############### MAX-KAPPA ## Function originally from A. Guisan (Unil-ECOSPAT, Switzerland)
############### Modified from L. Maiorano (Unil-ECOSPAT, Switzerland) and Olivier Broennimann
ecospat.max.kappa <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
FUN<-function(thresh){
xtab<-table(Pred >= thresh, Sp.occ)
return(ecospat.cohen.kappa(xtab)$kap )
}
threshold<-(1:100)/100
k<-sapply(threshold,FUN)
table<-data.frame(cbind(threshold,k))
max.Kappa <- max(table$k,na.rm = TRUE)
max.threshold <- table$threshold[which.max(table$k)][1]
return(list(table=table,max.Kappa=max.Kappa,max.threshold=max.threshold))
}
############### MAX-TSS ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
############### modyfying MAX-KAPPA from A. Guisan (Unil-ECOSPAT, Switzerland), modified by Olivier Broennimann and Flavien Collart
ecospat.max.tss <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
FUN<-function(thresh){
xtab<-table(Pred >= thresh, Sp.occ)
if(nrow(xtab)!=2){ ## if the model is bad the Pred values are not changing across the whole dataset
return(0)
}else{
a <- xtab[4]
b <- xtab[2]
c <- xtab[3]
d <- xtab[1]
se <- a/(a + c)
sp <- d/(b + d)
return(se + sp - 1)
}
}
threshold<-(1:100)/100
tss<-sapply(threshold,FUN)
table<-data.frame(cbind(threshold,tss))
max.TSS <- max(table$tss,na.rm = TRUE)
max.threshold <- table$threshold[which.max(table$tss)][1]
return(list(table=table,max.TSS=max.TSS,max.threshold=max.threshold))
}
################## PLOT-K ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
ecospat.plot.kappa <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
k <- 0.01
i <- 1
evak <- data.frame(k)
while (k < 1) {
a <- table(Pred >= k, Sp.occ)[4]
b <- table(Pred >= k, Sp.occ)[2]
c <- table(Pred >= k, Sp.occ)[3]
d <- table(Pred >= k, Sp.occ)[1]
N <- a + b + c + d
tab <- table(Pred >= k, Sp.occ)
evak[i, "THRESHOLD"] <- k
evak[i, "k"] <- ecospat.cohen.kappa(tab)$kap
# evak[i,'k'] <- ecospat.cohen.kappa(table(Pred >= k, Sp.occ))$kap
k <- k + 0.01
i <- i + 1
}
plot(evak$THRESHOLD, evak$k, type = "l", xlab = "THRESHOLD", ylab = "KAPPA")
}
################## PLOT-TSS ## Function originally from L. Maiorano (Unil-ECOSPAT, Switzerland)
ecospat.plot.tss <- function(Pred, Sp.occ) # Pred: vector of predicted probabilities Sp.occ: vector of binary observations
{
tss <- 0.01
i <- 1
evatss <- data.frame(tss)
while (tss < 1) {
a <- table(Pred >= tss, Sp.occ)[4]
b <- table(Pred >= tss, Sp.occ)[2]
c <- table(Pred >= tss, Sp.occ)[3]
d <- table(Pred >= tss, Sp.occ)[1]
N <- a + b + c + d
se <- a/(a + c)
sp <- d/(b + d)
evatss[i, "THRESHOLD"] <- tss
evatss[i, "tss"] <- se + sp - 1
tss <- tss + 0.01
i <- i + 1
}
plot(evatss$THRESHOLD, evatss$tss, type = "l", xlab = "THRESHOLD", ylab = "TRUE SKILL STATISTIC")
}
############################################################################################################################
ecospat.cohen.kappa <- function(xtab) {
# computes Cohen's kappa and variance estimates, 95% CI from a symmetric
# agreement table xtab return value is a list with elements kap and vark See
# Bishop et al. Disc. Mult. Analysis pp. 395-397 1975
# Originally coded by C. Randin, UNIL Adjusted by N.E.Zimmermann, WSL
ci <- vector(length = 2)
kap <- 0
vark <- 0
totn <- 0
if (nrow(xtab) != ncol(xtab)) {
# cat('\n freq. table not symmetric.\n')
result = list(kap = kap, vark = vark, totn = totn)
return(result)
}
# compute obs props:
totn <- sum(xtab)
# calc marginals, expected diag props.
obs <- xtab/totn
pi <- apply(xtab, 1, sum)/totn
pj <- apply(xtab, 2, sum)/totn
exp <- outer(pi, pj)
theta.2 <- sum(diag(exp))
theta.1 <- sum(diag(obs))
kap <- (theta.1 - theta.2)/(1 - theta.2)
theta.3 <- diag(obs) %*% (pi + pj)
theta.4 <- sum(obs * outer(pi, pj, FUN = "+")^2)
# cat('\n thetas 1,2,3,4 \n', theta.1, theta.2, theta.3, theta.4, '\n') calc
# the var est
vark <- (theta.1 * (1 - theta.1))/(1 - theta.2)^2
vark <- vark + (2 * (1 - theta.1) * (2 * theta.1 * theta.2 - theta.3))/(1 - theta.2)^3
vark <- vark + ((1 - theta.1)^2 * (theta.4 - 4 * theta.2^2))/(1 - theta.2)^4
vark <- as.vector(vark/totn)
ci[1] <- kap - 1.96 * sqrt(vark)
ci[2] <- kap + 1.96 * sqrt(vark)
result = list(kap = kap, vark = vark, totn = totn, ci = ci)
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.