dcar_proper posterior appearing biased with covariates in monte carlo study

[ConnorDonegan]

Hello,

I've been working on proper CAR models in Stan, and part of my process has been to compare results to Nimble. Since the results differ across software, I've been using simulated data to check for biases in my models, and all the results for Nimble are turning out biased. Specifically regression coefficients within a proper car model. I'm very new to Nimble, so I know there's a good chance that I'm just doing something wrong here. I'd be glad to hear if I've made a mistake.

I've created a reproducible example based on the Nimble ICAR example (https://r-nimble.org/examples). I simulate spatially auto-correlated log-rates (using sim_sar, the SAR model) with a covariate, then draw some count data using rpois(). The data is simulated like this:

x <- sim_sar(w = W, rho = rho.x, sigma = 0.2)
phi <- sim_sar(w = W, mu = alpha + x * beta, rho = rho.y, sigma = 0.2)
y <- rpois(n = n, lambda = E * exp(phi))

Then I fit the CAR models and collect the posterior mean of the regression coefficient on each iteration, and also collect results from lm(log(y/E) ~ x), which should also be unbiased.

library(nimble)
library(CARBayesdata, quietly = TRUE)
library(sp, quietly = TRUE)
library(spdep, quietly = TRUE)
Nsim <- 30
#' draw from simultaneous autoregressive (SAR) model
sim_sar <- function (m = 1, mu = rep(0, nrow(w)), w, rho, sigma = 1, ...) {
    stopifnot(inherits(w, "matrix"))
    stopifnot(ncol(w) == nrow(w))
    N <- nrow(w)
    if (missing(mu)) {
        mu <- rep(0, N)
    }
    I <- diag(N)
    S <- sigma^2 * solve( (I - rho * t(w)) %*% (I - rho * w) )
    x <- MASS::mvrnorm(n = m, mu = mu, Sigma = S, ...)
    return(x)
}

data(GGHB.IG)
data(respiratorydata)
respiratorydata_spatial <- merge(x = GGHB.IG, y = respiratorydata, by = "IG", all.x = FALSE)

W.nb <- poly2nb(respiratorydata_spatial, row.names =  rownames(respiratorydata_spatial@data))
nbInfo <- nb2WB(W.nb)
CM <- as.carCM(nbInfo$adj, nbInfo$weights, nbInfo$num)
nregions <- nrow(respiratorydata_spatial)

code <- nimbleCode({
    # priors
    alpha ~ dnorm(0, 10) #/#
    beta ~ dnorm(0, sd = 100)
    sigma ~ dunif(0, 100)   # prior for variance components based on Gelman (2006)
    prec <- 1 / sigma^2
    rho ~ dunif(0, 1)
    # latent process
    s[1:N] ~ dcar_proper(mu[1:N],
                         C = C[1:nC],
                         adj = adj[1:L],
                         num = num[1:N],
                         M = M[1:N],
                         tau = prec, 
                         gamma = rho)    
    # likelihood
    for(i in 1:N) {
        mu[i] <- alpha + beta * x[i]
        lambda[i] <- E[i] * exp(s[i])
        y[i] ~ dpois(lambda[i])
    }
})

#/# prep sim data
W <- spdep::nb2mat(spdep::poly2nb(respiratorydata_spatial), style = "W")
rho.x <- 0.75
rho.y <- 0.6
beta <- -1
alpha <-  0
n <- nrow(W)
E <- respiratorydata_spatial$expected
inits <- list(alpha = 0, beta = -1, prec = 6, s = rnorm(n))
data <- list()
constants <- list(N = nregions,
                  L = length(nbInfo$adj),
                  C = CM$C,
                  M = CM$M,
                  nC = length(CM$C),
                  adj = nbInfo$adj,
                  num = nbInfo$num, #x = x,
                  E = E
                  )

