mapMM: Colorful plots of predictions in two-dimensional space.

Description

These functions provide either a map of predicted response in analyzed locations, or a predicted surface. mapMM is a straightforward representation of the analysis of the data, while filled.mapMM uses interpolation to cope with the fact that all predictor variables may not be known in all locations on a fine spatial grid. map_ranef maps a single spatial random effect. These three functions takes an HLfit object as input. mapMM calls spaMMplot2D, which is similar but takes a more conventional (x,y,z) input.

Using filled.mapMM may involve questionable choices. Plotting a filled contour generally requires prediction in non-observed locations, where predictor variables used in the original data analysis may be missing. In that case, the original model formula cannot be used and an alternative model (controlled by the map.formula argument) must be used to interpolate (not smooth) the predicted values in observed locations (these predictions still resulting from the original analysis based on predictor variables). filled.mapMM always performs such interpolation (it does not allow one to provide values for the predictor variables). As a result (1) filled.mapMM will be slower than a mere plotting function, since it involves the analysis of spatial data; (2) the results may have little useful meaning if the effect of the original predictor variables is not correctly represented by this interpolation step. For example, prediction by interpolation may be biased in a way analogous to prediction of temperature in non-observed locations while ignoring effect of variation in altitude in such locations. Likewise, the variance argument of filled.mapMM allows one only to plot the prediction variance of its own interpolator, rather than that of the input object.

map_ranef is free of the limitations of filled.mapMM.

Usage

spaMMplot2D(x, y, z, xrange=range(x, finite = TRUE),
  yrange=range(y, finite = TRUE), margin=1/20, add.map= FALSE, 
  nlevels = 20, color.palette = spaMM.colors, map.asp=NULL,
  col = color.palette(length(levels) - 1), plot.title=NULL, plot.axes=NULL, 
  decorations=NULL, key.title=NULL, key.axes=NULL, xaxs = "i", 
  yaxs = "i", las = 1, axes = TRUE, frame.plot = axes, ...) 
mapMM(fitobject,Ztransf=NULL,coordinates,
  add.points,decorations=NULL,plot.title=NULL,plot.axes=NULL,envir=-3, ...)
filled.mapMM(
  fitobject, Ztransf = NULL, coordinates, xrange = NULL, yrange = NULL, 
  margin = 1/20, map.formula, phi = 1e-05, gridSteps = 41, 
  decorations = quote(points(pred[, coordinates], cex = 1, lwd = 2)),
  add.map = FALSE, axes = TRUE, plot.title = NULL, plot.axes = NULL, 
  map.asp = NULL, variance = NULL, var.contour.args = list(), 
  smoothObject = NULL, return.="smoothObject", ...)
map_ranef(fitobject, re.form, Ztransf=NULL, xrange = NULL, yrange = NULL, 
  margin = 1/20, gridSteps = 41, 
  decorations = quote(points(fitobject$data[, coordinates], cex = 1, lwd = 2)), 
  add.map = FALSE, axes = TRUE, plot.title=NULL, plot.axes=NULL, 
  map.asp = NULL, ...)

Value

filled.mapMM by default returns invisibly the fit object predicting the interpolated response surface; however, for any non-default return. argument (return.="raster" would be recommended to ensure future back-compatibility), it will return a raster of values as a list with elements x, y and z. map_ranef returns invisibly a 3-column matrix containing the spatial coordinates, and the predicted effect z on the linear predictor scale (which is also the scale of the plot, unless a Ztransf is used). mapMM returns invisibly a list with elements x, y and z. Plots are produced as side-effects.

Arguments

