SAEM.sclm {RcppCensSpatial}

ML estimation of spatial censored linear models via the SAEM algorithm

Description

It fits the left, right, or interval spatial censored linear model using the Stochastic Approximation EM (SAEM) algorithm. It provides estimates and standard errors of the parameters and supports missing values on the dependent variable.

Usage

SAEM.sclm(y, x, ci, lcl = NULL, ucl = NULL, coords, phi0, nugget0,
  type = "exponential", kappa = NULL, lower = c(0.01, 0.01),
  upper = c(30, 30), MaxIter = 300, M = 20, pc = 0.2, error = 1e-04,
  show_se = TRUE)

Arguments

y

vector of responses of length n.

x

design matrix of dimensions n\times q, where q is the number of fixed effects, including the intercept.

ci

vector of censoring indicators of length n. For each observation: 1 if censored/missing, 0 otherwise.

lcl, ucl

vectors of length n representing the lower and upper bounds of the interval, which contains the true value of the censored observation. Default =NULL, indicating no-censored data. For each observation: lcl=-Inf and ucl=c (left censoring); lcl=c and ucl=Inf (right censoring); and lcl and ucl must be finite for interval censoring. Moreover, missing data could be defined by setting lcl=-Inf and ucl=Inf.

coords

2D spatial coordinates of dimensions n\times 2.

phi0

initial value for the spatial scaling parameter.

nugget0

initial value for the nugget effect parameter.

type

type of spatial correlation function: 'exponential', 'gaussian', 'matern', and 'pow.exp' for exponential, gaussian, matérn, and power exponential, respectively.

kappa

parameter for some spatial correlation functions. See CovMat.

lower, upper

vectors of lower and upper bounds for the optimization method. If unspecified, the default is c(0.01,0.01) for lower and c(30,30) for upper.

MaxIter

maximum number of iterations of the SAEM algorithm. By default =300.

M

number of Monte Carlo samples for stochastic approximation. By default =20.

pc

percentage of initial iterations of the SAEM algorithm with no memory. It is recommended that 50<MaxIter*pc<100. By default =0.20.

error

maximum convergence error. By default =1e-4.

show_se

logical. It indicates if the standard errors should be estimated by default =TRUE.

Details

The spatial Gaussian model is given by

Y = X\beta + \xi,

where Y is the n\times 1 response vector, X is the n\times q design matrix, \beta is the q\times 1 vector of regression coefficients to be estimated, and \xi is the error term which is normally distributed with zero-mean and covariance matrix \Sigma=\sigma^2 R(\phi) + \tau^2 I_n. We assume that \Sigma is non-singular and X has full rank (Diggle and Ribeiro 2007).

The estimation process is performed via the SAEM algorithm, initially proposed by Delyon et al. (1999). The spatial censored (SAEM) algorithm was previously proposed by Lachos et al. (2017) and Ordoñez et al. (2018) and is available in the package CensSpatial. These packages differ in the random number generation and optimization procedure.

This model is also a particular case of the spatio-temporal model defined by Valeriano et al. (2021) when the number of temporal observations is equal to one. The computing codes of the spatio-temporal SAEM algorithm are available in the package StempCens.

Value

An object of class "sclm". Generic functions print and summary have methods to show the results of the fit. The function plot can extract convergence graphs for the parameter estimates.

Specifically, the following components are returned:

Theta

estimated parameters in all iterations, \theta = (\beta, \sigma^2, \phi, \tau^2).

theta

final estimation of \theta = (\beta, \sigma^2, \phi, \tau^2).

beta

estimated \beta.

sigma2

estimated \sigma^2.

phi

estimated \phi.

tau2

estimated \tau^2.

EY

stochastic approximation of the first conditional moment.

EYY

stochastic approximation of the second conditional moment.

SE

vector of standard errors of \theta = (\beta, \sigma^2, \phi, \tau^2).

InfMat

observed information matrix.

loglik

log-likelihood for the SAEM method.

AIC

Akaike information criterion.

BIC

Bayesian information criterion.

Iter

number of iterations needed to converge.

time

processing time.

call

RcppCensSpatial call that produced the object.

tab

table of estimates.

critFin

selection criteria.

range

effective range.

ncens

number of censored/missing observations.

MaxIter

maximum number of iterations for the SAEM algorithm.

Note

The SAEM final estimates correspond to the estimates obtained at the last iteration of the algorithm.

To fit a regression model for non-censored data, just set ci as a vector of zeros.

Author(s)

Katherine L. Valeriano, Alejandro Ordoñez, Christian E. Galarza, and Larissa A. Matos.

References

Delyon B, Lavielle M, Moulines E (1999). “Convergence of a stochastic approximation version of the EM algorithm.” The Annals of Statistics, 27(1), 94–128.

Diggle P, Ribeiro P (2007). Model-based Geostatistics. Springer.

Lachos VH, Matos LA, Barbosa TS, Garay AM, Dey DK (2017). “Influence diagnostics in spatial models with censored response.” Environmetrics, 28(7).

Ordoñez JA, Bandyopadhyay D, Lachos VH, Cabral CRB (2018). “Geostatistical estimation and prediction for censored responses.” Spatial Statistics, 23, 109–123. doi:10.1016/j.spasta.2017.12.001.

Valeriano KL, Lachos VH, Prates MO, Matos LA (2021). “Likelihood-based inference for spatiotemporal data with censored and missing responses.” Environmetrics, 32(3).

See Also

EM.sclm, MCEM.sclm, predict.sclm

Examples

# Example 1: 8% of right-censored observations
set.seed(1000)
n = 50   # Test with another values for n
coords = round(matrix(runif(2*n,0,15),n,2), 5)
x = cbind(rnorm(n), rnorm(n))
data = rCensSp(c(4,-2), 1, 3, 0.50, x, coords, "right", 0.08)

fit = SAEM.sclm(y=data$y, x=x, ci=data$ci, lcl=data$lcl, ucl=data$ucl,
                coords, phi0=2, nugget0=1, type="exponential", M=10,
                pc=0.18)
fit

# Example 2: censored and missing observations
set.seed(123)
n = 200
coords = round(matrix(runif(2*n,0,20),n,2), 5)
x = cbind(runif(n), rnorm(n), rexp(n))
data = rCensSp(c(1,4,-1), 2, 3, 0.50, x, coords, "left", 0.05, 0,
               "matern", 3)
data$y[c(10,120)] = NA
data$ci[c(10,120)] = 1
data$ucl[c(10,120)] = Inf

fit2 = SAEM.sclm(y=data$y, x=x, ci=data$ci, lcl=data$lcl, ucl=data$ucl,
                 coords, phi0=2, nugget0=1, type="matern", kappa=3,
                 M=10, pc=0.18)
fit2$tab
plot(fit2)