This page was generated from notebooks/8_spatial_multipliers.ipynb. Interactive online version:

Spatial Multipliers

Luc Anselin

(revised 09/17/2024)

Preliminaries

This notebook takes a closer look at the spatial pattern induced by spatial lag terms in a regression specification. The usual notion of average impacts is disassembled into its individual components, with a particular focus on its spatial distribution. Three broad types of models are considered: a linear SLX model, a model with a spatially lagged dependent variable (Wy), and nonlinear SLX models.

Prerequisites

This notebook includes material that relies more on familiarity with Python than the previous notebooks, specifically the manipulation of pandas data frames and geopandas geo data frames, as well as the interpretation of the helper functions. In addition, spatial weights functionality and choropleth mapping are used extensively. Both were covered in previous notebooks. Two specialized routines are included from spreg, i.e., i_multipliers (to compute the individual multipliers) and make_wnslx (to create the lag operator in nonlinear SLX models). It is also assumed that the sample data set Police is installed (if not, execute libpysal.examples.load_example("Police") first, with the module libpysal.examples imported).

Modules Needed

The usual imports include numpy, pandas, geopandas and matplotlib (for mapping). More specialized imports consist of examples.get_path and weights from libpysal, as well as i_multipliers (from spreg.sputils) and make_wnslx (from spreg.utils). To allow the complete set of observations (n=82) to be listed in the notebook, the pd.options.display.max_rows should be set to a value larger than 82. Also, as usual, numpy set_printoptions should be set to legacy="1.25".

import warnings
warnings.filterwarnings("ignore")
import os
os.environ['USE_PYGEOS'] = '0'

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from libpysal.examples import get_path
import libpysal.weights as weights
from spreg.sputils import i_multipliers
from spreg.utils import make_wnslx
np.set_printoptions(legacy="1.25")
pd.options.display.max_rows = 100

Functionality Used

• from numpy:

  • array

• from pandas/geopandas:

  • idxmax, idxmin

  • iloc

  • plot

  • concat

  • read_file

  • centroid

  • describe

• from libpysal.examples:

  • get_path

• from libpysal.weights:

  • Queen.from_dataframe

  • neighbors

  • cardinalities

  • KNN.from_dataframe

  • Kernel.from_dataframe

• from spreg:

  • sputils.i_multipliers

  • utils.make_wnslx

Helper Functions

Since many operations will be repeated almost verbatim for the different types of weights and multipliers, two helper functions are used. One, multmap, simplifies the mapping, with (very) limited customization. If additional customization is desired, it must be made in the function itself. The function is essentially a wrapper around the commands to make a quantile choropleth map for a geo data frame.

The second function, nbreffect, finds the minimum and maximum locations for a given type of multiplier effect and then calls multmap to create the associated map.

def multmap(dfm,column='EonNbrs',model='slx'):
    """
    Creates a quintile map (k=5) for the multipliers computed by means
    of i_multiplier

    Arguments
    ---------
    dfm       : merged geo data frame with original data and spatial multiplier
                data frame
    column    : type of multiplier to be mapped, default is EonNbrs
    model     : spatial model, default is slx (i.e., only weights matrix)

    Returns
    -------
    Draws map

    """
    ax = dfm.plot(
            column = column,
            scheme = 'Quantiles',
            k = 5,
            cmap = 'YlOrRd',
            edgecolor = "Black",
            linewidth = 0.2,
            figsize = (6,6),
            legend = True,
            legend_kwds={"loc":"center left","bbox_to_anchor":(1,0.5), "title": column}
    )
    newtitle = "Spatial Multipliers for " + model
    ax.set_title(newtitle)
    ax.set_axis_off()
    plt.show()
    return


def nbreffect(dfs,dfm,oid,mult='EonNbrs',model='slx'):
    """
    Computes the minimum and maximum value and location for a given type of multiplier
    and calls multmap to create a quintile map

    Arguments
    ---------
    dfs       : initial geo data frame
    dfm       : data frame with multipliers
    oid       : pandas series with ID variable
    multi     : type of multiplier, default is EonNbrs
    model     : spatial model, default is slx

    """
    if mult == 'Direct':
        coli = 1
    elif mult == 'EofNbrs':
        coli = 2
    elif mult == 'EonNbrs':
        coli = 3
    else:
        raise Exception("Invalid column")
    imax = int(dfm[[mult]].idxmax())
    print(f"The maximum is {dfm.iloc[imax,coli]:0.3f} with id {int(oid.iloc[imax])}")
    imin = int(dfm[[mult]].idxmin())
    print(f"The minimum is {dfm.iloc[imin,coli]:0.3f} with id {int(oid.iloc[imin])}")
    dfmap = pd.concat([dfs,dfm],axis=1)
    multmap(dfmap,column=mult,model=model)
    return

