[R-sig-Geo] spgwr and selecting aic for bandwidth configuration

Roger Bivand Roger.Bivand at nhh.no
Thu May 6 20:00:38 CEST 2010


On Sat, 10 Apr 2010, Christopher Cardinal wrote:

> I am trying to use gwr.sel to apply the AIC method for bandwidth selection
> in SPGWR. R version is 2.10.1 (used for compatibility with Rpy--however, the
> following problem occurs directly from within R). SPGWR version is 0.6-4.
> First I run global regression, then gwr()with default CV method for
> bandwidth selection with no problems. Then I run "aic" bandwidth selection
> method, which outputs bandwidth results. Then I run gwr(), which produces no
> output. Perhaps I'm using the wrong command? I can't find many examples of
> gwr() implemented with "aic". Please share your ideas if you have any.

Sorry not to have replied earlier. There is no summary() method for "gwr" 
objects, only a print() method. The output of gwr.sel() is a single 
number. It is impossible to say what is in data_bw_aic - it should be 
similar to the last number shown in the verbose output - but the 
optimisation in gwr.sel() does not seem to have converged - I don't think 
that you have paid attention to error messages, or at least not shown them 
in your posting. Without a reproducible example, it is impossible to help. 
Try to reproduce your problem with the columbus or georgia data sets 
shipped with the package.

Hope this helps,

Roger


