MCEM.sclm: Censored Spatial Model Estimation via MCEM Algorithm
In RcppCensSpatial: Spatial Estimation and Prediction for Censored/Missing Responses

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

View source: R/EstMCEMspatial_USER.R

This function returns the maximum likelihood (ML) estimates of the unknown parameters in Gaussian spatial models with censored/missing responses via the MCEM 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.

MCEM.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 = 500, nMin = 20, nMax = 5000, error = 1e-05,
  show.SE = TRUE)

`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 MCEM algorithm. By default `=500`.
`nMin`	initial sample size for Monte Carlo integration. By default `=20`.
`nMax`	maximum sample size for Monte Carlo integration. By default `=5000`.
`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`.

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 MCEM algorithm initially proposed by \insertCitewei1990monte;textualRcppCensSpatial. The Monte Carlo integration starts with a sample of size nMin; at each iteration, the sample size increases (nMax-nMin)/MaxIter, and at the last iteration, the sample size is nMax. The random observations are sampled through the slice sampling algorithm available in package relliptical.

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`	MC approximation of the first moment for the truncated normal distribution.
`EYY`	MC 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 MCEM method.
`AIC`	Akaike information criterion.
`BIC`	Bayesian information criterion.
`Iterations`	number of iterations needed to converge.
`ptime`	processing time.
`range`	the effective range.

The MCEM final estimates correspond to the mean of the estimates obtained at each iteration after deleting the half and applying a thinning of 3.

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.

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

\insertAllCited

EM.sclm, SAEM.sclm, predict.sclm

# Simulated example: censored and missing data
n = 50   # Test with another values for n
set.seed(1000)
coords = round(matrix(runif(2*n,0,15),n,2),5)
x = cbind(rnorm(n), rnorm(n))
data = rCensSp(c(2,-1),2,3,0.70,x,coords,"left",0.08,0,"matern",1)
data$yobs[20] = NA
data$cens[20] = 1; data$LI[20] = -Inf; data$LS[20] = Inf

fit = MCEM.sclm(y=data$yobs, x=data[,7:8], cens=data$cens, LI=data$LI,
             LS=data$LS, coords=data[,5:6], init.phi=2.50, init.nugget=0.75,
             type="matern", kappa=1, MaxIter=20, nMax=1000, error=1e-4)
print(fit)


# Application: TCDD concentration in Missouri
library(CensSpatial)
data("Missouri")
y = log(Missouri$V3)
cc = Missouri$V5
coord = cbind(Missouri$V1/100,Missouri$V2)
X = matrix(1,length(y),1)
LI = LS = y; LI[cc==1] = -Inf

fit2 = MCEM.sclm(y=y, x=X, cens=cc, LI=LI, LS=LS, coords=coord, init.phi=5,
              init.nugget=1, type="exponential", lower=c(1e-5,1e-5), upper=c(50,50),
              MaxIter=500, nMax=1000, error=1e-5)
summary(fit2)
plot(fit2)