# R/AuxiliaryFunctions_USER.R In RcppCensSpatial: Spatial Estimation and Prediction for Censored/Missing Responses

#### Documented in CovMatdist2Dmatrixpredict.sclm

#' Distance matrix computation
#'
#' It computes the Euclidean distance matrix for a set of coordinates.
#' @param coords 2D spatial coordinates of dimensions \eqn{n\times 2}.
#' @return An \eqn{n\times n} distance matrix.
#' @author Katherine L. Valeriano, Alejandro Ordoñez, Christian E. Galarza, and Larissa A. Matos.
#' @examples
#' n = 100
#' set.seed(1000)
#' x = round(runif(n,0,10), 5)     # X coordinate
#' y = round(runif(n,0,10), 5)     # Y coordinate
#' Mdist = dist2Dmatrix(cbind(x, y))

dist2Dmatrix = function(coords){

if (length(c(coords)) == 0) stop("coords must be specified")
coords = as.matrix(coords)
if (ncol(coords) != 2) stop("coords must contain 2 columns")
if (!all(c(is.finite(coords)))) stop ("coords must contain only finite values")

distancesM = crossdist(coords)
out.ST = (distancesM + t(distancesM))/2
return(out.ST)
}

#' Covariance matrix for spatial models
#'
#' It computes the spatial variance-covariance matrix considering exponential,
#' gaussian, matérn, or power exponential correlation function.
#'
#' @param phi spatial scaling parameter.
#' @param tau2 nugget effect parameter.
#' @param sig2 partial sill parameter.
#' @param coords 2D spatial coordinates of dimensions \eqn{n\times 2}.
#' @param type type of spatial correlation function: \code{'exponential'},
#' \code{'gaussian'}, \code{'matern'}, and \code{'pow.exp'} for exponential,
#' gaussian, matérn, and power exponential, respectively.
#' @param kappa parameter for some spatial correlation functions. For exponential
#' and gaussian \code{kappa=NULL}, for power exponential \code{0 < kappa <= 2},
#' and for matérn correlation function \code{kappa > 0}.
#'
#' @details The spatial covariance matrix is given by
#'
#' \eqn{\Sigma = [Cov(s_i, s_j )] = \sigma^2 R(\phi) + \tau^2 I_n},
#'
#' where \eqn{\sigma^2 > 0} is the partial sill, \eqn{\phi > 0} is the spatial scaling
#' parameter, \eqn{\tau^2 > 0} is known as the nugget effect in the geostatistical
#' framework, \eqn{R(\phi)} is the \eqn{n\times n} correlation matrix computed from a
#' correlation function, and \eqn{I_n} is the \eqn{n\times n} identity matrix.
#'
#' The spatial correlation functions available are:
#' \describe{
#'     \item{\strong{Exponential}:}{\eqn{Corr(d) = exp(-d/\phi)},}
#'
#'     \item{\strong{Gaussian}:}{\eqn{Corr(d) = exp(-(d/\phi)^2)},}
#'
#'     \item{\strong{Matérn}:}{\eqn{Corr(d) = \frac{1}{2^{(\kappa-1)}\Gamma(\kappa)}\left(\frac{d}{\phi}\right)^\kappa K_\kappa \left( \frac{d}{\phi} \right)},}
#'
#'     \item{\strong{Power exponential}:}{\eqn{Corr(d) = exp(-(d/\phi)^\kappa)},}
#'
#' where \eqn{d \geq 0} is the Euclidean distance between two observations,
#' \eqn{\Gamma(.)} is the gamma function, \eqn{\kappa} is the smoothness parameter,
#' and \eqn{K_\kappa(.)} is the modified Bessel function of the second kind of order
#' \eqn{\kappa}.}
#'
#' @return An \eqn{n\times n} spatial covariance matrix.
#'
#' @author Katherine L. Valeriano, Alejandro Ordoñez, Christian E. Galarza, and Larissa A. Matos.
#'
#'
#' @examples
#' set.seed(1000)
#' n = 20
#' coords = round(matrix(runif(2*n, 0, 10), n, 2), 5)
#' Cov = CovMat(phi=5, tau2=0.8, sig2=2, coords=coords, type="exponential")

