fitting Bernardinelli model with covariates

[Aminath Shausan]

Dear community, 

I am trying to fit a Bernardinelli model including several covariates. 

My code is below: 

formula  <- y ~ log(lag1)  + log(lag2) + season + 
  tas + dtr + RH +        
  tas_lag1 + dtr_lag1 + RH_lag1 +
  f(ID.area, model="iid", constr=TRUE, hyper=list(prec=list(prior=sdunif))) + 
  f(ID.area1, ID.month, model="iid", constr=TRUE, hyper=list(prec=list(prior=sdunif))) + 
  ID.month 

mdl <-  inla(formula,
             family = "nbinomial", data = data.req, E=popn,
             control.predictor=list(compute=TRUE, cdf=c(log(1)), link =1),
             control.compute=list(dic=TRUE, cpo=TRUE, waic=TRUE),  
             control.fixed = list(prec.intercept = 0.0001),
             control.inla=list(strategy=strategy)
)

where sdunif is Uniform(0, Inf) for prior for precision of tandom effects: 

sdunif="expression:  
  logdens=-log_precision/2;
  return(logdens)"

My questions: 

1. Do we need to use constr=TRUE for random effects? I have seen in some studies this is not used (https://inla.r-inla-download.org/r-inla.org/case-studies/Blangiardo-et-al-2012/Report-version-Oct-2012.pdf; https://www.paulamoraga.com/book-geospatial/sec-arealdataexamplest.html) . When do we use constraints and when do we omit it? 

2. In my data, the season variable is already factored. Do I rather factor it in the formula? 

3. At the moment, the precision for fixed effects parameters (except intercept) are given default fixed value  (0.001). Is it possible to assign the function sdunif to the  prior of precision for fixed effects? 

4. I am getting very large value for ID.area1 hyper parameter (mean = 1.14e+06, df= 1.98e+06). Can I check if this means variance of this parameter = (1/1.14e+06) =0.000000877 ?

Thank you

Aminath

[Helpdesk (Haavard Rue)]

Theory says its not needed. I think in practice it will help separating out
effects, as otherwise there is a strong confounding between the intercept and
the precision for the iid model (or parameters in the AR1, if that is the one),

try this

n=20
y=rnorm(n)

and compare

> r=inla(y~1 + f(idx,model="iid",constr=F), data=data.frame(y,idx=1:n),
family="normal", control.family=list(hyper=list(prec=list(initial=20,
fixed=TRUE))));

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.227 0.117 -0.456 -0.227 0.002 -0.227 0


> r=inla(y~1 + f(idx,model="iid",constr=T), data=data.frame(y,idx=1:n),
family="normal", control.family=list(hyper=list(prec=list(initial=20,
fixed=TRUE))));

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) -0.227 0 -0.227 -0.227 -0.227 -0.227 0

> 2. In my data, the season variable is already factored. Do I rather factor it
> in the formula? 

no, if 'x' is a factor, then it will expand automatically in the formual

> 3. At the moment, the precision for fixed effects parameters (except
> intercept) are given default fixed value  (0.001). Is it possible to assign
> the function sdunif to the  prior of precision for fixed effects? 

fixed effects as fixed precision, if you make it random, its more in the 'random
effect' category. I do not think you want that.

> 4. I am getting very large value for ID.area1 hyper parameter (mean =
> 1.14e+06, df= 1.98e+06). Can I check if this means variance of this parameter
> = (1/1.14e+06) =0.000000877 ?

yes, this is about right. you can always use 'inla.tmarginal' to transform
marginals from prec to variance,

like

> r=inla(y~1,data=data.frame(y))
> r$summary.hyperpar
mean sd 0.025quant
0.5quant 0.975quant mode
Precision for the Gaussian observations 0.7776451833 0.239891535 0.3817593116
0.753225591 1.314879932 0.7038111443
> m=inla.tmarginal(function(x) 1/x, r$marginals.hyperpar[[1]])
> inla.zmarginal(m)
Mean 1.41846
Stdev 0.476398
Quantile 0.025 0.761233
Quantile 0.25 1.08201
Quantile 0.5 1.32642
Quantile 0.75 1.64999
Quantile 0.975 2.60993


note, that E(1/x) != 1/E(x) , but the quantiles will transform in that way


>
> Thank you
> Aminath
> --
> 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/4f5ac774-6b1d-4cfb-8a85-8e05cd39ce6bn%40googlegroups.com
> .

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