inla stack and weight (multiple likelihood model)

[Alok Bohara]

Hi:

I am trying to reproduce a meb type model (measurement error) using the example data --seediling (Muff et al..),  In particular, I am trying to use the inla's stack option.  The model is of type: 

-w = -x + u in the second equation where weight.x is introduced in the original paper to account for the negative in front of x.  I think there is some problem in the use of the weight.x in my case below,  I get estimate of beta.x = -.0002  instead of ,34.  (By the way, when I do    w = x +u  and remove weight.x )and use w.red instead of -w.red etc.., I get the correct result for beta.x. I experimented with all sorts of ways to introduce weight.x but I could not get the correct result. I am curious if the  weight option is not supported by inla.stack ot if there is some trick I am missing below in my script. )  

###############################################################################

### Stack Version

####################################################################

#

stk_1 = inla.stack(data = list(Y = cbind(y, NA)),

A = list(1,1,1,1),

effects = list(

list(beta.0 = rep(1,n)),

list(beta.x = sh),

list(beta.z = z),

              list(gamma = 1:n)),

tag = "est_1"

)

stk_2 = inla.stack(data = list(Y = cbind(NA, -w.red)),  

A = list(1,1),

effects = list(

list(idx2.x = 1:s),

              list(weight2.x = rep(-1,s))),

tag = "est_2"

)

stk2 <- inla.stack(stk_1,stk_2)

formula2 <- Y ~  beta.0 - 1 +

  f(beta.x, copy = "idx2.x", values = 1:n,

    hyper = list(beta = list(param = prior.beta, fixed = FALSE))) +

  f(idx2.x, weight2.x, model = "iid",  values = 1:s,

    hyper = list(prec = list(initial = -15, fixed = TRUE))) +

  beta.z +

  f(gamma, model = "iid", values = 1:n,

    hyper = list(prec = list(initial = log(prec.tau), param = prior.tau)))

r.copy2 <- inla(formula2, data = inla.stack.data(stk2), 

               family = c("poisson", "gaussian"),

control.predictor=list(compute=TRUE, A = inla.stack.A(stk2)),

               control.family = list(

                 list(hyper = list()),

                 list(hyper = list(

                   prec = list(

                     initial=log(prec.u),

                     param = prior.prec.u,

                     fixed = FALSE)))),

               control.fixed = list(

                 mean.intercept = prior.beta[1],

                 prec.intercept = prior.beta[2],

                 mean = prior.beta[1],

                 prec = prior.beta[2])

)

r.copy2 <- inla.hyperpar(r.copy2)

summary(r.copy2)

plot(r.copy2)

^^^^^^^^^^^^^^^^^^^^^^^^^

The original version Nuff et al example using the data.joint option:

Y <- matrix(NA, n+s, 2)

Y[1:n, 1] <- y

Y[n+(1:s), 2] <-  -w.red   ## -w.red was the original model..  meb type

beta.0 <- c(rep(1, n), rep(NA, s))

beta.x <- c(sh, rep(NA, s))

idx.x <- c(rep(NA, n), (1:s))

weight.x <- c(rep(NA, n), -rep(1, s))

beta.z <- c(z, rep(NA, s))

gamma <- c(1:n, rep(NA, s))

data.joint <- data.frame(Y, beta.0, beta.x, idx.x,  beta.z, gamma)

formula <- Y ~  beta.0 - 1 +

  f(beta.x, copy = "idx.x", values = 1:n,

    hyper = list(beta = list(param = prior.beta, fixed = FALSE))) +

  f(idx.x, weight.x, model = "iid", values = 1:s,

    hyper = list(prec = list(initial = -15, fixed = TRUE))) +

  beta.z +

  f(gamma, model = "iid", values = 1:n,

    hyper = list(prec = list(initial = log(prec.tau), param = prior.tau)))

etc...

[Elias T Krainski]

Hi Alok,

On 13/12/15 15:29, Alok Bohara wrote:
> I am curious if the weight option is not supported by inla.stack ot
> if there is some trick I am missing below in my script. )

The inla.stack() just organize the data. By the way, you can simplify
the inla.stack arguments as. Instead list(1,1,1,1) you can have only
list(1) and consider one length list effects as well as they all have
the same projection.

I didn't understood which model you are trying to fit... For example,
how the 'sh' index set relates to the copied '1:n' index set.

Elias

[Alok Bohara]

Hi Elias

n index refers to repeated observation 1 1 1 1    2 2 2 2 ... .   15 15 15 15 , where 1, 2, ... 15 are station ids.  sh, on the other hand, goes from 1 2 3 4... 15 (station id).    The model I am trying to fit is

