R/georob_lgnpp.R

In georob: Robust Geostatistical Analysis of Spatial Data

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

#### lgnpp

lgnpp <-
  function(
    object,
    newdata, locations,
    is.block = FALSE, all.pred = NULL,
    extended.output = FALSE
  )
{

  ## Function post-processes predictions obtained from predict.georob
  ## for lognormal point and block kriging

  ## 2011-12-22 A. Papritz
  ## 2012-05-07 AP backtransformation of block predictions under permanence
  ## of lognormality assumption
  ## 2013-06-12 AP substituting [["x"]] for $x in all lists
  ## 2015-06-23 AP modifications for missing compontents 'lower' and 'upper'
  ## 2015-08-22 AP changes for nested variogram models
  ## 2019-12-13 AP correcting use of class() in if() and switch()
  ## 2020-02-14 AP sanity checks of arguments and for if() and switch()

  ## auxiliary function to backtransform the point predictions and
  ## optionally the mean squared errors and the prediction intervals

 f.backtrf <- function( object, btMb = 0., pred.only ){

#### -- auxiliary function

    ## object: dataframe with the prediction results of the
    ##         log-transformed predictions computed with extended.output = TRUE
    ## btMb  : vector with the terms t( coefficients ) M(B) coefficients for back-transformation
    ##         of block kriging results under assumption of permanence of lognormality
    ## pred.only: logical scalar controlling wether only predictions are computed

    t.result <- exp(
      object[, "pred"] + 0.5 * ( object[, "var.target"] + btMb - object[, "var.pred"] )
    )

    if( !pred.only ){
      t.mu <- exp( object[, "trend"] + 0.5 * object[, "var.target"] )
      t.result <- data.frame(
        lgn.pred = t.result,
        lgn.se = t.mu * sqrt(
          exp( object[, "var.target"] ) + exp( object[, "var.pred"] ) -
          2. * exp( object[, "cov.pred.target"] )
        )
      )
      if( all( c( "lower", "upper" ) %in% names( object ) ) ){
        t.result <- cbind( t.result,
          data.frame(
            lgn.lower = exp( object[, "lower"] ),
            lgn.upper = exp( object[, "upper"] )
          )
        )

      }
    }

    return( t.result )

  }

#### -- preparation

  ## check arguments

  if(!missing(newdata)) check.newdata(newdata)

  stopifnot(identical(length(is.block), 1L)        && is.logical(is.block))
  stopifnot(identical(length(extended.output), 1L) && is.logical(extended.output))

  #   if(!missing(newdata)) stopifnot(
  #     all(class(newdata) %in% c(
  #         "data.frame", "SpatialPointsDataFrame", "SpatialPixelsDataFrame",
  #         "SpatialGridDataFrame", "SpatialPolygonsDataFrame"
  #       )
  #     )
  #   )
  #   if(!missing(all.pred)) stopifnot(
  #     all(class(newdata) %in% c(
  #         "data.frame", "SpatialPointsDataFrame", "SpatialPixelsDataFrame",
  #         "SpatialGridDataFrame", "SpatialPolygonsDataFrame"
  #       )
  #     ) || (identical(length(all.pred), 1L) && is.integer(all.pred))
  #   )

  ## check whether back-transformation for object has been already done

  if( all( c( "lgn.pred", "lgn.se" ) %in% names( object ) ) ){
    return( object )
  }

  ## main body of function

  ## find out what object contains

  what <- switch(
    class( object )[1],
    "data.frame" =,
    "SpatialPointsDataFrame" =,
    "SpatialPixelsDataFrame" =,
    "SpatialGridDataFrame" = "point",
    "SpatialPolygonsDataFrame" = "block",
    "list" = switch(
      class( object[["pred"]] )[1],
      "data.frame" =,
      "SpatialPointsDataFrame" =,
      "SpatialPixelsDataFrame" =,
      "SpatialGridDataFrame" = "point",
      "SpatialPolygonsDataFrame" = "block",
      stop( "'object' not generated by predict method for class 'georob'" )
    ),
    stop( "'object' not generated by predict method for class 'georob'" )
  )

  if( what == "point" && ( is.block || !is.null( all.pred ) ) ) what <- "optimal"

  ## extract data frame with prediction results

  t.object <- switch(
    class( object )[1],
    "data.frame" = object,
    "SpatialPointsDataFrame" =,
    "SpatialPixelsDataFrame" =,
    "SpatialGridDataFrame" =,
    "SpatialPolygonsDataFrame" = object@data,
    "list" = switch(
      class( object[["pred"]] )[1],
      "data.frame" = object[["pred"]],
      "SpatialPointsDataFrame" =,
      "SpatialPixelsDataFrame" =,
      "SpatialGridDataFrame" =,
      "SpatialPolygonsDataFrame" = object[["pred"]]@data
    )
  )

  ## extract variogram and type of predictions

  variogram.object   <- attr( t.object, "variogram.object" )
  #   variogram.model   <- attr( t.object, "variogram.model" )
  #   param             <- attr( t.object, "param" )
  type              <- attr( t.object, "type" )

  if( is.null( variogram.object ) || is.null( type ) ) stop(
    "attributes 'variogram.object' or 'type' missing in 'object'"
  )

  if( any(
      sapply( variogram.object, function(x) x[["variogram.model"]])
      %in% control.georob()[["irf.models"]] )
  ) stop(
    "lognormal kriging requires a weakly stationary variogram model"
  )

  ## check whether object is complete

  if( is.null( t.object[["pred"]]) || is.null( t.object[["trend"]] ) ||
    is.null( t.object[["var.pred"]] ) || is.null( t.object[["cov.pred.target"]] ) ||
    is.null( t.object[["var.target"]] )
  ) stop( "some required items are missing, ",
    "re-run 'predict' with argument 'extended.output = TRUE'"
  )

  ## and now do the backtransformation for ...

  result <- switch(

    what,

#### -- backtransformation of point-predictions

    ## back-transform point predictions

    point = {

      t.result <- cbind(
        t.object,
        f.backtrf( t.object, pred.only = FALSE )
      )

      attr( t.result, "variogram.object" ) <- variogram.object
      #       attr( t.result, "variogram.model" ) <- variogram.model
      #       attr( t.result, "param" )           <- param
      attr( t.result, "type" )            <- type

      result <- switch(
        class( object )[1],
        "data.frame" = t.result,
        "SpatialPointsDataFrame" =,
        "SpatialPixelsDataFrame" =,
        "SpatialGridDataFrame" =,
        "SpatialPolygonsDataFrame" = {
          result <- object
          result@data <- t.result
          result
        },
        "list" = switch(
          class( object[["pred"]] )[1],
          "data.frame" = {
            result <- object
            result[["pred"]] <- t.result
            result
          },
          "SpatialPointsDataFrame" =,
          "SpatialPixelsDataFrame" =,
          "SpatialGridDataFrame" =,
          "SpatialPolygonsDataFrame" = {
            result <- object
            result[["pred"]]@data <- t.result
            result
          }
        )
      )

      result

    },

    block = {

#### -- back-transformation of block-predictions

      ## back-transform block predictions under the assumption of permanence
      ## of lognormality

      if( missing( newdata ) ) stop(
        "'newdata' must be provided for back-transformation under",
        " assumption of permance of lognormality "
      )

      ## extract terms, regression coefficients and locations formula

      tt           <- attr( t.object, "terms" )
      coefficients <- attr( t.object, "coefficients" )
      locations    <- attr( t.object, "locations" )

      if( is.null( tt ) || is.null( coefficients ) || is.null( locations ) ) stop(
        "attributes 'terms', 'coefficients' or 'locations' missing in object"
      )

      ## compute spatial covariance matrices of explanatory covariates

      if( !missing( locations ) ){
        locations <- as.formula(
          paste( deparse( locations ), "-1" ), env = parent.frame()
        )
      }

      ## extract modelframe and coordinates of prediction points

      Terms <- delete.response( tt )
      Terms.loc <- locations

      ## get the model frame for newdata

      mf.newdata <- switch(
        class( newdata )[1],
        "data.frame" = model.frame(
          Terms, newdata, na.action = na.pass, xlev = object[["xlevels"]]
        ),
        "SpatialPointsDataFrame" = model.frame(
          Terms, slot( newdata, "data" ), na.action = na.pass,
          xlev = object[["xlevels"]]
        ),
        "SpatialPixelsDataFrame" = model.frame(
          Terms, slot( newdata, "data" ), na.action = na.pass,
          xlev = object[["xlevels"]]
        ),
        "SpatialGridDataFrame" = model.frame(
          Terms, slot( newdata, "data" ), na.action = na.pass,
          xlev = object[["xlevels"]]
        ),
        "SpatialPolygonsDataFrame" = model.frame(
          Terms, slot( newdata, "data" ), na.action = na.pass,
          xlev = object[["xlevels"]]
        ),
        stop(
          "cannot construct model frame for class(newdata) ='",
          class( newdata )
        )
      )

      ## check whether variables that will be used to compute the
      ## predictions agree with those in object

      if( !is.null( cl <- attr(Terms, "dataClasses" ) ) )
        .checkMFClasses( cl, mf.newdata )

      ## get fixed effects design matrix for newdata

      pred.X <- model.matrix( Terms, mf.newdata, contrasts.arg = object[["contrasts"]] )


      #       ## deal with non-NULL offset
      #
      #       offset <- rep( 0, NROW(pred.X) )
      #       if( !is.null( off.num <- attr( tt, "offset" ) ) ){
      #         warning( "prediction with non-zero offset not yet debugged" )
      #         for( i in off.num ) {
      #           offset <- offset + eval( attr( tt, "variables" )[[i + 1L]], newdata )
      #         }
      #       }
      #       if( !is.null( object[["call"]][["offset"]] ) ){
      #         offset <- offset + eval( object[["call"]][["offset"]], newdata )
      #       }

      ## get matrix of coordinates of newdata for point kriging

      pred.coords <- switch(
        class( newdata )[1],
        "data.frame" = model.matrix(
          Terms.loc,
          model.frame(
            Terms.loc, newdata, na.action = na.pass
          )
        ),
        "SpatialPointsDataFrame" = model.matrix(
          Terms.loc,
          model.frame(
            Terms.loc, as.data.frame( coordinates( newdata ) ),
            na.action = na.pass
          )
        ),
        "SpatialPixelsDataFrame" = model.matrix(
          Terms.loc,
          model.frame(
            Terms.loc, as.data.frame( coordinates( newdata ) ),
            na.action = na.pass
          )
        ),
        "SpatialGridDataFrame" = model.matrix(
          Terms.loc,
          model.frame(
            Terms.loc, as.data.frame( coordinates( newdata ) ),
            na.action = na.pass
          )
        ),
        stop(
          "newdata must be of class 'data.frame', 'SpatialPointsDataFrame', 'SpatialPixelsDataFrame'\n",
          "  or 'SpatialGridDataFrame'"
        )
      )

      ## indices of polygons to which prediction points belong to

      t.ip <- over(
        t.points <- SpatialPoints( pred.coords ),
        t.polygons <- if( is.list( object ) ){
          geometry( object[["pred"]] )
        } else {
          geometry( object )
        }
      )

      ## compute t(coefficients) %*% Cov( covariates ) %*% coefficients
      ## cf Cressie, 2006, Eq. 18 & Appendix C

      btMb <- as.vector(
        tapply(
          1L:length( t.ip ),
          factor( t.ip ),
          function( i, x, b ){
            drop( b %*% cov( x[i, , drop = FALSE ] ) %*% b )
          },
          x = pred.X,
          b = coefficients
        )
      )

      ## compute back-transformed block predictions

      t.result <- cbind(
        t.object,
        f.backtrf( t.object, btMb, pred.only = FALSE )
      )

      attr( t.result, "variogram.object" ) <- variogram.object
      #       attr( t.result, "variogram.model" ) <- variogram.model
      #       attr( t.result, "param" )           <- param
      attr( t.result, "type" )            <- type
      attr( t.result, "terms" )           <- tt
      attr( t.result, "coefficients" )    <- coefficients
      attr( t.result, "locations" )       <- locations

      result <- object
      if( class( object )[1] == "SpatialPolygonsDataFrame" ){
        result@data <- t.result
      } else {
        result[["pred"]]@data <- t.result
      }

      result

    },

    optimal = {

#### -- back-transformation for a single block

      ## number of prediction items in object

      n.sample <- sum( complete.cases( t.object ) )

      if( is.null( all.pred ) ){

        ## object contains the predictions of all the points in the block

        n <- n.sample

      } else if( is.numeric( all.pred ) && length( all.pred == 1L ) ){

        ## complete predictions contains total number of points in block

        n <- all.pred

      } else if(
        is.data.frame( all.pred ) ||
        class( all.pred )[1] %in% c(
          "SpatialPointsDataFrame", "SpatialPixelsDataFrame", "SpatialGridDataFrame"
        )
      ){

        ## all.pred contains back-transformed point predictions of
        ## all the points in the block, extract the respective data frame

        all.pred <- switch(
          class( all.pred )[1],
          data.frame = all.pred,
          SpatialPointsDataFrame =,
          SpatialPixelsDataFrame =,
          SpatialGridDataFrame = all.pred@data
        )

        ## check whether back-transformation were done

        if( !all( c( "lgn.pred", "lgn.se", "lgn.lower", "lgn.upper" ) %in% names( all.pred ) ) ){
          all.pred <- f.backtrf( all.pred, pred.only = T )
        }

        ## check wether variogram models in all.pred and object match

        if( !all(
            c(
              identical( variogram.object, attr( all.pred, "variogram.object" ) ),
#               identical( variogram.model, attr( all.pred, "variogram.model" ) ),
#               identical( param, attr( all.pred, "param" ) ),
              identical( type, attr( all.pred, "type" ) )
            )
          )
        ) stop( "variogram or prediction type of 'object' and 'all.pred' differs" )

        n <- sum( complete.cases(all.pred ) )

      } else {

        stop(
          "'all.pred' does neither contain predictions results",
          "nor total number of points in block"
        )

      }

      if( n < n.sample ) stop( "nrow(all.pred) < nrow(object$pred)" )

      ## optimal lognormal backtransformation for block prediction

      if( is.null( object[["pred"]] ) || is.null( object[["mse.pred"]] ) ||
        is.null( object[["var.pred"]] ) || is.null( object[["cov.pred.target"]] )
      ) stop( "some required items are missing in 'object'\n",
        "re-run 'predict' with argument 'full.covmat = TRUE'"
      )

      ## back-transform trend and kriging predictions for the points in the block

      t.mu <- exp( t.object[, "trend"] + 0.5 * t.object[, "var.target"] )
      t.pred <- f.backtrf( t.object, pred.only = TRUE )

      ## covariance matrix of log(observations)

      t.cov.log.obs <- object[["mse.pred"]] - object[["var.pred"]] +
        object[["cov.pred.target"]] + t( object[["cov.pred.target"]] )


      ## covariance matrix of backtransformed prediction errors

      t.aux <- exp( t.cov.log.obs ) + exp( object[["var.pred"]] ) -
        exp( object[["cov.pred.target"]] ) - exp( t( object[["cov.pred.target"]] ) )

      t.cov.errors <- t( t.mu * t( t.mu * t.aux ) )

      ## compute weighted mean of diagonal and off-diagonal elements of
      ## covariance matrix of back-transformed point kriging predictions

      t.sum.diag <- sum( diag( t.cov.errors ), na.rm = TRUE )
      t.sum.off.diag <- sum( t.cov.errors, na.rm = TRUE ) - t.sum.diag

      if( !is.data.frame( all.pred ) ){

        t.mean <- mean( t.pred, na.rm = TRUE )

        t.mse <- (
          t.sum.diag / n.sample * n +
          t.sum.off.diag / ( n.sample * (n.sample-1L) ) * n * (n-1L)
        ) / n^2


      } else {

        t.mean <- mean( all.pred[["lgn.pred"]], na.rm = TRUE )

        t.sum.diag <- sum( all.pred[["lgn.se"]]^2, na.rm = TRUE )

        t.mse <- (
          t.sum.diag + t.sum.off.diag / ( n.sample * (n.sample-1L) ) * n * (n-1L)
        ) / n^2

      }

      t.result <- c( pred = t.mean, se = sqrt(t.mse) )

      if( extended.output ){
        attr( t.result, "mse.lgn.pred" ) <- t.cov.errors
      }

      t.result
    }
  )

  invisible( result )

}