plotting spatial residuals

[azi...@uw.edu]

Hi all - 

I'm interested in plotting my spatially smoothed residuals from a SPDE INLA model fit. 

My data is distributed randomly in space and to get my 'fitted values' I add in a datastack at the same locations as my original data but with NAs listed as the outcome. After running INLA, I can obtain my fitted values at my original data locations (and on my SPDE mesh) and calculate my residuals at my data locations. My question is, is there any fast or easy way to interpolate and plot these residual estimates (which are not on the SPDE mesh or on any grid) onto a regular grid over spatial my domain with the built-in INLA functions? Maybe something that automatically performs Kriging to a grid using the fitted INLA model? Pointers to a function or a tutorial that I've missed would be greatly appreciated.

much thanks,
Aaron

[Finn Lindgren]

On 20 April 2015 at 22:56, <azi...@uw.edu> wrote:
> My data is distributed randomly in space and to get my 'fitted values' I add in a datastack at the same locations as
> my original data but with NAs listed as the outcome.

The result of that should be (nearly) identical to the fitted values
given for the actual data values, since it's just the linear
predictor. See below for what you more realistically can use.

> my data locations. My question is, is there any fast or easy way to
> interpolate and plot these residual estimates (which are not on the SPDE
> mesh or on any grid) onto a regular grid over spatial my domain with the
> built-in INLA functions?

Since the residuals are defined only for the data locations, and are
essentially independent (in the perfect model case) it doesn't make
sense to interpolate them. For plotting residuals,
fields::quilt.plot() is a good option.

The field itself can however easily be interpolated onto a grid; see
inla.mesh.project() and inla.mesh.projector() (or more low-level by
directly multiplying with a matrix constructed with
inla.spde.make.A()).

Depending what you're going to do with the residuals, you may want to
consider leave-one-out crossvalidated residuals instead of the raw
residuals. These can, for Gaussian observation errors, be obtained
approximately by a
simple postprocessing (exact calculations are more tricky, and may
involve sampling from the posterior). See Section
2.3 of
http://people.bath.ac.uk/fl353/tmp/gmrf.pdf
for the needed calculation, where y_i - E_{(i)} are the
cross-validated residuals.
Also note that inla() can calculate a type of cross-validated residual
that also handles the non-Gaussian observation case; see
?control.inla, option "cpo".

> Maybe something that automatically performs Kriging
> to a grid using the fitted INLA model? Pointers to a function or a tutorial
> that I've missed would be greatly appreciated.

This would make sense in the case that the model did _not_ fully model
the data, so that there is "leftover correlation". But then that
leftover correlation would need to be modelled, so the "correct"
solution in the Bayesian setting would be to run a new model on the
original data, with a new model component capable of handling the
leftover correlation. In a stepwise residual modelling approach, one
would just take the residuals and estimate that extra model component
directly using only the residuals, which could be done with another
inla() model.

Finn

[azi...@uw.edu]

Thanks for the helpful and quick reply! 

best,

Aaron