E(y) = h(b0+b1*x + b2*z)  ,  where y is a poisson, x is a random term... 

-w = -x + u                        , w is gaussian 

The code for this model at the bottom of my earlier posting given in Muff et al works. Trying to estimate that model using the stack command is giving me the trouble.  

where x is measured with error and so w is a proxy.  y is observed 4 times for each station.  w is observed once for each of the 15 stations. Thus they have different length   (For simplicity, I aggregated the y data for 15 stations (n = sh) and estimated the model, but I was still not successful to handle the weight which was needed here because of the - in front of x.)

As for the (1), I agree.  I was just doing it following the Inla R book's code.  

Thanks. 

Alok 

[Elias T Krainski]

On 15/12/15 15:33, Alok Bohara wrote:
> E(y) = h(b0+b1*x + b2*z) , where y is a poisson, x is a random term...
> -w = -x + u , w is gaussian

How about using
w = x + u
instead? As "u" is a zero mean Gaussian error term.

Elias

[Alok Bohara]

I did that and deleted the weight variable and it worked.  But, I wanted to learn to use the weight option (within the stack) by estimating this model: -w = -x + u because there are other instances where weight variable/index need to be used to handle the negative in front of the index.  This is the case of the following 3-equation model Muff (et al.):   http://arxiv.org/pdf/1302.3065v2.pdf   (eq 7,8, and 9)   

E(y) =  h(b1*x  + b2*z)   poisson

0 = -x + a0 +  b3 *z +e    0 mean gaussian

w = x + u                         w is gaussian    

I could reproduce this model using their script provided in the appendix, but I wanted to use the stack option with the weight to reproduce this model.   Thanks. 

Best,

Alok 

[Elias T Krainski]

It shouldn't be a problem to handle weights with stack. See this example:

n = 10
stack <- inla.stack(
data=list(y=rpois(n, 1+1:n/n)),
effects=list(list(idx=1:n, weight=rep(-1,n))),
A=list(1))
dat <- inla.stack.data(stack)
str(dat)

r <- inla(y ~ f(idx, weight, model='rw1'),
family='poisson', data=dat)

cor(dat$y, r$summary.ran$idx$mean)

Elias

[Alok Bohara]

Thanks Elias.  I think the following trick worked:

A = list(1),

effects = list(

list(idx2.x = 1:s, weight2.x = rep(-1,s))),

Thanks a bunch for your patience!!!

[Alok Bohara]

I ran the 3-equation measurement model (blood pressure example -- Influence of systolic blood pressure on coronary heart disease,  Appendix Muff et al. 2013) one with the setup given in the paper (1st output below).  I then used the stack option and reestimated the same model (second output).  I used the same priors.   The results and coefficients are very similar (I am glad, finally to see stack option converge). But why are the marginal likelihood so different?  (-11688.31  VERSUS -17741.01) 

*********************************************

Time used:

 Pre-processing    Running inla Post-processing           Total 

         0.4461      18317.4020          0.3255      18318.1736 

Fixed effects:

           mean     sd 0.025quant 0.5quant 0.975quant    mode kld

beta.0  -2.3439 0.2699    -2.9026  -2.3338    -1.8419 -2.3133   0

beta.z   0.3838 0.2998    -0.1839   0.3761     0.9948  0.3606   0

alpha.0  0.0145 0.0186    -0.0221   0.0145     0.0511  0.0145   0

alpha.z -0.0196 0.0216    -0.0621  -0.0196     0.0229 -0.0196   0

Random effects:

Name      Model

 idx.x   IID model 

beta.x   Copy 

Model hyperparameters:

                                              mean     sd 0.025quant 0.5quant 0.975quant    mode

Precision for the Gaussian observations[2]  21.068 1.4261    18.4176   21.016     24.014  20.911

Precision for the Gaussian observations[3] 102.329 9.9572    83.7666  102.030    122.799 101.604

Beta for beta.x                              2.027 0.6076     0.8392    2.024      3.228   2.015

Expected number of effective parameters(std dev): 645.38(0.0171)

Number of equivalent replicates : 2.98 

Marginal log-Likelihood:  -11688.31 

***********************************************************

***********************************************************

Time used:

 Pre-processing    Running inla Post-processing           Total 

           0.39        82158.55            1.45        82160.39 

Fixed effects:

           mean     sd 0.025quant 0.5quant 0.975quant    mode kld

beta.0  -2.3439 0.2700    -2.9027  -2.3338    -1.8418 -2.3133   0

beta.z   0.3838 0.2998    -0.1840   0.3761     0.9949  0.3606   0