Data Input

The different multipliers are illustrated with the Police sample data set. A geo data frame is created which is used in the remainder to construct the spatial weights. The only input is the shape file:

• police.shp (shx,dbf): police expenditure data for 82 Mississippi counties

The variable FIPSNO is taken as the id variable and turned into a pandas Series fipsid. The variable itself is set to idvar. This makes it easy to use the same code for other data sets (replace infileshp, fipsid and idvar by the appropriate file name/variable name).

infileshp = "police.shp"   # input shape file
inpath = get_path(infileshp)
dfs = gpd.read_file(inpath)
print(dfs.columns)
fipsid = dfs[["FIPSNO"]]
idvar = "FIPSNO"

Coordinates

The construction of the sparse weights for the nonlinear SLX power and exponential functions requires a numpy array of X and Y coordinates as input. The Police data set does not include those variables explicitly, but they can be computed as the centroids of the polygons, using the centroid attribute of the geopandas data frame.

However, the make_wnslx function needs the coordinates as a numpy array, whereas the result of the centroid attribute is a geopandas geoseries with POINT geometries. This is turned into a numpy array by extracting the x and y coordinates and passing these to np.array. As it turns out, the result must be transposed, since as is, it becomes a 2 by 82 matrix and not the desired 82 by 2 matrix.

Note that the coordinates are decimal degrees (longitude-latitude), which will require the argument distance_metric = "Arc" in the make_wnslx function.

gcent = dfs.centroid
print(type(gcent))
print(gcent.head())
coords = np.array((gcent.x,gcent.y)).T
print(coords.shape)
coords[0:5,:]

Spatial Weights Characteristics

The simplest calculations are for the linear SLX model, where the spatial spillovers are limited to the neighbors as specified in the spatial weights matrix. With row-standardization, the effect of the neighbors always sums to one, so by constructions it is the same (=1) for all observations. Similarly, because the spatial weights have zero on the diagonal, there is no direct effect other than that already captured by coefficients of the X-variables. However, because the row-standardization introduces an asymmetry in the weights (but not in the contiguity structure), the effect on the neighbors of a change in X in a given location is not constant across all observations. This effect is the sum of the column elements associated with each observation. Whereas the mean of these sums is the same as the mean of the row sums and equals one, there remains considerable variation among the contributions of changes in a variable at a location on its neighbors. This is masked by using average effects.

The following three examples show the spatial distribution of the individual multipliers for queen contiguity, k-nearest neighbors and kernel weights. Queen contiguity is intrinsically symmetric, but becomes asymmetric after row-standardization. K-nearest neighbor weights tend not to be symmetric to begin with (k-nearest neighbor is not a symmetric relationship). Kernel weights that are based on a fixed distance band are symmetric. Since they are not row-standardized, the effect of the neighbors and effect on the neighbors is the same. However, when kernel weights are based on an adaptive bandwidth (such as k-nearest neighbors), they will be intrinsically asymmetric and the row and column elements for an observation will differ. In all cases, the average masks sometimes substantial spatial variation among observations.

Queen Contiguity Weights

As covered in a previous notebook, queen contiguity weights are constructed from the geo data frame by means of weights.Queen.from_dataframe, here using ids=idvar for the id-variable (which is "FIPSNO" in this notebook). The weights are row-standardized and then the individual multipliers are calculated using i_multipliers with model="slx" and id=fipsid (the pandas Series with FIPSNO created above). There is no need to specify a spatial coefficient since the default is coef=0.0. Some descriptive statistics are provided by describe().

The first column gives the ID-variable that was used and is not very meaningful. The second column describes the Direct effects, which are zero by construction. The last two columns pertain to the effect of neighbors (row sum, EofNbrs) and the effect on neighbors (column sum, EonNbrs).

wq = weights.Queen.from_dataframe(dfs,ids=idvar)
wq.transform = 'r'
queen = i_multipliers(wq,model='slx',id=fipsid)
queen.describe()

As expected, the mean of both effect of neighbors and effect on neighbors is the same and equals one. However, there is quite a bit of variation in the effect on neighbors, with a range from 0.560 to 1.667 and a standard deviation of 0.25. The full range of variability can be seen by printing out the complete data frame (to see all the rows, the pandas options.display.max_rows must be set to more than 82).

print(queen)

