Custom Prior for a Treatment Contrast

[Arlene Jiang]

Hi,

I am trying to use a custom shaped prior for the treatment effect (risk difference) using an rgeneric latent effect. I'm able to create a prior with the following code, but am unsure of the syntax to fix this to a treatment indicator (0/1) coefficient. Would you able to help with this? Thanks so much!

rgeneric.pearsonIV <- function(
    cmd = c("graph", "Q", "mu", "initial", "log.norm.const", "log.prior", "quit"),
    theta = NULL,
    m = 2, nu = -0.4,
    location = 0.01, scale = 0.08
) {
  interpret.theta <- function() list(value = as.numeric(theta))
 
  graph <- function() Matrix::Diagonal(1)
  Q     <- function() Matrix::Diagonal(1, x = 1e6)
  mu    <- function() numeric(0)
  log.norm.const <- function() numeric(0)
  initial <- function() rep(0, 1)
  quit <- function() invisible()
 
  log.prior <- function() {
    if (is.null(theta) || is.na(theta[1])) return(-Inf)
    logd <- PearsonDS::dpearsonIV(
      x = as.numeric(theta),
      m = m, nu = nu, location = location, scale = scale,
      log = TRUE
    )
    as.numeric(logd)
  }
 
  cmd <- match.arg(cmd)
  val <- switch(cmd,
                graph = graph(),
                Q = Q(),
                mu = mu(),
                initial = initial(),
                log.norm.const = log.norm.const(),
                log.prior = log.prior(),
                quit = quit())
  return(val)
}
model_obj <- inla.rgeneric.define(
  model=rgeneric.pearsonIV
)

# Run model with empty data
df_prior <- data.frame(y = NA, idx = 1)
result_prior <- inla(
  y ~ f(idx, model = model_obj),
  data = df_prior,
  family = "gaussian",
  control.family = list(hyper = list(prec = list(initial = log(1), fixed = T))), # Keeps precision constant
  control.compute = list(config=T)
)
summary(result_prior)

 

[Elias T. Krainski]

Hi, 

The "rgeneric" and the "cgeneric" are intended to implement new latent models, as the "cloglike” is for new likelihoods. As far as I can see, you want to do a new prior for a hyperparameter. Is that correct? If so, you can use either 'expression:', 'table:', or 'rprior:' as per the examples in inla.doc("rprior")

Best regards, 

Elias

--
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/d0164f57-a178-4f12-a033-3f88fc560fcan%40googlegroups.com.

[Arlene Jiang]

Hi,

Thanks for referring me to these resources. I'm hoping to set this custom prior onto a fixed effect, however, I have a non-normal prior distribution. From what I understand, using rprior would be setting a prior on the precision (or log-precision) rather than the fixed effect itself. Do you think coding this as a latent effect would be possible instead?

Thanks,

Arlene

[INLA help]

If you want a non Gaussian fixed effects, check model ‘clinear’ 

 

This will put the fixed effects into the hyperparameters , and u can use the tools Elias mentioned 

—

Haavard Rue

Helpdesk

he...@r-inla.org

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/f99cf974-4336-4ee6-82ec-374fca759987n%40googlegroups.com.

[Arlene Jiang]

Hi,

Thanks for pointing me to 'clinear'. I can only seem to apply my custom prior to the log precision, rather than the trt prior itself. What am I missing?

Thanks,

Arlene

> y <- NA;trt <- 1 > df <- data.frame(y=y, trt=trt) > rd_prior <- "expression: + m = 2; + nu = -0.4; + location = 5; + scale = 0.08; + log_dens = 2 * (-0.01082386-0.06468127*i); + log_dens = log_dens -m * log(1 + ((theta - location)/scale)^2); + log_dens = log_dens - nu * atan((theta-location)/scale); + return(log_dens); + " > model2 <- inla(y~f(trt, model="clinear"), + control.family=list(hyper=list(prec=list(prior=rd_prior))), + data=df, family="gaussian") > summary(model2) Time used: Pre = 0.305, Running = 0.435, Post = 0.114, Total = 0.854 Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld (Intercept) 0 3162.254 -6200.867 0 6200.867 0 0 Random effects: Name Model trt Constrained linear Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode Precision for the Gaussian observations 150.136 7.042 136.986 149.892 164.71 149.229 Beta for trt 0.988 0.294 0.409 0.988 1.57 0.987 Marginal log-Likelihood: 6.05 is computed Posterior summaries for the linear predictor and the fitted values are computed (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)') > inla.priors.used(model2) section=[family] tag=[INLA.Data1] component=[gaussian] theta1: parameter=[log precision] prior=[expression: m = 2; nu = -0.4; location = 5; scale = 0.08; log_dens = 2 * (-0.01082386-0.06468127*i); log_dens = log_dens -m * log(1 + ((theta - location)/scale)^2); log_dens = log_dens - nu * atan((theta-location)/scale); return(log_dens); ] param=[1e+00, 5e-05] section=[fixed] tag=[(Intercept)] component=[(Intercept)] beta: parameter=[(Intercept)] prior=[normal] param=[0, 0] section=[random] tag=[trt] component=[trt] theta1: parameter=[beta] prior=[normal] param=[ 1, 10] group.theta1: parameter=[logit correlation] prior=[normal] param=[0.0, 0.2]

[Helpdesk (Haavard Rue)]

the 'clinear' is implemented so that 'beta' is in the hyperpar, meaning you can
define your own general prior for it. you set low = -Inf and high = Inf

the 'precision=' argument; simply, do not touch

> > > https://groups.google.com/d/msgid/r-inla-discussion-group/f99cf974-4336-4ee6-82ec-374fca759987n%40googlegroups.com
> > > .

> --
> 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/42d409c9-2902-4ecd-83f9-0e8a2305073bn%40googlegroups.com
> .

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