R: ML estimation in spatio-temporal model with censored/missing...

EstStempCens {StempCens}

R Documentation

ML estimation in spatio-temporal model with censored/missing responses

Description

Return the maximum likelihood estimates of the unknown parameters of spatio-temporal model with censored/missing responses. The estimates are obtained using SAEM algorithm. The function also computes the observed information matrix using the method developed by Louis (1982). The types of censoring considered are left, right, interval or missing values.

Usage

EstStempCens(
  y,
  x,
  cc,
  time,
  coord,
  LI,
  LS,
  init.phi,
  init.rho,
  init.tau2,
  tau2.fixo = FALSE,
  type.Data = "balanced",
  method = "nlminb",
  kappa = 0,
  type.S = "exponential",
  IMatrix = TRUE,
  lower.lim = c(0.01, -0.99, 0.01),
  upper.lim = c(30, 0.99, 20),
  M = 20,
  perc = 0.25,
  MaxIter = 300,
  pc = 0.2,
  error = 1e-06
)

Arguments

`y`	a vector of responses.
`x`	a matrix or vector of covariates.
`cc`	a vector of censoring indicators. For each observation: `1` if censored/missing and `0` if non-censored/non-missing.
`time`	a vector of time.
`coord`	a matrix of coordinates of the spatial locations.
`LI`	lower limit of detection. For each observation: if non-censored/non-missing `=y`, if left-censored/missing `=-Inf` or `=LOD` if right/interval-censored.
`LS`	upper limit of detection. For each observation: if non-censored/non-missing `=y`, if right-censored/missing `=Inf` or `=LOD` if left/interval-censored.
`init.phi`	initial value of the spatial scaling parameter.
`init.rho`	initial value of the time scaling parameter.
`init.tau2`	initial value of the the nugget effect parameter.
`tau2.fixo`	`TRUE` or `FALSE`. Indicate if the nugget effect (`\tau^2`) parameter must be fixed. By default = `FALSE`.
`type.Data`	type of the data: '`balanced`' for balanced data and '`unbalanced`' for unbalanced data. By default = `balanced`.
`method`	optimization method used to estimate (`\phi`, `\rho` and `\tau^2`): '`optim`' for the function `optim` and '`nlminb`' for the function `nlminb`. By default = `nlminb`.
`kappa`	parameter for all spatial covariance functions. In the case of exponential, gaussian and spherical function `\kappa` is equal to zero. For the power exponential function `\kappa` is a number between 0 and 2. For the matern correlation function is upper than 0.
`type.S`	type of spatial correlation function: '`exponential`' for exponential, '`gaussian`' for gaussian, '`matern`' for matern, '`pow.exp`' for power exponential and '`spherical`' for spherical function, respectively. Default is `exponential` function.
`IMatrix`	`TRUE` or `FALSE`. Indicate if the observed information matrix will be computed. By default = `TRUE`.
`lower.lim`, `upper.lim`	vectors of lower and upper bounds for the optimization method. If unspecified, the default is `c(0.01,-0.99,0.01)` for the lower bound and `c(30,0.99,20)` for the upper bound if tau2.fixo=`FALSE`.
`M`	number of Monte Carlo samples for stochastic approximation. By default = `20`.
`perc`	percentage of burn-in on the Monte Carlo sample. By default = `0.25`.
`MaxIter`	the maximum number of iterations of the SAEM algorithm. By default = `300`.
`pc`	percentage of iterations of the SAEM algorithm with no memory. By default = `0.20`.
`error`	the convergence maximum error. By default = `1e-6`.

Details

The spatio-temporal Gaussian model is giving by:

Y(s_i,t_j)= \mu(s_i,t_j)+ Z(s_i,t_j) + \epsilon(s_i,t_j),

where the deterministic term \mu(s_i,t_j) and the stochastic terms Z(s_i,t_j), \epsilon(s_i,t_j) can depend on the observed spatio-temporal indexes for Y(s_i,t_j). We assume Z is normally distributed with zero-mean and covariance matrix \Sigma_z = \sigma^2 \Omega_{\phi\rho}, where \sigma^2 is the partial sill, \Omega_{\phi\rho} is the spatio-temporal correlation matrix,\phi and \rho are the spatial and time scaling parameters; \epsilon(s_i,t_j) is an independent and identically distributed measurement error with E[\epsilon(s_i,t_j)]=0, variance Var[\epsilon(s_i,t_j)]=\tau^2 (the nugget effect) and Cov[\epsilon(s_i,t_j), \epsilon(s_k,t_l)]=0 for all s_i =! s_k or t_j =! t_l.