The position and ID of the largest and smallest effects on neighbors can be found by means of idxmax and idxmin applied to the EonNbrs data series. This yields the county with FIPS 28131 (Stone) for the maximum and the county with FIPS 28045 (Hancock) for the minimum. The spatial variation (and location of minimum and maximum) is directly related to the graph-theoretic structure of the spatial weights.

imax = int(queen[['EonNbrs']].idxmax())
print(f'The maximum is {queen.iloc[imax,3]:0.3f} with id {int(fipsid.iloc[imax])}')
imin = int(queen[['EonNbrs']].idxmin())
print(f'The minimum is {queen.iloc[imin,3]:0.3f} with id {int(fipsid.iloc[imin])}')

To get a better insight into the effect of the asymmetry introduced by the row-standardization, the following lines of code extract the list of neighbors for county 28131, then checks their respective cardinalities from wq.cardinalities and computes the cumulative sum of the corresponding weights (the inverse of the cardinalities). Because many of the neighbors have themselves (much) fewer than seven neighbors, their contribution to the cumulative sum is more than 1/7, yielding an overall multiplier effect that can become larger than one. The opposite case is where the neighbors have more neighbors themselves, yielding a spatial weight less than 1/7 which contributes to a multiplier effect that can be smaller than one.

As an exercise, repeat this analysis for the county with the minimum multiplier of 0.560 by entering its fips code in wq.neighbors (county 28045 - this county only has three neighbors).

nbrs = wq.neighbors[28131]
print("neighbor list ",nbrs)
cards = wq.cardinalities
j = 0
for i in nbrs:
    u = cards[i]
    print(f"for {i}, the number of neighbors is {u}, with weight {1.0/u:0.3f}")
    j += 1.0/u
print(f"total effect on neighbors: {j:0.3f}")

The spatial distribution of the multiplier effects can be mapped. This functionality is wrapped in the functions listed at the top of the notebook, but first each step is spelled out in detail here. To begin, pd.concat is used to add the new dataframe to the existing spatial data frame. Because this operation will be repeated for each type of multiplier, the merged data frame is given a new name (otherwise, dfs would be accumulating all the multiplier effecs that have the same name). The multipliers are visualized in a quintile map (k=5) using a common set of customizations as covered in the mapping notebook.

The resulting map reveals considerable variation in the spatial pattern. From a policy perspective, the darkest counties are the ones with the most impact on their neighbors and thus they could be targeted for maximum effect place-based policies.

dfq = pd.concat([dfs,queen],axis=1)
ax = dfq.plot(
    column = 'EonNbrs',
    scheme = 'Quantiles',
    k = 5,
    cmap = 'YlOrRd',
    edgecolor = "Black",
    linewidth = 0.2,
    figsize = (6,6),
    legend = True,
    legend_kwds={"loc":"center left","bbox_to_anchor":(1,0.5), "title": 'EonNbrs'}
)
ax.set_title("Spatial Multipliers")
ax.set_axis_off()

Since the mapping will be repeated for each multiplier, the various settings are encapsulated in the multmap helper function listed at the top of the notebook. As a result, the map can now be obtained with a single line command.

multmap(dfq)

K-Nearest Neighbor Weights

K-nearest neighbor weights are constructed from the geo data frame by means of weights.KNN.from_dataframe. Since the median number of neighbors for queen contiguity was 5, k is set to 5 as well, to keep some degree of comparability. The multipliers are again obtained with i_multipliers passing the row-standardized weights, the model (slx) and the identifiers (idvar).

wk = weights.KNN.from_dataframe(dfs,k=5,ids=idvar)
wk.transform = 'r'
knn = i_multipliers(wk,model='slx',id=fipsid)
knn.describe()

Again, by construction, the direct effects are zero and the effect of the neighbors is 1, due to row-standardization. However, there is considerable variation in the effect on neighbors, again due to the asymmetry of the weights. The multipliers range from 0.4 to 1.8, a wider range than for queen contiguity.

Further insight into the spatial pattern of the multipliers is obtained by means of the helper function nbreffect, which combines identifying maximum and miminum with the mapping of multmap. The arguments are the original data frame, dfs, the multiplier data frame for k-nearest neighbor weights, knn, the id variable, fipsid, the type of multiplier, EonNbrs, and the model, slx. This gives all the desired results with a one line command.

nbreffect(dfs,knn,fipsid,mult='EonNbrs',model='slx')

The spatial pattern shows quite a few differences with that for the queen weights, highlighting the importance of the structure of the spatial weights. The maximum location is now county 28081 (Lee) and the minimum is county 28157 (Wilkinson). More specialized map comparison techniques (such as a co-location map) can be used to compare the two patterns, but this is not further pursued here.