fitobject: The return object of a corrHLfit or fitme call.
x,y,z: Three vectors of coordinates, with z being expectedly the response.
re.form: A model formula giving the single random effect term to plot, needed only if there are several spatial random effects in the fitted model. In that case, it must be formatted as . ~ <term>, as for the re.form argument of predict.HLfit.
Ztransf: A transformation of the predicted response, given as a function whose only required argument can be a one-column matrix. The name of this argument must be Z (not x), as is appropriate for use in do.call(Ztransf,list(Z=Zvalues)).
coordinates: The geographical coordinates. By default they are deduced from the model formula. For example if this formula is resp ~ 1 + Matern(1| x + y ) the default coordinates are c("x","y"). If this formula is resp ~ 1 + Matern(1| x + y + z ), the user must choose two of the three coordinates.
xrange: The x range of the plot (a vector of length 2); by default defined to cover all analyzed points.
yrange: The y range of the plot (a vector of length 2); by default defined to cover all analyzed points.
margin: This controls how far (in relative terms) the plot extends beyond the x and y ranges of the analyzed points, and is overriden by explicit xrange and yrange arguments.
map.formula: NULL, or a formula whose left-hand side is ignored. Provides the formula used for interpolation. If NULL, a default formula with the same spatial effect(s) as in the input fitobject is used.
phi: This controls the phi value assumed in the interpolation step. Ideally phi would be zero, but problems with numerically singular matrices may arise when phi is too small.
gridSteps: The number of levels of the grid of x and y values
variance: Either NULL, or the name of a component of variance of prediction by the interpolator to be plotted. Must name one of the components that can be returned by predict.HLfit. variance="predVar" is suitable for uncertainty in point prediction.
var.contour.args: A list of control parameters for rendering of prediction variances. See contour for possible arguments (except x, y, z and add).
add.map: Either a boolean or an explicit expression, enclosed in quote (see Examples). If TRUE, the map function from the maps package (which much therefore the loaded) is used to add a map from its default world database. xrange and yrange are used to select the area, so it is most convenient if the coordinates are longitude and latitude (in this order and in standard units). An explicit expression can also be used for further control.
levels: a set of levels which are used to partition the range of z. Must be strictly increasing (and finite). Areas with z values between consecutive levels are painted with the same color.
nlevels: if levels is not specified, the range of z, values is divided into *approximately* this many levels (a call to pretty determines the actual number of levels).
color.palette: a color palette function to be used to assign colors in the plot.
map.asp: the y/x aspect ratio of the 2D plot area (not of the full figure including the scale). By default, the scales for x and y are identical unless the x and y ranges are too different. Namely, the scales are identical if (plotted y range)/(plotted x range) is 1/4 < . < 4, and map.asp is 1 otherwise.
col: an explicit set of colors to be used in the plot. This argument overrides any palette function specification. There should be one less color than levels
plot.title: statements which add titles to the main plot. See Details for differences between functions.
plot.axes: statements which draw axes (and a box) on the main plot. See Details for differences between functions.
decorations: Either NULL or Additional graphic statements (points, polygon, etc.), enclosed in quote (the default value illustrates the latter syntax). .
add.points: Obsolete, use decorations instead.
envir: Controls the environment in which plot.title, plot.axes, and decorations are evaluated. mapMM calls spaMM2Dplot from where these graphic arguments are evaluated, and the default value -3 means that they are evaluated within the environment from where mapMM was called.
key.title: statements which add titles for the plot key.
key.axes: statements which draw axes on the plot key.
xaxs: the x axis style. The default is to use internal labeling.
yaxs: the y axis style. The default is to use internal labeling.
las: the style of labeling to be used. The default is to use horizontal labeling.
axes, frame.plot: logicals indicating if axes and a box should be drawn, as in plot.default.
smoothObject: Either NULL, or an object inheriting from class HLfit (hence, an object on which predict.HLfit can be called), predicting the response surface in any coordinates. See Details for typical usages.
return.: character string: see Value
...: further arguments passed to or from other methods. For mapMM, all such arguments are passed to spaMMplot2D; for spaMMplot2D, currently only additional graphical parameters passed to title() (see Details). For filled.mapMM and map_ranef, these parameters are those that can be passed to spaMM.filled.contour.

Details

