Precision parameters are very high. Is there an explanation?

[alyssa...@duke.edu]

I'm running a logistic model with iid random effects and have mean of the precision parameters on the random effects estimated at 18,000+. Instead of posting my own data, you can see this using the Seeds example (code as listed in Volume I). Is there some characteristic of the data/model that causes this? Is there some sort of transformation occurring? The shape of the posterior densities for the variance of the random effect are similar, am I missing something?:

WinBUGS output:

<---------------- Mean=0.267, sd=0.1471

then INLA output (code as listed in Volume I:

Fixed effects:

               mean     sd 0.025quant 0.5quant 0.975quant   kld

(Intercept) -0.5581 0.1261    -0.8076  -0.5573    -0.3130 0e+00

x1           0.1461 0.2233    -0.2933   0.1467     0.5823 0e+00

x2           1.3206 0.1776     0.9748   1.3197     1.6716 1e-04

x1:x2       -0.7793 0.3066    -1.3799  -0.7796    -0.1774 0e+00

Random effects:

Name      Model

 plate   IID model 

Model hyperparameters:

                    mean     sd       0.025quant 0.5quant 0.975quant

Precision for plate 18413.08 18280.69  1217.90   13003.80 66486.61  <----------------------------------------- Mean=18,413.08. sd=18280.69

[Håvard Rue]

On Fri, 2013-05-24 at 12:00 -0700, alyssa...@duke.edu wrote:
> I'm running a logistic model with iid random effects and have mean of
> the precision parameters on the random effects estimated at 18,000+.
> Instead of posting my own data, you can see this using the Seeds
> example (code as listed in Volume I). Is there some characteristic of
> the data/model that causes this? Is there some sort of transformation
> occurring? The shape of the posterior densities for the variance of
> the random effect are similar, am I missing something?:
> WinBUGS output:

> [seeds3]<---------------- Mean=0.267, sd=0.1471
>
>
>
> then INLA output (code as listed in Volume I:
> Fixed effects:
> mean sd 0.025quant 0.5quant 0.975quant kld
> (Intercept) -0.5581 0.1261 -0.8076 -0.5573 -0.3130 0e+00
> x1 0.1461 0.2233 -0.2933 0.1467 0.5823 0e+00
> x2 1.3206 0.1776 0.9748 1.3197 1.6716 1e-04
> x1:x2 -0.7793 0.3066 -1.3799 -0.7796 -0.1774 0e+00
>
>
> Random effects:
> Name Model
> plate IID model
>
>
> Model hyperparameters:
> mean sd 0.025quant 0.5quant 0.975quant
> Precision for plate 18413.08 18280.69 1217.90 13003.80 66486.61
> <----------------------------------------- Mean=18,413.08.
> sd=18280.69

is this a similar case to

http://www.r-inla.org/faq#TOC-INLA-seems-to-work-great-for-near-all-cases-but-are-there-cases-where-INLA-is-known-to-have-problems-

if you have an unique iid term for each observation, then there is an
experimental option to integrate out the iid term from the model....

Best
H


--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Mobile: +47-9260-0021 Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[Alyssa Platt]

There is not a unique observation for each iid term though there are some iid groups that contain only one member and the groups are small (20% of groups have 1 member, median number of members is 3 and highest is 15; this is in a sample of N=747). Perhaps groups are small enough that there is not enough smoothing strength? And I DID see that FAQ but I could not find the "discussion" where the problem was illustrated and/or a solution was proposed. Perhaps I was not looking in the right place?

[Elias Krainski]

The key, pointed on the link (#20) at http://www.r-inla.org/faq, is that
we have a strong prior sensitivity. So it need carefull on the prior
choise. For example, considering the reproducible example on that link,
see the results bellow using two different priors

res1 <- inla(y ~ 1 + f(idx, model="iid",
hyper=list(theta=list(param=c(1,1)))),
data =data.frame(y,idx), family = "binomial", Ntrials = 1)
res1$summary.hy

res2 <- inla(y ~ 1 + f(idx, model="iid",
hyper=list(theta=list(param=c(0.1,0.1)))),
data =data.frame(y,idx), family = "binomial", Ntrials = 1)
res2$summary.hy

we have different results...

[Mark James Adams]

I am running into a similar situation using individual data grouped by country, where I am interested in a spatial component as well as a country's deviation from its neighborhood average.

If I just run an IID model, the precision looks fine and translates into a variance that matches variance in raw country means:

Random effects:

Name      Model

 ISO   IID model 

Model hyperparameters:

                                        mean    sd      

Precision for the Gaussian observations  2.6035  0.0056  

Precision for ISO                       76.2270  8.4436   

However, when combined with the besag model the precision for the IID component is higher than I would expect

Random effects:

Name      Model

 ISO   IID model 

geo   Besags ICAR model 

Model hyperparameters:

                                        mean      sd        

Precision for the Gaussian observations    2.6034    0.0056    

Precision for ISO                        704.7466  384.7496   

Precision for geo                         53.6867   10.7937      

when I say 'higher than expected' what I mean is that the spatial averaging from besag looks like it is behaving well, yielding a spatially smoothed version of the raw data. However, if I subtract the spatial BLUPs (e.g., model$summary.random$geo$mean) from each country's raw average (Anomaly = Observed - Besag), they are is considerably more variable than the IID estimates

In the figure red = low, blue = high. All fill values are plotted on the same scale. If you look at, for example, China: its raw value is high (blue) but its spatial value is in the middle (white) as it is in between India (blue) and Mongolia and Russia (red). It's IID BLUP is just barely above average whereas the anomaly is more markedly so.

Group size shouldn't be a problem: median data points/country is 400 (range 6 - 10,000). The IID and anomaly values do correlate strongly  (r = .78) but the IID BLUPs have about 7 times less variance (too much shrinkage?). 

Does my model look like it is behaving correctly or is there anything that can be done with the priors for the IID component to give more reasonable looking estimates?

[Håvard Rue]

On Tue, 2014-04-29 at 05:57 -0700, Mark James Adams wrote:
> I am running into a similar situation using individual data grouped by
> country, where I am interested in a spatial component as well as a
> country's deviation from its neighborhood average.

Hi Mark,

I suggest to try to the new parametersisation of this model suggested
in


> http://arxiv.org/abs/1403.4630

and you find the code for the bym-case here:

http://www.r-inla.org/examples/case-studies/pc-priors-martins-et-al-2014


this allows controling the prior for the marginal variance.

Best,

[Mark James Adams]

On Tuesday, April 29, 2014 8:43:42 PM UTC+1, Havard Rue wrote:

On Tue, 2014-04-29 at 05:57 -0700, Mark James Adams wrote:
> I am running into a similar situation using individual data grouped by
> country, where I am interested in a spatial component as well as a
> country's deviation from its neighborhood average.

Hi Mark,

I suggest to try to the new parametersisation of this model 

http://www.r-inla.org/examples/case-studies/pc-priors-martins-et-al-2014

Hi Håvard

Thanks for the link. I will try this alternative. Can you clarify that when I have multiple measurements of each country I should also fit an IID effect (that is, the model='bym2' parameterization is is like 'besag' than 'bym')?

Best

Mark

[Håvard Rue]

On Thu, 2014-05-01 at 08:12 -0700, Mark James Adams wrote:

>
> Thanks for the link. I will try this alternative. Can you clarify that
> when I have multiple measurements of each country I should also fit an
> IID effect (that is, the model='bym2' parameterization is is like
> 'besag' than 'bym')?

the iid effect is part of this model. if you duplicated observations,
just duplicate the linear predictor and observe both.

Best

[Mark James Adams]

Hi Håvard

Thanks for the explanation. It makes sense now that the hyper parameter has two elements and thus represents two random effects.

I have run my model with the same set up as in the example code, i.e., defining Q using the graph and setting alpha, u, &c. The model runs for a while and then says

*** ERROR ***   table-prior returns NAN. Argument is 20.8466 but prior is defined on [-16,16] only

Here is the output from verbose=TRUE: https://gist.github.com/mja/ece0addcc30bf493c7d2

Thanks.

Best

Mark

[Håvard Rue]

Hi,

I made some changes. can you upgrade and retry?

formula2 <- agr ~ age + sex + age:sex +
f(geo, model='bym2', graph=border,
constr=TRUE, scale.model=TRUE,
hyper=list(phi=list(prior='pc',
param=c(phi.u, phi.alpha)),
prec=list(prior='pc.prec',
param=c(u, alpha))))

# Agreeableness
agr_geo <- inla(formula2, family='gaussian', data=CC, verbose=TRUE,
control.family = list(
hyper = list(
prec = list(
initial = -1))))



H

[Mark James Adams]

On Tuesday, May 6, 2014 8:22:11 PM UTC+1, Havard Rue wrote:

Hi,

I made some changes. can you upgrade and retry?

formula2 <- agr ~ age + sex + age:sex +
    f(geo, model='bym2', graph=border,
      constr=TRUE, scale.model=TRUE,
      hyper=list(phi=list(prior='pc',
                         param=c(phi.u, phi.alpha)),
              prec=list(prior='pc.prec',
                      param=c(u, alpha))))

model <- inla(formula2, family='gaussian', data=CC,  verbose=TRUE,

                control.family = list(
                        hyper = list(
                                prec = list(
                                        initial = -1))))

This runs for me now, thanks (INLA version 0.0-1399439934).

The "Random effects" and "Model hyperparamter" tables in the model summary list the precision for a single component. Is this the precision of besag + iid? model$summary.random$geo has twice as many rows as the number of levels in geo. Is it correct that the first N of these are the iid estimates and the second N are the besag ones (or the other way around)?

Best

Mark

[Håvard Rue]

I think inla.doc("bym2") should answer this; let me know if it
doesn't.

H

[Mark James Adams]

That makes perfect sense. Thanks for the link to the documentation.

[David Crow]

Silly question:  I'm trying to get a handle on how to interpret the precision estimates for latent effects.  They are in the scale of the original variable, right (that is, exponentiated from the internal, log-transformed representation of the parameters)?  

So, in Mark's example, the precision for the "geo" random effect, 53.6867, translates to a variance of 1/53.6867 = .0186.  Is that right?  

Thanks, and sorry to bother everyone with such a trivial question,

David

[Helpdesk]

they are in the same scale as the linear predictor. the default for the
poisson is

y ~ Poisson(exp(lin.predictor))

yes, so if an estimate for the precision is 50 then an estimate for the
variance is 1/50 = 0.02. you can compute the posterior the variance if
you want to get it 'right'.

be aware, that for intrinsic models, like default besag, rw1, rw2, this
does not correspond to the 'marginal precision' but a 'precision
parameter', as there is a buildt in scale in these models depending on
the dimension, locations and the graph. well, adding scale.model=TRUE
in the f(), will remove this effect and make the precision also the
marginal precision (in the sense that is defined with intrinsic
models). We recommend adding scale.model=TRUE by default.

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 on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/eb251466-abc0-44cb-9bdf-6d25a799e832%40googlegroups.com
> .

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

[Longwen Zhao]

I got a similar issue with precision. We are estimating covid-19 daily cases in Missouri in the US and we believe there is spatial heterogeneity, but the posterior of the precision of spatial random (suing besag)  telling us the mean of precision is around 20,000. Which means the variance is very small. 

After we added a population density as a predictor, the mean decreased from 20,000 to 19,000. But still way too high. 

Any idea or suggestion on what's happening here?

Thanks a lot

[Helpdesk]

using model='besag', then add also scale.model=TRUE, otherwise,
precision cannot be properly compared, as its model property depends on
the graph.

see
vignette("scale-model")

On Fri, 2021-03-19 at 13:25 -0700, Longwen Zhao wrote:
> I got a similar issue with precision. We are estimating covid-19 daily
> cases in Missouri in the US and we believe there is spatial
> heterogeneity, but the posterior of the precision of spatial random
> (suing besag)  telling us the mean of precision is around 20,000.
> Which
> means the variance is very small. 
> After we added a population density as a predictor, the mean decreased
> from 20,000 to 19,000. But still way too high. 
>
> Any idea or suggestion on what's happening here?
>
> Thanks a lot
>
> On Friday, May 24, 2013 at 2:00:20 PM UTC-5 alyssa...@duke.edu wrote:
> > I'm running a logistic model with iid random effects and have mean
> > of
> > the precision parameters on the random effects estimated at 18,000+.
> > Instead of posting my own data, you can see this using the Seeds
> > example (code as listed in Volume I). Is there some characteristic
> > of
> > the data/model that causes this? Is there some sort of
> > transformation
> > occurring? The shape of the posterior densities for the variance of
> > the random effect are similar, am I missing something?:
> > WinBUGS output:

> > [seeds3]<---------------- Mean=0.267, sd=0.1471
> >
> > then INLA output (code as listed in Volume I:
> > Fixed effects:
> >                mean     sd 0.025quant 0.5quant 0.975quant   kld
> > (Intercept) -0.5581 0.1261    -0.8076  -0.5573    -0.3130 0e+00
> > x1           0.1461 0.2233    -0.2933   0.1467     0.5823 0e+00
> > x2           1.3206 0.1776     0.9748   1.3197     1.6716 1e-04
> > x1:x2       -0.7793 0.3066    -1.3799  -0.7796    -0.1774 0e+00
> >
> > Random effects:
> > Name      Model
> >  plate   IID model 
> >
> > Model hyperparameters:
> >                     mean     sd       0.025quant 0.5quant 0.975quant
> > Precision for plate 18413.08 18280.69  1217.90   13003.80 66486.61
> >  <----------------------------------------- Mean=18,413.08.
> > sd=18280.69
> >
> >

> --
> 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 on the web, visit

> https://groups.google.com/d/msgid/r-inla-discussion-group/04e6b182-3e26-4317-81dc-d87192334d63n%40googlegroups.com
> .

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

[Helpdesk]

Hi

Seeds has only 21 observations, hence the prior for the precision used
would be cruial

also, with that new observations, the E() and Var() of precision ''does
not make much sense'' as its all in the tail... You would rather want to
do E() and Var() of log(precision) that is much much much more Gaussian
like, which you'll find under

results$internal. ...



let me know if this not answer your problem

H

On Fri, 2013-05-24 at 12:00 -0700, alyssa...@duke.edu wrote:

> I'm running a logistic model with iid random effects and have mean of
> the precision parameters on the random effects estimated at 18,000+.
> Instead of posting my own data, you can see this using the Seeds
> example
> (code as listed in Volume I). Is there some characteristic of the
> data/model that causes this? Is there some sort of transformation
> occurring? The shape of the posterior densities for the variance of
> the
> random effect are similar, am I missing something?:
> WinBUGS output:

> [seeds3]<---------------- Mean=0.267, sd=0.1471
>
> then INLA output (code as listed in Volume I:
> Fixed effects:
>                mean     sd 0.025quant 0.5quant 0.975quant   kld
> (Intercept) -0.5581 0.1261    -0.8076  -0.5573    -0.3130 0e+00
> x1           0.1461 0.2233    -0.2933   0.1467     0.5823 0e+00
> x2           1.3206 0.1776     0.9748   1.3197     1.6716 1e-04
> x1:x2       -0.7793 0.3066    -1.3799  -0.7796    -0.1774 0e+00
>
> Random effects:
> Name      Model
>  plate   IID model 
>
> Model hyperparameters:
>                     mean     sd       0.025quant 0.5quant 0.975quant
> Precision for plate 18413.08 18280.69  1217.90   13003.80 66486.61  <-

> ---------------------------------------- Mean=18,413.08. sd=18280.69
>
>
> --
> 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 post to this group, send an email to r-inla-discussion-
> gr...@googlegroups.com.
> Visit this group at
> http://groups.google.com/group/r-inla-discussion-group?hl=en-GB.
> For more options, visit https://groups.google.com/groups/opt_out.
>  
>  

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