Random effect logistic model prior specification problem

[Raisa Rashid]

sim <- 1000
n_clusters <- 10
n_obs <- 10
n <- n_clusters * n_obs
B0 <- -1
B1 <- 1.5
B2 <- 0.8
sigma2b <- 2.5
estimates <- matrix(NA, nrow = sim, ncol = 4)

#Half-t prior
HT.prior = "expression:
sigma = exp(-theta/2);
nu = 3;                        
sigma = exp(-theta/2);
nu = 3;
log_dens = 0-0.5 * log(nu * pi) - (-0.1207822);
log_dens = log_dens - 0.5 * (nu + 1) * log(1 + sigma * sigma);
log_dens = log_dens - log(2) - theta / 2;
return(log_dens);"

formula_HT <- Y ~ 1 + X1 + X2 +
  f(cluster, model = "iid", hyper = list(prec = list(prior = HT.prior)))

#Half-Normal prior
HN.prior = "expression:
tau0 = 0.001;                     # precision of Normal prior (small → diffuse)
sigma = exp(-theta/2);
log_dens = log(2) - 0.5 * log(2 * pi) + 0.5 * log(tau0);
log_dens = log_dens - 0.5 * tau0 * sigma^2;
log_dens = log_dens - log(2) - theta / 2;
return(log_dens);"
formula_HN <- Y ~ 1 + X1 + X2 +
  f(cluster, model = "iid", hyper = list(prec = list(prior = HN.prior)))

#Half-Cauchy prior
HC.prior = "expression:
sigma = exp(-theta/2);
gamma = 25;                       # scale parameter
log_dens = log(2) - log(pi) - log(gamma);
log_dens = log_dens - log(1 + (sigma / gamma)^2);
log_dens = log_dens - log(2) - theta / 2;
return(log_dens);"
formula_HC <- Y ~ 1 + X1 + X2 +
  f(cluster, model = "iid", hyper = list(prec = list(prior = HC.prior)))

set.seed(123)
for (i in 1:sim) {
  cluster_id <- rep(1:n_clusters, each = n_obs)
  bi <- rnorm(n_clusters, mean = 0, sd = sqrt(sigma2b))
  X1 <- rbinom(n_clusters * n_obs, 1, 0.5)
  X2 <- rnorm(n_clusters * n_obs, 0, 1)
  bi_expanded <- rep(bi, each = n_obs)
 
  eta <- B0 + B1 * X1 + B2 * X2 + bi_expanded
  p <- exp(eta) / (1 + exp(eta))
  Y <- rbinom(n, 1, p)
 
  data <- data.frame(cluster = cluster_id, Y = Y, X1 = X1, X2 = X2)
 
  model_inla <- inla(
    formula_HN,
    data = data,
    family = "binomial",
    control.predictor = list(compute = TRUE),verbose = FALSE,
    control.inla = list(strategy = "simplified.laplace"),
    control.compute = list(dic = TRUE, return.marginals.predictor = TRUE,config=FALSE)
  )
 
  sig_hatb <- 1 / model_inla$summary.hyperpar["Precision for cluster", "mean"]
  estimates[i, 1:3] <- model_inla$summary.fixed[, "mean"]
  estimates[i, 4] <- sig_hatb
}
summary(model_inla)

mean_estimates <- colMeans(estimates, na.rm = TRUE)
mean_estimates

I want to estimate random effect ligistic model parameters. However, for half-t prior, the code works and gives results but for both half normal and half cauchy, the code does no work. it gives the following error message:

Eval: Error while parsing expression:
----------
Message   :  "Unexpected end of expression at position 14"
Token     :    ""
Position  : 14
Error Code:     2

 *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1



Eval: Error while parsing expression:
----------
Message   :  "Unexpected end of expression at position 14"
Token     :    ""
Position  : 14
Error Code:     2
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  :
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.
The inla program failed and the maximum number of tries has been reached.

can u please help? is the prior specification okay or not ?

[Janet Van Niekerk]

Hi,

The comments in the expression terms are causing the issue - remove them and it works.

Remove  # precision of Normal prior (small → diffuse) and  # scale parameter

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/7fb1e86e-2d3a-4e12-b09e-b2e0496c5ed3n%40googlegroups.com.

[Helpdesk (Haavard Rue)]

I would recommend using 'rprior' instead of 'expression:', its just easier

inla.doc("rprior")

> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/7fb1e86e-2d3a-4e12-b09e-b2e0496c5ed3n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org