The smoothObject argument may be used to redraw a figure faster by recycling the predictor of the response surface returned invisibly by a previous call to filled.mapMM.

For smoothObject=NULL (the default), filled.mapMM interpolates the predicted response, with sometimes unpleasant effects. For example, if one interpolates probabilities, the result may not be within [0,1], and then (say) a logarithmic Ztransf may generate NaN values that would otherwise not occur. The smoothObject argument may be used to overcome the default behaviour, by providing an alternative predictor.

If you have values for all predictor variables in all locations of a fine spatial grid, filled.mapMM may not be a good choice, since it will ignore that information (see map.formula argument). Rather, one should use predict(<fitobject>,newdata= <all predictor variables >) to generate all predictions, and then either spaMM.filled.contour or some other raster functions.

The different functions are (currently) inconsistent among themselves in the way they handle the plot.title and plot.axes argument:

spaMM.filled.contour behaves like graphics::filled.contour, which (1) handles arguments which are calls such as title(.) or {axis(1);axis(2)}; (2) ignores ... arguments if plot.title is missing; and (3) draws axes by default when plot.axes is missing, given axes = TRUE.

By contrast, filled.mapMM handles arguments which are language expressions such as produced by quote(.) or substitute(.) (see Examples).

mapMM can handles language expressions, but also accepts at least some calls.

Examples

Run this code

data("blackcap")
bfit <- fitme(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
              fixed=list(lambda=0.5537,phi=1.376e-05,rho=0.0544740,nu=0.6286311))
mapMM(bfit,color.palette = function(n){spaMM.colors(n,redshift=1/2)},add.map=TRUE)
map_ranef(bfit) # providing argument re.form= . ~ Matern(1|longitude+latitude)

if (spaMM.getOption("example_maxtime")>1) {
  ## filled.mapMM takes a bit longer
  # showing 'add.map', 'nlevels', and contour lines for 'variance'
  filled.mapMM(bfit, nlevels=30, add.map=TRUE, plot.axes=quote({axis(1);axis(2)}),
             variance="respVar",
             plot.title=title(main="Inferred migration propensity of blackcaps",
                               xlab="longitude",ylab="latitude"))
                               
  ## Similar plots by ggplot2:
  if (FALSE) {
   library(rnaturalearth) # provides sea mask through 'ne_download' function
   library(ggplot2)
   library(sp)
  
   # sea mask 
   sea <- ne_download(scale = 10, type = 'ocean', category = "physical", returnclass = "sf")
   
   # Generation of data.frame for ggplot:
   rastr <- filled.mapMM(bfit, return.="raster")
   spdf <- data.frame(Long=rep(rastr$x, nc), Lat=rastr$y[gl(nr,nc)], z = as.vector(rastr$z))

   ggplot(spdf) + 
     geom_contour_filled(aes(Long,Lat,z=z), bins = 20) +
     guides(fill = "none") +
     geom_sf(data = sea, fill = "black") +
     coord_sf(ylim = range(rastr$y), xlim = range(rastr$x), expand = FALSE)
  }
                               
}

if (spaMM.getOption("example_maxtime")>3) {
 data("Loaloa")  
 lfit <- fitme(cbind(npos,ntot-npos)~elev1+elev2+elev3+elev4+maxNDVI1+seNDVI
                  +Matern(1|longitude+latitude), method="PQL", data=Loaloa,
                  family=binomial(), fixed=list(nu=0.5,rho=2.255197,lambda=1.075))   

 ## longer computation requiring interpolation of 197 points 
 ## Also illustrating effect of 'return.' argument
 res <- filled.mapMM(lfit,add.map=TRUE,plot.axes=quote({axis(1);axis(2)}),
    decorations=quote(points(pred[,coordinates],pch=15,cex=0.3)),
    return.="raster", # so that 'res' is a list representing a raster. 
    plot.title=title(main="Inferred prevalence, North Cameroon",
                     xlab="longitude",ylab="latitude"))
}