2 Distances
All of the above mentioned implementations employ the original econometric method that Conley introduced more than two decades ago. And at least with OLS as accompanying estimation strategy, there is no stochastic component to this standard error correction. Nonetheless, results differ across implementations. The root of these discrepancies is the heterogeneity in distance estimations. Computing spatial distances between observations is key to deriving Conley standard errors, and implementations vary in multiple aspects of this step. The disparities range from a slightly different Earth radius over alternative distance functions to unaligned definitions of space itself. Timothy Conley, Luis Calderon, and Leander Heldring do not restrict their functions to geographic space on Earth, but use a more abstract definition of space with an arbitrary number of dimensions and distance cutoffs delineating intervals along them. The flexibility of this abstract definition allows for a broader variety of applications, but comes with a number of drawbacks when applied to geographic space on Earth. It is somewhat restricted to planar projections, draws distance buffers in two-dimensional space as a rectangular box rather than a circle and faces a series of distortions and measurement errors arising from the planet’s spherical nature. The other packages and scripts focus on geographic distances on Earth, and largely rely on the Haversine equation. Haversine distances are not accurate down to the meter because the planet is not a perfect sphere. However, they are widely adopted in modern GIS software and are faster than more complex alternatives. Google’s S2 library, which the sf package now uses as default, employs that function - and conleyreg does so as well. The Haversine distance function handles spherical geometry and thereby does not apply to projected data. Unlike unprojected data, also referred to as longlat, projected data does not model locations as points on a sphere. It projects them to another spatial representation, usually onto a two-dimensional plain. Such transformation necessarily distorts the data in some regard. Hence, there are many different projections targeting specific purposes. A Mollweide equal-area projection generates equally large pixels. An azimuthal equidistant projection preserves the correct distances between the centroid and other points. Like the sf package, conleyreg uses Euclidean distances for projected data. Those computations do not correct for distortions introduced by the projection or the fact that the Earth is round and it is possible to leave a world map on one side and enter it on the other side. The latter would i.a. require calculating the optimal transition point for any projection, including the non-rectangular ones. Even though Haversine distances on longlat data are not accurate down to the meter, measurement error in distances computed from projected data can be much larger. Hence, it is recommended to use longlat data in conleyreg, unless you are certain that the selected projection gets the distances right.
By default, conleyreg uses Haversine distances on longlat input data and Euclidean distances on projected input data. You can change this behavior by specifying dist_comp = "spherical" or dist_comp = "planar", which refer to Haversine and Euclidean distances respectively. The functions then transform the data before computing distances. The resulting distance values are identical to those produced by sf::st_distance, unless you use sf’s old default (sf::sf_use_s2(FALSE)). In that case, sf employs the GEOS planar library. The difference between conleyreg’s internal distance functions and sf::st_distance is that the former are computationally highly optimized, making use of parallilization and the fact that only values below dist_cutoff have to be stored. You can still make conleyreg use the sf implementation via st_distance = TRUE. Though because of the lower efficiency, this only makes sense in specific applications - e.g. when you insist on using GEOS.
This package is built on memory efficiency, speed, and flexibility regarding applications and user preferences. Deriving distance matrices is the most time consuming task in calculating Conley standard errors. With a small data set of 5,000 observations, there are 25 million bilateral distances linking each point to all other points and itself. With a data set of 50,000 observations, the distance matrix contains 2.5 billion cells. conleyreg is highly optimized to run computations on even larger data at high speeds and efficient RAM usage. A considerable share of the internal code is written in C++, which is not just faster than R, but also allows to use more efficient data types than the ones available in R. While the exported functions conleyreg and dist_mat come with reasonable defaults, they offer a number of parameters with which you can tailor their behavior to your application. The subsequent paragraphs guide you through the options and explain what they imply in more detail.
2.1 Sparse Matrices
One feature that you can explicitly control in conleyreg is the use of sparse matrices. As the size of the distance matrix grows exponentially in the number of observations, it quickly surpasses your computer’s memory limit. A data set of 50,000 points implies 2.5 billion matrix cells. This is where sparse matrices come in handy. Instead of storing all cell values, they only store those that are not zero. Unless multiple points share a location, just \(1/N\) of all cells in a distance matrix - those along the diagonal - are zero. However, conleyreg does not directly store distances in those cells. The function transforms distances into weights based on the specified distance cutoff and kernel. Only cells with a distance below dist_cutoff are assigned a non-zero weight. If this radius is much smaller than the overall sample region, a vast majority of cells is set to zero - assuming an even distribution of points. In that case, a sparse matrix requires less RAM than a conventional dense matrix. Note that this advantage only holds with many zeros. Otherwise, sparse matrices can even be larger than their dense counterparts. We can illustrate the size differences using a global data set of 5,000 points combined with various distance cutoffs and the below described Bartlett kernel:
Dense and Sparse Matrix Sizes
The higher the distance cutoff is, the larger the sparse matrix becomes. At 15,000 km it is already substantially larger than the dense counterpart. This is admittedly not exactly what it looks like in conleyreg, as the package keeps the matrix by default in C++, while this plot compares dense and sparse matrices when exported to R. However, the general message still holds. Because sparse matrices are slower to manipulate than their dense counterpart and only offer a size advantage, if many distances are larger than dist_cutoff, conleyreg uses dense matrices by default. Specify sparse = TRUE to use sparse matrices instead.
Because sparse matrices’ compressed format with column-major ordering makes it slow to enter cell values one by one, conleyreg uses batch insertion by default when sparse = TRUE. With batch insertion, we loop through the \(N(N-1)/2\) unique bilateral distances (the matrix is symmetric - point A is as far from point B as point B is from A - and values along the diagonal - the distance between a point and itself - are zero), as we would do otherwise, but do not directly insert the values into the sparse matrix. Instead, we store the values and their indices in intermediate dense objects and then batch insert them into the sparse matrix all at once. That requires more RAM, but is much faster than direct insertion. Specify batch = FALSE to use direct insertion instead.
In the case of batch insertion, conleyreg allows to further play around with the trade-off between RAM efficiency and speed. The batch_ram_opt parameter controls the removal of the batch insertion’s intermediate objects. Permitted values are "none", "moderate" (default), and "heavy", where RAM efficiency increases and speed decreases from "none" to "heavy".
Speed of Element-Wise and Batch Insertion
The plot compares the computational speed of the three above described settings: element-wise insertion into a dense matrix, element-wise insertion into a sparse matrix, and batch insertion into a sparse matrix. In this case, the functions compute Haversine distances and convert them into weights based on the Bartlett kernel. Each derived matrix contains 100 million cells. Batch insertion uses batch_ram_opt = "moderate". All three functions are parallelized across eight CPU cores, but differ largely in their computational time. While the dense matrix function takes on average around two seconds to compute the 100 million cell matrix, the element-wise sparse matrix alternative takes up to 2.5 minutes to complete the exercise. If you choose to use sparse matrices nonetheless, it is adviseable to try batch insertion before utilizing the element-wise approach.
2.2 Data Types
When using conleyreg’s internal distance computations, i.e. st_distance = FALSE, it keeps the distance matrix or weights matrix in C++. Accordingly, the function can make use of data types that do not exist in R. With kernel = "uniform", the function stores values as short integers, each of which is two bytes in size. R, in contrast, only knows regular integers, measuring four bytes each. With kernel = "bartlett", you can choose between doubles (float = FALSE) and floats (float = TRUE). The former occupies eight bytes per value, as it would do in R, whereas the latter occupies four bytes per value. However, because floats store numbers at lower precision than doubles do, this choice can affect the results. In conclusion, you should only use floats, if your machine is unable to accomodate both a dense (sparse = FALSE) and a sparse matrix (sparse = TRUE). The above mentioned 2.5 billion matrix cells, thus, occupy 20GB when using doubles, 10GB when using floats, and 5GB when using short integers in a dense matrix. Sparse matrices can use the same data types, but the overall matrix size depends on how many values are below dist_cutoff.
Whereas conleyreg directly stores the weights, without storing the distances they are based on, the dist_mat function does return the actual distance values. And because these are exported to R, the float data type is not available. By specifying dist_round = TRUE, you can, however, round the distances to full kilometers. The resulting cell values then occupy four instead of eight bytes each.
2.3 Row-Wise Distances
As the above paragraphs illustrate, this package offers a number of options with which you can control distance matrices’ RAM requirements. And irrespective of whether you use sparse matrices, the float data type, or rounded distances, the functions loop over \(N(N-1)/2\) unique bilateral distances, given symmetry and zeros along the diagonal. The difference between \(N(N-1)/2\) and \(N^2\) can be billions of iterations in the case of larger data sets. Making use of this speed advantage implies computing the full distance matrix before inserting its values into the standard error corrections.
Dropping the \(N(N-1)/2\) technique, we can avoid storing the entire distance matrix. The spatial sandwich derivations do not make use of the full distance matrix at once anyway, but loop over its rows. By specifying rowwise = TRUE, you instruct the function to directly insert rows into the subsequent computations. The function then only has to store \(N \times\) ncores distance matrix cells at once. Though, because it discards rows after use, it has to derive \(N(N-1)\) distances - twice as many as in the default rowwise = FALSE application.
The rowwise parameter is one out of various options in conleyreg targeting a trade-off between computational speed and RAM efficiency. It can e.g. be useful when the machine’s RAM space is insufficient but the dist_cutoff is too high for sparse matrices to offer a size advantage.
2.4 Pre-Computed Distances
When estimating a model through conleyreg, it is usually the computationally best solution to let the function derive distances internally. Yet, there are a few cases in which you might want to pre-compute distance matrices. One example is an application where you estimate a number of models on the exact same data, e.g. trying a number of different distance cutoffs. Then, it is inefficient to have the function re-estimate the same distances in each specification. Instead, you can compute distances once and then pass them to conleyreg via the dist_mat parameter. While you can enter distances derived through functions outside this package, using conleyreg’s dist_mat function is the recommended way. It does not only make sure that the distances are stored in the right format, but is also computationally optimized and pre-computes other results that are constant across estimations on the same data. Employing functions outside the conleyreg package is mostly appropriate when requiring a distance function that the package does not include - e.g. Tobler’s hiking function, which is based on topography and expresses distances in terms of travel time.
In general, you must not modify the input data between the distance estimation and running conleyreg. That includes dropping observations or changing values of the unit, time, or coordinate variables. In cross-sectional settings, you must not re-order rows either. If you compute distances through a function other than dist_mat, there are a few additional issues to account for. (i) In the panel scenario, you must order observations by time and unit in ascending order. I.e. cells [1, 2] and [2, 1] of the distance matrix must refer to the distance between unit 1 and unit 2, cells [2, 3] and [3, 2] to the distance between unit 2 and unit 3 etc. The unit numbers in this example refer to their rank when they are sorted. (ii) dist_cutoff does not refer to kilometers, but to the units of the matrix. (iii) While in a balanced panel you only enter one matrix that is applied to all periods, you supply distances as a list of matrices in the unbalanced case. The matrices must be sorted, with the first list element containing the first period’s distance matrix etc. (iv) Zeros in sparse matrices are interpreted as values above the distance cutoff and NaN values are interpreted as zeros. (v) The matrices in the list must all be of the same type - all dense or all sparse. (vi) Distance matrices must only contain non-missing, finite numbers (and NaN in the case of sparse matrices).
Because the dist_mat function exports distances to R and does not keep them in C++, the float data type is not available. To reduce RAM requirements nonetheless, dist_mat offers a rounding option. With dist_round = TRUE, it rounds distances to full kilometers, storing them in a four byte integer rather than the default eight byte numeric format. Another way to reduce the output object’s size is to set a dist_cutoff combined with sparse = TRUE. Though note that the dist_cutoff used in subsequent conleyreg calls must not be larger than the dist_cutoff set in the dist_mat function.
