Short Course
Regularised Estimation Procedures in Spatiotemporal Statistics
Cyclists in the city centre of Helsinki (https://ecf.com/news-and-events/news/visionarycities-series-will-helsinki-be-next-cycling-capital)
# time lag k = 1, ...
k <- 1
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station 5 minutes ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
# time lag k = 1, ...
k <- 2
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station 10 minutes ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
# time lag k = 1, ...
k <- 3
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station 15 minutes ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
# time lag k = 1, ...
k <- 4
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station 20 minutes ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
# time lag k = 1, ...
k <- 12
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station one hour ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
# time lag k = 1, ...
k <- 24
plot(avg_bikes[-c(1:k)], avg_bikes[-c((length(avg_bikes)-k+1):length(avg_bikes))],
col = color_vector[-c(1:k)],
xlab = "Average number of bikes per station",
ylab = "Average number of bikes per station two hours ago",
axes = FALSE)
axis(1, seq(0, 10, by = 0.5))
axis(2, seq(0, 20, by = 0.5))
Expand for full code
Expand for full code
bikes_space <- bikes[which(months(times) == "July"), ]
k <- 3
t <- 1
knnW <- t(sapply(1:n, function(i) ifelse(dist_mat[i,] < sort(dist_mat[i,])[k+2] & dist_mat[i,] > 0, 1/k, 0))) # k+2, because of the diagonal zero entry and the strict inequality
plot(bikes_space[t,],
knnW %*% ifelse(is.na(bikes_space[t,]), mean(bikes_space[t,], na.rm = TRUE), bikes_space[t,]),
ylim = c(0, 60),
xlim = c(0, 60),
ylab = "Average bike count of k-nearest neighbours",
xlab = "Average bike count midnight",
col = color_vector[t])
x <- knnW %*% ifelse(is.na(bikes_space[t,]), mean(bikes_space[t,], na.rm = TRUE), bikes_space[t,])
ab <- lm(bikes_space[t,] ~ x)$coeff
text(50, 50, ab[2])
abline(a = ab[1], b = ab[2])
h <- 1 # index time point across the day
plot(1:100, M.I,
col = colors[h],
main = paste("Spatial Correlation (Moran's I) at ", h-1, ":00", sep = ""),
pch = ifelse(M.Ip < 0.05, 20, 1),
ylim = c(0, 0.8),
axes = FALSE,
xlab = "Spatial lag order (k-nearest neighbours)",
ylab = "Spatial correlation (Moran's I)")
axis(1, seq(-20, 120, by = 10))
axis(2, seq(-1, 1, by = 0.05))
h <- 8 # index time point across the day
plot(1:100, M.I,
col = colors[h],
main = paste("Spatial Correlation (Moran's I) at ", h-1, ":00", sep = ""),
pch = ifelse(M.Ip < 0.05, 20, 1),
ylim = c(0, 0.8),
axes = FALSE,
xlab = "Spatial lag order (k-nearest neighbours)",
ylab = "Spatial correlation (Moran's I)")
axis(1, seq(-20, 120, by = 10))
axis(2, seq(-1, 1, by = 0.05))
h <- 13 # index time point across the day
plot(1:100, M.I,
col = colors[h],
main = paste("Spatial Correlation (Moran's I) at ", h-1, ":00", sep = ""),
pch = ifelse(M.Ip < 0.05, 20, 1),
ylim = c(0, 0.8),
axes = FALSE,
xlab = "Spatial lag order (k-nearest neighbours)",
ylab = "Spatial correlation (Moran's I)")
axis(1, seq(-20, 120, by = 10))
axis(2, seq(-1, 1, by = 0.05))
load("sources/data_berlin.rda")
library("maptools")
library("spdep")
B <- data_berlin$map
IDs <- as.character(names(B))
n <- length(unique(IDs))
plot(B)
IDs4 <- substr(as.character(names(B)), 1, 4)
n_thin <- length(unique(IDs4))
if(!gpclibPermitStatus()){
gpclibPermit()
}[1] TRUEB_thin <- unionSpatialPolygons(B, IDs4)
coords <- array(, dim = c(length(B), 2))
for(i in 1:length(B)){
coords[i,] <- B@polygons[[i]]@labpt
}
cols <- colorRampPalette(c("white", "darkblue"))
# Inverse-distance matrix
W <- 1 / as.matrix(dist(coords))
diag(W) <- 0
image(1:n, 1:n, W, col = cols(10))
# Inverse-distance matrix (row-standardised)
W <- W / array(apply(W, 1, sum), dim = dim(W))
image(1:n, 1:n, W, col = cols(10))
# Contiguity matrix (row-standardised)
W <- nb2mat(poly2nb(B), style = "W")
image(1:n, 1:n, W, col = cols(10))
W_list <- poly2nb(B)
plot(B, border = "lightgrey")
plot(B_thin, add = TRUE, lwd = 1)
plot(W_list, coordinates(B), pch = 1, add = TRUE, col = "darkblue",
lwd = 2, # 10*apply(W, 1, max),
cex = 0.2*apply(W, 2, function(x) sum(x>0)))
Implementation in R
Implementation in Python
# Load the required libraries
from sklearn.linear_model import lasso_path
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from itertools import cycle
# Load the sample data
cars = pd.read_csv('sources/mtcars.csv', index_col = 0, header = 0)
lambdas = np.logspace(-4.5, 0.9, 100)
# Split the data into features and target
X = cars.drop('mpg', axis=1)
y = cars['mpg']
# Standardize data
X -= X.mean(axis=0)
X /= X.std(axis=0)
# Use lasso_path to compute a coefficient path
alphas_lasso, coef_path, _ = lasso_path(X, y, alphas = lambdas)
# Plot coefficients against lambda
plt.figure();
colors = cycle(["b", "r", "g", "c", "k"])
for coefs, c in zip(coef_path, colors):
l1 = plt.plot(alphas_lasso, coefs, c = c);
plt.xlabel("Log Lambda");
plt.ylabel("Coefficients");
plt.xscale('log');
plt.title("Lasso Coefficient Path");
plt.axis("tight");
plt.show()
3.4 Elastic Net
Convex combination of both penalties: \[\begin{equation*} \hat{\boldsymbol{\beta}} = \arg \min_{\boldsymbol{\beta} \in \mathbb{R}^{p}} \frac{1}{n} (\boldsymbol{Y} - \mathbf{X} \, \boldsymbol{\beta})^\prime (\boldsymbol{Y} - \mathbf{X} \, \boldsymbol{\beta}) + \lambda_1 ||\beta||_1 + \lambda_2 ||\beta||_2 \end{equation*}\] with regularisation parameter \(\lambda \geq 0\).
- Grid-search over two-dimensional space \((\lambda_1, \lambda_2)'\)
- If \(\lambda_1 \rightarrow 0\), the \(\hat{\boldsymbol{\beta}}(\lambda_1, \lambda_2) \rightarrow \hat{\boldsymbol{\beta}}_{\text{Ridge}}(\lambda_2)\)
- If \(\lambda_2 \rightarrow 0\), the \(\hat{\boldsymbol{\beta}}(\lambda_1, \lambda_2) \rightarrow \hat{\boldsymbol{\beta}}_{\text{Lasso}}(\lambda_1)\)
- Elastic nets can be reduced to the linear support vector machine (i.e., SVM solvers can be applied)
library(glmnet)
# Load the sample data
data(mtcars)
# Convert the data to matrices
X <- as.matrix(mtcars[, -1])
y <- as.matrix(mtcars[, 1])
X <- scale(X)
# Fit the LASSO model
fit <- glmnet(X, y, intercept = TRUE, alpha = 0.5)
# alpha elastic net mixing parameter
# Plot the lambda path
plot(fit, xvar = "lambda", label = TRUE)
# Load the required libraries
from sklearn.linear_model import LassoCV
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# Load the sample data
cars = pd.read_csv('sources/mtcars.csv', index_col = 0, header = 0)
lambdas = np.logspace(-4.5, 0.9, 100)
# Split the data into features and target
X = cars.drop('mpg', axis=1)
y = cars['mpg']
# Standardize data
X -= X.mean(axis=0)
X /= X.std(axis=0)
# Fit LASSO model with cross-validation and k = 10 folds
model = LassoCV(alphas = lambdas, cv = 10).fit(X, y);
# Plot average MAE across the lambda sequence
plt.semilogx(model.alphas_, model.mse_path_, ":");
plt.semilogx(model.alphas_, model.mse_path_.mean(axis=1), 'k',
label=r'Average across the folds', linewidth = 2);
plt.xlabel(r'$\lambda$');
plt.ylabel('Mean squared error');
plt.title('Mean squared error on each fold: Lasso');
plt.axis('tight');
plt.show()
4.2 Dimensionality Reduction Covariance Matrix
- Big-\(N\) Problem \(\leadsto\) clever methods to reduce the dimensionality are needed to apply geostatistical models to large data data sets
Methods to reduce the dimensionality of the covariance matrix:
- Sparse covariance matrix (i.e., covariance matrix resulting from \(C_\theta\) contains many zeros)
- covariance tapering (Furrer, Genton, and Nychka 2006)
- theoretical results on covariance tapering can be found in Stein (2013)
- covariance tapering for likelihood-based estimation procedures (see Kaufman, Schervish, and Nychka 2008; Furrer, Bachoc, and Du 2016)
Regularised estimation procedures (with a shrinkage target of zero)
We can typically find zeros if two observations across space/time are conditionally independent, but conditional independence is encoded in the precision matrix (i.e., inverse covariance matrix):
- Lasso procedure to induce zeros in the precision matrix (see Krock, Kleiber, and Becker 2021 for univariate processes; and Krock et al. 2021 for multivariate processes)
- Graphical Lasso
library(gstat)
library(glasso)
set.seed(5515)
# Simulation of geostatistical process
xy <- expand.grid(1:10, 1:10)
names(xy) <- c("x","y")
gridded(xy) = ~x+y
g.dummy <- gstat(formula = z ~ 1, dummy = TRUE, beta = 0, model = vgm(1, "Exp", 15), nmax = 100) # for speed -- 10 is too small!!
yy <- predict(g.dummy, xy, nsim = 500)
spplot(yy[1:10])
# graphical Lasso
x <- t(do.call(cbind, lapply(1:500, function(x) yy[x]@data)))
s <- var(x)
lambda.seq <- 2^seq(from = -2, to = 0.5, length = 50)
a <- glassopath(s, rholist = lambda.seq)
# Plot the share of zeros in the estimated precision matrix
par(mar = c(5, 4, 4, 4) + 0.3)
plot(log(a$rholist), apply(a$wi, 3, function(x) mean(x == 0)),
type = "b", pch = 20,
xlab = "Log Lambda", ylab = "Sparsity degree (precision matrix)")
par(new = TRUE)
plot(log(a$rholist), apply(a$w, 3, function(x) mean(x == 0)),
type = "b", , pch = 20, col = "darkblue", axes = FALSE, ylab = "", xlab = "")
axis(side=4, at = pretty(range(apply(a$w, 3, function(x) mean(x == 0)))), col.axis = "darkblue")
mtext("Sparsity degree (covariance matrix)", side = 4, line = 3, col = "darkblue")
Ahrens, Achim, and Arnab Bhattacharjee. 2015. “Two-Step Lasso Estimation of the Spatial Weights Matrix.” Econometrics 3 (1): 128–55.
Banerjee, Sudipto, Alan E Gelfand, Andrew O Finley, and Huiyan Sang. 2008. “Gaussian Predictive Process Models for Large Spatial Data Sets.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70 (4): 825–48.
Bhattacharjee, Arnab, and Chris Jensen-Butler. 2013. “Estimation of the Spatial Weights Matrix Under Structural Constraints.” Regional Science and Urban Economics 43 (4): 617–34.
Boonstra, Philip S, Bhramar Mukherjee, and Jeremy MG Taylor. 2015. “A Small-Sample Choice of the Tuning Parameter in Ridge Regression.” Journal Article. Statistica Sinica 25 (3): 1185.
Chang, Ya-Mei, Nan-Jung Hsu, and Hsin-Cheng Huang. 2010. “Semiparametric Estimation and Selection for Nonstationary Spatial Covariance Functions.” Journal of Computational and Graphical Statistics 19 (1): 117–39.
Cressie, Noel, and Gardar Johannesson. 2008. “Fixed Rank Kriging for Very Large Spatial Data Sets.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70 (1): 209–26.
Cressie, Noel, Tao Shi, and Emily L Kang. 2010. “Fixed Rank Filtering for Spatio-Temporal Data.” Journal of Computational and Graphical Statistics 19 (3): 724–45.
Cressie, Noel, and Christopher K Wikle. 2011. Statistics for Spatio-Temporal Data. Wiley.
Elhorst, J Paul. 2010. “Applied Spatial Econometrics: Raising the Bar.” Spatial Economic Analysis 5 (1): 9–28. https://doi.org/10.1080/17421770903541772.
Fisher, Ronald A. 1935. “The Design of Experiments.” Oliver and Boyd, Edinburgh.
Furrer, Reinhard, François Bachoc, and Juan Du. 2016. “Asymptotic Properties of Multivariate Tapering for Estimation and Prediction.” Journal of Multivariate Analysis 149: 177–91.
Furrer, Reinhard, and Thomas Bengtsson. 2007. “Estimation of High-Dimensional Prior and Posterior Covariance Matrices in Kalman Filter Variants.” Journal of Multivariate Analysis 98 (2): 227–55.
Furrer, Reinhard, Marc G Genton, and Douglas Nychka. 2006. “Covariance Tapering for Interpolation of Large Spatial Datasets.” Journal of Computational and Graphical Statistics 15 (3): 502–23.
Gibbons, Stephen, and Henry G Overman. 2012. “Mostly Pointless Spatial Econometrics?” Journal of Regional Science 52 (2): 172–91.
Gonella, Romina, Mathias Bourel, and Liliane Bel. 2022. “Facing Spatial Massive Data in Science and Society: Variable Selection for Spatial Models.” Spatial Statistics, 100627.
Hsiang, TC. 1975. “A Bayesian View on Ridge Regression.” Journal of the Royal Statistical Society: Series D (The Statistician) 24 (4): 267–68.
Hsu, Nan-Jung, Ya-Mei Chang, and Hsin-Cheng Huang. 2012. “A Group Lasso Approach for Non-Stationary Spatial–Temporal Covariance Estimation.” Environmetrics 23 (1): 12–23.
Jiang, Gaoxia, and Wenjian Wang. 2017. “Markov Cross-Validation for Time Series Model Evaluations.” Journal Article. Information Sciences 375: 219–33. https://doi.org/https://doi.org/10.1016/j.ins.2016.09.061.
Kang, Myeongjong, and Matthias Katzfuss. 2021. “Correlation-Based Sparse Inverse Cholesky Factorization for Fast Gaussian-Process Inference.” arXiv Preprint arXiv:2112.14591.
Kaufman, Cari G, Mark J Schervish, and Douglas W Nychka. 2008. “Covariance Tapering for Likelihood-Based Estimation in Large Spatial Data Sets.” Journal of the American Statistical Association 103 (484): 1545–55.
Kolaczyk, Eric D. 2010. Statistical Analysis of Network Data: Methods and Models. Springer.
Kolaczyk, Eric D, and Gábor Csárdi. 2014. Statistical Analysis of Network Data with R. Springer.
Kostov, Philip. 2010. “Model Boosting for Spatial Weighting Matrix Selection in Spatial Lag Models.” Environment and Planning B: Planning and Design 37 (3): 533–49.
Kostov, Phillip. 2013. “Spatial Weighting Matrix Selection in Spatial Lag Econometric Model.” Econometrics 1 (1): 20–30.
Krock, Mitchell, William Kleiber, and Stephen Becker. 2021. “Nonstationary Modeling with Sparsity for Spatial Data via the Basis Graphical Lasso.” Journal of Computational and Graphical Statistics 30 (2): 375–89.
Krock, Mitchell, William Kleiber, Dorit Hammerling, and Stephen Becker. 2021. “Modeling Massive Highly-Multivariate Nonstationary Spatial Data with the Basis Graphical Lasso.” arXiv Preprint arXiv:2101.02404.
Lam, Clifford, and Pedro CL Souza. 2020. “Estimation and Selection of Spatial Weight Matrix in a Spatial Lag Model.” Journal of Business & Economic Statistics 38 (3): 693–710.
Liu, Xuan, Jianbao Chen, and Suli Cheng. 2018. “A Penalized Quasi-Maximum Likelihood Method for Variable Selection in the Spatial Autoregressive Model.” Spatial Statistics 25: 86–104.
Manski, Charles F. 1993. “Identification of Endogenous Social Effects: The Reflection Problem.” The Review of Economic Studies 60 (3): 531–42.
Merk, Miryam S, and Philipp Otto. 2021. “Directional Spatial Autoregressive Dependence in the Conditional First-and Second-Order Moments.” Spatial Statistics 41: 100490.
———. 2022. “Estimation of the Spatial Weighting Matrix for Regular Lattice Data—an Adaptive Lasso Approach with Cross-Sectional Resampling.” Environmetrics 33 (1): e2705.
Otto, Philipp, and Rick Steinert. 2022. “Estimation of the Spatial Weighting Matrix for Spatiotemporal Data Under the Presence of Structural Breaks.” Journal of Computational and Graphical Statistics, no. just-accepted: 1–44.
Park, Trevor, and George Casella. 2008. “The Bayesian Lasso.” Journal of the American Statistical Association 103 (482): 681–86.
Rue, Håvard, Andrea Riebler, Sigrunn H Sørbye, Janine B Illian, Daniel P Simpson, and Finn K Lindgren. 2017. “Bayesian Computing with INLA: A Review.” Annual Review of Statistics and Its Application 4: 395–421.
Saha, Arkajyoti, Sumanta Basu, and Abhirup Datta. 2021. “Random Forests for Spatially Dependent Data.” Journal of the American Statistical Association, 1–19.
Schäfer, Florian, Matthias Katzfuss, and Houman Owhadi. 2021. “Sparse Cholesky Factorization by Kullback–Leibler Minimization.” SIAM Journal on Scientific Computing 43 (3): A2019–46.
Stein, Michael L. 2013. “Statistical Properties of Covariance Tapers.” Journal of Computational and Graphical Statistics 22 (4): 866–85.
Stein, Michael L, Zhiyi Chi, and Leah J Welty. 2004. “Approximating Likelihoods for Large Spatial Data Sets.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66 (2): 275–96.
Tobler, Waldo R. 1970. “A Computer Movie Simulating Urban Growth in the Detroit Region.” Economic Geography 46 (sup1): 234–40.
Van Erp, Sara, Daniel L Oberski, and Joris Mulder. 2019. “Shrinkage Priors for Bayesian Penalized Regression.” Journal of Mathematical Psychology 89: 31–50.
Vecchia, Aldo V. 1988. “Estimation and Model Identification for Continuous Spatial Processes.” Journal of the Royal Statistical Society: Series B (Methodological) 50 (2): 297–312.
Zhu, Jun, Hsin-Cheng Huang, and Perla E Reyes. 2010. “On Selection of Spatial Linear Models for Lattice Data.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 (3): 389–402.
Zhu, Zhengyuan, and Yufeng Liu. 2009. “Estimating Spatial Covariance Using Penalised Likelihood with Weighted \(L_1\) Penalty.” Journal of Nonparametric Statistics 21 (7): 925–42.