SAEMSCL {CensSpatial}

SAEM Algorithm estimation for censored spatial data.

Description

It estimates the parameters for a linear spatial model with censored observations

Usage

SAEMSCL(cc, y, cens.type="left", trend = "cte", LI = NULL, LS = NULL, x = NULL, coords,
kappa = 0, M = 20, perc = 0.25, MaxIter = 300, pc = 0.2, cov.model = "exponential",
fix.nugget = TRUE, nugget, inits.sigmae, inits.phi, search = FALSE, lower, upper)

Arguments

cc

(binary vector) indicator of censure (1: censored observation 0: observed).

y

(vector) corresponds to response variable.

cens.type

type of censure ("left":left or "right":right).

trend

linear trends options: "cte", "1st", "2nd" and "other", the three first are defined like in geoR, if trend="other", x (design matrix) must be defined.

LI

(vector) lower limit, if cens.type="both", LI must be provided, if cens.type="left" or "right" LI and LS are defined by the function through the indicator of censure cc.

LS

(vector) upper limit, if cens.type="both", LS must be provided, if cens.type="left" or "right" LI and LS are defined by the function through the indicator of censure cc.

x

design matrix.

coords

corresponds to the coordinates of the spatial data (2D coordinates).

kappa

value of kappa used in some covariance functions.

M

number of montecarlo samples for stochastic aproximation.

perc

percentage of burn-in on the Monte Carlo sample. Default=0.25.

MaxIter

maximum of iterations for the algorithm.

pc

percentage of initial iterations of the SAEM algorithm. (Default=0.2).

cov.model

covariance Structure (see, cov.spatial from geoR).

fix.nugget

(logical) indicates if the \tau^2 parameter must be fixed.

nugget

if fix.nugget=TRUE, the algorithm just estimates \beta, \sigma^2, and \phi, and fixed \tau^2 like nugget, else, \tau^2 is estimated and nugget corresponds to initial value for \tau^2.

inits.sigmae

corresponds to initial value for \sigma^2.

inits.phi

corresponds to initial value for \phi parameter.

search

(logical) this argument gives bounds where the optim routine can find the solution that maximizes the Maximum likelihood expectation. If search=F, the optim routine will try to search the solutions for maximization in all the domain for \phi and \tau^2 (if fix.nugget=FALSE). If search=TRUE, the optim routine search the solutions in a specific neighborhood. We recommended to use search=F (see details).

lower

(vector or numeric) lower bound from the optim solution. If fix.nugget=T, lower is numerical and corresponds to the lower bound for search the solution of the \phi parameter, if fix.nugget=FALSE lower is a vector and corresponds to the lower bounds for search the solution of \phi and \tau2 that maximizes the Maximum Likelihood Expectation (see details).

upper

(vector or numeric) upper bound from the optim solution. If fix.nugget==T, lower is numerical and corresponds to the lower bound for searching the solution of the phi parameter, if fix.nugget==F, lower is a vector and corresponds to the lower bounds for searching the solution for \phi and \tau^2 parameters that maximizes the Maximum Likelihood Expectation

Details

The estimation process was computed via SAEM algorithm initially proposed by Deylon et. al.(1999). This is a stochastic approximation of the EM procedure. This procedure circunvent the heavy computational time involved in the MCEM procedure necessary for estimating phi and tau2 parameters (when tau2 is not fixed) since there is not an analytical solution. The search interval was proposed because sometimes the maximization procedure used by optim function does not work for large intervals.

Value

beta

estimated \beta.

sigma2

estimated \sigma^2.

phi

estimated \phi.

nugget

estimated or fixed \tau^2.

Theta

estimated parameters in all iterations (\beta, \sigma^2, \phi) or (\beta, \sigma^2, \phi, \tau^2) if fix.nugget=F.

loglik

log likelihood for SAEM method.

AIC

Akaike information criteria.

BIC

Bayesian information criteria.

AICcorr

corrected AIC by the number of parameters.

X

design matrix.

Psi

estimated covariance matrix.

theta

final estimation of \theta=(\beta, \sigma^2, \phi) or \theta=(\beta, \sigma^2, \phi, \tau^2) if fix.nugget=F.

uy

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

uyy

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

cc

indicator of censure (0:observed, 1: censored).

type

covariance structure considered in the model.

kappa

\kappa parameter for some covariance structures.

coords

coordinates of the observed data.

iterations

number of iterations needed to convergence.

fitted

fitted values for the SAEM algortihm.

Author(s)

Alejandro Ordonez <<[e-mail]>>, Victor H. Lachos <<[e-mail]>> and Christian E. Galarza <<[e-mail]>>

Maintainer: Alejandro Ordonez <<[e-mail]>>

References

DELYON, B., LAVIELLE, M.,ANDMOULI NES, E. (1999). Convergence ofa stochastic approximation version of the EM algorithm.Annals of Statistic-s27, 1, 94-128.

Diggle, P. & Ribeiro, P. (2007). Model-Based Geostatistics. Springer Series in Statistics.

See Also

localinfmeas, derivQfun

Examples

n<-200 ### sample size for estimation.
n1=100 ### number of observation used in the prediction.

###simulated coordinates
r1=sample(seq(1,30,length=400),n+n1)
r2=sample(seq(1,30,length=400),n+n1)
coords=cbind(r1,r2)

coords1=coords[1:n,]

type="matern"
#xtot<-cbind(1,runif((n+n1)),runif((n+n1),2,3))
xtot=as.matrix(rep(1,(n+n1)))
xobs=xtot[1:n,]
beta=5
#beta=c(5,3,1)

###simulated data
obj=rspacens(cov.pars=c(3,.3,0),beta=beta,x=xtot,coords=coords,kappa=1.2,
cens=0.25,n=(n+n1),n1=n1,cov.model=type,cens.type="left")

data2=obj$datare
cc=obj$cc
y=obj$datare[,3]
coords=obj$datare[,1:2]
##initials values obtained from variofit.
cov.ini=c(0.13,0.86)

est=SAEMSCL(cc,y,cens.type="left",trend="cte",coords=coords,
kappa=1.2,M=15,perc=0.25,MaxIter=10,pc=0.2,cov.model=type,
fix.nugget=TRUE,nugget=0,inits.sigmae=cov.ini[1],
inits.phi=cov.ini[2],search=TRUE,lower=0.00001,upper=100)

summary(est)