Risks as hazardous (false positive) and safe (false negative) - R-SIG-Geo

Mon, Nov 23, 2009 7:47 AM #

(replying to the list)

Well, if you're willing to accept the indicator kriging values as
estimates of these probabilities, you're essentially done. Of course,
the kriged values can be outside [0, 1], so you need to deal with that.

Best regards,
--
Edzer

Zia Ahmed wrote:

Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de

Zia Ahmed

Mon, Nov 23, 2009 8:22 AM #

Thanks!  I may correct the values those are out side [0,1] using pmin 
and pmax function (?).
But  I am sill not clear how I will do this.
 Thanks again

Zia

Edzer Pebesma wrote:

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091123/083e46fc/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: zua3.vcf
Type: text/x-vcard
Size: 281 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20091123/083e46fc/attachment.vcf>

Zia Ahmed

Mon, Nov 23, 2009 6:26 PM #

/Dear Edzer,
/
I am trying to do  misclassification following way: could you please 
tell me know whether  I am  statistically correct or wrong!

Thanks again
Zia/

# Load data
#-------------
tala<-read.csv("Tala_data.csv",header=TRUE)
tala.grid<-read.csv("Tala_grid_2.csv",header=TRUE)
# load packages:
#----------------
library(sp)
library(gstat)
library(lattice)
library(geoR)
library(MASS)
library(car)
coordinates(tala) <- ~ x + y       # Observed data
coordinates(tala.grid) <-~ x + y   # Prediction locations

# Box-cox transformation; required pakage- car
tala$was.bc<-box.cox(tala$was, .47)  # Power (lambda)= .47

# Varigram modeling:
#---------------------
v.ok<-variogram(was.bc~1,data=tala)
plot(v.ok, pl=F, pch=20, cex=1, col="Black")
m.ok<-vgm(.06,"Exp",4000,0.02)
(m.ok.f<-fit.variogram(v.ok, m.ok))
attr(m.ok.f,"SSErr") # Sum of Squared Error (SSE)

# Ordinary Kriging
#------------------

ok.was<-krige(was.bc~1,tala,tala.grid, model=m.ok.f, nmax=50)

# Back-and Indicator- transformation of OK prediction:
#-----------------------------------------------------------------------------------------
#Power=0.47
k<-1/0.47                                        
ok.was$was.bt <-((ok.was$var1.pred *0.47+1)^k)

target <- box.cox(0.100,0.47)
ok.was$was.target <- (ok.was$var1.pred>= target)
ok.was$p.target <- pt((-target +ok.was$var1.pred)/sqrt(ok.was$var1.var),
    length(tala$was.bc))
summary(ok.was)
coordinates(ok.was)<-~x+y

# Probability of true indicator
#-----------------------------------------------

plot(coordinates(ok.was), asp = 1, col = ifelse(ok.was$was.target,
"grey", "yellow"),//"Easting (m)", ylab="Northing (m)",//  main = 
"Probability of TRUE indicator",
sub = "Actual indicator: yellow/grey = FALSE/TRUE")
grid()

#-----------------------------------------------------------------
# Probability-of-exceeding 200 ppb ground water As conc.
#--------------------------------------------------------------------
levelplot(p.target~x+y| was.target, main="  (b) Probability > 0.200 mg 
As/L",
              xlab="Easting (m)", ylab="Northing (m)",
              as.data.frame(ok.was), aspect = "iso",at = seq(0, 1, by = 
0.05),
              col.regions=topo.colors,                            
              panel = function(...) {
              panel.levelplot(...)
              panel.abline(h = 0:4*5000 + 545000, v= 0:4*5000 + 2650000,
              col = "light grey")
},
     )

/

Edzer Pebesma wrote:

Edzer Pebesma

Mon, Nov 23, 2009 11:47 PM #

Zia Ahmed wrote:

Zia, to do this, I would need your data, a lot of information about your
data, and a lot of time. And then, the answer would be a long one and
not one in the sense of "yes" or "no". I'm sorry, but you either
overestimate my capabilities, or those of statistical methods in
general, or both!

Thanks again
Zia/

