nk027 / bsreg Goto Github PK

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