Zoom: 18
Longitudinal Analysis of the Brumandhino Data
Code for Reproducing Analyses in Section 6
This document contains code and explanations to reproduce the results in Section 6. Most of the source code is not printed, however, it can be easily accessed by looking at the code blocks in the Quarto file.
Construction of UTM-distance ``Gram” Matrix
We provide the code for constructing the ``Gram” kernel matrix, where each entry is given by \(\exp[- \{ \delta_{\mathrm{UTM}} (i, \, j) / (2 \cdot 200 )\}^2 ]\), where \(\delta_{\mathrm{UTM}} (i, \, j)\) is the distance in UTM units between locations \(i\) and \(j\). We begin by plotting the geographical locations superimposed on a map with UTM coordinates on the \(\mathrm{x}\)- and \(\mathrm{y}\)-axis.
The resulting network is obtained by fitting the graphical_lasso to the ``Gram” matrix constructed above. We plot the network below.
UTM_coords <- matrix(c(X, Y), nrow = 98, ncol =2)
UTM_dist <- function(x1, x2) {
a = as.numeric(t(x1 - x2) %*% (x1 - x2))
return(sqrt(a))
}
UTM_dist_matrix <- matrix(rep(0.0, 98*98), nrow = 98, ncol = 98)
#UTM_dist_matrix <- vapply(1:98, function(i){UTM_dist_inner(UTM_coords[i, ])}, rep(0.0, 98))
for (i in 1:98) {
UTM_dist_matrix[i, ] = vapply(1:98, function(j){UTM_dist(UTM_coords[i, ], UTM_coords[j, ])}, 0.0)
}
View(UTM_dist_matrix)
mean_dist = median(UTM_dist_matrix)
UTM_kernel = exp(- (UTM_dist_matrix / (2 * 200))^2)
hist(UTM_dist_matrix[UTM_dist_matrix > 0], breaks = 50)
net_mat <- glasso(s = UTM_kernel, rho = 0.12, penalize.diagonal = FALSE)
View(net_mat$wi)
b =net_mat$wi
b[abs(b) > 0.0] = 1
diag(b) = 0
net_result = graph_from_adjacency_matrix(b, 'undirected')
V(net_result)$Xcoord <- X
V(net_result)$Ycoord <- Y
V(net_result)$color <- rep('black', 98)
node_labels <- vapply(1:98, function(i) {return(paste0('(', as.character(i), ')'))}, '')
group_legend_labels <- c(vapply(1:3, function(i) {c(paste0('Location: ', as.character(i), ' & Satellite: 1'),
paste0('Location: ', as.character(i), ' & Satellite: 2'))}, c('', '')))
location <- matrix(c(V(net_result)$Xcoord, V(net_result)$Ycoord), ncol = 2)
plot.igraph(net_result, vertex.label='', vertex.label.cex = .5, vertex.color = 'black', edge.width = 0.5, vertex.label.color = 'black',
vertex.shape = 'circle', vertex.size = 500, mark.col = c('#ADD5D7', '#B3DBB6', '#FDB06C', '#7871A8', '#AACAF5', '#F6C5E2'),
mark.groups = list(c(1:10), c(11:20), c(21:37), c(72:86), c(38:60), c(c(61:71), c(87:98))), mark.border = NA,
layout = location, rescale=FALSE, asp = 0, xlim = range(V(net_result)$Xcoord), ylim = range(V(net_result)$Ycoord))
title(main = "\n Brumadinho Dam Location Network", cex.main=3)
legend("bottomright",
legend = group_legend_labels,
pch = 16,
text.col = 'black',
pt.cex = 2,
cex = 2,
bty = 'o',
horiz = FALSE,
col = c('#ADD5D7', '#B3DBB6', '#FDB06C', '#7871A8', '#AACAF5', '#F6C5E2')
)
Model Order Selection and Exploratory Analysis
We begin by specifying model order for each group. There are six groups in total, one for each satellite and location pair.
locind_edt <- c(rep(1, 10), rep(2, 10), rep(3, 17), rep(5, 23), rep(6, 11), rep(4, 15), rep(6, 12))
bur_netGNAR <- igraphtoGNAR(net_result)
S1 <- as.matrix(GNARtoigraph(bur_netGNAR))
D = distances(graph_from_adjacency_matrix(S1))
D[!is.finite(D)] <- 0
rmax = max(D)
W_bru_data <- fastCorbit::weights_matrix_cpp(bur_netGNAR, rmax)
# Set community membership values
loc1_sat1 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==1, 1, 0))}, 0)
loc1_sat2 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==2, 1, 0))}, 0)
loc2_sat1 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==3, 1, 0))}, 0)
loc2_sat2 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==4, 1, 0))}, 0)
loc3_sat1 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==5, 1, 0))}, 0)
loc3_sat2 <- vapply(1:98, function(x) {return(ifelse(locind_edt[x]==6, 1, 0))}, 0)
bru_data_community_indicators <- list(loc1_sat1, loc1_sat2, loc2_sat1, loc2_sat2, loc3_sat1, loc3_sat2)
W_bru_data <- fastCorbit::weights_matrix_cpp(bur_netGNAR, rmax)
#bru_data_interactions <- list(c(2), c(1), c(4), c(3), c(6), c(5))
bru_data_interactions <- list(c(2:6), c(1, 3:6), c(1,2, 4:6), c(1, 2, 3, 5, 6), c(1, 2, 3, 4, 6), c(1:5))
alpha_order <- rep(1, 6)
beta_order <- list(c(2), c(2), c(2), c(2), c(3), c(3)) # one for each communityWe plot a histogram of the centred longitudinal data and analyse Corbit and R-Corbit plots for verifying model order selection of one-lag and stage-three \(r > 1\).
bru_vts <- vapply(1:98, function(x) {return(100.0 * MXD2[, x] - colMeans(100.0 * MXD2)[x])}, rep(0.0, 44))
hist(recursive_dmat_cols_cpp(lapply(1:43, function(j) {return(as.matrix(bru_vts[2:44, j]))})), breaks = 50, xlab = 'y', main ='Bru VTS')
nacf_raw_data <- corbit_plot_cpp(MXD2, bur_netGNAR, 10, 4, W_bru_data)
max(abs(nacf_raw_data))[1] 0.2165149corbit_plot_cpp(MXD2, bur_netGNAR, 10, 4, W_bru_data, partial = "yes")
# Community-wise decomposition for each vector time series
loc1_sat1_nts <- t(vapply(seq(1:44), function(x) {return(loc1_sat1 * MXD2[x, ])}, rep(0.0, 98)))
loc1_sat2_nts <- t(vapply(seq(1:44), function(x) {return(loc1_sat2 * MXD2[x, ])}, rep(0.0, 98)))
loc2_sat1_nts <- t(vapply(seq(1:44), function(x) {return(loc2_sat1 * MXD2[x, ])}, rep(0.0, 98)))
loc2_sat2_nts <- t(vapply(seq(1:44), function(x) {return(loc2_sat2 * MXD2[x, ])}, rep(0.0, 98)))
loc3_sat1_nts <- t(vapply(seq(1:44), function(x) {return(loc3_sat1 * MXD2[x, ])}, rep(0.0, 98)))
loc3_sat2_nts <- t(vapply(seq(1:44), function(x) {return(loc3_sat2 * MXD2[x, ])}, rep(0.0, 98)))
# R-Corbit plots
# NACF
# PNACF
r_corbit_plot_cpp(list(loc1_sat1_nts, loc1_sat2_nts, loc2_sat1_nts, loc2_sat2_nts, loc3_sat1_nts, loc3_sat2_nts), list(bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR),
max_lag = 10, max_stage = 4, list(W_bru_data, W_bru_data, W_bru_data, W_bru_data, W_bru_data, W_bru_data), c("loc1_sat1", "loc1_sat2", "loc2_sat1", "loc2_sat2", "loc3_sat1", "loc3_sat2"),
same_net = "no", partial = "no")
r_corbit_plot_cpp(list(loc1_sat1_nts, loc1_sat2_nts, loc2_sat1_nts, loc2_sat2_nts, loc3_sat1_nts, loc3_sat2_nts), list(bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR, bur_netGNAR),
max_lag = 5, max_stage = 4, list(W_bru_data, W_bru_data, W_bru_data, W_bru_data, W_bru_data, W_bru_data), c("loc1_sat1", "loc1_sat2", "loc2_sat1", "loc2_sat2", "loc3_sat1", "loc3_sat2"),
same_net = "no", partial = "yes")
Fitting the Dynamic GNAR Model
We proceed to fit each of the models in the paper and compute the corresponding diagnostic statistics. We fill the table one model at a time and plot a histogram of the residuals and their evolution as a linear function of time: \(\hat{u}_{i, t_j} = \hat{b}_{0, i} + \hat{b}_{1, i} \, t_j\) for \(i = 1, \dots, 98\) and \(j = 1, \dots 43\).
Dynamic GNAR Fit
bru_struct_GNARfit <- temp_trend_cv_GNARfit(longitudinal_network_data = 100.0 * MXD2, network = bur_netGNAR, alpha_order = alpha_order, beta_order = beta_order, weight_matrix = W_bru_data,
covariate_groups = bru_data_community_indicators, interaction_groups = bru_data_interactions, rmax = rmax, iter = 100)
bru_res_interaction_fit <- community_interaction_gnar_fit(bru_struct_GNARfit$`GNAR residuals`, bur_netGNAR, alpha_order, beta_order, W_bru_data, bru_data_community_indicators,
bru_data_interactions)
summary(bru_res_interaction_fit)
Call:
lm(formula = yvec ~ . + 0, data = data.frame(full_model, check.names = FALSE))
Residuals:
Min 1Q Median 3Q Max
-1.75527 -0.15516 -0.00624 0.13564 1.70022
Coefficients:
Estimate Std. Error t value Pr(>|t|)
alpha.1.1 -0.217627 0.100756 -2.160 0.030834 *
beta.1.1.1 -0.108612 0.423653 -0.256 0.797679
beta.1.2.1 -3.535296 2.337807 -1.512 0.130553
`gamma.1.1.1:2` -0.653016 0.379118 -1.722 0.085061 .
`gamma.1.2.1:2` 1.071690 1.007424 1.064 0.287485
`gamma.1.1.1:3` 3.120264 4.812901 0.648 0.516819
`gamma.1.2.1:3` 1.465068 0.874484 1.675 0.093941 .
`gamma.1.1.1:4` 1.305999 0.874310 1.494 0.135318
`gamma.1.2.1:4` 0.853828 0.523373 1.631 0.102883
`gamma.1.1.1:5` -1.539020 1.119571 -1.375 0.169314
`gamma.1.2.1:5` -0.223745 0.371615 -0.602 0.547148
`gamma.1.1.1:6` 0.254417 0.370718 0.686 0.492574
`gamma.1.2.1:6` -0.084756 0.223129 -0.380 0.704075
alpha.1.2 -0.138525 0.079776 -1.736 0.082563 .
beta.1.1.2 -0.149034 0.307010 -0.485 0.627393
beta.1.2.2 -1.569313 3.432669 -0.457 0.647573
`gamma.1.1.2:1` -0.608544 0.536524 -1.134 0.256762
`gamma.1.2.2:1` -1.082090 1.178886 -0.918 0.358729
`gamma.1.1.2:3` 0.122130 3.039305 0.040 0.967949
`gamma.1.2.2:3` 1.719757 0.991566 1.734 0.082924 .
`gamma.1.1.2:4` 0.596137 1.118079 0.533 0.593938
`gamma.1.2.2:4` 1.054231 0.598440 1.762 0.078206 .
`gamma.1.1.2:5` -0.674963 0.546612 -1.235 0.216971
`gamma.1.2.2:5` -0.191122 0.376921 -0.507 0.612139
`gamma.1.1.2:6` -0.911081 0.751328 -1.213 0.225342
`gamma.1.2.2:6` 0.345998 0.231953 1.492 0.135861
alpha.1.3 -0.126933 0.058611 -2.166 0.030393 *
beta.1.1.3 0.197151 0.206364 0.955 0.339454
beta.1.2.3 -0.203946 0.376315 -0.542 0.587878
`gamma.1.1.3:1` -1.756788 4.069953 -0.432 0.666020
`gamma.1.2.3:1` -0.355103 0.869257 -0.409 0.682918
`gamma.1.1.3:2` 1.180206 2.301672 0.513 0.608147
`gamma.1.2.3:2` 1.641600 0.786922 2.086 0.037031 *
`gamma.1.1.3:4` 0.173343 0.205245 0.845 0.398400
`gamma.1.2.3:4` 1.177946 0.319651 3.685 0.000232 ***
`gamma.1.1.3:5` -1.273202 0.752295 -1.692 0.090641 .
`gamma.1.2.3:5` -0.751895 0.389051 -1.933 0.053349 .
`gamma.1.1.3:6` -0.132016 0.891304 -0.148 0.882259
`gamma.1.2.3:6` -0.375712 0.259051 -1.450 0.147040
alpha.1.4 -0.170437 0.039352 -4.331 1.52e-05 ***
beta.1.1.4 0.028643 0.180515 0.159 0.873936
beta.1.2.4 0.631960 0.341846 1.849 0.064577 .
`gamma.1.1.4:1` -1.057984 1.014276 -1.043 0.296966
`gamma.1.2.4:1` -0.566792 0.784789 -0.722 0.470199
`gamma.1.1.4:2` -0.661051 1.507000 -0.439 0.660936
`gamma.1.2.4:2` 1.408716 0.772014 1.825 0.068114 .
`gamma.1.1.4:3` 0.328711 0.234939 1.399 0.161848
`gamma.1.2.4:3` -1.179271 0.394486 -2.989 0.002812 **
`gamma.1.1.4:5` -0.844534 1.650371 -0.512 0.608872
`gamma.1.2.4:5` 0.574224 0.427860 1.342 0.179643
`gamma.1.1.4:6` 0.158536 0.462504 0.343 0.731784
`gamma.1.2.4:6` 0.716913 0.281971 2.543 0.011042 *
alpha.1.5 -0.294074 0.032088 -9.165 < 2e-16 ***
beta.1.1.5 -0.310992 0.167189 -1.860 0.062939 .
beta.1.2.5 -0.192966 0.278280 -0.693 0.488083
beta.1.3.5 0.495713 0.598034 0.829 0.407207
`gamma.1.1.5:1` -0.187466 1.929726 -0.097 0.922615
`gamma.1.2.5:1` 0.039132 0.804035 0.049 0.961185
`gamma.1.3.5:1` -1.797161 0.772484 -2.326 0.020042 *
`gamma.1.1.5:2` 0.496286 0.991493 0.501 0.616719
`gamma.1.2.5:2` 0.042224 0.516183 0.082 0.934810
`gamma.1.3.5:2` 1.399409 0.592714 2.361 0.018271 *
`gamma.1.1.5:3` -0.313440 0.959233 -0.327 0.743865
`gamma.1.2.5:3` -0.786476 0.744779 -1.056 0.291037
`gamma.1.3.5:3` 0.774299 0.257573 3.006 0.002662 **
`gamma.1.1.5:4` -0.773193 1.719006 -0.450 0.652885
`gamma.1.2.5:4` 0.864081 0.636750 1.357 0.174850
`gamma.1.3.5:4` -0.039980 0.245149 -0.163 0.870461
`gamma.1.1.5:6` 0.235900 0.127487 1.850 0.064329 .
`gamma.1.2.5:6` -0.500892 0.186159 -2.691 0.007160 **
`gamma.1.3.5:6` 0.189470 0.368161 0.515 0.606833
alpha.1.6 -0.109822 0.023253 -4.723 2.40e-06 ***
beta.1.1.6 -0.030101 0.123951 -0.243 0.808137
beta.1.2.6 -0.089611 0.197085 -0.455 0.649362
beta.1.3.6 1.192041 0.836806 1.425 0.154374
`gamma.1.1.6:1` 0.471160 1.515372 0.311 0.755877
`gamma.1.2.6:1` 0.446841 0.847100 0.527 0.597878
`gamma.1.3.6:1` -1.335547 0.855636 -1.561 0.118628
`gamma.1.1.6:2` -1.902277 1.114270 -1.707 0.087861 .
`gamma.1.2.6:2` 0.780845 0.600941 1.299 0.193889
`gamma.1.3.6:2` 0.887610 0.591148 1.502 0.133303
`gamma.1.1.6:3` -2.079832 2.456166 -0.847 0.397167
`gamma.1.2.6:3` 1.018692 0.743198 1.371 0.170547
`gamma.1.3.6:3` 0.454536 0.247947 1.833 0.066845 .
`gamma.1.1.6:4` -1.159294 1.060624 -1.093 0.274444
`gamma.1.2.6:4` -0.941848 0.607556 -1.550 0.121165
`gamma.1.3.6:4` -0.797861 0.239023 -3.338 0.000851 ***
`gamma.1.1.6:5` 0.005798 0.174871 0.033 0.973553
`gamma.1.2.6:5` 0.068947 0.267800 0.257 0.796838
`gamma.1.3.6:5` 0.762579 0.606387 1.258 0.208616
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3078 on 4124 degrees of freedom
Multiple R-squared: 0.08334, Adjusted R-squared: 0.06334
F-statistic: 4.166 on 90 and 4124 DF, p-value: < 2.2e-16#print(bru_struct_GNARfit$`GNAR parameters`)
print(sum(bru_struct_GNARfit$`GNAR residuals`[2:44, ]^2))[1] 426.2146print(sum(bru_struct_GNARfit$`Residuals Vector`^2) / (98 * 43 - 34))[1] 0.09343159print(sum(abs(bru_struct_GNARfit$`Residuals Vector`)) / (98 * 43 - 34))[1] 0.213178print(sqrt(sum(bru_struct_GNARfit$`GNAR parameters`)^2))[1] 4.03221bru_res_interaction_fit_ts <- bru_struct_GNARfit$`Residuals Matrix`
plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(bru_res_interaction_fit_ts[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_res_interaction_fit_ts[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_res_interaction_fit_ts[, j])))
}
l1_dgnar = l1_norm
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 1.12833482922229hist(bru_struct_GNARfit$`Residuals Vector`, breaks = 50, xlab = expression(hat(u)), main = "Dynamic GNAR Fit Residuals", cex = 5)
print(GNAR_AIC(bru_struct_GNARfit$`Residuals Vector`, 98, 90))[1] -2233.23GNAR_BIC(bru_struct_GNARfit$`Residuals Vector`, 98, 90)[1] -2229.544GNAR_HQC(bru_struct_GNARfit$`Residuals Vector`, 98, 90)[1] -2231.87nacf_vals <- corbit_plot_cpp(bru_struct_GNARfit$`Residuals Matrix`, bur_netGNAR, 10, 4, W_bru_data)
cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.0381276412373079nacf_dgnar = max(abs(nacf_vals))Global-\(\alpha\) GNAR Fit
bru_global_fit <- global_gnar_fit(bru_vts, bur_netGNAR, 1, 3, W_bru_data)
summary(bru_global_fit)
Call:
lm(formula = yvec ~ . + 0, data = data.frame(model_data, check.names = FALSE))
Residuals:
Min 1Q Median 3Q Max
-2.44405 -0.24554 -0.02023 0.20241 2.87097
Coefficients:
Estimate Std. Error t value Pr(>|t|)
alpha.1.1 0.11887 0.01660 7.159 9.55e-13 ***
beta.1.1.1 0.35687 0.04305 8.290 < 2e-16 ***
beta.1.2.1 0.17732 0.06149 2.884 0.00395 **
beta.1.3.1 0.24779 0.05221 4.746 2.15e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3956 on 4210 degrees of freedom
Multiple R-squared: 0.2754, Adjusted R-squared: 0.2748
F-statistic: 400.1 on 4 and 4210 DF, p-value: < 2.2e-16bru_global_fit_ts = residuals_matrix(bru_global_fit$residuals, 98)
print(sum((bru_global_fit$residuals)^2) / (98 * 43 - 4))[1] 0.1564628print(sum(abs(bru_global_fit$residuals)))[1] 1231.983print(sqrt(sum((bru_global_fit$coefficients)^2)))[1] 0.4840794print(sum(abs(bru_global_fit$coefficients)))[1] 0.9008567plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(bru_global_fit_ts[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_global_fit_ts[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_global_fit_ts[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 17.0080743279115l1_gnar = l1_norm
hist(bru_global_fit$residuals, breaks = 50, xlab = expression(hat(u)), main = "Global-alpha GNAR Fit Residuals", cex = 5)
GNAR_AIC(bru_global_fit$residuals, 98, 3)[1] -2153.53GNAR_BIC(bru_global_fit$residuals, 98, 3)[1] -2153.407GNAR_HQC(bru_global_fit$residuals, 98, 3)[1] -2153.485nacf_vals <- corbit_plot_cpp(residuals_matrix(bru_global_fit$residuals, 98), bur_netGNAR, 10, 4, W_bru_data)
cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.0408472856557266nacf_gnar = max(abs(nacf_vals))Community-\(\alpha\) with \(\gamma\)-Interactions GNAR Fit
bru_interaction_fit <- community_interaction_gnar_fit(bru_vts, bur_netGNAR, alpha_order, beta_order, W_bru_data, bru_data_community_indicators, bru_data_interactions)
summary(bru_interaction_fit)
Call:
lm(formula = yvec ~ . + 0, data = data.frame(full_model, check.names = FALSE))
Residuals:
Min 1Q Median 3Q Max
-2.3948 -0.2355 -0.0114 0.1851 2.5294
Coefficients:
Estimate Std. Error t value Pr(>|t|)
alpha.1.1 0.137919 0.094243 1.463 0.143423
beta.1.1.1 0.455615 0.374545 1.216 0.223883
beta.1.2.1 -2.762270 2.307254 -1.197 0.231293
`gamma.1.1.1:2` -0.471485 0.434907 -1.084 0.278382
`gamma.1.2.1:2` 0.154660 1.120174 0.138 0.890193
`gamma.1.1.1:3` 3.173497 4.227396 0.751 0.452877
`gamma.1.2.1:3` 0.440050 0.826787 0.532 0.594588
`gamma.1.1.1:4` 0.330297 0.686769 0.481 0.630583
`gamma.1.2.1:4` 0.838799 0.440635 1.904 0.057030 .
`gamma.1.1.1:5` -0.359820 0.952880 -0.378 0.705737
`gamma.1.2.1:5` -0.531165 0.325001 -1.634 0.102262
`gamma.1.1.1:6` 0.517150 0.381004 1.357 0.174749
`gamma.1.2.1:6` 0.181722 0.249981 0.727 0.467303
alpha.1.2 -0.006341 0.091471 -0.069 0.944739
beta.1.1.2 -0.257374 0.363954 -0.707 0.479507
beta.1.2.2 -3.018918 4.155262 -0.727 0.467556
`gamma.1.1.2:1` -0.092411 0.474827 -0.195 0.845699
`gamma.1.2.2:1` 0.535467 0.946405 0.566 0.571567
`gamma.1.1.2:3` -0.099799 3.156950 -0.032 0.974783
`gamma.1.2.2:3` 0.928748 0.960089 0.967 0.333423
`gamma.1.1.2:4` -0.850276 1.062122 -0.801 0.423442
`gamma.1.2.2:4` 0.832720 0.478883 1.739 0.082130 .
`gamma.1.1.2:5` -0.721942 0.549017 -1.315 0.188593
`gamma.1.2.2:5` -0.141881 0.314771 -0.451 0.652199
`gamma.1.1.2:6` -0.607897 0.692440 -0.878 0.380046
`gamma.1.2.2:6` 0.505504 0.246473 2.051 0.040334 *
alpha.1.3 0.093470 0.062913 1.486 0.137436
beta.1.1.3 0.232056 0.209966 1.105 0.269135
beta.1.2.3 0.066101 0.332199 0.199 0.842288
`gamma.1.1.3:1` -0.563799 4.684016 -0.120 0.904199
`gamma.1.2.3:1` -0.472984 0.853449 -0.554 0.579470
`gamma.1.1.3:2` 1.256189 2.857574 0.440 0.660250
`gamma.1.2.3:2` 1.103262 0.925338 1.192 0.233220
`gamma.1.1.3:4` 0.340138 0.132843 2.560 0.010489 *
`gamma.1.2.3:4` 1.111661 0.282800 3.931 8.60e-05 ***
`gamma.1.1.3:5` 0.030928 0.662306 0.047 0.962757
`gamma.1.2.3:5` -0.089873 0.374894 -0.240 0.810552
`gamma.1.1.3:6` -0.048243 0.987336 -0.049 0.961031
`gamma.1.2.3:6` 0.043203 0.296094 0.146 0.883999
alpha.1.4 0.063552 0.041965 1.514 0.129995
beta.1.1.4 0.672563 0.156389 4.301 1.74e-05 ***
beta.1.2.4 1.012814 0.287723 3.520 0.000436 ***
`gamma.1.1.4:1` -0.613113 1.007544 -0.609 0.542875
`gamma.1.2.4:1` -0.448865 0.774918 -0.579 0.562458
`gamma.1.1.4:2` -0.510754 1.846253 -0.277 0.782068
`gamma.1.2.4:2` 1.674823 0.906552 1.847 0.064751 .
`gamma.1.1.4:3` 0.209073 0.220566 0.948 0.343239
`gamma.1.2.4:3` 0.101526 0.353841 0.287 0.774184
`gamma.1.1.4:5` -1.300640 1.133887 -1.147 0.251422
`gamma.1.2.4:5` 1.744354 0.426866 4.086 4.46e-05 ***
`gamma.1.1.4:6` 0.027681 0.552429 0.050 0.960038
`gamma.1.2.4:6` 1.399633 0.310008 4.515 6.51e-06 ***
alpha.1.5 0.056121 0.032184 1.744 0.081277 .
beta.1.1.5 0.111035 0.149465 0.743 0.457595
beta.1.2.5 0.281726 0.237412 1.187 0.235433
beta.1.3.5 0.869508 0.548096 1.586 0.112722
`gamma.1.1.5:1` 0.523725 1.917181 0.273 0.784733
`gamma.1.2.5:1` -2.368492 0.703591 -3.366 0.000769 ***
`gamma.1.3.5:1` -0.213908 0.663471 -0.322 0.747160
`gamma.1.1.5:2` 2.045783 1.130307 1.810 0.070379 .
`gamma.1.2.5:2` 1.845753 0.611823 3.017 0.002570 **
`gamma.1.3.5:2` -0.061496 0.701557 -0.088 0.930154
`gamma.1.1.5:3` -0.009463 0.681254 -0.014 0.988918
`gamma.1.2.5:3` -0.078138 0.593890 -0.132 0.895330
`gamma.1.3.5:3` 0.173481 0.219821 0.789 0.430045
`gamma.1.1.5:4` 1.330444 1.466120 0.907 0.364217
`gamma.1.2.5:4` 0.869355 0.467074 1.861 0.062776 .
`gamma.1.3.5:4` 0.870340 0.188891 4.608 4.20e-06 ***
`gamma.1.1.5:6` 0.217921 0.142883 1.525 0.127292
`gamma.1.2.5:6` -0.223395 0.203677 -1.097 0.272789
`gamma.1.3.5:6` 0.362349 0.404919 0.895 0.370911
alpha.1.6 0.086177 0.025979 3.317 0.000917 ***
beta.1.1.6 -0.313388 0.138324 -2.266 0.023527 *
beta.1.2.6 0.289238 0.216964 1.333 0.182567
beta.1.3.6 0.379524 0.975295 0.389 0.697195
`gamma.1.1.6:1` 0.129383 1.443176 0.090 0.928568
`gamma.1.2.6:1` 0.279547 0.741563 0.377 0.706215
`gamma.1.3.6:1` 0.220300 0.707421 0.311 0.755503
`gamma.1.1.6:2` -1.928096 1.234860 -1.561 0.118509
`gamma.1.2.6:2` 2.180667 0.697567 3.126 0.001784 **
`gamma.1.3.6:2` -0.697551 0.677962 -1.029 0.303590
`gamma.1.1.6:3` -0.870495 2.129925 -0.409 0.682783
`gamma.1.2.6:3` 1.100028 0.549566 2.002 0.045390 *
`gamma.1.3.6:3` 0.620095 0.216037 2.870 0.004121 **
`gamma.1.1.6:4` -3.083422 0.743012 -4.150 3.39e-05 ***
`gamma.1.2.6:4` 0.238954 0.517428 0.462 0.644241
`gamma.1.3.6:4` 0.427215 0.186291 2.293 0.021883 *
`gamma.1.1.6:5` 0.479457 0.152078 3.153 0.001629 **
`gamma.1.2.6:5` 0.593187 0.240943 2.462 0.013859 *
`gamma.1.3.6:5` 0.755636 0.545649 1.385 0.166176
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3804 on 4124 degrees of freedom
Multiple R-squared: 0.3438, Adjusted R-squared: 0.3294
F-statistic: 24 on 90 and 4124 DF, p-value: < 2.2e-16bru_interaction_fit_ts = residuals_matrix(bru_interaction_fit$residuals, 98)
sum((bru_interaction_fit$residuals)^2) / (43 * 98 - 34)[1] 0.1427297sum(abs(bru_interaction_fit$residuals))[1] 1161.581sqrt(sum((bru_interaction_fit$coefficients)^2))[1] 9.490969sum(abs(bru_interaction_fit$coefficients))[1] 61.36982plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(bru_interaction_fit_ts[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_interaction_fit_ts[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), bru_interaction_fit_ts[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 10.3038446906965l1_intgnar = l1_norm
hist(bru_interaction_fit$residuals, breaks = 50, xlab = expression(hat(u)), main = "Interactions GNAR Fit Residuals", cex = 5)
GNAR_AIC(bru_interaction_fit$residuals, 98, 0)[1] -2165.865GNAR_BIC(bru_interaction_fit$residuals, 98, 30)[1] -2163.241GNAR_HQC(bru_interaction_fit$residuals, 98, 30)[1] -2164.017nacf_vals <- corbit_plot_cpp(residuals_matrix(bru_interaction_fit$residuals, 98), bur_netGNAR, 10, 4, W_bru_data)
cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.0260958947662737nacf_intgnar = max(abs(nacf_vals))Smoother Fit: 98 Independent Splines
mu_hat_mat <- temporal_trend_spline_smoother_cv(spline_residuals = bru_vts[2:44, ], vts_dimension = 98, plag = 1)
mu_res_mat <- bru_vts[2:44, ] - mu_hat_mat[2:44, ]
mu_hat_res <- recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(mu_res_mat[x, ]))}))
sum(mu_hat_res^2) / (98 * 43 - 1)[1] 0.2386304sum(abs(mu_hat_res)) / (98 * 43 - 1)[1] 0.361332plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(mu_res_mat[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), mu_res_mat[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), mu_res_mat[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 43.768348150643l1_splines = l1_norm
hist(mu_hat_res, breaks = 50, xlab = expression(hat(u)), main = "98 Smoothers Fit Residuals", cex = 5)
GNAR_AIC(mu_hat_res, 98, 0)[1] -2114.173GNAR_BIC(mu_hat_res, 98, 30)[1] -2111.549GNAR_HQC(mu_hat_res, 98, 30)[1] -2112.325nacf_vals <- corbit_plot_cpp(mu_res_mat, bur_netGNAR, 10, 4, W_bru_data)
cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.190624125508691nacf_splines = max(abs(nacf_vals))Sparse VAR Fit
sparse_var_fit <- fitVAR(100 * MXD2, p = 1, nlambda = 100)
theta_mat = sparse_var_fit$A[[1]]
sum(abs(sparse_var_fit$residuals))[1] 1067.933sparse_var_fit$lambda[1] 0.0004273116sparse_var_res <- sparse_var_fit$residuals[2:44, ]
sqrt(sum(theta_mat^2))[1] 2.816858b0 = theta_mat
b0[abs(b0) > 0.0] = 1
cat(paste0("Number of non-zero coefficients: ", sum(rowSums(b0))))Number of non-zero coefficients: 696sum((sparse_var_res)^2) / (98 * 43 - sum(rowSums(b0)))[1] 0.1223716plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(sparse_var_res[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), sparse_var_res[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), sparse_var_res[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 9.57597861897959l1_spvar = l1_norm
hist(recursive_dmat_cols_cpp(lapply(1:98, function(j) {return(as.matrix(sparse_var_fit$residuals[2:44, j]))})), breaks = 50, xlab = 'u', main ='Sparse VAR Fit Res', cex = 5)
GNAR_AIC(recursive_dmat_cols_cpp(lapply(1:98, function(j) {return(as.matrix(sparse_var_fit$residuals[2:44, j]))})), 98, sum(rowSums(b0)))[1] -2159.198GNAR_BIC(recursive_dmat_cols_cpp(lapply(1:98, function(j) {return(as.matrix(sparse_var_fit$residuals[2:44, j]))})), 98, sum(rowSums(b0)))[1] -2130.692GNAR_HQC(recursive_dmat_cols_cpp(lapply(1:98, function(j) {return(as.matrix(sparse_var_fit$residuals[2:44, j]))})), 98, sum(rowSums(b0)))[1] -2148.686nacf_vals <- corbit_plot_cpp(sparse_var_fit$residuals[2:44, ], bur_netGNAR, 10, 4, W_bru_data)
max(abs(nacf_vals))[1] 0.02236734cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.0223673381766903VAR_mat_plot(1, 1, theta_mat)
nacf_var = max(abs(nacf_vals))Null Model
naive_mean <- bru_vts[2:44, ]
naive_mean_vec <- recursive_dmat_cols_cpp(lapply(2:44, function(x) {return(as.matrix(bru_vts[x, ]))}))
sum(naive_mean_vec^2) / (98 * 43 - 1)[1] 0.2157901sum(abs(naive_mean_vec)) / (98 * 43 - 1)[1] 0.3429551plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(naive_mean[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), naive_mean[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), naive_mean[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 43.978445628209l1_null = l1_norm
hist(recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(naive_mean[x, ]))})), breaks = 50, xlab = expression(hat(u)), main = "Centred Data", cex = 5)
nacf_vals <- corbit_plot_cpp(bru_vts, bur_netGNAR, 10, 4, W_bru_data)
max(abs(nacf_vals))[1] 0.2165149cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.216514916405674nacf_null = max(abs(nacf_vals))Naive Model
naive_res <- 100 * t(vapply(2:44, function(x) {MXD2[x, 1:98] - MXD2[x - 1, 1:98]}, rep(0, 98)))
naive_res_vec <- recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(naive_res[x, ]))}))
sum((naive_res)^2) / (98 * 43 - 1)[1] 0.2645339sum(abs(naive_res)) / (98 * 43 - 1)[1] 0.3812074plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(naive_res[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), naive_res[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), naive_res[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 3.80606171851404l1_naive = l1_norm
hist(recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(naive_res[x, ]))})), breaks = 50, xlab = expression(hat(u)), main = "Naive Fit Residuals", cex = 5)
nacf_vals <- corbit_plot_cpp(naive_res, bur_netGNAR, 10, 4, W_bru_data)
max(abs(nacf_vals))[1] 0.1772237cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.177223660548852nacf_naive = max(abs(nacf_vals))98 ARs Fits
ar_residuals <- vapply(1:98, function(x){forecast::auto.arima(bru_vts[1:44, x], d = 0, D = 0, max.p = 1, max.q = 0,
max.P = 0, max.Q = 0, stationary = TRUE, seasonal = FALSE, ic = 'bic',
allowmean = FALSE, allowdrift = FALSE, trace = FALSE)$residuals}, rep(0, 44))[2:44, ]Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo ar_coeffs <- vapply(1:98, function(x){forecast::auto.arima(bru_vts[1:44, x], d = 0, D = 0, max.p = 1, max.q = 0, max.P = 0, max.Q = 0, stationary = TRUE, seasonal = FALSE, ic = 'bic',
allowmean = FALSE, allowdrift = FALSE, trace = FALSE)$coef[1]}, 1)
ar_coeffs[is.na(ar_coeffs)] = 1
ar_res_vec <- recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(ar_residuals[x, ]))}))
sum(ar_res_vec^2) / (98 * 42)[1] 0.1799183sum(abs(ar_res_vec)) / (98 * 42)[1] 0.3144478plot(rep(0, 43), type = 'l', col = 'blue', ylim = c(-2, 2), ylab = expression(hat(u)), xlab = 't', cex = 5)
l1_norm = 0
for (j in 1:98) {
points(c(1:43), c(ar_residuals[, j]), pch = 16)
lines(c(1:43), vapply(1:43, function(i) {c(1, i) %*% least_squares_solver(cbind(rep(1, 43), c(1:43)), naive_res[, j])}, 0.0), col = orange_pal[j],
lty = 2)
l1_norm = l1_norm + sum(abs(least_squares_solver(cbind(rep(1, 43), c(1:43)), ar_residuals[, j])))
}
cat(paste0("Sum of residual linear trend coefficients: ", as.character(l1_norm)))Sum of residual linear trend coefficients: 25.5885199670834l1_ar = l1_norm
hist(recursive_dmat_cols_cpp(lapply(1:43, function(x) {return(as.matrix(ar_residuals[x, ]))})), breaks = 50, xlab = expression(hat(u)), main = "98 ARs Fit Residuals", cex = 5)
nacf_vals <- corbit_plot_cpp(ar_residuals, bur_netGNAR, 10, 4, W_bru_data)
max(abs(nacf_vals))[1] 0.08554123cat(paste0("Maximum NACF Value: ", max(abs(nacf_vals))))Maximum NACF Value: 0.0855412264746516nacf_ar = max(abs(nacf_vals))Model Comparisons
The table below reproduces Table 1 and summarizes model comparisons.
| MSE | MAE | NACF | Est. Theta Norm | Res. Trend Coeff. | |
|---|---|---|---|---|---|
| Dynamic GNAR | 9.0824 | 20.7229 | 0.04 | 4.03 | 1.13 |
| Global-a GNAR | 15.3188 | 28.6508 | 0.04 | 0.90 | 17.01 |
| g-Interactions GNAR | 13.8747 | 27.0135 | 0.03 | 13.73 | 10.30 |
| 98 Smoothers | 23.3802 | 35.4021 | 0.19 | 0.00 | 43.77 |
| Sparse VAR | 10.0117 | 23.7546 | 0.02 | 2.82 | 9.58 |
| 98 AR | 17.2219 | 30.0992 | 0.09 | 7.56 | 25.59 |
| Naive Model | 25.9182 | 37.3495 | 0.18 | 0.00 | 3.81 |
| Null Model | 21.1424 | 33.6016 | 0.22 | 0.00 | 43.98 |
Temporal Trend
We plot the fitted temporal trend, which also accounts for autocorrelation in the data. The plot shows the main effect of this group-wise decreasing trends.
integer(0)
integer(0)
integer(0)
integer(0)
integer(0)
integer(0)