MCEM.sclm: Censored Spatial Model Estimation via MCEM Algorithm

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

View source: R/EstMCEMspatial_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 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.

Usage

1
2
3
4
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)

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 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.

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 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.

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

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.

Note

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.

Author(s)

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

References

\insertAllCited

See Also

EM.sclm, SAEM.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: 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)

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