CovMat = function(phi, tau2, sig2, coords, type="exponential", kappa=NULL){

if (length(c(phi))>1 | !is.numeric(phi)) stop("phi must be specified")
if (phi <= 0) stop("The spatial parameter (phi) must be non-negative")

if (length(c(tau2))>1 | !is.numeric(tau2)) stop("tau2 must be specified")
if (tau2 <= 0) stop("The nugget effect (tau2) must be non-negative")

if (length(c(sig2))>1 | !is.numeric(sig2)) stop("sig2 must be specified")
if (sig2 <= 0) stop("The partial sill (sig2) must be non-negative")

if (length(c(coords)) == 0) stop("coords must be specified")
coords = as.matrix(coords)
if (ncol(coords) != 2) stop("coords must contain 2 columns")
if (!all(c(is.finite(coords)))) stop ("coords must contain only finite values")

if (is.null(type)) stop("type must be specified")
if (type!="matern" & type !="gaussian" & type != "pow.exp" & type != "exponential"){
stop("type should be one of matern, gaussian, pow.exp, exponential")}

if (type!="exponential" & type!="gaussian"){
if (is.null(kappa)) stop("kappa mut be provided for pow.exp and matern")
if (length(c(kappa))>1 | !is.numeric(kappa)) stop("kappa must be specified")
if (type=="pow.exp" & (kappa > 2| kappa<=0)) stop("kappa must be a real in (0,2]")
if (type=="matern" & kappa <= 0) stop("kappa must be a real number in (0,Inf)")
if (type=="matern" & is.infinite(kappa)) stop("kappa must be a real number in (0,Inf)")
} else { kappa = 0 }

dist = dist2Dmatrix(coords)
Covariance = varianceMat(phi, tau2, sig2, kappa, dist, type)$Sigma out.ST = (Covariance + t(Covariance))/2 return(out.ST) } #' Prediction in spatial models with censored/missing responses #' #' It performs spatial prediction in a set of new \code{S} spatial locations. #' #' @param object object of class \code{'sclm'} given as output of \code{\link{EM.sclm}}, #' \code{\link{MCEM.sclm}}, or \code{\link{SAEM.sclm}} function. #' @param locPre matrix of coordinates for which prediction is performed. #' @param xPre matrix of covariates for which prediction is performed. #' @param ... further arguments passed to or from other methods. #' #' @details This function predicts using the mean squared error (MSE) criterion, which #' takes the conditional expectation E(Y|X) as the best linear predictor. #' #' @return The function returns a list with: #' \item{coord}{matrix of coordinates.} #' \item{predValues}{predicted values.} #' \item{sdPred}{predicted standard deviations.} #' #' @author Katherine L. Valeriano, Alejandro Ordoñez, Christian E. Galarza, and Larissa A. Matos. #' #' @seealso \code{\link{EM.sclm}}, \code{\link{MCEM.sclm}}, \code{\link{SAEM.sclm}} #' #' @examples \donttest{set.seed(1000) #' n = 120 #' coords = round(matrix(runif(2*n,0,15),n,2), 5) #' x = cbind(rbinom(n,1,0.50), rnorm(n), rnorm(n)) #' data = rCensSp(c(1,4,-1), 2, 3, 0.50, x, coords, "left", 0.10, 20) #' #' ## Estimation #' data1 = data$Data
#'
#' # Estimation: EM algorithm
#' fit1 = EM.sclm(y=data1$y, x=data1$x, ci=data1$ci, lcl=data1$lcl,
#'                ucl=data1$ucl, coords=data1$coords, phi0=2.50, nugget0=1)
#'
#' # Estimation: SAEM algorithm
#' fit2 = SAEM.sclm(y=data1$y, x=data1$x, ci=data1$ci, lcl=data1$lcl,
#'                  ucl=data1$ucl, coords=data1$coords, phi0=2.50, nugget0=1)
#'
#' # Estimation: MCEM algorithm
#' fit3 = MCEM.sclm(y=data1$y, x=data1$x, ci=data1$ci, lcl=data1$lcl,
#'                  ucl=data1$ucl, coords=data1$coords, phi0=2.50, nugget0=1,
#'                  MaxIter=300)
#' cbind(fit1$theta, fit2$theta, fit3$theta) #' #' # Prediction #' data2 = data$TestData
#' pred1 = predict(fit1, data2$coords, data2$x)
#' pred2 = predict(fit2, data2$coords, data2$x)
#' pred3 = predict(fit3, data2$coords, data2$x)
#'
#' # Cross-validation
#' mean((data2$y - pred1$predValues)^2)
#' mean((data2$y - pred2$predValues)^2)
#' mean((data2$y - pred3$predValues)^2)}
#' @export
predict.sclm = function(object, locPre, xPre, ...){

if (is.null(object)) stop("object must be specified")

if (!is.null(locPre) & !is.null(xPre)){
locPre = as.matrix(locPre)
xPre = as.matrix(xPre)
if (!all(c(is.finite(locPre)))) stop("locPre must contain only finite values")
if (!all(c(is.finite(xPre)))) stop("xPre must contain only finite values")
if(ncol(xPre)!=length(c(object$beta))) stop("Non-conformable dimensions between xPred and beta") if (nrow(locPre)!=nrow(xPre) | ncol(locPre)!=2) stop("Non-conformable dimensions between locPre and xPre") } else { stop("locPre and xPre must be specified") } ypred = predict.new(object, xPre, locPre) out.ST = ypred return(out.ST) } #' @export summary.sclm = function(object, ...){ cat('----------------------------------------------------------------\n') cat(' Censored Linear Spatial Regression Model \n') cat('----------------------------------------------------------------\n') cat("Call:\n") print(object$call)
cat('\nEstimated parameters:\n')
print(object$tab) cat('\n') cat(paste('The effective range is',round(object$range,4),'units.\n'))
cat('\nModel selection criteria:\n')
print(object$critFin) cat('\nDetails:\n') cat('Number of censored/missing values:', object$ncens, '\n')
cat('Convergence reached?:', (object$Iter < object$MaxIter), '\n')
cat('Iterations:', object$Iter,'/',object$MaxIter, '\n')
cat('Processing time:', round(object$time,4), units(object$time), '\n')
}

#' @export
print.sclm = function(x, ...){
cat('----------------------------------------------------------------\n')
cat('           Spatial Censored Linear Regression Model   \n')
cat('----------------------------------------------------------------\n')
cat("Call:\n")
print(x$call) cat('\nEstimated parameters:\n') print(x$tab)
cat('\n')
cat(paste('The effective range is',round(x$range,4),'units.\n')) cat('\nModel selection criteria:\n') print(x$critFin)
cat('\nDetails:\n')
cat('Number of censored/missing values:', x$ncens,'\n') cat('Convergence reached?:',(x$Iter < x$MaxIter),'\n') cat('Iterations:',x$Iter,'/',x$MaxIter,'\n') cat('Processing time:',round(x$time,4),units(x\$time),'\n')
}

#' @export
plot.sclm = function(x, ...){
plot.convergence(x)
}


## Try the RcppCensSpatial package in your browser

Any scripts or data that you put into this service are public.

RcppCensSpatial documentation built on June 28, 2022, 1:07 a.m.