alpha.0  0.0145 0.0195    -0.0238   0.0145     0.0529  0.0145   0

alpha.z -0.0196 0.0224    -0.0636  -0.0196     0.0244 -0.0196   0

Random effects:

Name      Model

 idx2.x   IID model 

beta.x   Copy 

idx3.x   Copy 

Model hyperparameters:

                                              mean     sd 0.025quant 0.5quant 0.975quant    mode

Precision for the Gaussian observations[2]  21.084 1.4285    18.4121   21.037     24.025  20.946

Precision for the Gaussian observations[3] 102.331 9.9475    83.8315  102.008    122.801 101.513

Beta for beta.x                              2.027 0.6085     0.8376    2.024      3.231   2.015

Expected number of effective parameters(std dev): 645.38(0.0174)

Number of equivalent replicates : 2.98 

Marginal log-Likelihood:  -17741.01 

Posterior marginals for linear predictor and fitted values computed

******************************************************************************************************

[Håvard Rue]

This is now fixed in the newest testing version. thanks. 


> r1$mlik - r2a$mlik
                                                  [,1]
log marginal-likelihood (integration) -0.0004065311523
log marginal-likelihood (Gaussian)    -0.0001486655847
> inla.version()


INLA version ............: 0.0-1450401880
INLA date ...............: Fri 18 Dec 02:24:40 CET 2015
INLA hgid ...............: hgid: e18ca7638df2  date: Fri Dec 18
02:02:55 2015 +0100



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 post to this group, send email to r-inla-discussion-group@googlegr
> oups.com.
> Visit this group at https://groups.google.com/group/r-inla-discussion
> -group.
> For more options, visit https://groups.google.com/d/optout.

--
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

[Alok Bohara]

Thanks Havard.

On a related note:

Is there anyway we could compare models such as:  1) a poisson model of Y (with covariates x1 and x2) without the measurement error  VERSUS 2) the same poisson model but with measurement error in x1.  (I also want to add to it spatial effect besag versus spde/mesh and see which one gives me the most preferred specification.)  There will be two or three three likelihoods in the second model as opposed to just one in the first, so the total marginal log likelihood value will be naturally different.   Is there something else I could use in this Bayesian set-up to compare these (non-nested?) models?  

Alok

[Alok Bohara]

Am I doing  something  wrong here, where I get .00000 value for the a1 coefficient.   I allow zhi and u to follow spatial correlation (using spde/mesh), and I use copy and weight options.   I estimated this model with with and without the scale.model option.  Same problem...

require(INLA)

#source("http://www.math.ntnu.no/inla/givemeINLA-testing.R")

inla.upgrade(testing=TRUE)

inla.setOption(scale.model.default=TRUE)

#y2 = b0 + b1* zhi 

#0 = -zhi + a0 + a1*z1

#x2 = zhi + u 

stk_31 = inla.stack(data = list(Y = cbind(y2, NA, NA)),

A = list(A.est,1,1,1,1),

effects = list(

list(s1.field=spde$n.spde, weight1.s1 = rep(1,spde$n.spde)),

list(beta.01 = rep(1,n)),

list(beta.x1 = 1:n),  ## index for the errored var with b_x coeff

            list(beta.x11 = cx1),

        list(idx1.x1 = 1:n, weight1.x1 = rep(1,n))),

tag = "est_1"

)

stk_32 = inla.stack(data = list(Y = cbind(NA, zeros, NA)),

A = list(A.est,1,1,1),

effects = list(

list(s2.field=spde$n.spde,weight2.s2 = rep(-1,spde$n.spde)),

list(alpha.01 = rep(1,n)),

             list(alpha.z11 = cz1),

list(idx2.x1 = 1:n, weight2.x1 = rep(-1,n))),

tag = "est_2"

)

stk_33 = inla.stack(data = list(Y = cbind(NA, NA, x2)),

A = list(A.est,1),

effects = list(

list(u.field= spde$n.spde, s3.field=spde$n.spde,

                    weight3.s3 = rep(1,spde$n.spde)),

list(idx3.x1 = 1:n, weight3.x1 = rep(1,n))),

tag = "est_3"

)

stk_pool <-inla.stack(stk_31,stk_32,stk_33)

#Model with 2 equations (measurement err. r.e are spatial)

formula3b <- Y ~ -1 + 

f(s1.field, weight1.s1, copy="s3.field", fixed=FALSE, 

hyper=list(theta=list(param=c(-1, 10)))) +

f(s2.field, weight2.s2, copy= "s3.field") +   

  f(s3.field, weight3.s3, model = spde) +    ## , copy = "s2.field") + 

f(u.field, model=spde)  +              

