EstStempCens {StempCens} | R Documentation |
Return the maximum likelihood estimates of the unknown parameters of spatio-temporal model with censored/missing responses. The estimates are obtained using SAEM algorithm. Also, the function computes the observed information matrix using the method developed by Thomas (1982). The types of censoring considered are left, right or missing values.
EstStempCens(y, x, cc, tempo, coord, inits.phi, inits.rho, inits.tau2, tau2.fixo = FALSE, type.Data = "balanced", cens.type = "left", 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 = 10^-6)
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. |
tempo |
a vector of time. |
coord |
a matrix of coordinates of the spatial locations. |
inits.phi |
initial value of the spatial scaling parameter. |
inits.rho |
initial value of the time scaling parameter. |
inits.tau2 |
initial value of the the nugget effect parameter. |
tau2.fixo |
|
type.Data |
type of the data: ' |
cens.type |
type of censoring: ' |
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. |
type.S |
type of spatial covariance 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 aproximation. 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.2 |
error |
the convergence maximum error. By default = 10^-6. |
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 index for Y(s_{i},t_{j}). We assume Z is normally distributed with zero-mean and covariance matrix Σ_z = σ^2 Ω_{φρ}(s,t) where σ^2 is the model variance, φ 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 function 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 and 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 our spatio-temporal linear model as follows:
Y = X β + Z + ε,
Z ~ N(0,σ^2 Ω_{φρ}),
ε ~ N(0,τ^2 I).
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. Using standard geostatistical terms, τ^2 is the nugget effect (or the measurement error variance), σ^2 the sill, and Ω_{φρ}=[r_{ij}] is the m x m spatio-temporal correlation matrix with diagonal elements r_{ii}=1, for i=1,...,m.
The estimation process was computed via SAEM algorithm initially proposed by Deylon et. al.(1999).
The function returns an object of class Est.StempCens
which is a list given by:
Returns a list with all data components given in input.
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 ρ. |
PsiInv |
estimated Ψ^-1. |
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 A. L. Valeriano, Victor H. Lachos and Larissa A. Matos
# Initial parameter values beta <- c(-1,1.50); phi <- 5; rho <- 0.45; tau2 <- 0.80; sigma2 <- 1.5 # Simulating data n1 <- 9 # Number of spatial locations n2 <- 4 # Number of temporal index set.seed(700) x.coord <- round(runif(n1,0,10),9) # X coordinate y.coord <- round(runif(n1,0,10),9) # Y coordinate coordenadas <- 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,1)) # 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 H <- as.matrix(dist(coordenadas)) # Spatial distances Mt <- as.matrix(dist(time)) # Temporal distances Cov <- CovarianceM(phi,rho,tau2,sigma2,distSpa=H,disTemp=Mt,kappa=0,type.S="exponential") # Data require(mvtnorm) y <- as.vector(rmvnorm(1,mean=as.vector(media),sigma=Cov)) perc <- 0.2 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 # Estimation est_teste <- EstStempCens(y, x, cc, time2, coord2, inits.phi=3.5, inits.rho=0.5, inits.tau2=0.7, type.Data="balanced", cens.type="left", 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=25, pc=0.2, error = 10^-6) ## Not run: # New York data library(spTimer) library(sp) library(rgdal) # Transform coordinates station <- as.data.frame(NYdata[,2:3]) names(station) <- c("lon","lat") coordinates(station) <- ~lon+lat proj4string(station) <- CRS("+init=epsg:32116 +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") station2 <- spTransform(station,CRS("+proj=utm +zone=18 +north +datum=NAD83")) station <- as.data.frame(station2) names(station) <- c("x.Coord","y.Coord") NYdata <- cbind(NYdata,station) coord <- unique(station) EstEstimation <- c(1,2,3,5,6,7,8,9,11,12,13,14,15,16,17,18,19,20,22,23,24,26,27,28) EstEstimation <- coord[EstEstimation,] dataEstimation <- NYdata[(NYdata[,11]%in%EstEstimation[,1])&(NYdata[,12]%in%EstEstimation[,2]),] cc <- vector("numeric",length=nrow(dataEstimation)) cc[is.na(dataEstimation$o8hrmax)==1] <- 1 time <- rep(seq(1,62),nrow(EstEstimation)) dados <- cbind(dataEstimation,cc,time) names(dados) <- c("s.index","Longitude","Latitude","Year","Month","Day","o8hrmax","cMAXTMP", "WDSP","RH","x.Coord","y.Coord","censure","time") # Data y <- dados$o8hrmax cc <- dados$censure x <- as.matrix(cbind(dados$cMAXTMP,dados$WDSP,dados$RH)) tempo <- as.data.frame(dados$time) coord <- as.matrix(cbind(dados$x.Coord/1000,dados$y.Coord/1000)) # Coordinates in Km # Power exponential model for the spatial covariance structure est <- EstStempCens(y, x, cc, tempo, coord, inits.phi = 50, inits.rho = 0.30, inits.tau2 = 25, type.Data="balanced", cens.type="missing", method="nlminb", kappa=0.50, type.S="pow.exp", IMatrix=TRUE, lower.lim=c(0.01,-0.80,0.01), upper.lim=c(500,0.80,150), M=20, perc=0.25, MaxIter=500, pc=0.2, error = 10^-5) ## End(Not run)