SAEM.sclm: Censored Spatial Model Estimation via SAEM Algorithm

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/EstSAEMspatial_USER.R

Description

This function returns the maximum likelihood (ML) estimates of the unknown parameters in Gaussian spatial models with censored/missing responses via the SAEM algorithm. It supports left, right, interval, or missing values in the dependent variable. It also computes the observed information matrix using the method developed by \insertCitelouis1982finding;textualRcppCensSpatial.

Usage

1
2
3
4
SAEM.sclm(y, x, cens, LI, LS, coords, init.phi, init.nugget,
  type = "exponential", kappa = 0, lower = c(0.01, 0.01), upper = c(30,
  30), MaxIter = 300, M = 20, pc = 0.25, error = 1e-05,
  show.SE = TRUE)

Arguments

y

vector of responses.

x

design matrix.

cens

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

LI

lower limit of detection. For each observation: if non-censored =y, if left-censored/missing =-Inf, or =LOD if right/interval-censored.

LS

upper limit of detection. For each observation: if non-censored =y, if right-censored/missing =Inf, or =LOD if left/interval-censored.

coords

2D spatial coordinates.

init.phi

initial value for the spatial scaling parameter.

init.nugget

initial value for the nugget effect parameter.

type

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

kappa

parameter for all 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 iterations of the SAEM algorithm with no-memory. By default =0.25.

error

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

show.SE

TRUE or FALSE. It indicates if the standard errors should be estimated. By default =TRUE.

Details

The spatial Gaussian model is given by

Y = Xβ + ξ,

where Y is the n x 1 vector of response, X is the n x q design matrix, β is the q x 1 vector of regression coefficients to be estimated, and ξ is the error term which is normally distributed with zero-mean and covariance matrix Σ=σ^2 R(φ) + τ^2 I_n. We assume that Σ is non-singular and X has full rank \insertCitediggle2007springerRcppCensSpatial.

The estimation process was performed via the SAEM \insertCitedelyon1999convergenceRcppCensSpatial algorithm. The spatial SAEM algorithm was previously proposed by \insertCitelachos2017influence;textualRcppCensSpatial and \insertCiteordonez2018geostatistical;textualRcppCensSpatial and is available in package CensSpatial. The difference between this package to CensSpatial is that the random observations are sampled through the slice sampling algorithm available in package relliptical and the optimization procedure by the roptim package.

This model is also a particular case of the Spatio-temporal model defined by \insertCitevaleriano2021likelihood;textualRcppCensSpatial, 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

The function returns an object of class sclm which is a list given by:

Theta

estimated parameters in all iterations, θ = (β, σ^2, φ, τ^2).

theta

final estimation of θ = (β, σ^2, φ, τ^2).

beta

estimated β.

sigma2

estimated σ^2.

phi

estimated φ.

tau2

estimated τ^2.

EY

stochastic approximation of the first moment for the truncated normal distribution.

EYY

stochastic approximation of the second moment for the truncated normal distribution.

SE

vector of standard errors of θ = (β, σ^2, φ, τ^2).

InfMat

observed information matrix.

loglik

log-likelihood for the SAEM method.

AIC

Akaike information criterion.

BIC

Bayesian information criterion.

Iterations

number of iterations needed to converge.

ptime

processing time.

range

the effective range.

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 cens as a vector of zeros.

Functions print, summary, and plot work for objects of class sclm.

Author(s)

Katherine L. Valeriano, Alejandro Ordonez, Christian E. Galarza and Larissa A. Matos.

References

\insertAllCited

See Also

EM.sclm, MCEM.sclm, predict.sclm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Simulated example: 10% of right-censored observations
n = 50   # Test with another values for n
set.seed(1000)
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,-2),2,3,0.50,x,coords,"right",0.10,0,"matern",2)

fit = SAEM.sclm(y=data$yobs, x=data[,7:9], cens=data$cens, LI=data$LI,
             LS=data$LS, coords=data[,5:6], init.phi=2, init.nugget=1,
             type="matern", kappa=2, MaxIter = 20, error=1e-4)
summary(fit)


# Simulated example: censored and missing observations
n = 200
set.seed(123)
coords = round(matrix(runif(2*n,0,20),n,2),5)
x = cbind(1, rnorm(n), rexp(n))
data = rCensSp(c(1,4,-1),2,4,0.50,x,coords,"left",0.10,0,"exponential",0)
data$yobs[c(10,20)] = NA;   data$cens[c(10,20)] = 1
data$LI[c(10,20)] = -Inf;   data$LS[c(10,20)] = Inf

fit2 = SAEM.sclm(y=data$yobs, x=data[,7:9], cens=data$cens, LI=data$LI,
              LS=data$LS, coords=data[,5:6], init.phi=2, init.nugget=1,
              type="exponential", MaxIter = 300, error=1e-4)
fit2$theta  # Estimates
fit2$SE     # Standard error
fit2$InfMat # Information matrix
plot(fit2)

RcppCensSpatial documentation built on Sept. 21, 2021, 5:07 p.m.