Got an error when calculating the standard deviation for the RW1 component

[Jingjing]

Hi,

I ran an RW1 model in INLA. The model looks like 

form<-y~-1 + Intercept + income + f(year, model="rw1", hyper=list(theta = list(prior="pc.prec", param=c(1,0.01))), scale.model=TRUE)

I tried to get the standard deviation of the RW1 component using code: 

sd.rw1<- inla.rmarginal( 100000, inla.tmarginal(function(x) 1/sqrt(x), model$marginals.hyper$"Precision for Year")

but it gave an error saying "Warning in sqrt(x) : NaNs produced". After checking the output of "model$marginals.hyper$"Precision for Year", I noticed there were some negative values. This seems really strange. Would anyone be able to help me understand if there's an issue with my code?

Jingjing

[Helpdesk (Haavard Rue)]

the only failsafe option is to use the marginals in the internal representation,

> r=inla(y~1,data=data.frame(y=0))
> head(str(r$internal.marginals.hyperpar[[1]]))
num [1:43, 1:2] 4.75 4.95 5.13 5.81 6.44 ...
- attr(*, "hyperid")= chr "65001|INLA.Data1"
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:2] "x" "y"
NULL

the `internal.' prefix gives the marginals in the internal scale, which is
log(precision) in your example. for other models it is given in the
documentation. 

sampling is done this way

> m = r$internal.marginals.hyperpar[[1]]
> log.prec = inla.rmarginal(1000, m)
> stdev = 1/sqrt(exp(log.prec))
> summary(stdev)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.002589641 0.006071897 0.008338903 0.011762855 0.013853448 0.089429718

just transform the sample, no need to transform the density as this just gives
further approximation errors.


the deeper reason why you might get negative values in your example, is that the
internal representation for the density does not have constraints, like knowing
that x > 0. the internal representation does not have these contraints, so
you're fine. if your density has mass close to zero this might happen

best
H

> --
> 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/23b85d7a-e22e-448e-99f7-c2e4ad351391n%40googlegroups.com
> .

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

[Jingjing]

Thank you very much for your reply! I tried running the code and successfully obtained the results. 

I have a few follow-up questions:

1. Is internal.marginals.hyperpar giving the log(precision) for all types of models (e.g., Poisson, Binomial)? If my model includes multiple random effects, like RW2 and BYM2, can I obtained the log(precision) of a specific term simply by changing the index in internal.marginals.hyperpar[[]]?

2.  If I would like to get more information about the stdev, is it appropriate to use the following code for this purpose?

>summary_stdev<-c(median=median(stdev), sd=sd(stdev), CI_low=quantile(stdev, 0.025), CI_high=quantile(stdev, 0.975))

Best,

Jingjing

[Helpdesk (Haavard Rue)]

On Tue, 2026-02-03 at 07:02 -0800, Jingjing wrote:
> Thank you very much for your reply! I tried running the code and successfully
> obtained the results. 
>
> I have a few follow-up questions:
> 1. Is internal.marginals.hyperpar giving the log(precision) for all types of
> models (e.g., Poisson, Binomial)? If my model includes multiple random
> effects, like RW2 and BYM2, can I obtained the log(precision) of a specific
> term simply by changing the index in internal.marginals.hyperpar[[]]?

for models with a precision, you will find the marginal at 'some index'. check
the summary, as the order is the same, or the name of each one

>
> 2.  If I would like to get more information about the stdev, is it appropriate
> to use the following code for this purpose?
> > summary_stdev<-c(median=median(stdev), sd=sd(stdev), CI_low=quantile(stdev,
> > 0.025), CI_high=quantile(stdev, 0.975))

yes, something like this, if stdev is your simulated stdev's