res <- lapply(1:Nsim, function(i) {
    x <- sim_sar(w = W, rho = rho.x, sigma = 0.2)
    phi <- sim_sar(w = W, mu = alpha + x * beta, rho = rho.y, sigma = 0.2)
    y <- rpois(n = n, lambda = E * exp(phi))
    constants$x <- x 
    data$y <- y
    # par(mfrow=c(1,3)); hist(x); hist(y/E); plot(x, log(y/E))
    model <- nimbleModel(code, constants = constants, data = data, inits = inits)
    cModel <- compileNimble(model)
    conf <- configureMCMC(model, monitors = c('alpha', 'beta', 'sigma', 's'))
    MCMC <- buildMCMC(conf)
    cMCMC <- compileNimble(MCMC, project = cModel)
    samples <- runMCMC(cMCMC, niter = 15000, nburnin = 7500)
    f <- lm(log(y/E) ~ x)    
    c(nimble = mean(samples[,"beta"]), lm = coefficients(f)["x"])
}
)
r = do.call("rbind", res)
readr::write_csv(as.data.frame(r), "data/test-nimble-bias-mean-centered.csv")
#r <- readr::read_csv("data/test-nimble-bias-mean-centered.csv")

message("Sampling from log-linear Poisson model with beta = -1 (nsim = ",nrow(r), ")")
message("Sample mean of Nimble estimates: ", mean(r[,"nimble"]))
message("Sample mean of lm estimates: ", mean(r[,"lm.x"]))

par(mfrow = c(1,2))
hist(r[,"nimble"])
hist(r[,"lm.x"])

These are the results I've gotten:

> message("Sampling from log-linear Poisson model with beta = -1")
Sampling from log-linear Poisson model with beta = -1
> message("Sample mean of Nimble estimates: ", mean(r[,"nimble"]))
Sample mean of Nimble estimates: -0.801299816222076
> message("Sample mean of lm estimates: ", mean(r[,"lm.x"]))
Sample mean of lm estimates: -1.01871752123125

I also tried OpenBUGS, and found that results are different from Nimble but also appear biased. The amount of bias here is consistent with other results I've gotten. Obviously if I'm making a mistake here with Nimble, I could be making the same or similar mistake with OpenBUGS. However, when I do not use a covariate (intercept only model), the results from Nimble and Stan are identical.

Am I doing something wrong here? Thanks for your time!

> sessionInfo()
R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] spdep_1.1-8      sf_1.0-2         spData_0.3.10    sp_1.4-5        
[5] CARBayesdata_2.2 nimble_0.11.1   

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.7         compiler_4.1.1     pillar_1.6.2       LearnBayes_2.15.1 
 [5] class_7.3-19       shapefiles_0.7     tools_4.1.1        boot_1.3-28       
 [9] nlme_3.1-152       lifecycle_1.0.0    tibble_3.1.3       lattice_0.20-44   
[13] pkgconfig_2.0.3    rlang_0.4.11       Matrix_1.3-4       igraph_1.2.6      
[17] DBI_1.1.1          parallel_4.1.1     expm_0.999-6       e1071_1.7-8       
[21] coda_0.19-4        dplyr_1.0.7        raster_3.4-13      gtools_3.9.2      
[25] generics_0.1.0     vctrs_0.3.8        classInt_0.4-3     grid_4.1.1        
[29] tidyselect_1.1.1   glue_1.4.2         R6_2.5.0           fansi_0.5.0       
[33] foreign_0.8-81     gdata_2.18.0       deldir_0.2-10      purrr_0.3.4       
[37] magrittr_2.0.1     gmodels_2.18.1     splines_4.1.1      MASS_7.3-54       
[41] codetools_0.2-18   ellipsis_0.3.2     units_0.7-2        assertthat_0.2.1  
[45] KernSmooth_2.23-20 utf8_1.2.2         proxy_0.4-26       crayon_1.4.1     

The results:

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">

nimble | lm.x -- | -- -0.82232725857513 | -1.03677864295126 -0.876502927693985 | -0.974760506404128 -0.858686565413772 | -1.1000797449563 -0.786424043640384 | -0.929469727903681 -0.688338525708843 | -0.958121857059611 -0.646006763105514 | -0.764433142958986 -0.709161374581548 | -0.892886985978557 -0.835793984533608 | -1.14589608020523 -0.812737439486002 | -1.12577896910568 -0.86187317932852 | -1.0782901574334 -0.855261315967987 | -0.956724296180544 -0.889154185893735 | -1.02506647983045 -0.840048451024205 | -1.08381443465891 -0.80990753317543 | -1.00541709721234 -0.759248973226371 | -0.975112288165628 -0.858613957576768 | -1.14430364775781 -0.700019037284245 | -0.92014704006507 -0.859785451942923 | -1.12990489599549 -0.790357999019719 | -1.03179045854255 -0.745481917835726 | -1.02967259799018 -0.812743181991503 | -1.00116485232969 -0.796261269960158 | -0.972140732732058 -0.772647389490486 | -0.928329688788328 -0.859329771644722 | -1.01867771933016 -0.754810459171743 | -0.962827845998812 -0.831211743010605 | -1.13312305460461 -0.854018029248226 | -1.07643180240894 -0.816860284505026 | -1.05400433443646 -0.711645793700346 | -0.941255309575998 -0.823735678925053 | -1.16512124537676

