Currently, the package only keep track of the parameters, but evidently it would be good to have a nice way of computing partial effects (which differ for the supported models and can be computationally intensive). Beyond the full partial effects, the total, direct, and indirect effect summaries (by LeSage & Pace, 2009) would be useful.
library(Matrix)
library(spdep)
library(bsreg)
library(tidyverse)
# Simulated data
n = 1000
k = 3
sige = 0.5
x <- matrix(rnorm(n*2*k), n, 2*k)
t <- 5
betas <- rep(1, times = 2*k) * 3
epsilon <- rnorm(n) * sige
rho = .8
xc = rnorm(n)
yc = rnorm(n)
coord <- cbind(xc, yc)
nb.nn <- knn2nb(knearneigh(coord, k = 5, longlat = TRUE))
W.nn <- nb2mat(nb.nn, style = "W", zero.policy = TRUE)
Identity = diag(1, nrow = n, ncol = n)
A = (Identity - rho*W.nn)
y = solve(A, x%*%betas) + solve(A, epsilon)
p = as.integer(ncol(x)/2)
# Bayesian SAR regression
n_save = 10000
n_burn = 500
out_bsar <- bsar(y ~ x, W = W.nn,
n_save = n_save, n_burn = n_burn,
ldet_SAR = list(reps = t))
tab_bsar <- as_tibble(out_bsar$draws) %>%
transmute(model = "BSAR",
x1 = beta2,
x2 = beta3,
x3 = beta4,
x4 = beta5,
x5 = beta6,
x6 = beta7,
lambda = lambda_SAR,
)
bsar_coeff <- tab_bsar %>%
as_tibble() %>%
summarize_at(vars(x1:lambda), mean, na.rm = TRUE) %>%
round(., 3)
# Effects estimates
W <- Matrix(W.nn, sparse = TRUE)
sz <- dim(out_bsar$draws)[1]
sx <- dim(out_bsar$draws)[2] - 3
uiter <- 50
maxorderu <- 100
nobs <- as.numeric(dim(x)[1])
rv <- matrix(rnorm(nobs*uiter), nrow = nobs, ncol = uiter)
tracew = matrix(0, maxorderu, 1)
wjjju = rv
rv <- Matrix(rv, sparse = TRUE)
wjjju <- Matrix(wjjju, sparse = TRUE)
for (jjj in 1:maxorderu) {
wjjju = W %*% wjjju
dp = rv * wjjju
tracew[jjj] = mean(mean(dp))
}
traces <- tracew
traces[1] <- 0
W_sparse_cp <- tcrossprod(t(W)*W)
traces[2] <- sum(sum(W_sparse_cp))/nobs
trs <- c(1, traces)
ntrs <- length(trs)
trbig <- trs
ndraw <- sz
p <- dim(x)[2]
bdraws <- out_bsar$draws[,2:7]
pdraws <- as.matrix(out_bsar$draws[,9])
ree <- 0:(ntrs - 1)
rmat <- rep(0, ntrs)
total <- array(0, c(sz, sx, ntrs))
direct <- array(0, c(sz, sx, ntrs))
indirect <- array(0, c(sz, sx, ntrs))
for (i in 1:sz) {
rmat <- pdraws[i,1]^ree
for (j in 1:p) {
beta = matrix(bdraws[i,j])
total[i,j,] <- beta[1,1]*rmat
direct[i,j,] <- c(beta * trbig) * rmat
indirect[i,j,] <- total[i,j,] - direct[i,j,]
}
}
# Direct effects
direct_out = matrix(data = NA, nrow = p, ncol = 4)
direct_save = matrix(data = NA, nrow = ndraw, ncol = p)
for (i in 1:p) {
tmp = direct[,i,]
direct_mean <- apply(tmp, 1, mean)
direct_std <- apply(tmp, 1, sd)
direct_sum <- (apply(tmp, 1, sum))
direct_save[,i] <- direct_sum
bounds <- quantile(direct_sum, c(0.025, 0.975))
cmean <- mean(direct_sum)
smean <- sd(direct_sum)
direct_out[i,] <- matrix(c(cmean, smean,
bounds[1], bounds[2]))
}
colnames(direct_out) <- c("Mean", "SD", "Lower Bound", "Upper Bound")
rownames(direct_out) <- c("x1", "x2", "x3",
"x4", "x5", "x6")
direct_out
# Indirect effects
indirect_out = matrix(data = NA, nrow = p, ncol = 4)
indirect_save = matrix(data = NA, nrow = ndraw, ncol = p)
for (i in 1:p) {
tmp = indirect[,i,]
indirect_mean <- apply(tmp, 1, mean)
indirect_std <- apply(tmp, 1, sd)
indirect_sum <- (apply(tmp, 1, sum))
indirect_save[,i] <- indirect_sum
#bounds <- quantile(direct_sum, c(0.025, 0.975))
bounds <- quantile(indirect_sum, c(0.025, 0.975))
cmean <- mean(indirect_sum)
smean <- sd(indirect_sum)
indirect_out[i,] <- matrix(c(cmean, smean,
bounds[1], bounds[2]))
}
colnames(indirect_out) <- c("Mean", "SD", "Lower Bound", "Upper Bound")
rownames(indirect_out) <- c("x1", "x2", "x3",
"x4", "x5", "x6")
indirect_out
# Total effects
total_out = matrix(data = NA, nrow = p, ncol = 4)
total_save = matrix(data = NA, nrow = ndraw, ncol = p)
for (i in 1:p) {
tmp = total[,i,]
total_mean <- apply(tmp, 1, mean)
total_std <- apply(tmp, 1, sd)
total_sum <- (apply(tmp, 1, sum))
total_save[,i] <- total_sum
bounds <- quantile(total_sum, c(0.025, 0.975))
cmean <- mean(total_sum)
smean <- sd(total_sum)
total_out[i,] <- matrix(c(cmean, smean,
bounds[1], bounds[2]))
}
colnames(total_out) <- c("Mean", "SD", "Lower Bound", "Upper Bound")
rownames(total_out) <- c("x1", "x2", "x3",
"x4", "x5", "x6")
total_out