Skip to contents

Hierarchical and Geographically Weighted Regression

Source: R/hgwr.R
hgwr.Rd

A Hierarchical Linear Model (HLM) with group-level geographically weighted effects.

Usage

hgwr(
  formula,
  data,
  ...,
  bw = "CV",
  kernel = c("gaussian", "bisquared"),
  alpha = 0.01,
  eps_iter = 1e-06,
  eps_gradient = 1e-06,
  max_iters = 1e+06,
  max_retries = 1e+06,
  ml_type = c("D_Only", "D_Beta"),
  f_test = FALSE,
  verbose = 0
)

# S3 method for class 'sf'
hgwr(
  formula,
  data,
  ...,
  bw = "CV",
  kernel = c("gaussian", "bisquared"),
  alpha = 0.01,
  eps_iter = 1e-06,
  eps_gradient = 1e-06,
  max_iters = 1e+06,
  max_retries = 1e+06,
  ml_type = c("D_Only", "D_Beta"),
  f_test = FALSE,
  verbose = 0
)

# S3 method for class 'data.frame'
hgwr(
  formula,
  data,
  ...,
  coords,
  bw = "CV",
  kernel = c("gaussian", "bisquared"),
  alpha = 0.01,
  eps_iter = 1e-06,
  eps_gradient = 1e-06,
  max_iters = 1e+06,
  max_retries = 1e+06,
  ml_type = c("D_Only", "D_Beta"),
  f_test = FALSE,
  verbose = 0
)

hgwr_fit(
  formula,
  data,
  coords,
  bw = c("CV", "AIC"),
  kernel = c("gaussian", "bisquared"),
  alpha = 0.01,
  eps_iter = 1e-06,
  eps_gradient = 1e-06,
  max_iters = 1e+06,
  max_retries = 1e+06,
  ml_type = c("D_Only", "D_Beta"),
  f_test = FALSE,
  verbose = 0
)

Arguments

formula

A formula. Its structure is similar to lmer function in lme4 package. Models can be specified with the following form:

response ~ L(local.fixed) + global.fixed + (random | group)

For more information, please see the formula subsection in details.

data

The data.

...

Further arguments for the specified type of data.

bw

A numeric value. It is the value of bandwidth or "CV". In this stage this function only support adaptive bandwidth. And its unit must be the number of nearest neighbours. If "CV" is specified, the algorithm will automatically select an optimized bandwidth value.

kernel

A character value. It specify which kernel function is used in GWR part. Possible values are

gaussian

Gaussian kernel function \(k(d)=\exp\left(-\frac{d^2}{b^2}\right)\)

bisquared

Bi-squared kernel function. If \(d<b\) then \(k(d)=\left(1-\frac{d^2}{b^2}\right)^2\) else \(k(d)=0\)

alpha

A numeric value. It is the size of the first trial step in maximum likelihood algorithm.

eps_iter

A numeric value. Terminate threshold of back-fitting.

eps_gradient

A numeric value. Terminate threshold of maximum likelihood algorithm.

max_iters

An integer value. The maximum of iteration.

max_retries

An integer value. If the algorithm tends to be diverge, it stops automatically after trying max_retires times.

ml_type

An integer value. Represent which maximum likelihood algorithm is used. Possible values are:

D_Only

Only \(D\) is specified by maximum likelihood.

D_Beta

Both \(D\) and \(beta\) is specified by maximum likelihood.

f_test

A logical value. Determine whether to do F test on GLSW effects. If f_test=TURE, there will be a f_test item in the returned object showing the F test for each GLSW effect.

verbose

An integer value. Determine the log level. Possible values are:

0

no log is printed.

1

only logs in back-fitting are printed.

2

all logs are printed.

coords

A 2-column matrix. It consists of coordinates for each group.

Value

A list describing the model with following fields.

gamma

Coefficients of group-level spatially weighted effects.

beta

Coefficients of fixed effects.

mu

Coefficients of sample-level random effects.

D

Variance-covariance matrix of sample-level random effects.

sigma

Variance of errors.

effects

A list including names of all effects.

call

Calling of this function.

frame

The DataFrame object sent to this call.

frame.parsed

Variables extracted from the data.

groups

Unique group labels extracted from the data.

f_test

A list of F test for GLSW effects. Only exists when f_test=TRUE. Each item contains the F value, degrees of freedom in the numerator, degrees of freedom in the denominator, and \(p\) value of \(F>F_\alpha\).

