Variance/covariance of predictions in joint gaussian model

[sbh...@gmail.com]

I've fit a joint spatial model with two gaussian likelihoods based on the example in Ch 3.1 of Advanced Spatial Modeling by Krainski. When I summarize samples from the posterior, the mean values line up with what I know about the data, but the standard deviations are much lower and the covariance of the two responses does not line up with the data. My data (2200 rows) has multiple observations for many locations (up to 500 observations), so when I say the SD seems low I'm comparing the SD of the data at a given location to the SD of the posterior sample at that location.

Things I've attempted to no avail:

- Widening the spde prior - no change

- Widening the pc.prec prior on the error term - no change

- Making the mean.linear and prec.linear terms on the intercepts less informative, which made the mean predictions a lot worse

Included below is my code, any feedback would be appreciated, thanks!

spde <- inla.spde2.pcmatern(Mesh1, alpha = 2,

                            prior.range=c(3, 0.05),

                            prior.sigma=c(6, 0.05)

)

s_index_horz <- inla.spde.make.index(name="spatial.field.horz",

                                     n.spde=spde$n.spde)

stk.horz <- inla.stack(data=list(y=cbind(data$horz_break,NA)),

                       A=list(A),

                       effects=list(c(s_index_horz, list(intercept1=1))),

                       tag="horz")

s_index_induced <- inla.spde.make.index(name="spatial.field.induced",

                                        n.spde=spde$n.spde)

stk.induced <- inla.stack(data=list(y=cbind(NA,data$induced_vert_break)),

                          A=list(A),

                          effects=list(c(s_index_induced, list(intercept2=1, s12=1:spde$n.spde))),

                          tag="induced")

stk <- inla.stack(stk.horz, stk.induced)

horz_mean <- mean(data$horz_break)

horz_prec <- 1/sd(data$horz_break)^2

induced_mean <- mean(data$induced_vert_break)

induced_prec <- 1/sd(data$induced_vert_break)^2

hyper <- list(theta=list(prior='normal', param=c(0,10)))

formula <- y ~ -1 + 

  f(intercept1, model='linear', prec.linear=horz_prec, mean.linear=horz_mean) +

  f(intercept2, model='linear', prec.linear=induced_prec, mean.linear=induced_mean) +

  f(spatial.field.horz, model = spde) + 

  f(spatial.field.induced, model = spde) +

  f(s12, copy="spatial.field.horz", fixed=FALSE, hyper=hyper)

hyper.eps <- list(hyper = list(theta = list(prior = 'pc.prec', param = c(1, 0.01))))

int.strategy <- 'auto'

fit <- inla(formula, family=c('gaussian', 'gaussian'),

            data = inla.stack.data(stk), 

            control.predictor = list(A=inla.stack.A(stk)),

            control.inla=list(strategy='adaptive',

                              int.strategy=int.strategy),

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

            verbose=TRUE

)

# Prediction

n=1000

samples = inla.posterior.sample(n=n, fit, num.threads=6)

contents=fit$misc$configs$contents

## intercept1

idx_intercept1=which(contents$tag=="intercept1")

idx_intercept1=contents$start[idx_intercept1]:(contents$start[idx_intercept1]+contents$length[idx_intercept1]-1)

intercept1 = lapply(samples, function(x) x$latent[idx_intercept1])

intercept1 = matrix(unlist(intercept1), ncol = length(idx_intercept1), byrow = T)

## intercept2

idx_intercept2=which(contents$tag=="intercept2")

idx_intercept2=contents$start[idx_intercept2]:(contents$start[idx_intercept2]+contents$length[idx_intercept2]-1)

intercept2 = lapply(samples, function(x) x$latent[idx_intercept2])

intercept2 = matrix(unlist(intercept2), ncol = length(idx_intercept2), byrow = T)

## spatial.field.horz

idx_spat.horz=which(contents$tag=="spatial.field.horz")

