diff --git a/.RData b/.RData new file mode 100644 index 0000000..0f0066c Binary files /dev/null and b/.RData differ diff --git a/DESCRIPTION b/DESCRIPTION index f5bc1e3..3246b10 100644 --- a/DESCRIPTION +++ b/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, diff --git a/NAMESPACE b/NAMESPACE index 913b8f0..0b942be 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index b115a46..fc407c2 100644 --- a/NEWS.md +++ b/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. diff --git a/R/RcppExports.R b/R/RcppExports.R index 9f5c35e..bffa37b 100644 --- a/R/RcppExports.R +++ b/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) } diff --git a/R/hsar.R b/R/hsar.R index 06a742a..05455cd 100644 --- a/R/hsar.R +++ b/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) } diff --git a/R/mcmc_impacts.R b/R/mcmc_impacts.R new file mode 100644 index 0000000..4f3b957 --- /dev/null +++ b/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) +} diff --git a/R/sar.R b/R/sar.R index cfaed48..29fd513 100644 --- a/R/sar.R +++ b/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) } diff --git a/R/summary.hsar.R b/R/summary.hsar.R index 9448b78..7a53c88 100644 --- a/R/summary.hsar.R +++ b/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))) ) diff --git a/_pkgdown.yml b/_pkgdown.yml index 85cdfff..f25c270 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -21,6 +21,7 @@ reference: contents: - hsar - sar + - mcmc_impacts - title: Data contents: diff --git a/man/hsar.Rd b/man/hsar.Rd index 84d4d26..ae8ba64 100644 --- a/man/hsar.Rd +++ b/man/hsar.Rd @@ -10,6 +10,7 @@ hsar( W = NULL, M = NULL, Delta, + Durbin = FALSE, burnin = 5000, Nsim = 10000, thinning = 1, @@ -31,6 +32,8 @@ less than 1, M is also row-normalised before running the HSAR model. As with W, \item{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.} +\item{Durbin}{\code{logical}. Estimate Durbin model (i.e. include spatial lags of \code{X} as predictors)? Default \code{FALSE}.} + \item{burnin}{The number of MCMC samples to discard as the burnin period.} \item{Nsim}{The total number of MCMC samples to generate.} @@ -45,14 +48,19 @@ A \code{list}. \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.} diff --git a/man/mcmc_impacts.Rd b/man/mcmc_impacts.Rd new file mode 100644 index 0000000..b917ae2 --- /dev/null +++ b/man/mcmc_impacts.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mcmc_impacts.R +\name{mcmc_impacts} +\alias{mcmc_impacts} +\title{Estimate impact per MCMC draw of coefficients} +\usage{ +mcmc_impacts(object, p = 5) +} +\arguments{ +\item{object}{an \code{mcmc_hsar} or \code{mcmc_hsar_lambda_0} or \code{mcmc_sar} object} + +\item{p}{\code{numeric}. Maximum weight matrix power. Default \code{5}.} +} +\value{ +a list of three matrices named \code{direct}, \code{indirect}, and \code{total} giving impact estimates for each MCMC draw +} +\description{ +Estimate impact per MCMC draw of coefficients +} diff --git a/man/sar.Rd b/man/sar.Rd index d7401b7..f877185 100644 --- a/man/sar.Rd +++ b/man/sar.Rd @@ -8,6 +8,7 @@ sar( formula, data = NULL, W, + Durbin = FALSE, burnin = 5000, Nsim = 10000, thinning = 1, @@ -30,6 +31,8 @@ time. More specifically, W should be a column-oriented numeric sparse matrices o \code{Matrix} package. The converion between a dense numeric matrix and a sparse numeric matrix is made quite convenient through the \code{Matrix} library.} +\item{Durbin}{\code{logical}. Estimate Durbin model (i.e. include spatial lags of \code{X} as predictors)? Default \code{FALSE}.} + \item{burnin}{The number of McMC samples to discard as the burnin period.} \item{Nsim}{The total number of McMC samples to generate.} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dbd2c3f..5e4852c 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,8 +12,8 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // hsar_cpp_arma -List hsar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat M, arma::sp_mat Z, arma::mat detval, arma::mat detvalM, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float lambda_start, float sigma2e_start, float sigma2u_start, arma::vec betas_start); -RcppExport SEXP _HSAR_hsar_cpp_arma(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP MSEXP, SEXP ZSEXP, SEXP detvalSEXP, SEXP detvalMSEXP, SEXP UnumSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP lambda_startSEXP, SEXP sigma2e_startSEXP, SEXP sigma2u_startSEXP, SEXP betas_startSEXP) { +List hsar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat M, arma::sp_mat Z, arma::mat detval, arma::mat detvalM, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float lambda_start, float sigma2e_start, float sigma2u_start, arma::vec betas_start, bool Durbin); +RcppExport SEXP _HSAR_hsar_cpp_arma(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP MSEXP, SEXP ZSEXP, SEXP detvalSEXP, SEXP detvalMSEXP, SEXP UnumSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP lambda_startSEXP, SEXP sigma2e_startSEXP, SEXP sigma2u_startSEXP, SEXP betas_startSEXP, SEXP DurbinSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -33,13 +33,14 @@ BEGIN_RCPP Rcpp::traits::input_parameter< float >::type sigma2e_start(sigma2e_startSEXP); Rcpp::traits::input_parameter< float >::type sigma2u_start(sigma2u_startSEXP); Rcpp::traits::input_parameter< arma::vec >::type betas_start(betas_startSEXP); - rcpp_result_gen = Rcpp::wrap(hsar_cpp_arma(X, y, W, M, Z, detval, detvalM, Unum, burnin, Nsim, thinning, rho_start, lambda_start, sigma2e_start, sigma2u_start, betas_start)); + Rcpp::traits::input_parameter< bool >::type Durbin(DurbinSEXP); + rcpp_result_gen = Rcpp::wrap(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)); return rcpp_result_gen; END_RCPP } // hsar_cpp_arma_lambda_0 -List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat Z, arma::mat detval, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float sigma2e_start, float sigma2u_start, arma::vec betas_start); -RcppExport SEXP _HSAR_hsar_cpp_arma_lambda_0(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP ZSEXP, SEXP detvalSEXP, SEXP UnumSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP sigma2e_startSEXP, SEXP sigma2u_startSEXP, SEXP betas_startSEXP) { +List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat Z, arma::mat detval, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float sigma2e_start, float sigma2u_start, arma::vec betas_start, bool Durbin); +RcppExport SEXP _HSAR_hsar_cpp_arma_lambda_0(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP ZSEXP, SEXP detvalSEXP, SEXP UnumSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP sigma2e_startSEXP, SEXP sigma2u_startSEXP, SEXP betas_startSEXP, SEXP DurbinSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -56,7 +57,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< float >::type sigma2e_start(sigma2e_startSEXP); Rcpp::traits::input_parameter< float >::type sigma2u_start(sigma2u_startSEXP); Rcpp::traits::input_parameter< arma::vec >::type betas_start(betas_startSEXP); - rcpp_result_gen = Rcpp::wrap(hsar_cpp_arma_lambda_0(X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start)); + Rcpp::traits::input_parameter< bool >::type Durbin(DurbinSEXP); + rcpp_result_gen = Rcpp::wrap(hsar_cpp_arma_lambda_0(X, y, W, Z, detval, Unum, burnin, Nsim, thinning, rho_start, sigma2e_start, sigma2u_start, betas_start, Durbin)); return rcpp_result_gen; END_RCPP } @@ -83,9 +85,36 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// impact_cpp_arma +List impact_cpp_arma(arma::mat betas, arma::vec rhos, arma::sp_mat W); +RcppExport SEXP _HSAR_impact_cpp_arma(SEXP betasSEXP, SEXP rhosSEXP, SEXP WSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type betas(betasSEXP); + Rcpp::traits::input_parameter< arma::vec >::type rhos(rhosSEXP); + Rcpp::traits::input_parameter< arma::sp_mat >::type W(WSEXP); + rcpp_result_gen = Rcpp::wrap(impact_cpp_arma(betas, rhos, W)); + return rcpp_result_gen; +END_RCPP +} +// impact_Durbin_cpp_arma +List impact_Durbin_cpp_arma(arma::mat betas, arma::mat thetas, arma::vec rhos, arma::sp_mat W); +RcppExport SEXP _HSAR_impact_Durbin_cpp_arma(SEXP betasSEXP, SEXP thetasSEXP, SEXP rhosSEXP, SEXP WSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::mat >::type betas(betasSEXP); + Rcpp::traits::input_parameter< arma::mat >::type thetas(thetasSEXP); + Rcpp::traits::input_parameter< arma::vec >::type rhos(rhosSEXP); + Rcpp::traits::input_parameter< arma::sp_mat >::type W(WSEXP); + rcpp_result_gen = Rcpp::wrap(impact_Durbin_cpp_arma(betas, thetas, rhos, W)); + return rcpp_result_gen; +END_RCPP +} // sar_cpp_arma -List sar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::mat detval, int burnin, int Nsim, int thinning, float rho_start, float sigma2e_start, arma::vec betas_start); -RcppExport SEXP _HSAR_sar_cpp_arma(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP detvalSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP sigma2e_startSEXP, SEXP betas_startSEXP) { +List sar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::mat detval, int burnin, int Nsim, int thinning, float rho_start, float sigma2e_start, arma::vec betas_start, bool Durbin); +RcppExport SEXP _HSAR_sar_cpp_arma(SEXP XSEXP, SEXP ySEXP, SEXP WSEXP, SEXP detvalSEXP, SEXP burninSEXP, SEXP NsimSEXP, SEXP thinningSEXP, SEXP rho_startSEXP, SEXP sigma2e_startSEXP, SEXP betas_startSEXP, SEXP DurbinSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -99,16 +128,19 @@ BEGIN_RCPP Rcpp::traits::input_parameter< float >::type rho_start(rho_startSEXP); Rcpp::traits::input_parameter< float >::type sigma2e_start(sigma2e_startSEXP); Rcpp::traits::input_parameter< arma::vec >::type betas_start(betas_startSEXP); - rcpp_result_gen = Rcpp::wrap(sar_cpp_arma(X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start)); + Rcpp::traits::input_parameter< bool >::type Durbin(DurbinSEXP); + rcpp_result_gen = Rcpp::wrap(sar_cpp_arma(X, y, W, detval, burnin, Nsim, thinning, rho_start, sigma2e_start, betas_start, Durbin)); return rcpp_result_gen; END_RCPP } static const R_CallMethodDef CallEntries[] = { - {"_HSAR_hsar_cpp_arma", (DL_FUNC) &_HSAR_hsar_cpp_arma, 16}, - {"_HSAR_hsar_cpp_arma_lambda_0", (DL_FUNC) &_HSAR_hsar_cpp_arma_lambda_0, 13}, + {"_HSAR_hsar_cpp_arma", (DL_FUNC) &_HSAR_hsar_cpp_arma, 17}, + {"_HSAR_hsar_cpp_arma_lambda_0", (DL_FUNC) &_HSAR_hsar_cpp_arma_lambda_0, 14}, {"_HSAR_hsar_cpp_arma_rho_0", (DL_FUNC) &_HSAR_hsar_cpp_arma_rho_0, 13}, - {"_HSAR_sar_cpp_arma", (DL_FUNC) &_HSAR_sar_cpp_arma, 10}, + {"_HSAR_impact_cpp_arma", (DL_FUNC) &_HSAR_impact_cpp_arma, 3}, + {"_HSAR_impact_Durbin_cpp_arma", (DL_FUNC) &_HSAR_impact_Durbin_cpp_arma, 4}, + {"_HSAR_sar_cpp_arma", (DL_FUNC) &_HSAR_sar_cpp_arma, 11}, {NULL, NULL, 0} }; diff --git a/src/diagnostics.cpp b/src/diagnostics.cpp index 9c0acba..677eb41 100644 --- a/src/diagnostics.cpp +++ b/src/diagnostics.cpp @@ -1,6 +1,8 @@ #include "diagnostics.h" // Calculate the summary measures of the impacts (using aproximation for inv(SW) ) +// TODO 2025-05-27: pass in SW as an argument to make mcmc_impacts faster, I guess +// TODO 2025-05-27: consider increasing number of powers of rho * W for accuracy void diagnostic_impacts(const mat& betas, double rho, const sp_mat& W, mat& direct, mat& indirect, mat& total) { @@ -17,6 +19,31 @@ void diagnostic_impacts(const mat& betas, double rho, const sp_mat& W, indirect = total - direct; } +// 2025-05-22: as a non-C++ programmer I don't usually think about scope this way -GA +void diagnostic_impacts_Durbin(const mat& betas, const mat& thetas, double rho, const sp_mat& W, + mat& direct, mat& indirect, mat& total) { + + int n = W.n_rows; + + sp_mat I_sp = speye(n,n); + + sp_mat SW = I_sp+rho*W+pow(rho,2)*(W*W)+pow(rho,3)*(W*W*W)+pow(rho,4)*(W*W*W*W)+pow(rho,5)*(W*W*W*W*W); + + vec d( SW.diag() ); + + direct = sum(d)/n * betas; + total = accu(SW)/n * betas; + // indirect = total - direct; + + SW = W * SW; + + vec t( SW.diag() ); + + direct += sum(t)/n * thetas; + total += accu(SW)/n * thetas; + indirect = total - direct; +} + // Estimate Deviance Information Criterion (DIC) void diagnostic_dic_pd(vec log_likelihood_post_samples,double log_likelihood_mean_theta, double& DIC, double& pd) { diff --git a/src/diagnostics.h b/src/diagnostics.h index 8735a22..15bb6c2 100644 --- a/src/diagnostics.h +++ b/src/diagnostics.h @@ -8,6 +8,9 @@ using namespace arma; void diagnostic_impacts(const mat& betas, double rho, const sp_mat& W, mat& direct, mat& indirect, mat& total); +void diagnostic_impacts_Durbin(const mat& betas, const mat& thetas, double rho, const sp_mat& W, + mat& direct, mat& indirect, mat& total); + void diagnostic_dic_pd(vec log_likelihood_post_samples, double log_likelihood_mean_theta, double& DIC, double& pd); diff --git a/src/hsar.cpp b/src/hsar.cpp index 1a2019f..544bbb7 100644 --- a/src/hsar.cpp +++ b/src/hsar.cpp @@ -70,7 +70,8 @@ List hsar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat M, arma::sp_mat Z, arma::mat detval, arma::mat detvalM, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float lambda_start, - float sigma2e_start, float sigma2u_start, arma::vec betas_start){ + float sigma2e_start, float sigma2u_start, arma::vec betas_start, + bool Durbin){ //Starting the MCMC SET UP //arma_rng::set_seed(124); @@ -263,20 +264,39 @@ List hsar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat M, double r2 = diagnostic_Rsquared(y, y_hat_hsar(X, mean( Betas ), mean( rho ), W, Z ,mean( Us ) )); mat direct, indirect, total; - diagnostic_impacts( mean( Betas ), mean( rho ),W , - direct, indirect, total); + if(Durbin) { + int half_p = (p - 1) / 2; + + mat means = mean( Betas ); + + mat mean_betas = means.cols(1, half_p); + mat mean_thetas = means.cols(half_p + 1, p - 1); + + diagnostic_impacts_Durbin(mean_betas, mean_thetas, mean(rho), W, + direct, indirect, total); + } else { + mat mean_betas = mean( Betas ); + + diagnostic_impacts( mean_betas.cols(1, p - 1), mean( rho ),W , // 2025-05-27: don't include impact of `(Intercept)` + direct, indirect, total); + } return List::create(Named("cbetas")= Betas, Named("Mbetas")= mean( Betas ), Named("SDbetas") = stddev( Betas ), + Named("crho") = rho, // 2025-05-06: why not return the draws of all the things Named("Mrho")= mean( rho ), Named("SDrho") = stddev( rho ), + Named("clambda") = lambda, Named("Mlambda")= mean( lambda ), Named("SDlambda") = stddev( lambda ), + Named("csigma2e")= sigma2e, Named("Msigma2e")= mean( sigma2e ), Named("SDsigma2e") = stddev( sigma2e ), + Named("csigma2u")= mean( sigma2u ), Named("Msigma2u")= mean( sigma2u ), Named("SDsigma2u") = stddev( sigma2u ), + Named("cus")= Us, Named("Mus")= mean( Us ), Named("SDus") = stddev( Us ), Named("DIC") = dic, @@ -285,7 +305,8 @@ List hsar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat M, Named("R_Squared") = r2, Named("impact_direct") = direct, Named("impact_indirect") = indirect, - Named("impact_total") = total); + Named("impact_total") = total, + Named("W") = W); // 2025-05-27: return W so we can take samples of impacts later } diff --git a/src/hsar_lambda_0.cpp b/src/hsar_lambda_0.cpp index 4cc5cba..eb85634 100644 --- a/src/hsar_lambda_0.cpp +++ b/src/hsar_lambda_0.cpp @@ -71,7 +71,8 @@ List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, arma::sp_mat Z, arma::mat detval, arma::vec Unum, int burnin, int Nsim, int thinning, float rho_start, float sigma2e_start, - float sigma2u_start, arma::vec betas_start){ + float sigma2u_start, arma::vec betas_start, + bool Durbin){ //Starting the MCMC SET UP //arma_rng::set_seed(124); @@ -120,6 +121,8 @@ List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, mat invT0 = inv(T0); mat T0M0 = invT0 * M0; + + // throw std::logic_error("fixed"); mat tX = trans(X); mat Zfull(Z); mat tZ = trans(Zfull); @@ -246,17 +249,36 @@ List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, double r2 = diagnostic_Rsquared(y, y_hat_hsar_lambda_0(X, mean( Betas ), mean( rho ), W, Z ,mean( Us ) )); mat direct, indirect, total; - diagnostic_impacts( mean( Betas ), mean( rho ),W , direct, indirect, total); + if(Durbin) { + int half_p = (p - 1) / 2; + + mat means = mean( Betas ); + + mat mean_betas = means.cols(1, half_p); + mat mean_thetas = means.cols(half_p + 1, p - 1); + + diagnostic_impacts_Durbin(mean_betas, mean_thetas, mean(rho), W, + direct, indirect, total); + } else { + mat mean_betas = mean( Betas ); + + diagnostic_impacts( mean_betas.cols(1, p - 1), mean( rho ),W , // 2025-05-27: don't include impact of `(Intercept)` + direct, indirect, total); + } return List ::create(Named("cbetas")= Betas, Named("Mbetas")= mean( Betas ), Named("SDbetas") = stddev( Betas ), + Named("crho") = rho, Named("Mrho")= mean( rho ), Named("SDrho") = stddev( rho ), + Named("csigma2e") = sigma2e, Named("Msigma2e")= mean( sigma2e ), Named("SDsigma2e") = stddev( sigma2e ), + Named("csigma2u") = sigma2u, Named("Msigma2u")= mean( sigma2u ), Named("SDsigma2u") = stddev( sigma2u ), + Named("cus") = Us, Named("Mus")= mean( Us ), Named("SDus") = stddev( Us ), Named("DIC") = dic, @@ -265,7 +287,8 @@ List hsar_cpp_arma_lambda_0(arma::mat X, arma::vec y, arma::sp_mat W, Named("R_Squared") = r2, Named("impact_direct") = direct, Named("impact_idirect") = indirect, - Named("impact_total") = total); + Named("impact_total") = total, + Named("W") = W); // 2025-05-27: return W so we can take samples of impacts later } diff --git a/src/hsar_rho_0.cpp b/src/hsar_rho_0.cpp index a87bd3d..6a83faa 100644 --- a/src/hsar_rho_0.cpp +++ b/src/hsar_rho_0.cpp @@ -214,12 +214,16 @@ List hsar_cpp_arma_rho_0(arma::mat X, arma::vec y, arma::sp_mat M, return List::create(Named("cbetas")= Betas, Named("Mbetas")= mean( Betas ), Named("SDbetas") = stddev( Betas ), + Named("clambda") = lambda, Named("Mlambda")= mean( lambda ), Named("SDlambda") = stddev( lambda ), + Named("csigma2e") = sigma2e, Named("Msigma2e")= mean( sigma2e ), Named("SDsigma2e") = stddev( sigma2e ), + Named("csigma2u") = sigma2u, Named("Msigma2u")= mean( sigma2u ), Named("SDsigma2u") = stddev( sigma2u ), + Named("cus") = Us, Named("Mus")= mean( Us ), Named("SDus") = stddev( Us ), Named("DIC") = dic, diff --git a/src/impact.cpp b/src/impact.cpp new file mode 100644 index 0000000..71d9f74 --- /dev/null +++ b/src/impact.cpp @@ -0,0 +1,51 @@ +// #include // 2025-05-27: I'm not an RcppArmadillo coder -GA +#include +#include +#include "diagnostics.h" +#include "gibbs_method.h" + +// [[Rcpp::depends(RcppArmadillo)]] + +using namespace Rcpp; +using namespace arma; +using namespace std; + +// 2025-05-27: try as best we can to define funcs for taking impact estimates sample by sample + +// [[Rcpp::export]] + +List impact_cpp_arma(arma::mat betas, arma::vec rhos, arma::sp_mat W) { + mat direct, indirect, total; + mat directs, indirects, totals; + + for(int i = 0; i < rhos.n_elem; i++) { + diagnostic_impacts(betas.row(i), rhos[i], W, direct, indirect, total); + + directs.insert_rows(directs.n_rows, direct); + indirects.insert_rows(indirects.n_rows, indirect); + totals.insert_rows(totals.n_rows, total); + } + + return List::create(Named("direct") = directs, + Named("indirect") = indirects, + Named("total") = totals); +} + +// [[Rcpp::export]] + +List impact_Durbin_cpp_arma(arma::mat betas, arma::mat thetas, arma::vec rhos, arma::sp_mat W) { + mat direct, indirect, total; + mat directs, indirects, totals; + + for(int i = 0; i < rhos.n_elem; i++) { + diagnostic_impacts_Durbin(betas.row(i), thetas.row(i), rhos[i], W, direct, indirect, total); + + directs.insert_rows(directs.n_rows, direct); + indirects.insert_rows(indirects.n_rows, indirect); + totals.insert_rows(totals.n_rows, total); + } + + return List::create(Named("direct") = directs, + Named("indirect") = indirects, + Named("total") = totals); +} diff --git a/src/sar.cpp b/src/sar.cpp index be88179..812b46f 100644 --- a/src/sar.cpp +++ b/src/sar.cpp @@ -98,7 +98,10 @@ double SAR_draw_rho(const mat& detval, const mat& e0e0, List sar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::mat detval, int burnin, int Nsim, int thinning, - float rho_start, float sigma2e_start, arma::vec betas_start){ + float rho_start, float sigma2e_start, arma::vec betas_start, + bool Durbin){ + + // if(Durbin) {throw std::runtime_error("here");} //Starting the MCMC SET UP @@ -218,13 +221,30 @@ List sar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::mat detval, double r2 = diagnostic_Rsquared(y, y_hat_sar(X, mean( Betas ), mean( rho ),W )); mat direct, indirect, total; - diagnostic_impacts( mean( Betas ), mean( rho ),W , direct, indirect, total); + if(Durbin) { + int half_p = (p - 1) / 2; + + mat means = mean( Betas ); + + mat mean_betas = means.cols(1, half_p); + mat mean_thetas = means.cols(half_p + 1, p - 1); + + diagnostic_impacts_Durbin(mean_betas, mean_thetas, mean(rho), W, + direct, indirect, total); + } else { + mat mean_betas = mean( Betas ); + + diagnostic_impacts( mean_betas.cols(1, p - 1), mean( rho ),W , // 2025-05-27: don't include impact of `(Intercept)` + direct, indirect, total); + } return List::create(Named("cbetas")= Betas, Named("Mbetas")= mean( Betas ), Named("SDbetas") = stddev( Betas ), + Named("crho") = rho, Named("Mrho")= mean( rho ), Named("SDrho") = stddev( rho ), + Named("csigma2e") = sigma2e, Named("Msigma2e")= mean( sigma2e ), Named("SDsigma2e") = stddev( sigma2e ), Named("DIC") = dic, @@ -233,7 +253,8 @@ List sar_cpp_arma(arma::mat X, arma::vec y, arma::sp_mat W, arma::mat detval, Named("R_Squared") = r2, Named("impact_direct") = direct, Named("impact_indirect") = indirect, - Named("impact_total") = total); + Named("impact_total") = total, + Named("W") = W); // 2025-05-27: return W so we can take samples of impacts later }