Note how the quantile map for knn only shows four categories, even though k was set to 5. This is due to the large number of ties, which results in the two bottom categories being collapsed. This is a common problem for quantile maps when there is insufficient variation in the variable considered (here, all multiples of 0.2).

Kernel Weights

Kernel weights are illustrated for a triangular function (the default) computed for a variable bandwidth (fixed=False) determined by the 10 k-nearest neighbors. With diagonal=True, the diagonal elements all equal one. The function i_multiplier is called with the kernel weights (wkern), model="kernel" and the same id (fipsid) as before.

The result data frame is kern.

wkern = weights.Kernel.from_dataframe(dfs,k=10,ids=idvar,fixed=False,diagonal=True)
kern = i_multipliers(wkern,model='kernel',id=fipsid)
kern.describe()

The direct effects are all equal to one by construction. Both the effect of the neighbors (row sum - without the diagonal) and effect on the neighbors (column sum - without the diagonal) show substantial variation, even though their means are the same. The effect of the neighbors ranges from 1.789 to 3.833, whereas the effect on the neighbors goes from 1.562 to 3.890, a slightly larger range.

The location of the extrema and the map for the individual multipliers can again be obtained by means of nbreffect. Since this now has to be done for EofNbrs as well as EonNbrs, both are put in a list and the results are produced by a simple loop, minimizing code repetition.

