In timcdlucas/INLAutils: Utility Functions for 'INLA'

INLAutils

A package containing utility functions for the R-INLA package.

There's a fair bit of overlap with inlabru.

Installation

To install, first install INLA.

install.packages("INLA", repos="https://www.math.ntnu.no/inla/R/stable")

then install INLAutils

# From github
library(devtools)
install_github('timcdlucas/INLAutils')

# Load packages
library(INLA)
library(INLAutils)

# Now actually load the packages. This chunk is hidden

library(INLA)
library(INLAutils)

Overview

Plotting

I find the the plot function in INLA annoying and I like ggplot2. So INLAutils provides an autoplot method for INLA objects.

      data(Epil)
      ##Define the model
      formula = y ~ Trt + Age + V4 +
               f(Ind, model="iid") + f(rand,model="iid")
      result = inla(formula, family="poisson", data = Epil, control.predictor = list(compute = TRUE))

      p <- autoplot(result)

Because these are ggplot2 objects, we can easily modify them.

  # Find data names with names(p[[1]]$data)
  p[[1]] + 
    geom_line(aes(colour = var), size = 1.3) +
    palettetown::scale_colour_poke(pokemon = 'Oddish', spread = 4)

There is an autoplot method for INLA SPDE meshes.

    m = 100
    points = matrix(runif(m * 2), m, 2)
    mesh = inla.mesh.create.helper(
      points = points,
      cutoff = 0.05,
      offset = c(0.1, 0.4),
      max.edge = c(0.05, 0.5))

    autoplot(mesh)

There are functions for plotting more diagnostic plots.

 data(Epil)
 observed <- Epil[1:30, 'y']
 Epil <- rbind(Epil, Epil[1:30, ])
 Epil[1:30, 'y'] <- NA
 ## make centered covariates
 formula = y ~ Trt + Age + V4 +
          f(Ind, model="iid") + f(rand,model="iid")
 result = inla(formula, family="poisson", data = Epil,
               control.predictor = list(compute = TRUE, link = 1))
 ggplot_inla_residuals(result, observed, binwidth = 0.1)
 ggplot_inla_residuals2(result, observed, se = FALSE)

Finally there is a function for combining shapefiles, rasters (or INLA projections) and meshes. For more fine grained control the geoms defined in inlabru might be useful.

# Create inla projector
n <- 20
loc <- matrix(runif(n*2), n, 2)
mesh <- inla.mesh.create(loc, refine=list(max.edge=0.05))
projector <- inla.mesh.projector(mesh)

field <- cos(mesh$loc[,1]*2*pi*3)*sin(mesh$loc[,2]*2*pi*7)
projection <- inla.mesh.project(projector, field)

# And a shape file
crds <- loc[chull(loc), ]
SPls <- SpatialPolygons(list(Polygons(list(Polygon(crds)), ID = 'a')))

# plot
ggplot_projection_shapefile(projection, projector, SPls, mesh)

Analysis

There are some helper functions for general analyses.

INLAstep runs stepwise variable selection with INLA.

  data(Epil)
  stack <- inla.stack(data = list(y = Epil$y),
                      A = list(1),
                      effects = list(data.frame(Intercept = 1, Epil[3:5])))

  result <- INLAstep(fam1 = "poisson", 
                     Epil,
                     in_stack = stack,
                     invariant = "0 + Intercept",
                     direction = 'backwards',
                     include = 3:5,
                     y = 'y',
                     y2 = 'y',
                     powerl = 1,
                     inter = 1,
                     thresh = 2)

  result$best_formula
  autoplot(result$best_model, which = 1)

makeGAM helps create a function object for fitting GAMs with INLA.

 data(Epil)
 formula <- makeGAM('Age', invariant = '', linear = c('Age', 'Trt', 'V4'), returnstring = FALSE)
 formula
 result = inla(formula, family="poisson", data = Epil)

Spatial leave-one-out cross-validation (sloo-cv)

The package INLAutils provides an approach to run sloo-cv for INLA objects.

# generate a dataframe and INLA output for the function
set.seed(10)
coords <- data.frame(long = c(rnorm(70), rnorm(30, 3)), lat = rnorm(100))
x <- data.frame(x1 = rnorm(100), x2 = rnorm(100))# x1 no relat., x2 pos. relat.
y <- x$x2 * 2 + rnorm(100)
dataf1 <- sp::SpatialPointsDataFrame(coords = coords, data = data.frame(y = y, x))
mesh <- INLA::inla.mesh.2d(loc = sp::coordinates(dataf1), max.edge = c(3, 3),cutoff = 1.3)
spde <- INLA::inla.spde2.matern(mesh, alpha=2)#SPDE model is defined
A <- INLA::inla.spde.make.A(mesh, loc = sp::coordinates(dataf1))#projector matrix
dataframe <- data.frame(dataf1) #generate dataframe with response and covariate
modform<-stats::as.formula(paste('y ~ -1+ x1 + x2 + y.intercept + f(spatial.field, model=spde)'))
stk <- INLA::inla.stack(data = list(y=dataframe$y), 
                        A = list(A, 1),
                        effects = list(list(spatial.field=1:spde$n.spde),
                        list(y.intercept = rep(1, length(dataframe$y)),
                             covariate = dataframe[c(-1)])), 
                        tag='est')
out <- INLA::inla(modform, family='normal',Ntrials = 1, data=INLA::inla.stack.data(stk, spde=spde),
                  control.predictor = list(A =INLA::inla.stack.A(stk),link=1),
                  control.compute = list( config=TRUE),control.inla = list(int.strategy='eb'))
out.field <- INLA::inla.spde2.result(out,'spatial.field', spde, do.transf = TRUE)
range.out <- INLA::inla.emarginal(function(x) x, out.field$marginals.range.nominal[[1]])

# parameters for the SLOO process
ss <- 20 # sample size to process (number of SLOO runs)
rad <- min(range.out, max(dist(coords)) / 4) # define the radius of the spatial buffer surrounding the removed point. Make sure it isn't bigger than 25% of the study area (Le Rest 2014)
modform <- y ~ -1+ y.intercept + x1 + x2 + f(spatial.field, model=spde)
alpha <- 0.05 # rmse and mae confidence intervals (1-alpha)

# run the function
cv <- inlasloo(dataframe = dataframe, 
               long = 'long', lat = 'lat',
               y = 'y', ss = ss, 
               rad = rad, modform = modform,
               mesh = mesh, family = 'normal',
               mae = TRUE)