Linear models¶

An overview of linear and logistic regression models.

import geopandas as gpd
import matplotlib.pyplot as plt
from geodatasets import get_path
from sklearn import metrics

from gwlearn.linear_model import GWLinearRegression, GWLogisticRegression

Linear regression¶

Linear regression is a classical example of geographically weighted regression models. In this example, you can predict the number of suicides based on other population data in the Guerry dataset.

gdf = gpd.read_file(get_path("geoda.guerry"))

It is a relatively small dataset covering Frech data from 1800s.

gdf.plot().set_axis_off()

Specify the model and fit the data.

adaptive = GWLinearRegression(bandwidth=25, fixed=False, kernel="tricube")
adaptive.fit(
    gdf[["Crm_prp", "Litercy", "Donatns", "Lottery"]],
    gdf["Suicids"],
    geometry=gdf.representative_point(),
)

GWLinearRegression(bandwidth=25, kernel='tricube')

You can get a number of outputs from the fitted model. For example, focal predictions (a value is predicted on an observations using the local model fitted around it).

adaptive.pred_

0     65929.297258
1     14289.540382
2     60668.901578
3     25191.170857
4     12777.621045
          ...     
80    41304.669559
81    19766.026593
82    29125.591134
83    25717.441407
84    19381.885833
Length: 85, dtype: float64

This allows you to measure model performance metrics like $R^2$. This value would be comparable to the R2 reported by mgwr.

metrics.r2_score(gdf["Suicids"], adaptive.pred_)

0.7030520942210708

Local version of $R^2$ is accessible directly.

adaptive.local_r2_

0     0.655274
1     0.564506
2     0.616975
3     0.690529
4     0.704666
        ...   
80    0.444105
81    0.554521
82    0.646275
83    0.174557
84    0.215379
Length: 85, dtype: float64

Any of the values linked to individual observations can then be mapped.

gdf.plot(adaptive.local_r2_, legend=True).set_axis_off()

Similarly, you can access residuals.

gdf.plot(adaptive.resid_, legend=True).set_axis_off()

Each local model then have its own local coefficients.

adaptive.local_coef_

Again, this can be explored visually.

f, axs = plt.subplots(2, 2, figsize=(12, 10))

for column, ax in zip(adaptive.local_coef_.columns, axs.flat):
    gdf.plot(adaptive.local_coef_[column], legend=True, ax=ax)
    ax.set_title(column)
    ax.set_axis_off()

Alongside local coefficients, you can retireve a local intercept value.

gdf.plot(adaptive.local_intercept_, legend=True).set_axis_off()

For details on prediction, see the Prediction guide.

Logistic regression¶

Due to the nature of the principle in which gwlearn extends scikit-learn, you can also fit the local version of logistic regression for classification tasks. However, there’s a caveat. Given it is unlikely that all categories are present in all local models, fitting a non-binary would lead to inconsistent local models. Hence gwlearn currently supports only binary regression, akin to binomial model in mgwr.

Let’s illustrate it by prediciting whether the number of suicides is over or under median.

y = gdf["Suicids"] > gdf["Suicids"].median()

The API then looks the same, following scikit-learn’s model.

binary = GWLogisticRegression(bandwidth=25, fixed=False, max_iter=500)
binary.fit(
    gdf[["Crm_prp", "Litercy", "Donatns", "Lottery"]],
    y,
    geometry=gdf.representative_point(),
)

GWLogisticRegression(bandwidth=25)

Given this is a classification task, we can retrieve probabilites from focal predictions.

binary.proba_

Or their binary counterparts (based on the maximum value).

binary.pred_

0      True
1      <NA>
2      True
3     False
4     False
      ...  
80     True
81    False
82     True
83    False
84    False
Length: 85, dtype: boolean

The regions with missing values are not fitted due to extreme class imbalance in their local context (less than 20% of local y values are of a minority class). This is a specific behavior linked to classification models in gwlearn. Rather than simply failing, invariant or extremely imbalanced local y does not make the model to fail, it just records the location as skipped. That also means that the location is not used in prediction and for certain locations, there is no possible prediction.

You can check the proportion of local models that exist using prediction_rate_.

binary.prediction_rate_

np.float64(0.8705882352941177)

This shows that 87% of focal geometries have their respective models, the rest does not.

For more on dealing with class imbalance, see the imbalance guide.

The focal metrics (based on a prediction on a focal geometry using a single local model fitted around that geometry) can be measured using the two outputs shown above.

na_mask = binary.pred_.notna()

metrics.accuracy_score(y[na_mask], binary.pred_[na_mask])

0.9594594594594594

However, you can also extract all data from all local models and use those to measure the performance of the model, rather than using only focal geometries.

The data pooled from all local models are accessible as y_pooled_ and pred_pooled_.

metrics.accuracy_score(binary.y_pooled_, binary.pred_pooled_)

0.7864864864864864

Pooled metrics are typically showing worse results as the distance decay of local observation weights in individual models are playing against a good prediction of values far from the focal point.

Alternatively, you can measure performance metrics per each local model. That is relying on the same pooled data but split to their respective parent local models.

local_accuracy = binary.local_metric(metrics.accuracy_score)
local_accuracy

array([0.64,  nan, 0.84, 0.72, 0.6 , 0.84,  nan, 0.92, 0.84, 0.92,  nan,
       0.8 , 0.96, 0.88, 0.76, 0.84, 0.8 , 0.84, 0.64, 0.6 , 0.8 , 0.88,
       0.56, 0.76, 0.96, 0.88, 0.52, 0.92, 0.8 , 0.8 , 0.92, 0.92, 0.72,
       0.72, 0.6 , 0.72, 0.64, 0.92, 0.64, 0.72, 0.88, 0.76, 0.84,  nan,
       0.88, 0.96, 0.68, 0.8 ,  nan, 0.68, 0.8 , 0.84, 0.88, 0.56, 0.8 ,
       0.8 ,  nan,  nan, 0.88, 1.  , 0.76, 0.92, 0.76, 0.88, 0.76, 0.56,
       0.8 , 0.64, 0.72, 0.72,  nan, 0.96,  nan, 0.88, 0.8 , 0.96,  nan,
        nan, 0.72, 0.76, 0.72, 0.68, 0.76, 0.72, 0.84])

gdf.plot(
    local_accuracy, legend=True, missing_kwds=dict(color="lightgray")
).set_axis_off()

The local coefficients are available as in the regression counterpart.

binary.local_coef_

Again, this can be explored visually.

f, axs = plt.subplots(2, 2, figsize=(12, 10))

for column, ax in zip(binary.local_coef_.columns, axs.flat):
    gdf.plot(
        binary.local_coef_[column],
        legend=True,
        ax=ax,
        missing_kwds=dict(color="lightgray"),
    )
    ax.set_title(column)
    ax.set_axis_off()

For details on prediction, see the Prediction guide.