[paciorek]

hi Connor, you have

alpha ~ dnorm(0, 10)

which is using 10 as the precision by default. So one possibility is that you are overly constraining the intercept and that then affects the beta estimates.

If that doesn't change things, I'm happy to look more closely, though if I did, I'd ask for your Stan code as well for comparison.

(Note that in our example online we use the improper CAR and hence don't have an intercept in the example code.)

[ConnorDonegan]

That would do it, thanks! I'll check it out and close this if that's it

[ConnorDonegan]

Thanks; I've switched to alpha ~ dnorm(0, sd = 100), but unfortunately that doesn't make a difference here.

I can provide some Stan code if its really helpful for you but I think that might just raise more questions since its still being tested. I can use Stan's multivariate normal distribution to build a CAR model, but it'll be slow. I know I opened by explaining that I'm comparing Nimble to Stan, but, to be more clear, I didn't intend to raise an issue about Nimble giving a different posterior than Stan (that was just me being chatty and giving background on why I'm doing this test); the issue is really that there appears to be a bias in the Nimble results relative to what we should expect from this Monte Carlo study. Its still possible that its due to user error though.

[paciorek]

So with your Stan code, you are getting a posterior that doesn't seem to be biased? Is that right? And is the Stan model equivalent to the nimble model?

I'm asking because it could be that the posterior gives biased estimates of the parameters rather than that nimble is giving a biased estimate of the posterior. So having some evidence that the posterior is not biased would be helpful before digging deeper. If both WinBUGS and nimble are giving similar 'bias', then it's probably that the posterior is a biased estimate of the true parameters.

[ConnorDonegan]

I added the Stan model to the code above to confirm those things for you (and for me again). I gave them equivalent priors. In the output I record, the lm and Stan results are centered on the correct value, and the Stan results have less variation than lm (lower MSE, eyeballing it). The distribution of Nimble results is pushed over towards zero by 0.2.

To code the CAR model in Stan, I'm using the multi_normal_prec_lpdf, which uses the precision matrix to calculate the log pdf of the multivariate normal. So the precision matrix is: M^(-1) %*% (I - rho * C), where M contains the conditional variances. The diag_pre_multiply Stan function creates a diagonal matrix from the vector M_inv then does matrix multiplication with the second argument:

matrix[n, n] Sigma_inv = diag_pre_multiply(1/sigma^2 * M_inv, (I - rho * C));

I don't draw a large number of samples from the Stan model because that's enough to get effective sample sizes with small monte carlo standard errors for the mean of beta. (1/2 are discarded as warmup by default.)

Other than adding Stan, I changed the code slightly so that it writes the results to disk as it goes.

library(nimble)
library(CARBayesdata, quietly = TRUE)
library(sp, quietly = TRUE)
library(spdep, quietly = TRUE)
library(rstan)
Nsim <- 30
#' draw from simultaneous autoregressive (SAR) model
sim_sar <- function (m = 1, mu = rep(0, nrow(w)), w, rho, sigma = 1, ...) {
    stopifnot(inherits(w, "matrix"))
    stopifnot(ncol(w) == nrow(w))
    N <- nrow(w)
    if (missing(mu)) {
        mu <- rep(0, N)
    }
    I <- diag(N)
    S <- sigma^2 * solve( (I - rho * t(w)) %*% (I - rho * w) )
    x <- MASS::mvrnorm(n = m, mu = mu, Sigma = S, ...)
    return(x)
}

nimble_code <- nimbleCode({
    # priors
    alpha ~ dnorm(0, sd = 10) #/#
    beta ~ dnorm(0, sd = 10)
    sigma ~ dgamma(2, 2)   
    prec <- 1 / sigma^2
    rho ~ dunif(0, 1)
    # latent process
    s[1:N] ~ dcar_proper(mu[1:N],
                         C = C[1:nC],
                         adj = adj[1:L],
                         num = num[1:N],
                         M = M[1:N],
                         tau = prec, 
                         gamma = rho)    
    # likelihood
    for(i in 1:N) {
        mu[i] <- alpha + beta * x[i]
        lambda[i] <- E[i] * exp(s[i])
        y[i] ~ dpois(lambda[i])
    }
})

stan_code <- "
data {
  // mortality data
  int n;
  int y[n];
  vector[n] x;  
  vector[n] E;
  // CAR data
  vector[n] M_inv;
  matrix[n, n] C;
}

transformed data {
  matrix[n, n] I = diag_matrix(rep_vector(1, n));
}

parameters {
  real alpha;
  real beta;  
  vector[n] phi;
  real<lower=0> sigma;
  real<lower=0, upper=1> rho;
}

model {
  matrix[n, n] Sigma_inv = diag_pre_multiply(1/sigma^2 * M_inv, (I - rho * C));
  vector[n] mu_y = E .* exp(phi);
  vector[n] mu_phi = alpha + x * beta;
  target += poisson_lpmf(y | mu_y);
  target += multi_normal_prec_lpdf(phi | mu_phi, Sigma_inv);
  target += gamma_lpdf(sigma | 2, 2);
  target += normal_lpdf(alpha | 0, 10);
  target += normal_lpdf(beta | 0, 10);  
  // implicit uniform(0,1) on rho
}

"
stan_car <- rstan::stan_model(model_code = stan_code)

# the data
data(GGHB.IG)
data(respiratorydata)
respiratorydata_spatial <- merge(x = GGHB.IG, y = respiratorydata, by = "IG", all.x = FALSE)
E <- respiratorydata_spatial$expected
N <- length(E)

# Nimble data
W.nb <- poly2nb(respiratorydata_spatial, row.names =  rownames(respiratorydata_spatial@data))
nbInfo <- nb2WB(W.nb)
CM <- as.carCM(nbInfo$adj, nbInfo$weights, nbInfo$num)
constants <- list(N = N,
                  L = length(nbInfo$adj),
                  C = CM$C,
                  M = CM$M,
                  nC = length(CM$C),
                  adj = nbInfo$adj,
                  num = nbInfo$num, 
                  E = E
                  )
nimble_data <- list()
inits <- list(alpha = rnorm(1, 0, 1), beta = rnorm(1, 0, 1), prec = 4, s = rnorm(N))

# Stan data
C <- spdep::nb2mat(spdep::poly2nb(respiratorydata_spatial), style = "W")
stan_data <- list(
    n = nrow(C),
    E = respiratorydata_spatial$expected,
    M_inv = 1 / CM$M,
    C = C
    )

# Sim data
rho.x <- 0.75
rho.y <- 0.6
beta <- -1
alpha <-  0

res <- lapply(1:Nsim, function(i) {
    # draw data
    x <- sim_sar(w = C, rho = rho.x, sigma = 0.2)
    phi <- sim_sar(w = C, mu = alpha + x * beta, rho = rho.y, sigma = 0.2)
    y <- rpois(n = N, lambda = E * exp(phi))
    stan_data$x <- constants$x <- x 
    stan_data$y <- nimble_data$y <- y
    # par(mfrow=c(1,3)); hist(x); hist(y/E); plot(x, log(y/E))
    # sample Nimble
    model <- nimbleModel(nimble_code, constants = constants, data = nimble_data, inits = inits)
    cModel <- compileNimble(model)
    conf <- configureMCMC(model, monitors = c('alpha', 'beta', 'sigma', 's'))
    MCMC <- buildMCMC(conf)
    cMCMC <- compileNimble(MCMC, project = cModel)
    nimble.samples <- runMCMC(cMCMC, niter = 15000, nburnin = 7500)
    nimble.res <- mean(nimble.samples[,"beta"])

    # sample Stan
    S <- rstan::sampling(stan_car, data = stan_data, chains = 1, iter = 1e3, pars = "beta")
    stan.res <- mean(as.matrix(S, pars = "beta"))

    # fit lm
    fit.lm <- lm(log(y/E) ~ x)
    lm.res <- coefficients(fit.lm)["x"]

    # write results to disk
    res <- matrix(c(nimble = nimble.res, stan = stan.res, lm = lm.res), nrow = 1)
    write.table(as.data.frame(res),
                file = "sim-study-res.txt",
                append = TRUE,
                row.names = FALSE,
                col.names = FALSE
                )
    return (res)
}
)

##results = do.call("rbind", res)
results <- read.table("sim-study-res.txt", header = FALSE)
names(results) <- c("nimble", "stan", "lm")

message("Sampling from log-linear Poisson model with beta = -1 (nsim = ",nrow(results), ")")
message("Sample mean of Nimble estimates: ", mean(results[,"nimble"]))
message("Sample mean of Stan estimates: ", mean(results[,"stan"]))
message("Sample mean of lm estimates: ", mean(results[,"lm"]))

par(mfrow = c(1,3))
hist(results[,1])
hist(results[,2])
hist(results[,3])

Results:

> message("Sampling from log-linear Poisson model with beta = -1 (nsim = ",nrow(results), ")")
Sampling from log-linear Poisson model with beta = -1 (nsim = 34)
> message("Sample mean of Nimble estimates: ", mean(results[,"nimble"]))
Sample mean of Nimble estimates: -0.784268517370144
> message("Sample mean of Stan estimates: ", mean(results[,"stan"]))
Sample mean of Stan estimates: -0.999127910957048
> message("Sample mean of lm estimates: ", mean(results[,"lm"]))
Sample mean of lm estimates: -1.02010261324937

-0.869454220869161 -1.11441602278424 -1.10393272798944 -0.84119625465296 -1.11071508243597 -1.12079492982895 -0.800128392734161 -0.988316639931817 -1.06258154392354 -0.776819840463203 -0.938352608682106 -0.986763487133298 -0.792577668122557 -1.02208182896783 -0.979232853530182 -0.976453894222416 -1.19158695522661 -1.12225311022952 -0.740169374376679 -0.954207779781052 -1.07013744151004 -0.873376679996589 -1.0648919424246 -1.04558067401789 -0.724048339785209 -1.0080065846679 -1.00920131855944 -0.807474986904927 -0.996162062869756 -0.960728447973934 -0.73303293118834 -0.940479203284062 -0.994227107364423 -0.854968287498536 -0.985146857240806 -0.966415127323065 -0.776805521594512 -0.971800018628652 -1.01006883960142 -0.828408968073413 -1.0570795949291 -1.18683690272438 -0.79767688270586 -1.00572965600818 -1.00901297746585 -0.687822171652225 -0.902619838933142 -0.895493864956022 -0.772841193385266 -0.968033844390289 -0.936703309471773 -0.838221630553068 -1.05991583300346 -1.15261358649034 -0.816355511186022 -1.00201745889098 -0.948745638150729 -0.789788879631268 -1.10169418202533 -1.11972871406528 -0.767592460964726 -0.976893644419625 -0.900424789914287 -0.769887361464003 -0.982048088264885 -1.07134848049032 -0.837946145938048 -1.06072124415297 -1.11658243355564 -0.808884993812324 -1.00108613803788 -0.972081210171968 -0.710107378303388 -0.873966797825023 -0.796884219728618 -0.756694078375703 -1.01565550444563 -0.936535262319068 -0.812711014652591 -1.00041841936907 -1.13257559911405 -0.679599909649094 -0.883041246176913 -0.871967739230857 -0.757009048237822 -0.971855632286032 -1.05090172762642 -0.677675951799145 -0.876468795237584 -0.938243921290235 -0.664339061520318 -0.881208093005099 -0.95184423510215 -0.646853416369382 -0.821342599368114 -0.948134618193292 -0.813473539665214 -1.05308365560774 -1.01175369500405 -0.864733600236752 -1.18930511923717 -1.30315831642805

[paciorek]

I don't understand your model code. You have:

lambda[i] <- E[i] * exp(s[i])
y[i] ~ dpois(lambda[i])

which means that y is not a function of mu (and therefore not a function of the regression coefficients). Is this just a typo in the code in this GH issue or the actual model code you are using?

[ConnorDonegan]

It's written as intended. Here it is re-written using the same notation as the Stan model and the simulated data, and re-arranged for clarity:

nimble_code <- nimbleCode({
    alpha ~ dnorm(0, sd = 10) 
    beta ~ dnorm(0, sd = 10)
    sigma ~ dgamma(2, 2)   
    prec <- 1 / sigma^2
    rho ~ dunif(0, 1)
    for(i in 1:N) {
        phi_mu[i] <- alpha + beta * x[i]
    }
    phi[1:N] ~ dcar_proper(phi_mu[1:N],
                         C = C[1:nC],
                         adj = adj[1:L],
                         num = num[1:N],
                         M = M[1:N],
                         tau = prec, 
                         gamma = rho)    
  for (i in 1:N) {
        lambda[i] <- E[i] * exp(phi[i])
        y[i] ~ dpois(lambda[i])
}
})

The model for the log-rates, phi, is multivariate normal with mean of mu = alpha + x * beta.

This could be re-parameterized so that phi is centered on zero (pseudo code; call it phi_zc):

lambda[i] = E[i] * exp(alpha + beta * x[i] + phi_zc[i])
phi_zc[1:N] ~ dcar_proper(zero_mean[1:N], ...)

[paciorek]

Ah, right, I was reading too quickly and didn't see that mu was the mean of s.

So I tried some things and if one has the mean of phi be zero with betadirectly as part of the Poisson mean (as indicated in the previous comment), then one gets results similar to Stan. Similarly if one sets up a Nimble model with the covariance matrix calculated explicitly (as in Stan) with a dmnorm and uses AF_slice, one gets result similar to Stan. And again, using the ICAR (dcar_normal), we get results like Stan. So the outlying result seems to be for dcar_proper with beta directly in the mean for s, namely a "centered" parameterization.

Also, if I fix s at the posterior mean (for phi) from Stan, the samples for beta from the original (centered) CAR (proper) are similar to Stan's. But if I fix s at the posterior mean from the original CAR, then samples for beta are similar to the original CAR with s as part of the MCMC. And yet the two posterior means for phi and s are pretty similar.

[ConnorDonegan]

Oh that's interesting. I thought that might be it, but my testing must have been screwed up by the prior I had on the intercept (presumably, that's why the zero-centered OpenBUGS model failed so terribly). With the same priors as Stan (i.e., fixing my prior on the intercept), OpenBUGS does give identical results as Stan, using the covariate-centered car_proper specification.

In Stan, the zero-centered models tend to sample really poorly. Stan will throw a bunch of divergent transition warnings at you.

[paciorek]

@ConnorDonegan just to make sure I am understanding: with the covariate-centered proper CAR, openBUGS and Stan give the same results, but NIMBLE gives different results. Right? (We don't have easy access to a Windows/WinBUGS machine amongst the core NIMBLE developers...)

[paciorek]

Ok, I think I see the problem, and it looks like a bug affecting cases where the mean of the CAR is not the same for all locations.

In CAR_proper_evaluateDensity$getMean, we calculate the mean for a location (conditional on the neighbors) without accounting for the fact that the marginal mean at each location differs (this is also reflected in the conditional density in Section 9.2.1.2 of the manual):

mean <- mu + gamma * sum(neighborCs[1,1:numNeighbors] * (neighborValues - mu))

The last mu should be a vector of values specific to the neighborValues. So we need to grab the values of model$getParam(targetDCAR, 'mu') corresponding to the neighbors.

@danielturek please see what you think.

[ConnorDonegan]

Yeah, that's correct @paciorek. Here are some results using the same monte carlo study design, but a different county/population structure. The vertical line on the histograms show the correct value of beta, the lines on the scatter plots are abline(a = 0, b = 1):

After some digging around, I was able to install OpenBUGS on ubuntu and run it with R2OpenBUGS.

[danielturek]

@ConnorDonegan @paciorek

I've looked at this carefully. The original implementation of NIMBLE's dcar_proper sampler did indeed contain a bug, which arose when the mean vector mu was non-constant.

The unfortunate fact is that the original dcar_proper sampler implementation correctly followed the equations as presented in the GeoBUGS User Manual, Version 1.2 (Sept 2004), specifically that of the conditional specification of the proper CAR distribution, which contains a typo relating to the mean term. Thanks for noticing this oddity @paciorek.

https://www.mrc-bsu.cam.ac.uk/wp-content/uploads/geobugs12manual.pdf (Appendix 1, section on "Conditional Specification")

I've fixed this on the branch fix_car_proper_mean of the nimble-dev/nimble repository.

Using this fix, I've checked results with Stan on two different models, including that of @ConnorDonegan. The posterior results of NIMBLE, after this fix, now agree with those of Stan.

I'm opening a PR with this fix.