> Thanks.
>
>> library(spgwr)
> Loading required package: sp
> Loading required package: maptools
> Loading required package: foreign
> Loading required package: lattice
> NOTE: default kernel and CV criteria changed
> see help pages for details
>> library(foreign)
>> data <-
> read.dbf("C:/GIS/Hunter/Thesis/GISData/Parcel/Albany/AlbCnty/Polygon2004/test2.dbf")
>> data <- subset(data, GOOGDIST > 0)
>> data <- subset(data, ACRES > .5)
>> data_lm <- lm(SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH, data =
> data)
>> data_lm_sum <- summary(data_lm)
>> tmp<-print(data_lm_sum)
>
> Call:
> lm(formula = SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH,
>    data = data)
>
> Residuals:
>    Min      1Q  Median      3Q     Max
> -890883  -92357  -66457   76176 1550181
>
> Coefficients:
>              Estimate Std. Error t value Pr(>|t|)
> (Intercept)  8.997e+04  2.566e+04   3.507 0.000483 ***
> ACRES       -1.999e+03  3.915e+02  -5.106 4.26e-07 ***
> CUR_TOT_A    1.993e+00  9.443e-02  21.104  < 2e-16 ***
> HHINC       -7.636e-02  5.363e-01  -0.142 0.886810
> PERGTBACH    9.287e+03  5.020e+04   0.185 0.853272
> ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>
> Residual standard error: 186700 on 689 degrees of freedom
> Multiple R-squared: 0.4061,     Adjusted R-squared: 0.4027
> F-statistic: 117.8 on 4 and 689 DF,  p-value: < 2.2e-16
>
>> data_bw <- gwr.sel(SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH ,
> data = data, coords = cbind(data$X, data$Y))
> Bandwidth: 79378.23 CV score: 2.508302e+13
> Bandwidth: 128308.4 CV score: 2.549703e+13
> Bandwidth: 49137.69 CV score: 2.444211e+13
> Bandwidth: 30448.02 CV score: 2.369592e+13
> Bandwidth: 18897.16 CV score: 2.35517e+13
> Bandwidth: 17792.42 CV score: 2.365362e+13
> Bandwidth: 23918.27 CV score: 2.340727e+13
> Bandwidth: 26412.41 CV score: 2.348778e+13
> Bandwidth: 23178.31 CV score: 2.339931e+13
> Bandwidth: 22740.56 CV score: 2.339894e+13
> Bandwidth: 22908.52 CV score: 2.339868e+13
> Bandwidth: 22910.54 CV score: 2.339868e+13
> Bandwidth: 22909.77 CV score: 2.339868e+13
> Bandwidth: 22909.78 CV score: 2.339868e+13
> Bandwidth: 22909.78 CV score: 2.339868e+13
> Bandwidth: 22909.78 CV score: 2.339868e+13
> Bandwidth: 22909.78 CV score: 2.339868e+13
>> data_gauss <- gwr(SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH, data
> = data, coords = cbind(data$X, data$Y), bandwidth = data_bw, hatmatrix =
> TRUE)
>> data_gauss_sum <- summary(data_gauss)
>> print (data_gauss_sum)
>          Length Class                  Mode
> SDF            1 SpatialPointsDataFrame S4
> lhat      481636 -none-                 numeric
> lm            11 -none-                 list
> results       14 -none-                 list
> bandwidth      1 -none-                 numeric
> adapt          0 -none-                 NULL
> hatmatrix      1 -none-                 logical
> gweight        1 -none-                 character
> this.call      6 -none-                 call
>> print (data_gauss)
> Call:
> gwr(formula = SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH,
>    data = data, coords = cbind(data$X, data$Y), bandwidth = data_bw,
>    hatmatrix = TRUE)
> Kernel function: gwr.Gauss
> Fixed bandwidth: 22909.78
> Summary of GWR coefficient estimates:
>                   Min.    1st Qu.     Median    3rd Qu.       Max.
> Global
> X.Intercept. -2.883e+05 -4.166e+04  1.135e+05  1.655e+05  2.834e+05
> 89969.4089
> ACRES        -3.259e+03 -2.678e+03 -1.572e+03  9.080e+02  8.611e+03
> -1998.9124
> CUR_TOT_A     2.233e-02  1.261e+00  1.628e+00  2.049e+00  2.380e+00
> 1.9928
> HHINC        -5.563e+00 -1.219e+00  3.129e-01  2.753e+00  8.023e+00
> -0.0764
> PERGTBACH    -3.023e+05 -9.237e+04 -5.303e+04 -7.264e+03  1.014e+05
> 9287.3433
> Number of data points: 694
> Effective number of parameters (residual: 2traceS - traceS'S): 29.87522
> Effective degrees of freedom (residual: 2traceS - traceS'S): 664.1248
> Sigma (residual: 2traceS - traceS'S): 170145.2
> Effective number of parameters (model: traceS): 21.99972
> Effective degrees of freedom (model: traceS): 672.0003
> Sigma (model: traceS): 169145.3
> Sigma (ML): 166442.7
> AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 18704.23
> AIC (GWR p. 96, eq. 4.22): 18678.59
> Residual sum of squares: 1.922601e+13
>
>> data_bw_aic <- gwr.sel(SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH
> , data = data, coords = cbind(data$X, data$Y), method="aic")
> Bandwidth: 79378.23 AIC: 18795.23
> Bandwidth: 128308.4 AIC: 18811.03
> Bandwidth: 49137.69 AIC: 18768.93
> Bandwidth: 30448.02 AIC: 18733.32
> Bandwidth: 18897.16 AIC: 18678.52
> Bandwidth: 11758.34 AIC: 18593.90
>> data_bw_aic_sum <- summary(data_bw_aic)
>> data_bw_aic_sum
>> print(data_bw_aic)
>> data_gauss_aic <- gwr(SALE_PRICE ~ ACRES + CUR_TOT_A + HHINC + PERGTBACH,
> data = data, coords = cbind(data$X, data$Y), bandwidth = data_bw_aic,
> hatmatrix = TRUE)
>> print(data_gauss_aic)
>> summary(data_gauss_aic)
>>
>
>
> data looks like this:
>        UID YEAR        X       Y      LAT      LONG       ACRES SLOPE
> CUR_TOT_A     FLOODPER
> 1      6482 2004 698260.3 1446285 42.80053 -73.73233   0.7309970     1
> 40000 0.0000000000
> 5      8554 2004 695351.2 1438002 42.77801 -73.74342   0.8890091     1
> 7400 0.0000000000
> 7     15413 2004 696033.0 1431745 42.76088 -73.74108   0.5565335     1
> 50000 0.0000000000
> 10     8220 2004 680872.2 1435917 42.77268 -73.79744   0.8693349     1
> 27200 0.0000000000
> 11     8221 2004 681023.1 1436053 42.77303 -73.79688   0.8018674     1
> 25900 0.0000000000
> 12     8792 2004 691488.3 1439555 42.78245 -73.75765   0.8236566     1
> 42000 0.0000000000
> 13     8793 2004 691556.9 1439692 42.78277 -73.75746   0.8610513     1
> 42000 0.0000000000
> 14    14411 2004 682812.3 1426310 42.74625 -73.79048   1.9767706     1
> 150000 0.0000000000
> 15    14656 2004 688672.5 1432607 42.76348 -73.76833   3.5036509     1
> 65700 0.0000000000
> 16    15103 2004 686513.4 1428931 42.75343 -73.77653   4.5957896     1
> 352800 0.0000000000
> 19    89916 2004 702825.7 1439163 42.78107 -73.71557   5.4003465     1
> 23200 0.0000000000
> 22    93898 2004 704736.3 1432229 42.76197 -73.70865   0.9085830     1
> 8400 0.0000000000
> 25    10109 2004 657502.9 1427039 42.74891 -73.88458   3.1770620     1
> 53600 0.0000000000
> 26    10296 2004 657568.5 1427583 42.75031 -73.88438   0.9898759     1
> 39200 0.5419490988
>
>                                                         GOOGADDR GOOGDIST
> PROP_CLASS SALE_PRICE HHINC
> 1                       1251 New Loudon Rd, Cohoes, NY 12047, USA
> 21597        311      42500 25964
> 5                             22 Meadow St, Cohoes, NY 12047, USA
> 19195        311      31500 25964
> 7                            44 Graffin Dr, Latham, NY 12110, USA
> 18764        311     396500 25964
> 10                      4119-4195 River Rd, Latham, NY 12110, USA
> 18941        311      71454 16222
> 11                      4119-4195 River Rd, Latham, NY 12110, USA
> 18941        311      71454 16222
> 12                   112 Dunsbach Ferry Rd, Cohoes, NY 12047, USA
> 19164        311      65000 16222
> 13                   114 Dunsbach Ferry Rd, Cohoes, NY 12047, USA
> 19208        311      68000 16222
> 14                             155 Wade Rd, Latham, NY 12110, USA
> 13289        330     100000 16222
> 15                      100 Sparrowbush Rd, Latham, NY 12110, USA
> 16458        311     275000 16222
> 16                           4 Stanley Cir, Latham, NY 12110, USA
> 15605        330     307000 16222
> 19                            28 Edward St, Cohoes, NY 12047, USA
> 21813        311     114685 51214
> 22                            8 Oxford Cir, Cohoes, NY 12047, USA
> 22499        311     342551 51214
> 25                           72 Morris Rd, Colonie, NY 12205, USA
> 14924        330     190000 19865
> 26                       79 Morris Rd, Schenectady, NY 12304, USA
> 14929        330      18000 19865
>
>      PERGTBACH PERLTBACH PERCAPINC
> 1    0.11521926 0.8847807     14998
> 5    0.11521926 0.8847807     14998
> 7    0.11521926 0.8847807     14998
> 10   0.12627425 0.8737257     12147
> 11   0.12627425 0.8737257     12147
> 12   0.12627425 0.8737257     12147
> 13   0.12627425 0.8737257     12147
> 14   0.12627425 0.8737257     12147
> 15   0.12627425 0.8737257     12147
> 16   0.12627425 0.8737257     12147
> 19   0.37195452 0.6280455     25915
> 22   0.37195452 0.6280455     25915
> 25   0.16253444 0.8374656     10992
> 26   0.16253444 0.8374656     10992
>
> 	[[alternative HTML version deleted]]
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no



More information about the R-sig-Geo mailing list