beta.01 + beta.x11 +  alpha.01 + alpha.z11 ## alpha.01 + alpha.z11 

precprior <- list(theta=list(param=c(1, 0.1)))

r_mesh <- inla(formula3b, family=c("gaussian","gaussian","gaussian"),

               data=inla.stack.data(stk_pool),  

               control.predictor=list(compute=TRUE, A=inla.stack.A(stk_pool)),

               control.family=list(list(hyper=precprior), list(hyper=precprior),list(hyper=precprior))  )

summary(r_mesh)

Call:

c("inla(formula = formula3b, family = c(\"gaussian\", \"gaussian\", ",  "    \"gaussian\"), data = inla.stack.data(stk_pool), control.predictor = list(compute = TRUE, ",  "    A = inla.stack.A(stk_pool)), control.family = list(list(hyper = precprior), ",  "    list(hyper = precprior), list(hyper = precprior)))")

Time used:

 Pre-processing    Running inla Post-processing           Total 

         2.2969          7.1167          0.8348         10.2485 

Fixed effects:

                mean     sd 0.025quant 0.5quant 0.975quant    mode kld

beta.01      16.6612 1.3807    13.8080  16.7393    18.9168 17.8443   0

beta.x11      1.6195 0.2311     1.1619   1.6195     2.0764  1.6196   0

alpha.01     -3.8777 0.6849    -5.2148  -3.8845    -2.5024 -3.8982   0

I(alpha.z11)  0.0000 0.0142    -0.0282   0.0000     0.0281  0.0000   0

Random effects:

Name      Model

 s3.field   SPDE2 model 

s1.field   Copy 

s2.field   Copy 

Model hyperparameters:

                                              mean      sd 0.025quant 0.5quant

Precision for the Gaussian observations     0.6914  0.2188     0.3442   0.6655

Precision for the Gaussian observations[2] 99.3691 31.4423    49.4677  95.6421

Precision for the Gaussian observations[3]  0.1156  0.0357     0.0586   0.1115

Theta1 for s3.field                        -3.5886  1.5549    -6.3335  -3.7087

Theta2 for s3.field                         1.4720  1.2773    -1.2671   1.5703

Beta for s1.field                          -1.0069  0.3160    -1.6263  -1.0071

                                           0.975quant    mode

Precision for the Gaussian observations        1.1938  0.6164

Precision for the Gaussian observations[2]   171.5494 88.5983

Precision for the Gaussian observations[3]     0.1973  0.1037

Theta1 for s3.field                           -0.2428 -4.0196

Theta2 for s3.field                            3.7437  1.8313

Beta for s1.field                             -0.3851 -1.0078

Expected number of effective parameters(std dev): 5.01(0.047)

Number of equivalent replicates : 11.98 

Marginal log-Likelihood:  -3121.08 

Posterior marginals for linear predictor and fitted values computed

[Håvard Rue]

On Wed, 2015-12-23 at 09:23 -0800, Alok Bohara wrote:
> Am I doing  something  wrong here, where I get .00000 value for the
> a1 coefficient.   I allow zhi and u to follow spatial correlation
> (using spde/mesh), and I use copy and weight options.   I estimated
> this model with with and without the scale.model option.  Same
> problem...
>

Hi

maybe try with a simulated dataset to check if you're doing something
'wrong'. you might in a such check, assume known noise variance.


Best,
H

[Alok Bohara]

Hi Havard:

Yes, I did use the simulated data too (similar to INLA book page 261/262...).  I considered two types of spatial structures: 1)  zhi is allowed to have spde/mesh spatial err structure  (see my codes above), but I get a1 = -.000001,--a0 is usually something reasonable.  I dropped a0 , but the a1 was still .00001.  So, I was wondering if my codes above have some errors.  2)  The second option, however, works and gives me something reasonable for a1, where I considered zhi as an unstructured iid, and allowed "e" and/or "u" as spatially correlated... 

#y2 = b0 + b1* zhi 

#0 = -zhi + a0 + a1*z1 +  e

#x2 = zhi + u  

Thanks.

Alok

[Elias T Krainski]

Hi Alok,

On 28/12/15 15:19, Alok Bohara wrote:
>
> Yes, I did use the simulated data too (similar to INLA book page
> 261/262...). I considered two types of spatial structures: 1) zhi is
> allowed to have spde/mesh spatial err structure (see my codes above),
> but I get a1 = -.000001,--a0 is usually something reasonable. I
> dropped a0 , but the a1 was still .00001. So, I was wondering if my
> codes above have some errors. 2) The second option, however, works
> and gives me something reasonable for a1, where I considered zhi as an
> unstructured iid, and allowed "e" and/or "u" as spatially correlated...
>
> #y2 = b0 + b1* zhi
> #0 = -zhi + a0 + a1*z1 + e
> #x2 = zhi + u

