Durbin = TRUE and more posteriors

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: HSAR
 Title: Hierarchical Spatial Autoregressive Model
-Version: 0.6.1
+Version: 0.7.3
 Authors@R: 
     c(
       person(given = "Guanpeng", family = "Dong",
@@ -38,6 +38,7 @@ Suggests:
     rmarkdown,
     sdsfun,
     sf,
+    stringr,
     tidyverse
 LinkingTo: 
     Rcpp,

NAMESPACE

@@ -9,5 +9,6 @@ S3method(summary,mcmc_hsar_lambda_0)
 S3method(summary,mcmc_hsar_rho_0)
 S3method(summary,mcmc_sar)
 export(hsar)
+export(mcmc_impacts)
 export(sar)
 useDynLib(HSAR, .registration = TRUE)

NEWS.md

@@ -1,3 +1,7 @@
+# HSAR 0.7.0
+
+* Add `Durbin` arg to `hsar`
+
 # HSAR 0.6.0
 
-* Resolve the CRAN check issues for `HASR` and resubmit it to CRAN.
+* Resolve the CRAN check issues for `HSAR` and resubmit it to CRAN.

R/RcppExports.R

@@ -1,19 +1,27 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-hsar_cpp_arma <- function(X, y, W, M, Z, detval, detvalM, Unum, burnin, Nsim, thinning, rho_start, lambda_start, sigma2e_start, sigma2u_start, betas_start) {
-    .Call(`_HSAR_hsar_cpp_arma`, X, y, W, M, Z, detval, detvalM, Unum, burnin, Nsim, thinning, rho_start, lambda_start, sigma2e_start, sigma2u_start, betas_start)
+hsar_cpp_arma <- function(X, y, W, M, Z, detval, detvalM, Unum, burnin, Nsim, thinning, rho_start, lambda_start, sigma2e_start, sigma2u_start, betas_start, Durbin) {
+    .Call(`_HSAR_hsar_cpp_arma`, X, y, W, M, Z, detval, detvalM, Unum, burnin, Nsim, thinning, rho_start, lambda_start, sigma2e_start, sigma2u_start, betas_start, Durbin)
 }
 
-hsar_cpp_arma_lambda_0 <- function(X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start) {
-    .Call(`_HSAR_hsar_cpp_arma_lambda_0`, X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start)
+hsar_cpp_arma_lambda_0 <- function(X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start, Durbin) {
+    .Call(`_HSAR_hsar_cpp_arma_lambda_0`, X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start, Durbin)
 }
 
 hsar_cpp_arma_rho_0 <- function(X, y, M, Z, detvalM, Unum, burnin, Nsim, thinning, lambda_start, sigma2e_start, sigma2u_start, betas_start) {
     .Call(`_HSAR_hsar_cpp_arma_rho_0`, X, y, M, Z, detvalM, Unum, burnin, Nsim, thinning, lambda_start, sigma2e_start, sigma2u_start, betas_start)
 }
 
-sar_cpp_arma <- function(X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start) {
-    .Call(`_HSAR_sar_cpp_arma`, X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start)
+impact_cpp_arma <- function(betas, rhos, W) {
+    .Call(`_HSAR_impact_cpp_arma`, betas, rhos, W)
+}
+
+impact_Durbin_cpp_arma <- function(betas, thetas, rhos, W) {
+    .Call(`_HSAR_impact_Durbin_cpp_arma`, betas, thetas, rhos, W)
+}
+
+sar_cpp_arma <- function(X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start, Durbin) {
+    .Call(`_HSAR_sar_cpp_arma`, X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start, Durbin)
 }
 

R/hsar.R

@@ -34,6 +34,7 @@
 #' less than 1, M is also row-normalised before running the HSAR model. As with W, M should also be a column-oriented numeric sparse matrices.
 #' @param Delta The N by J random effect design matrix that links the J by 1 higher-level random effect vector back to the N by 1 response variable under investigation. It is simply how lower-level units are grouped into each high-level units with columns of the matrix being each higher-level units. As with W and M, \eqn{\delta} should also be a
 #' column-oriented numeric sparse matrices.
+#' @param Durbin `logical`. Estimate Durbin model (i.e. include spatial lags of `X` as predictors)? Default `FALSE`.
 #' @param burnin The number of MCMC samples to discard as the burnin period.
 #' @param Nsim The total number of MCMC samples to generate.
 #' @param thinning MCMC thinning factor.
@@ -44,14 +45,19 @@
 #'  \item{cbetas}{A matrix with the MCMC samples of the draws for the coefficients.}
 #'  \item{Mbetas}{A vector of estimated mean values of regression coefficients. }
 #'  \item{SDbetas}{The standard deviations of estimated regression coefficients.}
+#'  \item{crho}{A vector with the MCMC samples of the draws for the lower-level spatial autoregressive parameter.}
 #'  \item{Mrho}{The estimated mean of the lower-level spatial autoregressive parameter \eqn{\rho}.}
 #'  \item{SDrho}{The standard deviation of the estimated lower-level spatial autoregressive parameter.}
-#'  \item{Mlamda}{The estimated mean of the higher-level spatial autoregressive parameter \eqn{\lambda}.}
+#'  \item{clambda}{A vector with the MCMC samples of the draws for the higher-level spatial autoregressive parameter.}
+#'  \item{Mlambda}{The estimated mean of the higher-level spatial autoregressive parameter \eqn{\lambda}.}
 #'  \item{SDlambda}{The standard deviation of the estimated higher-level spatial autoregressive parameter.}
+#'  \item{csigma2e}{A vector with the MCMC samples of the draws for the lower-level variance parameter.}
 #'  \item{Msigma2e}{The estimated mean of the lower-level variance parameter \eqn{\sigma^2_e}.}
 #'  \item{SDsigma2e}{The standard deviation of the estimated lower-level variance parameter \eqn{\sigma^{2}_{e} }.}
+#'  \item{csigma2u}{A vector with the MCMC samples of the draws for the higher-level variance parameter.}
 #'  \item{Msigma2u}{The estimated mean of the higher-level variance parameter \eqn{\sigma^2_u}.}
 #'  \item{SDsigma2u}{The standard deviation of the estimated higher-level variance parameter \eqn{\sigma^2_u}.}
+#'  \item{cus}{A matrix with the MCMC samples of the draws of \eqn{\theta}.}
 #'  \item{Mus}{Mean values of \eqn{\theta} }
 #'  \item{SDus}{Standard deviation of \eqn{\theta} }
 #'  \item{DIC}{The deviance information criterion (DIC) of the fitted model.}
@@ -128,7 +134,7 @@
 #'   palette <- RColorBrewer::brewer.pal(4, "Blues")
 #'   plot(Beijingdistricts,col=palette[groups],border="grey")
 #' }
-hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
+hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta, Durbin = FALSE,
                  burnin=5000, Nsim=10000, thinning=1, parameters.start = NULL) {
 
     ## check input data and formula
@@ -147,7 +153,21 @@ hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
     if( !is.null(W) ) check_matrix_dimensions(W,n,'Wrong dimensions for matrix W' )
     if( !is.null(M) ) check_matrix_dimensions(M,p,'Wrong dimensions for matrix M' )
 
-    Unum <- apply(Delta,2,sum)
+    if(length(Durbin) != 1)           {stop("`Durbin` must be either TRUE or FALSE")}
+    if(!(Durbin %in% c(TRUE, FALSE))) {stop("`Durbin` must be either TRUE or FALSE")}
+
+    Xlabel <- colnames(X)
+    Xnames <- setdiff(Xlabel, "(Intercept)")
+
+    # 2025-05-20: attach lag X if `Durbin == TRUE`
+    # TODO 2025-05-20: is it OK to treat factor dummies the same as any other column?
+    lag_X <- if(Durbin) {W %*% X[, Xnames]}
+    lag_X <- if(Durbin) {as.matrix(lag_X)}
+    lag_X <- if(Durbin) {`colnames<-`(lag_X, paste0("lag_", Xnames))}
+
+    X <- if(Durbin) {cbind(X, lag_X)} else {X}
+
+    Unum <- apply(Delta,2,sum) # TODO 2025-05-20: eventually allow for multiple membership somehow
 
     #start parameters
     if (! is.null(parameters.start)){
@@ -159,14 +179,15 @@ hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
         betas <- parameters.start$betas
         if (dim(X)[2]!= length(betas) ) stop("Starting values for Betas have got wrong dimension", call. = FALSE)
       }
-      else betas <- stats::coef(stats::lm(formula,data))
+      else betas <- stats::coef(stats::lm.fit(X, y)) # 2025-07-15: fixed this long ago in `sar.R` but not here -GA
     }
     else{
       rho <- 0.5
       lambda <- 0.5
       sigma2e <- 2.0
       sigma2u <- 2.0
-      betas <- stats::coef(stats::lm(formula,data))
+      # betas <- stats::coef(stats::lm(formula,data)) # 2025-07-03: omits betas for lag X vars if Durbin = TRUE -GA
+      betas <- stats::lm.fit(X, y)$coefficients # 2025-07-03: I thought I tried this before, though -GA # 2025-07-15: see above -GA
     }
 
     ## Call various models
@@ -181,7 +202,7 @@ hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
     # Special case where lamda =0 ; independent regional effect
     if ( is.null(M)){
       detval <- lndet_imrw(W)
-      result <- hsar_cpp_arma_lambda_0(X, y, W, Delta, detval, Unum, burnin, Nsim, thinning, rho, sigma2e, sigma2u, betas)
+      result <- hsar_cpp_arma_lambda_0(X, y, W, Delta, detval, Unum, burnin, Nsim, thinning, rho, sigma2e, sigma2u, betas, Durbin)
         #.Call("HSAR_hsar_cpp_arma_lambda_0", PACKAGE = 'HSAR', X, y, W, Delta, detval, Unum,
          #             burnin, Nsim, thinning, rho, sigma2e, sigma2u, betas)
       class(result) <- "mcmc_hsar_lambda_0"
@@ -190,7 +211,7 @@ hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
     if ( (!is.null(M)) & (!is.null(W))){
       detval <- lndet_imrw(W)
       detvalM <- lndet_imrw(M)
-      result <- hsar_cpp_arma(X, y, W, M, Delta, detval, detvalM, Unum, burnin, Nsim, thinning, rho, lambda, sigma2e, sigma2u, betas)
+      result <- hsar_cpp_arma(X, y, W, M, Delta, detval, detvalM, Unum, burnin, Nsim, thinning, rho, lambda, sigma2e, sigma2u, betas, Durbin)
         #.Call("HSAR_hsar_cpp_arma", PACKAGE = 'HSAR', X, y, W, M, Delta, detval, detvalM, Unum,
          #             burnin, Nsim, thinning, rho, lambda, sigma2e, sigma2u, betas)
       class(result) <- "mcmc_hsar"
@@ -200,9 +221,16 @@ hsar <- function(formula, data = NULL, W=NULL, M=NULL, Delta,
     result$Mbetas<-put_labels_to_coefficients(result$Mbetas, colnames(X))
     result$SDbetas<-put_labels_to_coefficients(result$SDbetas, colnames(X))
 
-    result$labels <- colnames(X)
+    if(!is.null(W)) { # 2025-05-27: name the impacts here and avoid complications arising from `Durbin == TRUE` downstream
+      result$impact_direct   <- put_labels_to_coefficients(result$impact_direct, Xnames)
+      result$impact_indirect <- put_labels_to_coefficients(result$impact_direct, Xnames)
+      result$impact_total    <- put_labels_to_coefficients(result$impact_direct, Xnames)
+    }
+
+    result$labels <- Xlabel
     result$call <- match.call()
     result$formula <- formula
+    result$Durbin <- Durbin
 
     return(result)
 }

R/mcmc_impacts.R

@@ -0,0 +1,59 @@
+#' Estimate impact per MCMC draw of coefficients
+#' @param object an `mcmc_hsar` or `mcmc_hsar_lambda_0` or `mcmc_sar` object
+#' @param p `numeric`. Maximum weight matrix power. Default `5`.
+#' @return a list of three matrices named `direct`, `indirect`, and `total` giving impact estimates for each MCMC draw
+#' @export
+mcmc_impacts <- function(object, p = 5) {
+  betas <- object$cbetas
+  betas <- as.matrix(betas)
+  betas <- betas[, setdiff(colnames(betas), "(Intercept)")]
+
+  if(object$Durbin) {
+    thetas <- betas[,  grepl("^lag_", colnames(betas)), drop = FALSE]
+    betas  <- betas[, !grepl("^lag_", colnames(betas)), drop = FALSE]
+    thetas <- `colnames<-`(thetas, stringr::str_remove(colnames(thetas), "^lag_"))
+  }
+
+  # take first p powers of wates
+  wates <- object$W
+  n     <- nrow(wates)
+  eyes  <- Matrix::Diagonal(n)
+  mates <- list(eyes)
+  w8pow <- eyes
+  i     <- 0
+  while(i < p) {
+    w8pow <- w8pow %*% wates
+    mates <- c(mates, list(w8pow))
+    i <- i + 1
+  }
+
+  # first p powers of each sample of rho, divided by n
+  rhopo <- outer(as.numeric(object$crho), 0:p, `^`) / n
+
+  diags <- lapply(mates, Matrix::diag)
+  disum <- sapply(diags, sum)
+  dirho <- rhopo %*% disum
+  direc <- betas * as.numeric(dirho)
+
+  # m8sum <- sapply(mates, sum) # these are all equal tho?
+  m8sum <- rep(n, p + 1)
+  m8rho <- rhopo %*% m8sum
+  total <- betas * as.numeric(m8rho)
+
+  # also impacts of lag vars
+  if(object$Durbin) {
+    mates <- lapply(mates, `%*%`, y = wates)
+
+    liags <- lapply(mates, Matrix::diag)
+    lisum <- sapply(liags, sum)
+    lirho <- rhopo %*% lisum
+    lirec <- thetas * as.numeric(lirho)
+    direc <- direc + lirec
+
+    lotal <- thetas * as.numeric(m8rho)
+    total <- total + lotal
+  }
+
+
+  list(direct = direc, indirect = total - direc, total = total)
+}

R/sar.R

@@ -19,6 +19,7 @@
 #' time. More specifically, W should be a column-oriented numeric sparse matrices of a `dgCMatrix` class defined in the
 #' `Matrix` package. The converion between a dense numeric matrix and a sparse numeric matrix is made quite convenient through
 #' the `Matrix` library.
+#' @param Durbin `logical`. Estimate Durbin model (i.e. include spatial lags of `X` as predictors)? Default `FALSE`.
 #' @param burnin The number of McMC samples to discard as the burnin period.
 #' @param Nsim The total number of McMC samples to generate.
 #' @param thinning MCMC thinning factor.
@@ -70,7 +71,7 @@
 #'            parameters.start=pars)
 #' summary(res)
 #' }
-sar <- function(formula, data = NULL, W,
+sar <- function(formula, data = NULL, W, Durbin = FALSE,
                 burnin = 5000, Nsim = 10000, thinning = 1,
                 parameters.start = NULL) {
 
@@ -82,10 +83,24 @@ sar <- function(formula, data = NULL, W,
     if (any(is.na(y))) stop("NAs in dependent variable", call. = FALSE)
     if (any(is.na(X))) stop("NAs in independent variable", call. = FALSE)
 
+    if(length(Durbin) != 1)           {stop("`Durbin` must be either TRUE or FALSE")}
+    if(!(Durbin %in% c(TRUE, FALSE))) {stop("`Durbin` must be either TRUE or FALSE")}
+
     n <- nrow(X)
 
     check_matrix_dimensions(W,n,'Wrong dimensions for matrix W' )
 
+    Xlabel <- colnames(X)
+    Xnames <- setdiff(Xlabel, "(Intercept)")
+
+    # 2025-05-20: attach lag X if `Durbin == TRUE`
+    # TODO 2025-05-20: is it OK to treat factor dummies the same as any other column?
+    lag_X <- if(Durbin) {W %*% X[, Xnames]}
+    lag_X <- if(Durbin) {as.matrix(lag_X)}
+    lag_X <- if(Durbin) {`colnames<-`(lag_X, paste0("lag_", Xnames))}
+
+    X <- if(Durbin) {cbind(X, lag_X)} else {X}
+
     detval <- lndet_imrw(W)
 
     #start parameters
@@ -101,10 +116,10 @@ sar <- function(formula, data = NULL, W,
     else{
       rho<-0.5
       sigma2e <-2.0
-      betas <- stats::coef(stats::lm(formula,data))
+      betas <- if(Durbin) {stats::coef(stats::lm.fit(X, y))} else {stats::coef(stats::lm(formula,data))}
     }
 
-    result <- sar_cpp_arma(X, y, W, detval, burnin, Nsim, thinning, rho, sigma2e, betas )
+    result <- sar_cpp_arma(X, y, W, detval, burnin, Nsim, thinning, rho, sigma2e, betas, Durbin )
       #.Call("HSAR_sar_cpp_arma", PACKAGE = 'HSAR', X, y, W, detval,
        #            burnin, Nsim, thinning, rho, sigma2e, betas )
 
@@ -114,9 +129,10 @@ sar <- function(formula, data = NULL, W,
     result$Mbetas<-put_labels_to_coefficients(result$Mbetas, colnames(X))
     result$SDbetas<-put_labels_to_coefficients(result$SDbetas, colnames(X))
 
-    result$labels <- colnames(X)
+    result$labels <- Xlabel
     result$call <-match.call()
     result$formula <- formula
+    result$Durbin <- Durbin
 
     return(result)
 }

R/summary.hsar.R

@@ -25,7 +25,7 @@ summary.mcmc_hsar <- function(object, ...)
   cat("\n Impacts:\n")
   df <- as.data.frame(cbind( t(x$impact_direct), t(x$impact_indirect), t(x$impact_total)))
   names(df)<-c("direct","indirect","total")
-  row.names(df)<- x$labels
+  row.names(df)<- setdiff(x$labels, "(Intercept)") # 2025-05-27: or just set the labels correctly here, it's fine
 
   print(df)
 
@@ -64,7 +64,7 @@ summary.mcmc_sar <- function(object, ...)
   cat("\n Impacts:\n")
   df <- as.data.frame(cbind( t(x$impact_direct), t(x$impact_indirect), t(x$impact_total)))
   names(df)<-c("direct","indirect","total")
-  row.names(df)<- x$labels
+  row.names(df)<- setdiff(x$labels, "(Intercept)") # 2025-05-27: or just set the labels correctly here, it's fine
 
   print( df )
 
@@ -132,6 +132,11 @@ summary.mcmc_hsar_lambda_0 <- function(object, ...)
   cat("Log likelihood:", x$Log_Likelihood, "\n")
   cat("Pseudo R squared:", x$R_Squared, "\n")
 
+  cat("\n Impacts:\n")
+  df <- as.data.frame(cbind( t(x$impact_direct), t(x$impact_indirect), t(x$impact_total)))
+  names(df)<-c("direct","indirect","total")
+  row.names(df)<- setdiff(x$labels, "(Intercept)") # 2025-05-27: or just set the labels correctly here, it's fine
+
   cat("\n Quantiles:\n")
   v <- c(0.05, 0.25, 0.5, 0.75, 0.95)
   print( t(sapply( x$cbetas, function(i) stats::quantile(i, v))) )

_pkgdown.yml

@@ -21,6 +21,7 @@ reference:
   contents:
   - hsar
   - sar
+  - mcmc_impacts
 
 - title: Data
   contents: