StempCens (version 1.1.0)

EstStempCens: 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 louis1982finding;textualStempCens. 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.

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.

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 delyon1999convergence;textualStempCens.

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)
# }
# NOT RUN {
# }