## 1006Re: AI-GEOSTATS: kriging estimates in S+

Expand Messages
• May 9, 2003
Vanessa,

you are right. I was confused, as MASS has a free library
called "spatial", but you meant the commercial add-on module
"S+SpatialStats", which is also called "module spatial" by
the vendor of S-Plus.

It seems that S+SpatialStats, when predicting at point locations,
does as if the locations are shifted a very small amount; I guess
this is saying the same thing that Isobel meantioned earlier (see
my sample below). When predicting without a nugget effect, the
differences are also too large, IMO (see x2 example); I do not
have an explanation for that.

I also tried the examples with the (my) gstat S library, which can
what you'd expect; again see below.

Thanks for sharing this with us.
--
Edzer

Here's the output of the S-Plus session:

S-PLUS : Copyright (c) 1988, 2002 Insightful Corp.
S : Copyright Lucent Technologies, Inc.
Version 6.1.2 Release 2 for Linux 2.2.12 : 2002
Working data will be in .Data
> module(spatial)
> x<-as.numeric(1:10)
> y<-as.numeric(1:10)
> z<-as.numeric(1:10)
> d<-data.frame(x=x,y=y,z=z)
> d
x y z
1 1 1 1
2 2 2 2
3 3 3 3
4 4 4 4
5 5 5 5
6 6 6 6
7 7 7 7
8 8 8 8
9 9 9 9
10 10 10 10
> krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=1,nugget=1)
Call:
krige(formula = z ~ loc(x, y), data = d, covfun = gauss.cov, range = 3, sill =
1, nugget = 1)

Coefficients:
constant
5.5

Number of observations: 10
> x<-krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=1,nugget=1)
> predict(x, d)
x y fit se.fit
1 1 1 2.827990 1.206120
2 2 2 2.778034 1.155678
3 3 3 3.432520 1.152567
4 4 4 4.299386 1.153139
5 5 5 5.113279 1.153538
6 6 6 5.886721 1.153538
7 7 7 6.700614 1.153139
8 8 8 7.567480 1.152567
9 9 9 8.221966 1.155678
10 10 10 8.172010 1.206120
> d2<-data.frame(x=y+1e-7,y=y+1e-7)
> d2
x y
1 1 1
2 2 2
3 3 3
4 4 4
5 5 5
6 6 6
7 7 7
8 8 8
9 9 9
10 10 10
> predict(x, d2)
x y fit se.fit
1 1 1 2.827990 1.206120
2 2 2 2.778034 1.155678
3 3 3 3.432520 1.152567
4 4 4 4.299386 1.153139
5 5 5 5.113279 1.153538
6 6 6 5.886721 1.153538
7 7 7 6.700614 1.153139
8 8 8 7.567480 1.152567
9 9 9 8.221966 1.155678
10 10 10 8.172010 1.206120
> x2<-krige(z~loc(x,y),data=d,covfun=gauss.cov,range=3,sill=2,nugget=0)
> predict(x2, d)
x y fit se.fit
1 1 1 1.000004 0.001272791
2 2 2 1.999994 0.001272785
3 3 3 3.000009 0.001272774
4 4 4 3.999990 0.001272763
5 5 5 5.000010 0.001272757
6 6 6 5.999990 0.001272757
7 7 7 7.000010 0.001272763
8 8 8 7.999991 0.001272774
9 9 9 9.000006 0.001272785
10 10 10 9.999996 0.001272791
>
# Now using the gstat library:
> library(gstat)
> krige(z~1,~x+y,d,d,vgm(1,"Gau",3,1))
[using ordinary kriging]
x y var1.pred var1.var
1 1 1 1 0.000000e+00
2 2 2 2 0.000000e+00
3 3 3 3 5.185705e-33
4 4 4 4 2.074282e-32
5 5 5 5 5.185705e-33
6 6 6 6 2.074282e-32
7 7 7 7 0.000000e+00
8 8 8 8 5.185705e-33
9 9 9 9 0.000000e+00
10 10 10 10 0.000000e+00
> krige(z~1,~x+y,d,d,vgm(2,"Gau",3,0))
[using ordinary kriging]
x y var1.pred var1.var
1 1 1 1 7.106437e-33
2 2 2 2 0.000000e+00
3 3 3 3 4.440892e-16
4 4 4 4 2.220446e-16
5 5 5 5 2.220446e-16
6 6 6 6 4.440892e-16
7 7 7 7 2.220446e-16
8 8 8 8 2.842575e-32
9 9 9 9 2.842575e-32
10 10 10 10 0.000000e+00
>

--
* To post a message to the list, send it to ai-geostats@...
* As a general service to the users, please remember to post a summary of any useful responses to your questions.
* To unsubscribe, send an email to majordomo@... with no subject and "unsubscribe ai-geostats" followed by "end" on the next line in the message body. DO NOT SEND Subscribe/Unsubscribe requests to the list
* Support to the list is provided at http://www.ai-geostats.org
• Show all 7 messages in this topic