inlabru-org / inlabru

inlabru
https://inlabru-org.github.io/inlabru/
73 stars 21 forks source link

Adding 'scale' option in the 'like()' function #230

Closed osupplisson closed 3 months ago

osupplisson commented 3 months ago

Hi,

I was trying to fit a Gaussian approximation to the betabinomial model using the 'betabinomialna' likelihood from INLA (https://inla.r-inla-download.org/r-inla.org/doc/likelihood/betabinomialna.pdf).

It appeared that the scale argument from the inla() call is currently not a valid option for the like() function. I believe it should be the case as this is a parameter that can change across likelihood (like E and N) ?

Best, Olivier

(ps: I think, the same thing holds for the strata argument)

Please, find a small reprex hereafter:

library(inlabru)

#Fake data
n = 1000
rho = 0.1
z = rnorm(n, sd=0.4)
Ntrials = sample(1:20, n, replace=TRUE)
eta = 1 + z
p = exp(eta)/(1+exp(eta))
s = runif(n)
m = Ntrials * p
v = Ntrials * p * (1.0 - p) * (1.0 + s * (Ntrials - 1) * rho)
y = rnorm(n, mean = m, sd = sqrt(v))

#INLA
formula = y ~ 1 + z
data = data.frame(y, z, s)
r = inla(formula, data = data, scale = s,
         family = "betabinomialna", Ntrials=Ntrials, verbose = TRUE,
         control.inla = list(strategy = "adaptive"))
summary(r)
r$.arg$scale

#BRU
cmp <- ~intercept(1) + linear(z,model = "linear",mean.linear = 0,prec.linear = 1)
lik1 <- like(
  data = data,
  Ntrials = Ntrials,
  ###Should be a likelihood arguement 
  scale = s,
  family = "betabinomialna",
  formula = y ~ intercept + z
)

#Error in like(data = data, Ntrials = Ntrials, scale = s, family = "betabinomialna",  : unused argument (scale = s)

#Trying without the argument
lik1 <- like(
  data = data,
  Ntrials = Ntrials,
  family = "betabinomialna",
  formula = y ~ intercept + z
)

rr <- bru(components = cmp,
    lik1)
rr$.args$scale
#NULL
finnlindgren commented 3 months ago

Hi, this should now be fixed in development version 2.10.1.9001. Also note a small correction to your component definition:

cmp <- ~intercept(1) + z(z,model = "linear",mean.linear = 0,prec.linear = 1)
lik1 <- like(
  data = data,
  Ntrials = Ntrials,
  ###Should be a likelihood arguement 
  scale = s,
  family = "betabinomialna",
  formula = y ~ intercept + z
)