idx_spat.horz=contents$start[idx_spat.horz]:(contents$start[idx_spat.horz]+contents$length[idx_spat.horz]-1)

spat.horz = lapply(samples, function(x) x$latent[idx_spat.horz])

spat.horz = matrix(unlist(spat.horz), ncol = length(idx_spat.horz), byrow = T)

## spatial.field.induced

idx_spat.induced=which(contents$tag=="spatial.field.induced")

idx_spat.induced=contents$start[idx_spat.induced]:(contents$start[idx_spat.induced]+contents$length[idx_spat.induced]-1)

spat.induced = lapply(samples, function(x) x$latent[idx_spat.induced])

spat.induced = matrix(unlist(spat.induced), ncol = length(idx_spat.induced), byrow = T)

# s12

idx_spat.horz.induced=which(contents$tag=="s12")

idx_spat.horz.induced=contents$start[idx_spat.horz.induced]:(contents$start[idx_spat.horz.induced]+contents$length[idx_spat.horz.induced]-1)

spat.horz.induced = lapply(samples, function(x) x$latent[idx_spat.horz.induced])

spat.horz.induced = matrix(unlist(spat.horz.induced), ncol = length(idx_spat.horz.induced), byrow = T)

A_i <- inla.spde.make.A(Mesh1, loc=as.matrix(loc.domain))

pred.horz <- A_i %*% t(spat.horz)

pred.horz <- intercept1 + t(pred.horz)

pred_horz_mean <- apply(pred.horz, 2, mean)

# Much lower than expected

pred_horz_sd <- apply(pred.horz,2,sd)

pred.induced <- A_i %*% t(spat.induced)

pred.horz.induced <- A_i %*% t(spat.horz.induced)

pred.induced <- intercept2 + t(pred.induced) + t(pred.horz.induced)

pred_induced_mean <- apply(pred.induced, 2, mean)

# Much lower than expected

pred_induced_sd <- apply(pred.induced, 2, sd)

# This looks very different than the covariance of the data

cov(data.frame(x=pred.horz[,1],y=pred.induced[,1]))

[Finn Lindgren]

Hi,

the posterior standard deviation of the linear predictor should not match the observation (residual) standard deviation, especially not when you have multiple observations at some locations.

To do data-level prediction you need to include the observation model in the prediction, not just the "linear predictor".

As an extreme illustration, consider a model with only a single parameter (the intercept). If we ignore the effect of the priors, then the posterior variance of that would be roughly sigma^2/n, where sigma is the observation variance, and n is the number of observations. For full models the relationship won't be that simple.

You can do a cheap approximation of the data level standard deviation by adding the estimated posterior observation standard deviation to the posterior standard deviation of the linear predictor.

(another approximation would be to add the two variances instead. The std.dev. sum is appropriate if the two quantities are positively correlated in the joint posterior distribution)

Finn

--
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/6bf5bbb0-8d49-4411-9476-9c0152a4d622n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Finn Lindgren]

Addition:

Since you're using inla.posterior.sample, you can simply add the information derived directly from the sampled observation model precisions; that will take care of any correlation, since they are sampled from the joint distribution.

Finn

[Helpdesk]

You might have a different correlation structure in the latent than in
the data itself, but I trust you when you say the results are
'weird'/'surprising'.

We need to check this out, so we need to rerun it here. can you provide
data and full R script? you can use the private email he...@r-inla.org

if you have repeated observations, (with the same predictor), you might
benefit from a new likelihood, 'agaussian', see

inla.doc("agaussian")

if its not there, you'll need to update to a recent testing version.

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/6bf5bbb0-8d49-4411-9476-9c0152a4d622n%40googlegroups.com
> .

--

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

[Helpdesk]

aaah, I see Finn have responded to this one, hence the first part of my
answer might not be relevant.

the 'agaussian' might still be, though.

[S. Brian Huey]

Thanks for the advice. I was able to use the observation model to capture the variance/covariance in the data. I'll look into the gaussian likelihood as that seems relevant to my problem. Cheers,

Brian