I didn't understood the a1*z1 term... Is it zhi or z1 in the last
equation? What you want different from the model in page 261 of the INLA
book?

Elias

[Alok Bohara]

z1 is a covariate...    This model is very similar to 4.1: Classical Measurement 3-equation model described on page 9: in Bayesian analysis of measurement error models using INLA by Muff, Riebler,  Rue....  except here I am adding spatial to the measurement variable like the one on page 261 (INLA book). 

[Elias T Krainski]

I have added an example on this in the draft SPDE tutorial

[Alok Bohara]

Thanks Elias.  I looked at your SPDE Tutorial example and estimated it. A very interesting set-up which seems different than the one given in Muff et al page 10. Some observations below:

1. In the third equation  stk.o,  your set-up has 0 = a.c + b.w*W -W + m, where W is an observed proxy for unknown "c".  I noticed that the b.w*W is treated as a fixed part  and estimated by introducing W in the formula, whereas "-W" is treated as an iid index. I wonder why the same variable W is treated differently. 

2.  I added an extra covariate  Z in   C  = a.c + b.w * W + m + b.z*Z  and estimated b.z.  No problem.  

3.  In your first equation, y = a.y + b.c*C + x +e, the spatial correlation is introduced in x, whereas C is treated as iid index. I  wonder, if the two could be switched -- i.e., C to be treated as a spatial variable and x as an iid index. This will be very similar to what's presented on page 261.  My attempt to treat C (in b.c*C) as a Marten GF especially with an extra covariate Z is giving me a lot of problem:  

That is,  could we estimate the following model?

y = a.y + b.c*C + x + e  (where x is iid index to capture extra variation, C is Marten GF)

0 = -C + a.c + a.z * Z    (Z is an observed covariate)

W = C + u   (where u could be be an iid index)

This specification is similar to Muff page 10 with the flavor of page 261 in INLA book.

Best,

Alok

[Elias T Krainski]

On 31/12/15 15:44, Alok Bohara wrote:
> 1. In the third equation stk.o, your set-up has 0 = a.c + b.w*W -W +
> m, where W is an observed proxy for unknown "c". I noticed that the
> b.w*W is treated as a fixed part and estimated by introducing W in
> the formula, whereas "-W" is treated as an iid index. I wonder why the
> same variable W is treated differently.

you may think -c instead -w... since they are not the same think. I'll
fix it.

>
> 2. I added an extra covariate Z in C = a.c + b.w * W + m + b.z*Z
> and estimated b.z. No problem.

good

>
> 3. In your first equation, y = a.y + b.c*C + x +e, the spatial
> correlation is introduced in x, whereas C is treated as iid index. I
> wonder, if the two could be switched -- i.e., C to be treated as a
> spatial variable and x as an iid index. This will be very similar to
> what's presented on page 261. My attempt to treat C (in b.c*C) as a
> Marten GF especially with an extra covariate Z is giving me a lot of
> problem:
>
> That is, could we estimate the following model?
>
> y = a.y + b.c*C + x + e (where x is iid index to capture extra
> variation, C is Marten GF)
> 0 = -C + a.c + a.z * Z (Z is an observed covariate)
> W = C + u (where u could be be an iid index)
>

I think that the measurement error here 'u' has to be in the second
equation to have
C = a.c + a.z*Z + u
otherwise it does not makes sense (to me).

Elias

[Alok Bohara]

Thanks.    The model I picked was from Muff et al. page 10, and yes there was a typo on mine.   Yes, I added e2 like Muff. 

Similar to Muff et al. traditional measurement error model  with a slight twist (page 10):

y = a.y + b.c*C + x + e1  (where x is iid index to capture extra iid  variation, C is assumed Marten GF) 
0 = -C + a.c + a.z * Z  + e2   (Z is an observed covariate, e2 is distributed as Normal) 
W = C + u   (where u could be be an iid index to capture some extra variation) 

(I was able to estimate above model with x as a Marten GF, and C as an iid.  I could not make it to work when the two were switched.)

As for the new  model in SPDE tutorial with -W or -C , I have the following observation:

y = a.y + b.c * C + x + e

C = a.c + b.w*W + m

0 = a.c + b.w*W + m - C    (replacing -W with -C as per your feedback)

But, the eq 3 is just eq 2 written differently.  Is it not just a 2-equation system?   Or am I missing  something here?

Thanks for your patience and help...

[Alok Bohara]

Just curious to see if someone has any idea about the question below.

Alok