effects = ['EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,kern,fipsid,mult=eff,model='kernel')

The resulting respective maximum and minimum locations are quite distinct, as are the associated spatial patterns. A closer examination of the spatial (mis)match is left as an exercise.

Spatial Lag

The spatial spillovers in a model that includes a spatially lagged dependent variable Wy are determined by the structure of the inverse matrix \((I - \rho W)^{-1}\). This is illustrated here for the queen contiguity weights. Similar results can be obtained for the k-nearest neighbor weights.

The i_multipliers function needs as arguments the spatial weights (wq), a value for the spatial autoregressive coefficient (coef=0.5), the model type (lag) and the id. Descriptive statistics are provided by describe.

lag5 = i_multipliers(wq,coef=0.5,model='lag',id=fipsid)
lag5.describe()

The mean of Direct, 1.066, is what is typically referred to as ADI or average direct impact, whereas the mean of either EofNbrs or EonNbrs, 0.934, is the AII or average indirect impact. Since \(\rho = 0.5\), the average total impact (ATI) is \(1.0 / (1.0 - \rho)\), or 2.0. Clearly this equals the sum of 1.066 and 0.934.

The averages mask considerable individual variation. The range of variation is fairly small for the direct effects - from 1.046 to 1.096 - and the effect of neighbors - from 0.904 to 0.954 -, but is quite substantial for the effect on neighbors - from 0.532 to 1.482.

Details are again provided by means of the nbreffects function, now with model = "lag".

effects = ['Direct','EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,lag5,fipsid,mult=eff,model='lag')

The spatial pattern of the direct effects and effects from neighbors are the exact opposite. This reflects how the total effect is split between these two parts: larger direct effects imply smaller indirect effects and vice versa. On the other hand, the spatial pattern for the effect on neighbors is quite distinct. It shows great similarity (but is not identical) to the spatial pattern for the queen weights as such (SLX).

Interestingly, the magnitude of the spatial autoregressive coefficient, while it affects the mean impact measures, does not affect the spatial pattern of the multipliers, which remains exactly the same. This is because the spatial pattern is determined by the network structure in the weights and the coefficient is just a scaling factor.

This is illustrated by using 0.3 as the spatial autoregressive coefficient.

lag3 = i_multipliers(wq,coef=0.3,model='lag',id=fipsid)
lag3.describe()

effects = ['Direct','EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,lag3,fipsid,mult=eff,model='lag')

Negative Exponential Distance Function

The multiplier effects for a nonlinear SLX model can be computed in the same way as before, from the diagonals, row sums and column sums of the parameterized weights matrix \(W(\alpha)\). Two cases are considered: a negative exponential distance transformation and an inverse power distance transformation. The input for the estimation of this model, in spreg.NSLX, is a sparse CSR matrix of distance measures that are standardized with respect to the bandwidth for each observation (fixed or variable). This ensures a well-behaved distance metric and a distance decay process.

For the negative exponential transformation in the nonlinear SLX model, the transformation is \(e^{-\alpha z_{ij}}\), where \(z_{ij} = d_{ij} / d_{bw}\) for \(d_{ij} \le d_{bw}\), and 0 otherwise, with \(d_{bw}\) as the bandwidth.

The input is a sparse CSR matrix, constructed by means of make_wnslx with the parameter exponential. As input it takes a numpy array of coordinates (see above), a tuple of parameters setting the number of nearest neighbors (here 10), whether the bandwidth is adaptive (np.inf) or fixed (a value), and the type of transformation (exponential). Since the coordinates are lat-lon decimal degrees, distance_metric = "Arc" to compute great circle distances. The result is a sparse CSR matrix. Its contents can be shown by means of toarray, which turns it into a regular (full) numpy array. This is just a view of the array and does not affect its sparse representation.

Since the transformed distance input is the fraction of the bandwidth, it results in smaller values closer to the origin and a value of 1.0 for the farthest nearest neighbor.

wexp = make_wnslx(coords,params=(10,np.inf,"exponential"),leafsize=30,distance_metric='Arc')
print(type(wexp))
wexp.toarray()[0:1,:]

Some descriptive statistics and the spatial pattern are obtained in the same way as before, using i_multipliers and nbreffect. The parameter is set to coef = 2.0 and the model is exponential.

exp2 = i_multipliers(wexp,coef=2.0,model='exponential',id=fipsid)
exp2.describe()

The mean multiplier of 2.648 ranges from 2.117 to 3.315 for the effect of neighbors, and from 1.368 to 3.871 for the effect on neighbors. The spatial distribution of the multipliers is again distinct.

effects = ['EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,exp2,fipsid,mult=eff,model='exponential')

Again, a different value of the coefficient affects the mean multiplier, but the spatial distribution remains the same. For example, with coef = 1.5, the mean multiplier becomes 3.636, reflecting the less steep distance decay that results from the smaller parameter.

This is important for the interpretation of the model. Unlike the other spatial models, a larger value of the coefficient in the nonlinear SLX models results in a smaller multiplier.

exp15 = i_multipliers(wexp,coef=1.5,model='exponential',id=fipsid)
exp15.describe()

effects = ['EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,exp15,fipsid,mult=eff,model='exponential')

Inverse Distance Power Function

For the inverse distance power function, the transformation is \(1.0 / d_{ij}^{\alpha}\). This is implemented as \(z_{ij}^{\alpha}\), where \(z_{ij}\) is guaranteed to be less than one, as \(z_{ij} = 1 - d_{ij}/d_{bw}\) for \(d_{ij} \le d_{bw}\), and zero otherwise. Hence, the distance measure for this transformation corresponds to a triangular kernel for the given bandwidth, with the diagonals zeroed out. The fraction of the bandwidth decreases as the distance increases, eventually yielding zero for the value of k.

Again, the corresponding weights are computed by means of the make_wnslx function (from spreg.utils). As input it takes a numpy array of coordinates (see above), a tuple of parameters setting the number of nearest neighbors (here 10), whether the bandwidth is variable (np.inf) or fixed (a value), and the type of transformation (power). Since the coordinates are lat-lon decimal degrees, distance_metric = "Arc" to compute great circle distances. The result is a sparse CSR matrix. Its contents can be shown by means of toarray.

winv = make_wnslx(coords,params=(10,np.inf,"power"),distance_metric='Arc')
print(type(winv))
winv.toarray()[0:1,:]

The sparse array winv is the spatial weights argument to the function i_multipliers. A coefficient of 2.0 is used as coef and the model is specified as power. The id is the same as before. Some summary characteristics are shown by means of describe.

pow2 = i_multipliers(winv,coef=2.0,model='power',id=fipsid)
pow2.describe()

The mean multiplier effect is 1.279, but it ranges from 0.628 to 2.124 for the effect of neighbors and from 0.709 to 1.917 for the effect on neighbors. Again, the mean masks considerable spatial variation. The direct effects are zero by construction.

The spatial pattern of the individual multipliers can be investigated in the same way as for the other multipliers by means of the nbreffects function.

effects = ['EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,pow2,fipsid,mult=eff,model='power')

Again, the spatial pattern of multipliers is not affected by the spatial parameter, but solely determined by the graph structure implied by the spatial weights before applying the transformation. This is illustrated with coef = 1.5. The mean effect of 1.902 is different and larger than for coef=2.0. This is due to the slower distance decay with the smaller coefficient, which gives larger weight to neighbors further away. The spatial pattern is unchanged.

pow15 = i_multipliers(winv,coef=1.5,model='power',id=fipsid)
pow15.describe()

effects = ['EofNbrs','EonNbrs']
for eff in effects:
    nbreffect(dfs,pow15,fipsid,mult=eff,model='power')

Practice

The similarities and differences between the different spatial layouts can be investigated more closely. In addition, the correlation between the multiplier vectors, respective locations of extrema and other characteristics can be examined. Of course, the same type of analysis can be investigated for other spatial weights and/or another data set.