Details

Effect Specification in Formula

In the HGWR model, there are three types of effects specified by the formula argument:

Group-level spatially weighted (GLSW, aka. local fixed) effects

Effects wrapped by functional symbol L.

Sample-level random (SLR) effects

Effects specified outside the functional symbol L but to the left of symbol |.

Fixed effects

Other effects

For example, the following formula in the example of this function below is written as

y ~ L(g1 + g2) + x1 + (z1 | group)

where g1 and g2 are GLSW effects, x1 is the fixed effects, and z1 is the SLR effects grouped by the group indicator group. Note that SLR effects can only be specified once!

Functions

  • hgwr_fit(): Fit a HGWR model

Examples

data(multisampling)
hgwr(formula = y ~ L(g1 + g2) + x1 + (z1 | group),
     data = multisampling$data,
     coords = multisampling$coords,
     bw = 10)
#> Hierarchical and geographically weighted regression model
#> =========================================================
#> Formula: y ~ L(g1 + g2) + x1 + (z1 | group)
#>  Method: Back-fitting and Maximum likelihood
#>    Data: multisampling$data
#> 
#> Fixed Effects
#> -------------
#>  Intercept        x1 
#>   2.844509  0.966059 
#> 
#> Group-level Spatially Weighted Effects
#> --------------------------------------
#> Bandwidth: 10 (nearest neighbours)
#> 
#> Coefficient estimates:
#>  Coefficient        Min  1st Quartile     Median  3rd Quartile        Max 
#>    Intercept  -1.430985     -1.350693  -1.140097     -0.911629  -0.727237 
#>           g1   5.969397      6.328338   7.136143      7.447083   8.481628 
#>           g2  -0.683021      0.043467   0.887209      1.174541   1.521918 
#> 
#> Sample-level Random Effects
#> ---------------------------
#>    Groups       Name  Std.Dev.      Corr 
#>     group  Intercept  1.921255           
#>                   z1  1.921255  0.000000 
#>  Residual             1.921255           
#>    Groups       Name  Std.Dev.      Corr 
#>     group  Intercept  1.921255           
#>                   z1  1.921255  0.000000 
#>  Residual             1.921255           
#> 
#> Other Information
#> -----------------
#> Number of Obs: 484
#>        Groups: group , 16

mod_Ftest <-hgwr(
 formula = y ~ L(g1 + g2) + x1 + (z1 | group),
 data = multisampling$data,
 coords = multisampling$coords,
 bw = 10
)
summary(mod_Ftest)
#> Hierarchical and geographically weighted regression model
#> =========================================================
#> Formula: y ~ L(g1 + g2) + x1 + (z1 | group)
#>  Method: Back-fitting and Maximum likelihood
#>    Data: multisampling$data
#> 
#> Parameter Estimates
#> -------------------
#> Fixed effects:
#>             Estimated   Sd. Err      t.val  Pr(>|t|)      
#>  Intercept   2.844509  0.254491  11.177234  0.000000  *** 
#>         x1   0.966059  0.046763  20.658467  0.000000  *** 
#> 
#> Bandwidth: 10 (nearest neighbours)
#> 
#> GLSW effects:
#>             Mean Est.  Mean Sd.   ***    **      *      . 
#>  Intercept  -1.115285  0.557561  0.0%  0.0%  56.2%  18.8% 
#>         g1   6.976060  3.208388  0.0%  6.2%  62.5%  31.2% 
#>         g2   0.663998  3.249016  0.0%  0.0%   0.0%   0.0% 
#> 
#> SLR effects:
#>    Groups       Name      Mean  Std.Dev.      Corr 
#>     group  Intercept  0.000000  1.921255           
#>                   z1  0.005422  1.921255  0.000000 
#>  Residual             0.233641  1.921255           
#>    Groups       Name      Mean  Std.Dev.      Corr 
#>     group  Intercept  0.000000  1.921255           
#>                   z1  0.005422  1.921255  0.000000 
#>  Residual             0.233641  1.921255           
#> 
#> 
#> Diagnostics
#> -----------
#>  rsquared  0.641698 
#>    logLik       NaN 
#>       AIC       NaN 
#> 
#> Scaled Residuals
#> ----------------
#>        Min         1Q    Median        3Q       Max 
#>  -5.222809  -0.982000  0.161400  1.480594  5.562956 
#> 
#> Other Information
#> -----------------
#> Number of Obs: 484
#>        Groups: group , 16