Uncertainty_assessment.R

library(sp)
library(lattice)

data(meuse)
coordinates(meuse)=~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

# Cross validation
loo <- krige.cv(log(zinc) ~ sqrt(dist), meuse, meuse.grid, 
                model = vgm(.149, "Sph", 700, .0674), 
                nmax = 40, nmin = 20, maxdist = 1000)
loo <- na.omit(loo)


# Does the value fall in the right PI ? (Goovaert 2001, equ. 13)
q <- seq(from=0.01, to=0.99, by=0.01)
m <- rep(NA,99)
a=0
ch <- rep(NA,length(loo$observed))

for (j in q){
  
  a=a+1
  min_q <- qnorm((1-j)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  max_q <- qnorm((1+j)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  for (i in 1:length(loo$observed)){
    if(loo$observed[i]>=min_q[i] & loo$observed[i]<=max_q[i]){ch[i]=1}else{ch[i]=0}
  }
  m[a] <- mean(ch, na.rm=T)
  
}

plot(q,m,xlim=c(0,1), ylim=c(0,1))
abline(0,1)


# absolute area between the 1:1 line and the curve
FUN <- function(x){
  min_q <- qnorm((1-x)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  max_q <- qnorm((1+x)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  for (i in 1:length(loo$observed)){
    if(loo$observed[i]>=min_q[i] & loo$observed[i]<=max_q[i]){ch[i]=1}else{ch[i]=0}
  }
  m <- mean(na.omit(ch))
  abs(x-m)
  
}
abs.dev <-  integrate(f=Vectorize(FUN), lower = 0, upper=1, subdivisions = 1000)$value
abs.dev 

# proportion of the deviation due to overestimation of the uncertainty
FUN <- function(x){
  min_q <- qnorm((1-x)/2,  mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  max_q <- qnorm((1+x)/2,  mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  for (i in 1:length(loo$observed)){
    if(loo$observed[i]>=min_q[i] & loo$observed[i]<=max_q[i]){ch[i]=1}else{ch[i]=0}
  }
  m <- mean(na.omit(ch))
  if(m>x){abs(x-m)}else{0}
  
}
over.unc <-  integrate(f=Vectorize(FUN), lower = 0, upper=1, subdivisions = 100)$value
over <- (over.unc*100)/abs.dev 


# proportion of the deviation due to overestimation of the uncertainty
FUN <- function(x){
  min_q <- qnorm((1-x)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  max_q <- qnorm((1+x)/2, mean=loo$var1.pred, sd=sqrt(loo$var1.var))
  for (i in 1:length(loo$observed)){
    if(loo$observed[i]>=min_q[i] & loo$observed[i]<=max_q[i]){ch[i]=1}else{ch[i]=0}
  }
  m <- mean(na.omit(ch))
  if(m<x){abs(x-m)}else{0}
  
}
under.unc <-  integrate(f=Vectorize(FUN), lower = 0, upper=1, subdivisions = 100)$value
under <- (under.unc*100)/abs.dev

# over- and under-estimation must sum to 100
over+under