Geographically weighted regression (GWR) has become popular in recent years to deal with spatial autocorrelation and heterogeneity in forestry and ecological data. However, researchers have realized that GWR has some limitations, such as correlated model coefficients across study areas, strong influence of outliers, weak data problem, etc. In this study, we applied Bayesian geographically weighted regression (BGWR) and a robust BGWR (rBGWR) to model the relationship between tree crown and diameter using observed tree data and simulated data to investigate model fitting and performance in order to overcome some limitations of GWR. Our results indicated that, for observed tree data, the rBGWR estimated tree crown more accurate than both BGWR and GWR. For the simulated data, 74.1% of the estimated slope coefficients by rBGWR and 73.4% of the estimated slope coefficients by BGWR were not significantly different (α = 0.05) from the corresponding simulated slope coefficients. The estimation of model coefficients by rBGWR was not sensitive to outliers, but the coefficient estimation by GWR was strongly affected by those outliers. The majority of the coefficient estimates by rBGWR and BGWR for weak observations were not significantly (α = 0.05) different from the simulated coefficients. Therefore, BGWR (including rBGWR) may be a better alternative to overcome some limitations of GWR. In addition, both BGWR and rBGWR were more powerful than GWR to detect the spatial areas with non-constant variance or spatial outliers.
Most tree or plant communities display some degree of spatial structure in the form of geographical patchiness or gradients. Trees or plants at close distances interact to each other (
When spatial autocorrelation and heterogeneity exist in forestry and ecological data, the independence and homogeneity assumptions of traditional statistical methods,
In contrast, local models, such as geographically weighted regression (GWR), fit a regression relationship for each spatial location using the neighbors within a given bandwidth (
However, researchers have realized that GWR has several limitations: (i) lack of a unifying approach to estimate local coefficients; (ii) incorrect estimates of dispersion of model coefficients, as GWR utilizes the same sample data observations (with different spatial weights) to produce a sequence of parameter estimates for all locations in the space. Consequently, the estimated model coefficients are dependent across spatial locations, and the correct expression for the variance of local coefficients could be non-linear; (iii) the influence of outliers
To overcome the limitations of GWR, a Bayesian approach was proposed to estimate the local regression coefficients, namely Bayesian geographically weighted regression (BGWR -
However, BGWR requires the use of Gibbs sampler, a Monte Carlo simulation, to estimate its parameters. The Gibbs sampler (
To date, limited studies have used the BGWR methods.
A linear regression relationship between variables can be expressed as (
where
Now suppose that we have a set of locations
where
where
where
The spatial weight matrix
The underlying model of Bayesian geographically weighted regression (BGWR) is the same as that of GWR in
where
The parameter smoothing function and assumption of heterogeneous variance ensure that BGWR captures spatial heterogeneity and autocorrelation. On the other hand, these assumptions add complexity to the modeling process. Thus, a Bayesian approach is required to estimate the mod-el coefficients. The likelihood of BGWR (
where
Using
where
Choosing a prior distribution or prior information for the hyper-parameters is an important issue in Bayesian inference, as it may be sensitive or robust to the choice of the prior distribution (
(1) The conditional prior distribution of scale factor
where
Note that the value of
(2) The conditional prior distribution of
where
(3) The conditional prior distribution of
where
which provides greater flexibility to
Given the likelihood function and the conditional prior distributions of BGWR hyper-parameters, the posterior distribution of βÌ(
We briefly discuss the observed data and simulated data, crown area regression model, model fitting process for each method (GWR, BGWR, and rBGWR), and the evaluation of the three methods.
The observed data used in this study were a subset of the stem map data of a mature, second growth, uneven-aged softwood stand located near Sault Ste. Marie, Ontario, Canada (
Based on the relationship between tree dbh and crown in the data, we chose the following linear regression model (
where
To objectively verify the performance of BGWR, a simulated tree data set with geographic coordinates, tree dbh and crown area, and known local regression coefficients was generated for each tree. The procedure of generating the simulated data is as follows:
(1) The AMORPHYS software (
(2) The generated tree data were used to fit
(3) Given the estimated association between β0 and β1, the bivariate Gaussian copula was used to randomly simulate the two coefficients (β0 and β1) in a unit scale [0, 1], and then transform these coefficients back to the original scale and assign the coefficients to the spatial locations generated in step (1).
(4) Create 5 outliers by artificially making the β1 coefficient very large (
(5) Given the randomly generated regression coefficients (β0 and β1) and tree dbh generated in step (1), tree crown area was computed using
Three modeling methods, GWR, BGWR, and rBGWR, were used to fit the linear relationship between ln(crown) and ln(dbh) (
The BGWR and rBGWR models were fitted using bgwr.m and bgwrv.m Matlab functions available in the Econometric Toolbox (
Posterior distributions of model parameters for each model were sampled using the Gibbs sampler. A single long chain of MCMC was used to generate samples of each coefficient. The convergence of sampled coefficients was diagnosed using the method of
Using the post convergence values of model coefficients, the Deviance Information Criterion (DIC -
The local coefficients estimated by the three methods (GWR, BGWR, and rBGWR) were summarized using descriptive statistics. For the Ontario softwood data, the two regression coefficients estimated by each method were interpolated using inverse distance weighting method and contour maps of these coefficients were generated for each method. Further, the model residuals from each method were assessed across a range of the explanatory variable (dbh).
The local regression coefficients were also evaluated for spatial heterogeneity using intra-block variance (
where
For the simulated data, the “true” local coefficients of each observation were known. Hence, the 95% of credible limits for each coefficient estimated by BGWR and rBGWR were compared with the corresponding known coefficients for each location. If the 95% credible limits of the posterior samples of a coefficient estimated for a particular location included the known coefficient, the coefficient estimated by the method for the location was declared to be unbiased. The percentages of unbiased intercept, slope, and both (intercept and slope) coefficients for all locations by BGWR and rBGWR were calculated. However, this assessment was not used for GWR because the estimation of variance of local coefficients was incorrect (
Model fitting and performance were evaluated for special cases such as outliers and weak data (any observation with insufficient neighbors within its effective bandwidth) by comparing the estimated coefficients and predicted crown area by each modeling method. For the Ontario softwood data, the comparison was only for weak observations (no obvious outliers in the data), while for simulated data, the modeling fitting and performance were evaluated for both outliers and weak observations in terms of the 95% credible limits interval of the corresponding BGWR and rBGWR coefficients, and the estimation bias of model coefficients and tree crown area.
The identification of the “best” BGWR and rBGWR models for different priors of the hyper-parameters
Interestingly, the BGWR and rBGWR models had consistently smaller DIC values for
The two local regression coefficients were used to generate contour maps across the study area for the three modeling methods (
Similar to the results obtained for the observed data, the posterior mean of
A summary of descriptive statistics on the local regression coefficients estimated by the three modeling methods are given in
The 95% credible limits of each model coefficient for both BGWR and rBGWR were used to assess the accuracy of estimates of the local model coefficients. If the 95% credible limit of a coefficient for a particular location included the “known” simulated coefficient, the coefficient estimated by that method for that location was declared to be unbiased. The results indicated that, out of 1230 comparisons, 853 (69.3%) of β0 and 903 (73.4%) of β1 by BGWR, and 864 (70.2%) of β0 and 912 (74.1%) of β1 by rBGWR, respectively, were unbiased. When the two coefficients were compared simultaneously, the rBGWR method (839 or 68.2%) was slightly better than the BGWR method (826 or 67.2%).
To assess the sensitivity of GWR and rBGWR to outliers in the data, two outliers (Tree ID 141 and 541) were selected to compare the local coefficients by the two methods.
In this study, we established a relationship between tree crown and dbh for observed and simulated data using GWR, BGWR and rBGWR, and demonstrated the advantages of Bayesian approach to improve the model fitting and performance of GWR. Firstly, BGWR and rBGWR not only provided better estimation of local regression coefficients, but also correct estimation of the variances of the coefficients because Bayesian procedures (
However, the BGWR methods need to be used with caution. Given their nature, the local coefficients estimated by the BGWR methods may not be used to make prediction outside the study area. Attention should be given to assess potential model misspecification, because the spatial heterogeneities of the relationships between variables might have been produced as a result of model misspecification (
One of the challenges on modeling spatial effects (
Another challenge to the BGWR methods is related to the convergence of parameters. In this study, longer MCMC chains were used to sample the posterior distributions of local coefficients. While sampling posterior of coefficients for nine BGWR models (with a combination of the priors on variance and scale factor) and for three rBGWR models (with three priors on variance), the slope and intercept coefficients estimated by these models converged at the desired level of precision (quartile = 0.025, precision = 0.005 and probability = 0.95) for most of the observations within the 4.000 runs after the few thousands initial burn-in and storing one sample out of 10 samples. However, for a few observations, those that fall at the edge of the study area required almost 15.000 iterations (at the maximum for the BGWR method with prior on scale factor at the smallest level) to converge. The assessment of histogram of coefficients and its movement across iterations also suggested the desirability of longer runs for the latter observations. One of the reasons for the convergence problem may be the correlation between the parameters within the parameter set (
Despite these limitations, BGWR methods have been recommended as an effective tool for investigating spatial heterogeneity of covariate effects across a study area. Future studies on developing an optimization method as a function of both bandwidth and covariates are recommended. The current approach of bandwidth optimization (
The relationship between tree crown area and diameter at breast height was modeled using GWR, BGWR, and rBGWR models. Both observed and simulated tree data were used to compare model fitting and performance. For the observed Ontario softwood stand, BGWR and rBGWR produced more accurate parameter estimates and better model prediction than GWR. Further, rBGWR produced the largest localized spatial variability for both coefficients, followed by BGWR, while GWR had the smallest localized spatial variability for both coefficients. Thus, the contour maps of local regression coefficients show more “hot” or “cold” spots of the two coefficients estimated by GWR than those obtained by BGWR and rBGWR.
For the simulated data, the percentage of the unbiased parameter estimation for the two local coefficients by rBGWR was slightly higher than those by BGWR. However, the credible limits cannot be computed for GWR due to the lack of an accurate expression of the variance of the model coefficients in GWR. Further, the coefficients estimates using the rBGWR method were not affected by outliers, while these estimates by GWR were “contaminated” by the outliers, resulting in erroneous parameter estimates. Similarly, for the weak data, the majority of the coefficients estimated by the rBGWR and BGWR methods were not different from the simulated “known” coefficients. Therefore, BGWR and rBGWR are better alternatives to overcome the limitations of GWR and are more powerful than GWR to detect the spatial areas with non-constant variance or spatial outliers.
This work was supported by the National Natural Science Foundation of China, proj. no. 31400491: “Three dimensional individual tree crown delineation based on multiple remotely sensed data sources and tree competition mechanism”, and by the Fundamental Research Funds for the Central Universities, proj. no. 2572017CA02.
Map of tree locations for Ontario softwood data. Circles are proportional to the diameter at breast height of trees.
Map of tree locations for simulated data. Circles are proportional to the diameter at breast height of trees.
Average model prediction errors by GWR, BGWR, and rBGWR (Ontario softwood data).
Contour maps of local regression coefficients. (a) GWR - β0; (b) BGWR - β0; (c) rBGWR - β0; (d) GWR - β1; (e) BGWR - β1; and (f) rBGWR - β1.
Intra-block variances (
The outlier and neighboring observations within the bandwidth. The two regression lines of GWR and rBGWR were generated using the estimated slope and intercept coefficients. The response variable (crown area) and explanatory variable (dbh) are weighted by proximity to outlier observations and expressed in logarithmic scale.
Descriptive statistics of observed tree data (Ontario softwood stand). Moran’s
Variable | n | Mean | SD | Min | Max | Moran’s |
SH% | Distance(m) |
---|---|---|---|---|---|---|---|---|
dbh (cm) | 1698 | 18.68 | 8.47 | 11.43 | 84.07 | 0.087 | 91.64 | 2.2(0.3-11.0) |
Crown area (m2) | 1698 | 7.79 | 7.80 | 0.93 | 90.48 | 0.115 | 80.62 | - |
Descriptive statistics of simulated tree data (n=1230).
Variable | Mean | SD | Min | Max | Moran’s |
SH% |
---|---|---|---|---|---|---|
DBH (cm) | 18.23 | 9.44 | 7.43 | 80.45 | 0.03 | 86.40 |
Crown area (m2) | 12.94 | 7.00 | 4.61 | 74.91 | 0.05 | 86.90 |
Intercept β0 | 0.201 | 0.143 | -0.260 | 0.580 | 0.01 | 87.80 |
Slope β1 | 0.807 | 0.130 | 0.600 | 1.233 | 0.01 | 92.60 |
Nearest neighbor distance | 1.4 | 1.26 | 0.15 | 6.47 | - | - |
Model fitting statistics of nine BGWR models and three rBGWR models for the Ontario softwood data.
Method | Model | Prior Meanof |
Scale Factor |
DIC | D(β) | pD |
---|---|---|---|---|---|---|
BGWR | 1 | 4 | 10 | 4183.1 | 893.9 | 3289.2 |
2 | 4 | 25 | 4149.6 | 899.7 | 3249.9 | |
3 | 4 | 50 | 4115.8 | 904.0 | 3211.8 | |
4 | 5 | 10 | 3891.4 | 1039.6 | 2851.8 | |
5 | 5 | 25 | 3852.6 | 1042.2 | 2810.3 | |
6 | 5 | 50 | 3826.7 | 1044.2 | 2782.5 | |
7 | 6 | 10 | 4516.4 | 1160.0 | 3356.4 | |
8 | 6 | 25 | 4500.0 | 1160.3 | 3339.7 | |
9 | 6 | 50 | 3751.1 | 1161.1 | 2590.1 | |
rBGWR | 1 | 4 | - | 4620.2 | 1279.2 | 3341.1 |
2 | 5 | - | 4193.0 | 1307.0 | 2886.0 | |
3 | 6 | - | 5273.4 | 1323.7 | 3949.7 |
Summary of local regression coefficients by the three modeling methods (Ontario softwood data).
Coefficient | Method | n | Mean | Median | SD | Skewness | Kurtosis | Min | Max |
---|---|---|---|---|---|---|---|---|---|
Intercept β0 | GWR | 1698 | -2.258 | -2.293 | 1.245 | 0.541 | 4.845 | -6.087 | 3.672 |
BGWR | 1698 | -2.250 | -2.225 | 1.407 | 0.231 | 4.877 | -8.395 | 5.098 | |
rBGWR | 1698 | -2.242 | -2.212 | 1.480 | 0.227 | 4.942 | -8.665 | 5.384 | |
Slope β1 | GWR | 1698 | 1.426 | 1.440 | 0.424 | -0.596 | 4.845 | -0.736 | 2.774 |
BGWR | 1698 | 1.428 | 1.425 | 0.482 | -0.252 | 4.877 | -1.113 | 3.485 | |
rBGWR | 1698 | 1.426 | 1.421 | 0.507 | -0.255 | 4.942 | -1.216 | 3.576 |
Model fitting statistics for BGWR and rBGWR to the simulated data. BGWR models were fitted with scale factor = 50%.
Method | Model | Prior Meanof |
DIC | D( |
pD | |||
---|---|---|---|---|---|---|---|---|
Min | Max | STD | ||||||
BGWR | 1 | 4 | 1567.24 | 632.04 | 935.21 | 0.9949 | 1.0212 | 0.0022 |
2 | 5 | 1587.61 | 699.24 | 888.36 | 0.9957 | 1.0204 | 0.0019 | |
3 | 6 | 1651.40 | 750.07 | 875.07 | 0.9964 | 1.0196 | 0.0017 | |
rBGWR | 4 | 4 | 1550.99 | 746.42 | 804.57 | 1.0400 | 30.8000 | 0.8891 |
5 | 5 | 1568.63 | 764.23 | 804.40 | 1.0300 | 14.7800 | 0.4192 | |
6 | 6 | 1583.89 | 776.21 | 807.68 | 1.0000 | 8.3100 | 0.2276 |
Summary of local regression coefficients from GWR, BGWR and rBGWR methods for the simulated data. Q1, Q3, P2.5 and P97.5 are lower quartile, upper quartile, 2.5th percentile, and 97.5th percentile, respectively. (SD): standard deviation.
Coefficient | Method | n | Mean | Median | Q1 | Q3 | P2.5 | P97.5 | Min | Max | SD |
---|---|---|---|---|---|---|---|---|---|---|---|
Intercept β0 | Simulated | 1230 | 0.201 | 0.198 | 0.096 | 0.312 | -0.08 | 0.478 | -0.260 | 0.58 | 0.143 |
GWR | 1230 | 0.231 | 0.265 | -0.124 | 0.685 | -1.526 | 1.683 | -7.656 | 4.801 | 0.917 | |
BGWR | 1230 | 0.300 | 0.302 | -0.131 | 0.788 | -1.436 | 2.002 | -7.702 | 3.303 | 0.847 | |
rBGWR | 1230 | 0.298 | 0.302 | -0.16 | 0.812 | -1.456 | 2.035 | -7.805 | 3.462 | 0.87 | |
Slope β1 | Simulated | 1230 | 0.811 | 0.800 | 0.696 | 0.903 | 0.610 | 1.079 | 0.600 | 2.248 | 0.154 |
GWR | 1230 | 0.797 | 0.781 | 0.622 | 0.939 | 0.278 | 1.418 | -0.448 | 3.613 | 0.334 | |
BGWR | 1230 | 0.766 | 0.765 | 0.583 | 0.941 | 0.158 | 1.401 | -0.349 | 3.272 | 0.313 | |
rBGWR | 1230 | 0.766 | 0.764 | 0.576 | 0.948 | 0.143 | 1.423 | -0.430 | 3.304 | 0.322 |
Simulated and geographically weighted regression (GWR) estimated coefficients, and Bayesian geographically weighted regression (BGWR) and robust BGWR (rBGWR) estimated 95% credible limits of coefficients for weak observations in simulated data. (‡): indicates the 95% of credible limits of BGWR and rBGWR coefficients that do not contain simulated coefficients.
Coefficients | Tree ID | Location | Simulated | GWR | BGWR | rBGWR | |||
---|---|---|---|---|---|---|---|---|---|
Easting | Northing | 95% Credible limits | 95% Credible limits | ||||||
Intercept | 133 | 75.1 | 15.5 | -0.023 | -0.392 | -0.661 ‡ | -0.094 ‡ | -0.553 ‡ | -0.179 ‡ |
174 | 37.5 | 21.9 | 0.014 | 1.556 | -0.544 | 2.607 | -0.438 | 2.583 | |
265 | 59.9 | 31.8 | 0.080 | -0.708 | -6.598 | 5.564 | -4.230 | 3.138 | |
620 | 32.7 | 74.7 | 0.020 | -0.880 | -2.913 | 1.048 | -2.623 | 0.619 | |
649 | 85.2 | 76.6 | 0.340 | -0.040 | -0.189 ‡ | 0.132 ‡ | -0.190 ‡ | 0.130 ‡ | |
924 | 142.6 | 23.1 | 0.135 | -0.120 | -1.103 | 0.871 | -0.918 | 0.757 | |
925 | 137.5 | 21.9 | 0.399 | -1.681 | -4.007 ‡ | -0.286 ‡ | -4.034 ‡ | -0.484 ‡ | |
1132 | 132.7 | 74.7 | 0.119 | 0.130 | -0.802 | 0.791 | -0.738 | 0.679 | |
Slope | 133 | 75.1 | 15.5 | 0.925 | 1.045 | 0.945 ‡ | 1.136 ‡ | 0.974 ‡ | 1.096 ‡ |
174 | 37.5 | 21.9 | 0.971 | 0.341 | -0.082 | 1.191 | -0.070 | 1.153 | |
265 | 59.9 | 31.8 | 0.951 | 1.235 | -1.083 | 3.409 | -0.174 | 2.543 | |
620 | 32.7 | 74.7 | 1.034 | 1.301 | 0.653 | 1.995 | 0.827 | 1.903 | |
649 | 85.2 | 76.6 | 0.700 | 0.864 | 0.793 ‡ | 0.929 ‡ | 0.792 ‡ | 0.930 ‡ | |
924 | 142.6 | 23.1 | 0.853 | 0.941 | 0.584 | 1.313 | 0.620 | 1.248 | |
925 | 137.5 | 21.9 | 0.649 | 1.496 | 0.931 ‡ | 2.433 ‡ | 1.008 ‡ | 2.444 ‡ | |
1132 | 132.7 | 74.7 | 0.837 | 0.824 | 0.601 | 1.140 | 0.650 | 1.119 |