In particular, we define \mu(s_i,t_j), the mean of the stochastic process as

\mu(s_i,t_j)=\sum_{k=1}^{p} x_k(s_i,t_j)\beta_k,

where x_1(s_i,t_j),..., x_p(s_i,t_j) are known functions of (s_i,t_j), and \beta_1,...,\beta_p are unknown parameters to be estimated. Equivalently, in matrix notation, we have the spatio-temporal linear model as follows:

Y = X \beta + Z + \epsilon,

Z ~ N(0,\sigma^2 \Omega_{\phi\rho}),

\epsilon ~ N(0,\tau^2 I_m).

Therefore the spatio-temporal process, Y, has normal distribution with mean E[Y]=X\beta and variance \Sigma=\sigma^2\Omega_{\phi\rho}+\tau^2 I_m. We assume that \Sigma is non-singular and X has full rank.

The estimation process was computed via SAEM algorithm initially proposed by Delyon et al. (1999).

Value

The function returns an object of class Est.StempCens which is a list given by:

m.data: Returns a list with all data components given in input.
m.results: A list given by:

`theta`	final estimation of `\theta = (\beta, \sigma^2, \tau^2, \phi, \rho)`.
`Theta`	estimated parameters in all iterations, `\theta = (\beta, \sigma^2, \tau^2, \phi, \rho)`.
`beta`	estimated `\beta`.
`sigma2`	estimated `\sigma^2`.
`tau2`	estimated `\tau^2`.
`phi`	estimated `\phi`.
`rho`	estimated `\rho`.
`Eff.range`	estimated effective range.
`PsiInv`	estimated `\Psi^-1`, where `\Psi=\Sigma/\sigma^2`.
`Cov`	estimated `\Sigma`.
`SAEMy`	stochastic approximation of the first moment for the truncated normal distribution.
`SAEMyy`	stochastic approximation of the second moment for the truncated normal distribution.
`Hessian`	Hessian matrix, the negative of the conditional expected second derivative matrix given the observed values.
`Louis`	the observed information matrix using the Louis' method.
`loglik`	log likelihood for SAEM method.
`AIC`	Akaike information criteria.
`BIC`	Bayesian information criteria.
`AICcorr`	corrected AIC by the number of parameters.
`iteration`	number of iterations needed to convergence.

Author(s)

Katherine L. Valeriano, Victor H. Lachos and Larissa A. Matos

Examples

## Not run: 
# Simulating data
# Initial parameter values
beta <- c(-1,1.50)
phi <- 5;     rho <- 0.45
tau2 <- 0.80; sigma2 <- 1.5
n1 <- 5    # Number of spatial locations
n2 <- 5    # Number of temporal index
set.seed(1000)
x.coord <- round(runif(n1,0,10),9)    # X coordinate
y.coord <- round(runif(n1,0,10),9)    # Y coordinate
coord  <- cbind(x.coord,y.coord)      # Cartesian coordinates without repetitions
coord2 <- cbind(rep(x.coord,each=n2),rep(y.coord,each=n2)) # Cartesian coordinates with repetitions
time <- as.matrix(seq(1,n2))          # Time index without repetitions
time2 <- as.matrix(rep(time,n1))      # Time index with repetitions
x1 <- rexp(n1*n2,2)
x2 <- rnorm(n1*n2,2,1)
x  <- cbind(x1,x2)
media <- x%*%beta
# Covariance matrix
Ms  <- as.matrix(dist(coord))   # Spatial distances
Mt  <- as.matrix(dist(time))    # Temporal distances
Cov <- CovarianceM(phi,rho,tau2,sigma2,Ms,Mt,1.5,"matern")
# Data
require(mvtnorm)
y <- as.vector(rmvnorm(1,mean=as.vector(media),sigma=Cov))
perc <- 0.20
aa <- sort(y); bb <- aa[1:(perc*n1*n2)]; cutof <- bb[perc*n1*n2]
cc <- matrix(1,(n1*n2),1)*(y<=cutof)
y[cc==1] <- cutof
LI <- y; LI[cc==1] <- -Inf    # Left-censored
LS <- y

# Estimation
est_teste <- EstStempCens(y, x, cc, time2, coord2, LI, LS, init.phi=3.5,
                 init.rho=0.5, init.tau2=0.7,tau2.fixo=FALSE, kappa=1.5,
                 type.S="matern", IMatrix=TRUE, M=20, perc=0.25,
                 MaxIter=300, pc=0.2)
## End(Not run)

[Package StempCens version 1.1.0 Index]