SAEMSCL {CensSpatial} R Documentation

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 <<ordonezjosealejandro@gmail.com>>, Victor H. Lachos <<hlachos@ime.unicamp.br>> and Christian E. Galarza <<cgalarza88@gmail.com>>

Maintainer: Alejandro Ordonez <<ordonezjosealejandro@gmail.com>>

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.

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)



[Package CensSpatial version 3.6 Index]