[R-sig-Geo] Moving window correlation- Further
Jaime R. Garcia M.
jaime.garcia at uni-bonn.de
Wed Dec 9 22:27:53 CET 2009
On Wed, 2009-12-09 at 11:46 -0500, Zia Ahmed wrote:
> Thank you so much! Before running my data, I have tested your code
> using "meuse" data. It works fine. But it produce following output
> like below from V1 to V3103. I do not understand what does the values
> in rows mean (1 to 4).
> Here is the code that I used. I have questions How do I plot
> r-values with coordinates. Is it possible to get p-value for each
> estimate?
>
> Thanks again
> Zia
>
> # R code:
> #----------
> library(gstat)
> library(sp)
> library(spdep)
>
> data(meuse)
> data(meuse.grid)
> coordinates(meuse) <- ~x + y
> coordinates(meuse.grid) <- ~x + y
> str(meuse.grid)
>
> # GRID 1:
> v1 <- variogram(log(zinc) ~ 1, meuse)
> x1<-fit.variogram(v1, vgm(1, "Sph", 800, 1))
> G1 <- krige(zinc ~ 1, meuse, meuse.grid, x1, nmax = 30)
> gridded(G1)<-TRUE
> # GRID-2:
> v2 <- variogram(log(lead) ~ 1, meuse)
> x2<-fit.variogram(v2, vgm(.1, "Sph", 1000, .6))
> G2 <- krige(zinc ~ 1, meuse, meuse.grid, x2, nmax = 30)
> gridded(G2)<-TRUE
names(G1)
names(G2)
Note that G1 and G2 consist each of two columns, the prediction
(var1.pred) and the variance of predictions (var1.var)
> # Mooving Windows Correlation
> # Define the size of the moving window,
> dist = sqrt(2*((40*3)^2))
> # Create a list of neighbours for the distance criteria
> nb <- dnearneigh(coordinates(G1),0, dist)
> v = sapply(nb, function(i) cor(G1 at data[i,], G2 at data[i,]))
>
if you run the line above as it is that means you are calculating the
correlation of all possible combination of the two G1 columns and the
two G2 columns. That is, 1) prediction G1 - prediction G2, 2) prediction
G1 - variance G2, 3) variance G1 - prediction G2 and 4) variance G1 -
variance G2, in that order.
If you want to leave the code like this you can do then
class(v) # in this case a matrix
v = as.data.frame(t(v))
coordinates(v) = coordinates(G1)
gridded(v)=TRUE
spplot(v)
or you can calculate at once the correlation between both prediction:
v = sapply(nb, function(i) cor(G1 at data[i,1], G2 at data[i,1]))
class(v) # in this case is a numeric vector
v= as.data.frame(v)
coordinates(v) = coordinates(G1)
gridded(v) = TRUE
spplot(v)
Now you are interested in the p-value of each correlation!!!
Then I would calculate everything at one (correlation between
predictions and the p-value of each correlation)
v = sapply(nb, function(i) cor.test(G1 at data[i,1], G2 at data[i,1]))
# correlation coefficient:
c = unlist(v["estimate",])
c = as.data.frame(v)
coordinates(c) = coordinates(G1)
gridded(c) = TRUE
spplot(c)
# p-value
p = unlist(v["p.value",])
p = as.data.frame(v)
coordinates(p) = coordinates(G1)
gridded(p) = TRUE
spplot(p)
HTH,
Jaime -R.
> When I run following codes it shows error.
>
> coordinates(v) = coordinates(G1)
>
> Error in function (classes, fdef, mtable) :
> unable to find an inherited method for function "coordinates<-", for
> signature "matrix"
>
> gridded(v) = TRUE
>
> Error in `gridded<-`(`*tmp*`, value = TRUE) :
> gridded<- only works for SpatialPoints[DataFrame] or
> SpatialGrid[DataFrame]
>
>
> # Output #--------
>
> V1 V2 V3 V4 V5 V6
> V7
> 1 0.9996752 0.9996896 0.9997172 0.99971074 0.9996085 0.9996784
> 0.99972521
> 2 -0.3807374 -0.3576051 -0.2260482 0.05467783 -0.2882067 -0.2124191
> -0.08551635
> 3 -0.3943537 -0.3750882 -0.2435564 0.03511977 -0.3081105 -0.2300869
> -0.10426607
> 4 0.9999324 0.9999349 0.9999442 0.99994815 0.9999047 0.9999104
> 0.99992347
>
>
>
> Jaime R. Garcia M. wrote:
> > On Tue, 2009-12-08 at 21:54 +0100, Jaime R. Garcia M. wrote:
> >
> > > On Tue, 2009-12-08 at 14:32 -0500, Zia Ahmed wrote:
> > >
> > > > I have two grid maps (100 m cell size). I want to a create another
> > > > raster gird of correlation coefficient ( r value) of two maps using
> > > > moving windows approach. Has someone any idea how to do it in R?
> > > >
> > > > Thanks
> > > >
> > > > Zia
> > > >
> > > > _______________________________________________
> > > > R-sig-Geo mailing list
> > > > R-sig-Geo at stat.math.ethz.ch
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > > >
> > > Dear Zia,
> > >
> > > Some time ago I did exactly that and the script was based on a earlier
> > > post by Prof. Roger Bivand in this list. (Sorry I couldn't find the link
> > > anymore)
> > >
> > > Suppose "G1" and "G2" are the two grid maps.
> > > I assume they are "SpatialGridDataFrames" and have exactly the same
> > > extent and number of columns and rows. Also that you have libraries (sp)
> > > and (spdep) loaded
> > >
> > > # Define the size of the moving window, for example (7*7) so you have
> > > enough neighbors for reliable correlation estimation:
> > >
> > > dist = sqrt(2*((100*3)^2))
> > >
> > > # Create a list of neighbours for the distance criteria
> > > nb <- dnearneigh(coordinates(G1),0, dist)
> > >
> > > # And then calculate correlations
> > > v <- sapply(nb, function(i) cor(G1 at data, G2 at data))
> > >
> >
> > # Small correction
> >
> > v = sapply(nb, function(i) cor(G1 at data[i,], G2 at data[i,]))
> >
> > # And in this case the result should be a numeric vector
> >
> > coordinates(v) = coordinates(G1)
> > gridded(v) = TRUE
> > fullgrid(v) = TRUE
> > _______________________________________________
> >
> > > R-sig-Geo mailing list
> > > R-sig-Geo at stat.math.ethz.ch
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> >
> >
> >
More information about the R-sig-Geo
mailing list