Specify Prior for control.fixed

[hon...@uw.edu]

I have the following poisson model y ~ Possion(beta0+beta1x). Now I want to specify the normal prior for beta0 and beta1, however, , the results are different depending on whether I use list() in mean=list(-0.02) and prec=list(1/0.10^2). As Blangiardo & Cameletti (2015), they should be the same.Which result should I believe? Thank you!

mod.pois1 <- inla(newy ~ 1 + newx1, 

                  family = 'poisson',

                  E = newE,

                  data = lungradon1,

                  control.fixed = list(mean=-0.02,    # mean=list(-0.02),

                                       prec=1/0.10^2, # prec=list(1/0.10^2),

                                       mean.intercept=0,

                                       prec.intercept=1/0.21^2)

)

summary(mod.pois1)

[hon...@uw.edu]

The full code is as follows:

 

lung <- read.table("http://faculty.washington.edu/jonno/book/MNlung.txt", header=TRUE, sep="\t")

radon <- read.table("http://faculty.washington.edu/jonno/book/MNradon.txt", header=TRUE)

Obs <- apply(cbind(lung[,3], lung[,5]), 1, sum)

Exp <- apply(cbind(lung[,4], lung[,6]), 1, sum)

rad.avg <- rep(0, length(lung$X))

for(i in 1:length(lung$X)) {

  rad.avg[i] <- mean(radon[radon$county==i,2])

}

x <- rad.avg

rad.avg[26]<-0

rad.avg[63]<-0

x[26] <- NA

x[63] <- NA

newy <- Obs[is.na(x)==F]

newx <- x[is.na(x)==F]

newE <- Exp[is.na(x)==F]

newx1 <- scale(newx, center = T, scale = F)

lungradon1 <- list(newy = newy, newx1 = newx1, newE = newE)

library(INLA)

[INLA help]

On Sat, 2016-11-05 at 14:13 -0700, hon...@uw.edu wrote:
> I have the following poisson model y ~ Possion(beta0+beta1x). Now I
> want to specify the normal prior for beta0 and beta1, however, , the
> results are different depending on whether I
> use list() in mean=list(-0.02) and prec=list(1/0.10^2). As Blangiardo
> & Cameletti (2015), they should be the same.Which result should I
> believe? Thank you!

Hi,

if you use a  mean=list() it must be named. [the question is if inla()
should throw an error if you don't...], so like

inla(y~1+a+b, control.fixed = list(mean = list(a=-1111)), 
data = data.frame(y=0:1, a=1:2, b=4:3))

the  named elements of the list, like  'a' and 'b', refer to the name
of the variables in the formula, while 'default' is referred to as the
default value, if not mentioned in the list. if you pass an NULL name
(ie no name), then it is simply ignored. 



you can check the priors used either reading the output using
verbose=T, or by result$logfile, or check the result object all.hyper
and then 'fixed', like

> inla(y~1+a+b, control.fixed = list(mean = list(a=-1111)),  data =
data.frame(y=0:1, a=1:2, b=4:3))$all.hyper$fixed

[[1]]
[[1]]$label
[1] "(Intercept)"

[[1]]$prior.mean
[1] 0

[[1]]$prior.prec
[1] 0


[[2]]
[[2]]$label
[1] "a"

[[2]]$prior.mean
[1] -1111

[[2]]$prior.prec
[1] 0.001


[[3]]
[[3]]$label
[1] "b"

[[3]]$prior.mean
[1] 0

[[3]]$prior.prec
[1] 0.001



hope this helps.
H

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

[hon...@uw.edu]

Hi Håvard,

Thank you for quick reply! In control.fixed, the default type of prior is Normal(Gaussian), right? How can I change the type of prior? 

Thanks again!

-Hongjian

[INLA help]

On Sun, 2016-11-06 at 10:56 -0800, hon...@uw.edu wrote:
> Hi Håvard,
>
> Thank you for quick reply! In control.fixed, the default type of
> prior is Normal(Gaussian), right? How can I change the type of
> prior? 
> 

Hi

you cannot, as the latent field must be Gaussian. you can, for a few of
them, throw the regression coefficient into the hyperparameters and
then you're more free. if the mean is zero, you can instead of

y ~ ... + x + 

use

y ~ ... + f(beta.x, x, model="iid") + ...

where beta.x = rep(1,n)

then you can control the prior of the variance to beta.x, which,
integrated out, gives the non-gaussian prior for beta.x.

otherwise, you can use "rgeneric" to test it out; see the doc for
examples. 

Best