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.
1 2 3 4 |
y |
vector of responses. |
x |
design matrix. |
cens |
vector of censoring indicators. For each observation: |
LI |
lower limit of detection. For each observation: if non-censored |
LS |
upper limit of detection. For each observation: if non-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: ' |
kappa |
parameter for all spatial correlation functions. See |
lower, upper |
vectors of lower and upper bounds for the optimization method. If unspecified, the default is
|
MaxIter |
maximum number of iterations of the MCEM algorithm. By default |
nMin |
initial sample size for Monte Carlo integration. By default |
nMax |
maximum sample size for Monte Carlo integration. By default |
error |
maximum convergence error. By default |
show.SE |
|
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.
EM.sclm
, SAEM.sclm
, predict.sclm
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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.