# Load data
#-------------
tala<-read.csv("Tala_data.csv",header=TRUE)
tala.grid<-read.csv("Tala_grid_2.csv",header=TRUE)
# load packages:
#----------------
library(sp)
library(gstat)
library(lattice)
library(geoR)
library(MASS)
library(car)
coordinates(tala) <- ~ x + y       # Observed data
coordinates(tala.grid) <-~ x + y   # Prediction locations

# Box-cox transformation; required pakage- car
tala$was.bc<-box.cox(tala$was, .47)  # Power (lambda)= .47

# Varigram modeling:
#---------------------
v.ok<-variogram(was.bc~1,data=tala)
plot(v.ok, pl=F, pch=20, cex=1, col="Black")
m.ok<-vgm(.06,"Exp",4000,0.02)
(m.ok.f<-fit.variogram(v.ok, m.ok))
attr(m.ok.f,"SSErr") # Sum of Squared Error (SSE)

# Ordinary Kriging
#------------------

ok.was<-krige(was.bc~1,tala,tala.grid, model=m.ok.f, nmax=50)

# Back-and Indicator- transformation of OK prediction:
#-----------------------------------------------------------------------------------------

#Power=0.47
k<-1/0.47                                        ok.was$was.bt
<-((ok.was$var1.pred *0.47+1)^k)

target <- box.cox(0.100,0.47)
ok.was$was.target <- (ok.was$var1.pred>= target)
ok.was$p.target <- pt((-target +ok.was$var1.pred)/sqrt(ok.was$var1.var),
   length(tala$was.bc))
summary(ok.was)
coordinates(ok.was)<-~x+y

# Probability of true indicator
#-----------------------------------------------

plot(coordinates(ok.was), asp = 1, col = ifelse(ok.was$was.target,
"grey", "yellow"),//"Easting (m)", ylab="Northing (m)",//  main =
"Probability of TRUE indicator",
sub = "Actual indicator: yellow/grey = FALSE/TRUE")
grid()

#-----------------------------------------------------------------
# Probability-of-exceeding 200 ppb ground water As conc.
#--------------------------------------------------------------------
levelplot(p.target~x+y| was.target, main="  (b) Probability > 0.200 mg
As/L",
             xlab="Easting (m)", ylab="Northing (m)",
             as.data.frame(ok.was), aspect = "iso",at = seq(0, 1, by =
0.05),
             col.regions=topo.colors,                           
             panel = function(...) {
             panel.levelplot(...)
             panel.abline(h = 0:4*5000 + 545000, v= 0:4*5000 + 2650000,
             col = "light grey")
},
    )

/
Edzer Pebesma wrote:

(replying to the list)

Well, if you're willing to accept the indicator kriging values as
estimates of these probabilities, you're essentially done. Of course,
the kriged values can be outside [0, 1], so you need to deal with that.

Best regards,
-- 
Edzer

Zia Ahmed wrote:

Dear Edzer,

Thanks for your replay. I have done IK of water arsenic using
quantiles of the data. I want to do   miss classification  of my risk
that  proposed  by  Goovaerts.
Two misclassification  of risks are proposed by Goovaerts (1997) based
on ccdf model F(*u*;z_k |(n)) of Z(*u*)?Z_k .


1.       The risk ?(*u*) of incorrectly classified a location u as
hazardous (false positive):

?(*u*) = Prob{Z(*u*) ?z_c |z*_L (u) > Z_c , (n)}

       = F(*u*;z_c |(n))

Where z*_L (u) is the estimator.

2.       The risk ?(u)  wrongly classified a location u as safe (false
negative)

?(*u*) = Prob{Z(*u*) > z_c |z*_L (u) ? Z_c , (n)}

      = 1- F(*u*;z_c |(n))





Ref: Goovaerts, P (1997). Geostatistics for Natural Resource
Evaluation, pp: 259-368



Zia

Edzer Pebesma wrote:

Dear Zia,

If you mean by Goovaerts (1997) his 483 page book "Geostatistics for
Natural Resource Evaluation", then please give us a bit more precise
detail on what you want (which page? which equation?)
-- 
Edzer

Zia Ahmed wrote:

Hi all,

Is it possible in gstat to get misclassification  of risks as
hazardous (false positive) and safe (false negative) as describe by
Goovaerts (1997)  from parametric or non-parametric probability
mapping?.  help will appreciated.

Thanks

Zia

Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de