Issues with setting priors

[anmiri]

Thank you, first of all, for developing this very useful package for spatial analysis.

I had been using it successfully for a while, but occasionally would get a strange error when setting my priors as described in the documentation. I have five covariates, all of which I want to assign a normal(0,1) prior:

prior = list() 
prior$beta <- geostan::normal(c(0, 1), c(0, 1), c(0,1), c(0,1), c(0,1))
prior$intercept <- normal(0,3)

Upon which I get the following error:

Error in normal(c(0, 1), c(0, 1), c(0, 1), c(0, 1), c(0, 1)) : 
  unused arguments (c(0, 1), c(0, 1))

I specified that R should use the norma() function from geostan in order to avoid any possible interference from other packages. However, I'm still at a loss as to how to fix this; I'm not able to set priors in other ways such as prior$beta = c(normal(0,1), normal(0,1)) as this results in a different error. I also can't set them by setting a single normal(0,1) prior for beta, as this results in other errors with the priors being underspecified.

Any help on this would be greatly appreciated!

[ConnorDonegan]

Hi, @anmiri! Glad you've found it useful.

Try this:

K <- 5
prior <- list()
prior$beta <- normal(
                     location = rep(0, times = K),
                     scale = rep(1, times = K)
                     )
prior$intercept <- normal(0, 3)
print(prior)

Does this work for you?

The function is looking for a vector of location parameters and then a vector of scale parameters. Instead of a series of paired values, provide each parameter independently. I'll add this example to the documentation.

Thanks for asking. Let me know if you have any other issues or if there are any changes or additions you think could improve your experience with geostan.

[anmiri]

Hi @ConnorDonegan, thank you so much for your quick response.

I tried setting the prior as you specified. It works, but when I run the chains, I come across another error that I've seen in this context before:

Error in dimnames(x)$parameters : 
  $ operator is invalid for atomic vectors

Do you have any idea why this might occur? I currently have only sf, spdep, tidyverse loaded in addition to geostan.

[ConnorDonegan]

Alright. When I use it with k=2 and the example model from here the documentation I don't get an error. Can you show your code, from setting priors through the model call?

[ConnorDonegan]

As in (standard normal priors for two covariates):

library(geostan)
K <- 2
prior <- list()
prior$beta <- normal(
                     location = rep(0, times = K),
                     scale = rep(1, times = K)
                     )
prior$intercept <- normal(0, 3)
print(prior)

fit <- stan_glm(log(rate.male) ~ I(income/1e3) + college,
                data = georgia,
                prior = prior,
                prior_only = TRUE,
                chains = 2,
                iter = 600) 
print(fit)

[ConnorDonegan]

Hi Miri, any update on this? Is it working for you, or is this error still coming up?

I haven't seen that error before; based on the error message, I'm curious if its an intercept-only model or if you're using a formula that tries to drop the intercept (e.g., ~ 0 + x).

[anmiri]

Hi @ConnorDonegan, Apologies for my delayed response on this.

It's a model with an intercept (~ 1 + x1). Since I can't share my actual data, I came up with a MWE based on the georgia dataset. However, the error is gone now, so I unfortunately can't illustrate what caused it - I will let you know if it re-emerges and try to find out why it was happening in the first place, but I think it was to do with my specification of the priors in the list. It would be helpful to put the way in which you specify the same or different priors for multiple covariates in the documentation, I think.

Feel free to close the issue and thanks again for your help.

data(georgia)

K <- 3
prior <- list()
prior$beta <- normal(
  location = rep(0, times = K),
  scale = rep(1, times = K)
)
prior$intercept <- normal(0, 3)

data = drop_na(georgia)
W <- shape2mat(data, "B")
cp = prep_car_data(W)

fit <- stan_car(log(deaths.female) ~ 1 + college + log(income) + log(pop.at.risk.female),
         data = data,
         car_parts = cp,
         chains = 2,
         iter = 600
         )

[ConnorDonegan]

Alright, glad it went away at least. Let me know if it recurs or any other issues come up. I'll add this case to the documentation of the priors so its more clear. Thanks for your help!

I realize your example code isn't exactly what you're working with and it didn't need to be written exactly correct and so on; particularly because I see this issue often, I'll just note (to potentially avoid any issues in the future) to remember to use the Poisson family = poisson() and the offset wrapper for population at risk/denominator data + offset(log(pop.at.risk)) when the outcome is a count.