Description Usage Arguments Details Value Author(s) Examples
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 \insertCitelouis1982finding;textualStempCens. The types of censoring considered are left, right, interval or missing values.
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 | 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
)
|
y |
a vector of responses. |
x |
a matrix or vector of covariates. |
cc |
a vector of censoring indicators. For each observation: |
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 |
LS |
upper limit of detection. For each observation: if non-censored/non-missing |
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 |
|
type.Data |
type of the data: ' |
method |
optimization method used to estimate (φ, ρ and τ^2):
' |
kappa |
parameter for all spatial covariance functions. In the case of exponential, gaussian and spherical function κ is equal to zero. For the power exponential function κ is a number between 0 and 2. For the matern correlation function is upper than 0. |
type.S |
type of spatial correlation function: ' |
IMatrix |
|
lower.lim, upper.lim |
vectors of lower and upper bounds for the optimization method.
If unspecified, the default is |
M |
number of Monte Carlo samples for stochastic approximation. By default = |
perc |
percentage of burn-in on the Monte Carlo sample. By default = |
MaxIter |
the maximum number of iterations of the SAEM algorithm. By default = |
pc |
percentage of iterations of the SAEM algorithm with no memory. By default = |
error |
the convergence maximum error. By default = |
The spatio-temporal Gaussian model is giving by:
Y(s_i,t_j)= μ(s_i,t_j)+ Z(s_i,t_j) + ε(s_i,t_j),
where the deterministic term μ(s_i,t_j) and the stochastic terms Z(s_i,t_j), ε(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 Σ_z = σ^2 Ω_{φρ}, where σ^2 is the partial sill, Ω_{φρ} is the spatio-temporal correlation matrix,φ and ρ are the spatial and time scaling parameters; ε(s_i,t_j) is an independent and identically distributed measurement error with E[ε(s_i,t_j)]=0, variance Var[ε(s_i,t_j)]=τ^2 (the nugget effect) and Cov[ε(s_i,t_j), ε(s_k,t_l)]=0 for all s_i =! s_k or t_j =! t_l.
In particular, we define μ(s_i,t_j), the mean of the stochastic process as
μ(s_i,t_j)=∑_{k=1}^{p} x_k(s_i,t_j)β_k,
where x_1(s_i,t_j),..., x_p(s_i,t_j) are known functions of (s_i,t_j), and β_1,...,β_p are unknown parameters to be estimated. Equivalently, in matrix notation, we have the spatio-temporal linear model as follows:
Y = X β + Z + ε,
Z ~ N(0,σ^2 Ω_{φρ}),
ε ~ N(0,τ^2 I_m).
Therefore the spatio-temporal process, Y, has normal distribution with mean E[Y]=Xβ and variance Σ=σ^2Ω_{φρ}+τ^2 I_m. We assume that Σ is non-singular and X has full rank.
The estimation process was computed via SAEM algorithm initially proposed by \insertCitedelyon1999convergence;textualStempCens.
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 θ = (β, σ^2, τ^2, φ, ρ). |
Theta |
estimated parameters in all iterations, θ = (β, σ^2, τ^2, φ, ρ). |
beta |
estimated β. |
sigma2 |
estimated σ^2. |
tau2 |
estimated τ^2. |
phi |
estimated φ. |
rho |
estimated ρ. |
Eff.range |
estimated effective range. |
PsiInv |
estimated Ψ^-1, where Ψ=Σ/σ^2. |
Cov |
estimated Σ. |
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. |
Katherine L. Valeriano, Victor H. Lachos and Larissa A. Matos
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 30 31 32 33 34 35 36 37 38 39 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.