Geographically weighted regression
Description
The function implements the basic geographically weighted regression approach to exploring spatial non-stationarity for given global bandwidth and chosen weighting scheme.
Usage
gwr(formula, data=list(), coords, bandwidth, gweight=gwr.Gauss,
adapt=NULL, hatmatrix = FALSE, fit.points, longlat=NULL,
se.fit=FALSE, weights, cl=NULL, predictions = FALSE,
fittedGWRobject = NULL, se.fit.CCT = TRUE, use_snow=FALSE)
## S3 method for class 'gwr'
print(x, ...)
Arguments
formula |
regression model formula as in lm |
data |
model data frame, or SpatialPointsDataFrame or SpatialPolygonsDataFrame as defined in package sp |
coords |
matrix of coordinates of points representing the spatial positions of the observations; may be omitted if the object passed through the data argument is from package sp |
bandwidth |
bandwidth used in the weighting function, possibly calculated by gwr.sel |
gweight |
geographical weighting function, at present gwr.Gauss() default, or gwr.gauss(), the previous default orgwr.bisquare() |
adapt |
either NULL (default) or a proportion between 0 and 1 of observations to include in weighting scheme (k-nearest neighbours) |
hatmatrix |
if TRUE, return the hatmatrix as a component of the result, ignored if fit.points given |
fit.points |
an object containing the coordinates of fit points; often an object from package sp; if missing, the coordinates given through the data argument object, or the coords argument are used |
longlat |
TRUE if point coordinates are longitude-latitude decimal degrees, in which case distances are measured in kilometers; if x is a SpatialPoints object, the value is taken from the object itself |
se.fit |
if TRUE, return local coefficient standard errors – if hatmatrix is TRUE and no fit.points are given, two effective degrees of freedom sigmas will be used to generate alternative coefficient standard errors |
weights |
case weights used as in weighted least squares, beware of scaling issues, probably unsafe |
cl |
if NULL, ignored, otherwise cl must be an object describing a “cluster” created using makeCluster in the parallelpackage. The cluster will then be used to hand off the calculation of local coefficients to cluster nodes, if fit points have been given as an argument, and hatmatrix=FALSE |
predictions |
default FALSE; if TRUE and no fit points given, return GW fitted values at data points, if fit points given and are a Spatial*DataFrame object containing the RHS variables in the formula, return GW predictions at the fit points |
fittedGWRobject |
a fitted gwr object with a hatmatrix (optional), if given, and if fit.points are given and if se.fit is TRUE, two effective degrees of freedom sigmas will be used to generate alternative coefficient standard errors |
se.fit.CCT |
default TRUE, compute local coefficient standard errors using formula (2.14), p. 55, in the GWR book |
use_snow |
default FALSE, use parallel; if “cl” is given, and use_snow=TRUE, use snow instead of parallel to access more types of cluster |
x |
an object of class “gwr” returned by the gwr function |
... |
arguments to be passed to other functions |
Details
The function applies the weighting function in turn to each of the observations, or fit points if given, calculating a weighted regression for each point. The results may be explored to see if coefficient values vary over space. The local coefficient estimates may be made on a multi-node cluster using the cl argument to pass through a parallel cluster. The function will then divide the fit points (which must be given separately) between the clusters for fitting. Note that each node will need to have the “spgwr” package present, so initiating by clusterEvalQ(cl, library(spgwr)) may save a little time per node. The function clears the global environment on the node of objects sent. Using two nodes reduces timings to a little over half the time for a single node.
The section of the examples code now includes two simulation scenarios, showing how important it is to check that mapped pattern in local coefficients is actually there, rather than being an artefact.
Value
A list of class “gwr”:
SDF |
a SpatialPointsDataFrame (may be gridded) or SpatialPolygonsDataFrame object (see package “sp”) with fit.points, weights, GWR coefficient estimates, R-squared, and coefficient standard errors in its “data” slot. |
lhat |
Leung et al. L matrix |
lm |
Ordinary least squares global regression on the same model formula, as returned by lm.wfit(). |
bandwidth |
the bandwidth used. |
this.call |
the function call used. |
Author(s)
Roger Bivand Roger.Bivand@nhh.no
References
Fotheringham, A.S., Brunsdon, C., and Charlton, M.E., 2002, Geographically Weighted Regression, Chichester: Wiley; Paez A, Farber S, Wheeler D, 2011, “A simulation-based study of geographically weighted regression as a method for investigating spatially varying relationships”, Environment and Planning A 43(12) 2992-3010; http://www.nuim.ie/ncg/GWR/index.htm
See Also
gwr.sel, gwr.gauss, gwr.bisquare
Examples
data(columbus)
col.lm <- lm(crime ~ income + housing, data=columbus)
summary(col.lm)
col.bw <- gwr.sel(crime ~ income + housing, data=columbus,
coords=cbind(columbus$x, columbus$y))
col.gauss <- gwr(crime ~ income + housing, data=columbus,
coords=cbind(columbus$x, columbus$y), bandwidth=col.bw, hatmatrix=TRUE)
col.gauss
col.d <- gwr.sel(crime ~ income + housing, data=columbus,
coords=cbind(columbus$x, columbus$y), gweight=gwr.bisquare)
col.bisq <- gwr(crime ~ income + housing, data=columbus,
coords=cbind(columbus$x, columbus$y), bandwidth=col.d,
gweight=gwr.bisquare, hatmatrix=TRUE)
col.bisq
data(georgia)
g.adapt.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
PctPov + PctBlack, data=gSRDF, adapt=TRUE)
res.adpt <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
PctBlack, data=gSRDF, adapt=g.adapt.gauss)
res.adpt
pairs(as(res.adpt$SDF, "data.frame")[,2:8], pch=".")
brks <- c(-0.25, 0, 0.01, 0.025, 0.075)
cols <- grey(5:2/6)
plot(res.adpt$SDF, col=cols[findInterval(res.adpt$SDF$PctBlack, brks,
all.inside=TRUE)])
# simulation scenario with patterned dependent variable
set.seed(1)
X0 <- runif(nrow(gSRDF)*3)
X1 <- matrix(sample(X0), ncol=3)
X1 <- prcomp(X1, center=FALSE, scale.=FALSE)$x
gSRDF$X1 <- X1[,1]
gSRDF$X2 <- X1[,2]
gSRDF$X3 <- X1[,3]
bw <- gwr.sel(PctBach ~ X1 + X2 + X3, data=gSRDF, verbose=FALSE)
out <- gwr(PctBach ~ X1 + X2 + X3, data=gSRDF, bandwidth=bw, hatmatrix=TRUE)
out
spplot(gSRDF, "PctBach", col.regions=grey.colors(20))
spplot(gSRDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
# pattern in the local coefficients
spplot(out$SDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
# but no "significant" pattern
spplot(out$SDF, c("X1_se", "X2_se", "X3_se"), col.regions=grey.colors(20))
out$SDF$X1_t <- out$SDF$X1/out$SDF$X1_se
out$SDF$X2_t <- out$SDF$X2/out$SDF$X2_se
out$SDF$X3_t <- out$SDF$X3/out$SDF$X3_se
spplot(out$SDF, c("X1_t", "X2_t", "X3_t"), col.regions=grey.colors(20))
# simulation scenario with random dependent variable
yrn <- rnorm(nrow(gSRDF))
gSRDF$yrn <- sample(yrn)
bw <- gwr.sel(yrn ~ X1 + X2 + X3, data=gSRDF, verbose=FALSE)
# bandwidth selection maxes out at 620 km, equal to upper bound
# of line search
out <- gwr(yrn ~ X1 + X2 + X3, data=gSRDF, bandwidth=bw, hatmatrix=TRUE)
out
spplot(gSRDF, "yrn", col.regions=grey.colors(20))
spplot(gSRDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
# pattern in the local coefficients
spplot(out$SDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
# but no "significant" pattern
spplot(out$SDF, c("X1_se", "X2_se", "X3_se"), col.regions=grey.colors(20))
out$SDF$X1_t <- out$SDF$X1/out$SDF$X1_se
out$SDF$X2_t <- out$SDF$X2/out$SDF$X2_se
out$SDF$X3_t <- out$SDF$X3/out$SDF$X3_se
spplot(out$SDF, c("X1_t", "X2_t", "X3_t"), col.regions=grey.colors(20))
# end of simulations
data(meuse)
coordinates(meuse) <- c("x", "y")
meuse$ffreq <- factor(meuse$ffreq)
data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid$ffreq <- factor(meuse.grid$ffreq)
gridded(meuse.grid) <- TRUE
xx <- gwr(cadmium ~ dist, meuse, bandwidth = 228, hatmatrix=TRUE)
xx
x <- gwr(cadmium ~ dist, meuse, bandwidth = 228, fit.points = meuse.grid,
predict=TRUE, se.fit=TRUE, fittedGWRobject=xx)
x
spplot(x$SDF, "pred")
spplot(x$SDF, "pred.se")
## Not run:
g.bw.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
PctPov + PctBlack, data=gSRDF)
res.bw <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
PctBlack, data=gSRDF, bandwidth=g.bw.gauss)
res.bw
pairs(as(res.bw$SDF, "data.frame")[,2:8], pch=".")
plot(res.bw$SDF, col=cols[findInterval(res.bw$SDF$PctBlack, brks,
all.inside=TRUE)])
g.bw.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
PctPov + PctBlack, data=gSRDF, longlat=TRUE)
data(gSRouter)
SG <- GE_SpatialGrid(gSRouter, maxPixels = 100)
SPxMASK0 <- over(SG$SG, gSRouter)
SGDF <- SpatialGridDataFrame(slot(SG$SG, "grid"),
data=data.frame(SPxMASK0=SPxMASK0),
proj4string=CRS(proj4string(gSRouter)))
SPxDF <- as(SGDF, "SpatialPixelsDataFrame")
res.bw <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
PctBlack, data=gSRDF, bandwidth=g.bw.gauss, fit.points=SPxDF,
longlat=TRUE)
res.bw
res.bw$timings
spplot(res.bw$SDF, "PctBlack")
cl <- makeCluster(detectCores())
res.bwc <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
PctBlack, data=gSRDF, bandwidth=g.bw.gauss, fit.points=SPxDF,
longlat=TRUE, cl=cl)
res.bwc
res.bwc$timings
stopCluster(cl)
## End(Not run)
Results
R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(spgwr)
Loading required package: sp
Loading required package: maptools
Loading required package: foreign
Loading required package: grid
Loading required package: lattice
Checking rgeos availability: TRUE
NOTE: This package does not constitute approval of GWR
as a method of spatial analysis; see example(gwr)
> png(filename="/home/oogasawa/snapshot/RGM3/R_CC/result/spgwr/gwr.Rd_%03d_medium.png", width=480, height=480)
> ### Name: gwr
> ### Title: Geographically weighted regression
> ### Aliases: gwr print.gwr
> ### Keywords: spatial
>
> ### ** Examples
>
> data(columbus)
> col.lm <- lm(crime ~ income + housing, data=columbus)
> summary(col.lm)
Call:
lm(formula = crime ~ income + housing, data = columbus)
Residuals:
Min 1Q Median 3Q Max
-34.418 -6.388 -1.580 9.052 28.649
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.6189 4.7355 14.490 < 2e-16 ***
income -1.5973 0.3341 -4.780 1.83e-05 ***
housing -0.2739 0.1032 -2.654 0.0109 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.43 on 46 degrees of freedom
Multiple R-squared: 0.5524, Adjusted R-squared: 0.5329
F-statistic: 28.39 on 2 and 46 DF, p-value: 9.341e-09
> col.bw <- gwr.sel(crime ~ income + housing, data=columbus,
+ coords=cbind(columbus$x, columbus$y))
Bandwidth: 12.65221 CV score: 7432.171
Bandwidth: 20.45127 CV score: 7462.669
Bandwidth: 7.832129 CV score: 7323.494
Bandwidth: 4.853154 CV score: 7307.484
Bandwidth: 5.122256 CV score: 7322.64
Bandwidth: 3.012046 CV score: 6461.665
Bandwidth: 1.874178 CV score: 6473.183
Bandwidth: 2.475225 CV score: 6109.765
Bandwidth: 2.447682 CV score: 6098.241
Bandwidth: 2.228623 CV score: 6063.971
Bandwidth: 2.264529 CV score: 6060.644
Bandwidth: 2.280636 CV score: 6060.52
Bandwidth: 2.274941 CV score: 6060.473
Bandwidth: 2.275073 CV score: 6060.473
Bandwidth: 2.275032 CV score: 6060.473
Bandwidth: 2.274991 CV score: 6060.473
Bandwidth: 2.275032 CV score: 6060.473
> col.gauss <- gwr(crime ~ income + housing, data=columbus,
+ coords=cbind(columbus$x, columbus$y), bandwidth=col.bw, hatmatrix=TRUE)
> col.gauss
Call:
gwr(formula = crime ~ income + housing, data = columbus, coords = cbind(columbus$x,
columbus$y), bandwidth = col.bw, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 2.275032
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 23.23000 54.13000 63.90000 68.76000 80.90000 68.6189
income -3.13100 -1.91300 -0.98440 -0.36860 1.29100 -1.5973
housing -1.05300 -0.37670 -0.09739 0.03006 0.79460 -0.2739
Number of data points: 49
Effective number of parameters (residual: 2traceS - traceS'S): 29.61664
Effective degrees of freedom (residual: 2traceS - traceS'S): 19.38336
Sigma (residual: 2traceS - traceS'S): 8.027396
Effective number of parameters (model: traceS): 23.92826
Effective degrees of freedom (model: traceS): 25.07174
Sigma (model: traceS): 7.058251
Sigma (ML): 5.048836
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 403.6193
AIC (GWR p. 96, eq. 4.22): 321.6617
Residual sum of squares: 1249.046
Quasi-global R2: 0.9070521
> col.d <- gwr.sel(crime ~ income + housing, data=columbus,
+ coords=cbind(columbus$x, columbus$y), gweight=gwr.bisquare)
Bandwidth: 12.65221 CV score: 8180.49
Bandwidth: 20.45127 CV score: 7552.802
Bandwidth: 25.27135 CV score: 7508.191
Bandwidth: 23.68116 CV score: 7519.828
Bandwidth: 28.25033 CV score: 7491.815
Bandwidth: 30.09144 CV score: 7486.638
Bandwidth: 31.69367 CV score: 7483.628
Bandwidth: 31.08167 CV score: 7484.671
Bandwidth: 32.21954 CV score: 7482.811
Bandwidth: 32.54454 CV score: 7482.336
Bandwidth: 32.74541 CV score: 7482.053
Bandwidth: 32.86955 CV score: 7481.881
Bandwidth: 32.94627 CV score: 7481.777
Bandwidth: 32.99369 CV score: 7481.712
Bandwidth: 33.023 CV score: 7481.673
Bandwidth: 33.04111 CV score: 7481.649
Bandwidth: 33.0523 CV score: 7481.634
Bandwidth: 33.05922 CV score: 7481.624
Bandwidth: 33.06349 CV score: 7481.619
Bandwidth: 33.06614 CV score: 7481.615
Bandwidth: 33.06777 CV score: 7481.613
Bandwidth: 33.06878 CV score: 7481.612
Bandwidth: 33.0694 CV score: 7481.611
Bandwidth: 33.06979 CV score: 7481.61
Bandwidth: 33.07003 CV score: 7481.61
Bandwidth: 33.07017 CV score: 7481.61
Bandwidth: 33.07027 CV score: 7481.61
Bandwidth: 33.07032 CV score: 7481.609
Bandwidth: 33.07036 CV score: 7481.609
Bandwidth: 33.07036 CV score: 7481.609
Warning message:
In gwr.sel(crime ~ income + housing, data = columbus, coords = cbind(columbus$x, :
Bandwidth converged to upper bound:33.0704127582345
> col.bisq <- gwr(crime ~ income + housing, data=columbus,
+ coords=cbind(columbus$x, columbus$y), bandwidth=col.d,
+ gweight=gwr.bisquare, hatmatrix=TRUE)
> col.bisq
Call:
gwr(formula = crime ~ income + housing, data = columbus, coords = cbind(columbus$x,
columbus$y), bandwidth = col.d, gweight = gwr.bisquare, hatmatrix = TRUE)
Kernel function: gwr.bisquare
Fixed bandwidth: 33.07036
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 68.3100 68.9600 69.2600 69.5700 71.6400 68.6189
income -1.7270 -1.6480 -1.6140 -1.5690 -1.4680 -1.5973
housing -0.3331 -0.3052 -0.2806 -0.2535 -0.1883 -0.2739
Number of data points: 49
Effective number of parameters (residual: 2traceS - traceS'S): 4.604881
Effective degrees of freedom (residual: 2traceS - traceS'S): 44.39512
Sigma (residual: 2traceS - traceS'S): 11.37646
Effective number of parameters (model: traceS): 3.917653
Effective degrees of freedom (model: traceS): 45.08235
Sigma (model: traceS): 11.28942
Sigma (ML): 10.82871
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 383.698
AIC (GWR p. 96, eq. 4.22): 376.4294
Residual sum of squares: 5745.791
Quasi-global R2: 0.5724264
> data(georgia)
> g.adapt.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
+ PctPov + PctBlack, data=gSRDF, adapt=TRUE)
Adaptive q: 0.381966 CV score: 2048.736
Adaptive q: 0.618034 CV score: 2028.616
Adaptive q: 0.763932 CV score: 2041.588
Adaptive q: 0.5934745 CV score: 2027.578
Adaptive q: 0.5284474 CV score: 2026.646
Adaptive q: 0.537984 CV score: 2026.588
Adaptive q: 0.5414424 CV score: 2026.575
Adaptive q: 0.5613169 CV score: 2026.795
Adaptive q: 0.5507625 CV score: 2026.522
Adaptive q: 0.5547939 CV score: 2026.587
Adaptive q: 0.5478453 CV score: 2026.5
Adaptive q: 0.5474828 CV score: 2026.5
Adaptive q: 0.5474097 CV score: 2026.5
Adaptive q: 0.547369 CV score: 2026.5
Adaptive q: 0.5474097 CV score: 2026.5
> res.adpt <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
+ PctBlack, data=gSRDF, adapt=g.adapt.gauss)
> res.adpt
Call:
gwr(formula = PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
PctPov + PctBlack, data = gSRDF, adapt = g.adapt.gauss)
Kernel function: gwr.Gauss
Adaptive quantile: 0.5474097 (about 87 of 159)
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 1.358e+01 1.409e+01 1.522e+01 1.570e+01 1.597e+01 14.7766
TotPop90 1.870e-05 2.053e-05 2.279e-05 2.482e-05 2.586e-05 0.0000
PctRural -5.253e-02 -4.829e-02 -4.191e-02 -3.600e-02 -3.206e-02 -0.0439
PctEld -1.586e-01 -1.335e-01 -1.146e-01 -9.678e-02 -5.908e-02 -0.0619
PctFB 8.164e-01 9.940e-01 1.354e+00 1.742e+00 1.958e+00 1.2562
PctPov -1.838e-01 -1.643e-01 -1.256e-01 -9.503e-02 -6.538e-02 -0.1554
PctBlack -1.194e-02 -1.639e-04 1.160e-02 2.760e-02 3.895e-02 0.0219
> pairs(as(res.adpt$SDF, "data.frame")[,2:8], pch=".")
> brks <- c(-0.25, 0, 0.01, 0.025, 0.075)
> cols <- grey(5:2/6)
> plot(res.adpt$SDF, col=cols[findInterval(res.adpt$SDF$PctBlack, brks,
+ all.inside=TRUE)])
> ## No test:
> # simulation scenario with patterned dependent variable
> set.seed(1)
> X0 <- runif(nrow(gSRDF)*3)
> X1 <- matrix(sample(X0), ncol=3)
> X1 <- prcomp(X1, center=FALSE, scale.=FALSE)$x
> gSRDF$X1 <- X1[,1]
> gSRDF$X2 <- X1[,2]
> gSRDF$X3 <- X1[,3]
> bw <- gwr.sel(PctBach ~ X1 + X2 + X3, data=gSRDF, verbose=FALSE)
> out <- gwr(PctBach ~ X1 + X2 + X3, data=gSRDF, bandwidth=bw, hatmatrix=TRUE)
> out
Call:
gwr(formula = PctBach ~ X1 + X2 + X3, data = gSRDF, bandwidth = bw,
hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 110.5668
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 9.6300 11.3700 12.1300 12.4600 16.5800 11.9185
X1 -3.0180 -0.9243 0.7687 2.9440 8.0570 1.1207
X2 -1.9000 -0.5299 0.3944 1.2410 2.6240 0.1101
X3 -1.8910 -1.2040 -0.8930 -0.6000 2.7620 -0.6651
Number of data points: 159
Effective number of parameters (residual: 2traceS - traceS'S): 18.6728
Effective degrees of freedom (residual: 2traceS - traceS'S): 140.3272
Sigma (residual: 2traceS - traceS'S): 5.61679
Effective number of parameters (model: traceS): 13.51085
Effective degrees of freedom (model: traceS): 145.4892
Sigma (model: traceS): 5.516248
Sigma (ML): 5.276677
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1012.31
AIC (GWR p. 96, eq. 4.22): 993.6616
Residual sum of squares: 4427.088
Quasi-global R2: 0.1366961
> spplot(gSRDF, "PctBach", col.regions=grey.colors(20))
> spplot(gSRDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
> # pattern in the local coefficients
> spplot(out$SDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
> # but no "significant" pattern
> spplot(out$SDF, c("X1_se", "X2_se", "X3_se"), col.regions=grey.colors(20))
> out$SDF$X1_t <- out$SDF$X1/out$SDF$X1_se
> out$SDF$X2_t <- out$SDF$X2/out$SDF$X2_se
> out$SDF$X3_t <- out$SDF$X3/out$SDF$X3_se
> spplot(out$SDF, c("X1_t", "X2_t", "X3_t"), col.regions=grey.colors(20))
> # simulation scenario with random dependent variable
> yrn <- rnorm(nrow(gSRDF))
> gSRDF$yrn <- sample(yrn)
> bw <- gwr.sel(yrn ~ X1 + X2 + X3, data=gSRDF, verbose=FALSE)
Warning message:
In gwr.sel(yrn ~ X1 + X2 + X3, data = gSRDF, verbose = FALSE) :
Bandwidth converged to upper bound:620.141464185404
> # bandwidth selection maxes out at 620 km, equal to upper bound
> # of line search
> out <- gwr(yrn ~ X1 + X2 + X3, data=gSRDF, bandwidth=bw, hatmatrix=TRUE)
> out
Call:
gwr(formula = yrn ~ X1 + X2 + X3, data = gSRDF, bandwidth = bw,
hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 620.1414
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. -0.062010 -0.040040 -0.023470 -0.008171 0.011440 -0.0206
X1 -0.037920 -0.013900 0.001424 0.016200 0.037410 0.0038
X2 -0.007911 -0.002395 0.002073 0.007433 0.015880 0.0111
X3 -0.019610 0.001542 0.012040 0.025520 0.050420 0.0087
Number of data points: 159
Effective number of parameters (residual: 2traceS - traceS'S): 4.54916
Effective degrees of freedom (residual: 2traceS - traceS'S): 154.4508
Sigma (residual: 2traceS - traceS'S): 1.079237
Effective number of parameters (model: traceS): 4.281352
Effective degrees of freedom (model: traceS): 154.7186
Sigma (model: traceS): 1.078303
Sigma (ML): 1.063686
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 481.8531
AIC (GWR p. 96, eq. 4.22): 475.1373
Residual sum of squares: 179.8971
Quasi-global R2: 0.003021284
> spplot(gSRDF, "yrn", col.regions=grey.colors(20))
> spplot(gSRDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
> # pattern in the local coefficients
> spplot(out$SDF, c("X1", "X2", "X3"), col.regions=grey.colors(20))
> # but no "significant" pattern
> spplot(out$SDF, c("X1_se", "X2_se", "X3_se"), col.regions=grey.colors(20))
> out$SDF$X1_t <- out$SDF$X1/out$SDF$X1_se
> out$SDF$X2_t <- out$SDF$X2/out$SDF$X2_se
> out$SDF$X3_t <- out$SDF$X3/out$SDF$X3_se
> spplot(out$SDF, c("X1_t", "X2_t", "X3_t"), col.regions=grey.colors(20))
> # end of simulations
> ## End(No test)
> ## No test:
> data(meuse)
> coordinates(meuse) <- c("x", "y")
> meuse$ffreq <- factor(meuse$ffreq)
> data(meuse.grid)
> coordinates(meuse.grid) <- c("x", "y")
> meuse.grid$ffreq <- factor(meuse.grid$ffreq)
> gridded(meuse.grid) <- TRUE
> xx <- gwr(cadmium ~ dist, meuse, bandwidth = 228, hatmatrix=TRUE)
> xx
Call:
gwr(formula = cadmium ~ dist, data = meuse, bandwidth = 228,
hatmatrix = TRUE)
Kernel function: gwr.Gauss
Fixed bandwidth: 228
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 1.290 3.750 5.646 7.260 12.950 5.8795
dist -36.230 -20.570 -12.510 -6.306 12.200 -10.9730
Number of data points: 155
Effective number of parameters (residual: 2traceS - traceS'S): 33.66762
Effective degrees of freedom (residual: 2traceS - traceS'S): 121.3324
Sigma (residual: 2traceS - traceS'S): 2.172593
Effective number of parameters (model: traceS): 25.01465
Effective degrees of freedom (model: traceS): 129.9854
Sigma (model: traceS): 2.099034
Sigma (ML): 1.92221
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 705.4599
AIC (GWR p. 96, eq. 4.22): 667.4631
Residual sum of squares: 572.7084
Quasi-global R2: 0.7004953
> x <- gwr(cadmium ~ dist, meuse, bandwidth = 228, fit.points = meuse.grid,
+ predict=TRUE, se.fit=TRUE, fittedGWRobject=xx)
> x
Call:
gwr(formula = cadmium ~ dist, data = meuse, bandwidth = 228,
fit.points = meuse.grid, se.fit = TRUE, predictions = TRUE,
fittedGWRobject = xx)
Kernel function: gwr.Gauss
Fixed bandwidth: 228
Summary of GWR coefficient estimates:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. 1.267 3.176 4.742 6.763 14.150 5.8795
dist -39.280 -16.460 -8.957 -4.693 21.180 -10.9730
> spplot(x$SDF, "pred")
> spplot(x$SDF, "pred.se")
> ## End(No test)
> ## Not run:
> ##D g.bw.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
> ##D PctPov + PctBlack, data=gSRDF)
> ##D res.bw <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
> ##D PctBlack, data=gSRDF, bandwidth=g.bw.gauss)
> ##D res.bw
> ##D pairs(as(res.bw$SDF, "data.frame")[,2:8], pch=".")
> ##D plot(res.bw$SDF, col=cols[findInterval(res.bw$SDF$PctBlack, brks,
> ##D all.inside=TRUE)])
> ##D g.bw.gauss <- gwr.sel(PctBach ~ TotPop90 + PctRural + PctEld + PctFB +
> ##D PctPov + PctBlack, data=gSRDF, longlat=TRUE)
> ##D data(gSRouter)
> ##D SG <- GE_SpatialGrid(gSRouter, maxPixels = 100)
> ##D SPxMASK0 <- over(SG$SG, gSRouter)
> ##D SGDF <- SpatialGridDataFrame(slot(SG$SG, "grid"),
> ##D data=data.frame(SPxMASK0=SPxMASK0),
> ##D proj4string=CRS(proj4string(gSRouter)))
> ##D SPxDF <- as(SGDF, "SpatialPixelsDataFrame")
> ##D res.bw <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
> ##D PctBlack, data=gSRDF, bandwidth=g.bw.gauss, fit.points=SPxDF,
> ##D longlat=TRUE)
> ##D res.bw
> ##D res.bw$timings
> ##D spplot(res.bw$SDF, "PctBlack")
> ##D cl <- makeCluster(detectCores())
> ##D res.bwc <- gwr(PctBach ~ TotPop90 + PctRural + PctEld + PctFB + PctPov +
> ##D PctBlack, data=gSRDF, bandwidth=g.bw.gauss, fit.points=SPxDF,
> ##D longlat=TRUE, cl=cl)
> ##D res.bwc
> ##D res.bwc$timings
> ##D stopCluster(cl)
> ## End(Not run)
>
>
>
>
>
> dev.off()
null device
1
>
|
