The goal of StempCens
is to estimates the parameters of a censored or
missing data in spatio-temporal models using the SAEM algorithm (Delyon,
Lavielle, and Moulines 1999). This algorithm is a stochastic
approximation of the widely used EM algorithm and an important tool for
models in which the E-step does not have an analytic form. Besides the
expressions obtained to estimate the parameters to the proposed model,
we include the calculations for the observed information matrix using
the method developed by Louis (1982). To examine the performance of the
fitted model, case-deletion measure are provided (see also Cook 1977;
Zhu et al. 2001). Moreover, it computes the spatio-temporal covariance
matrix and the effective range for an isotropic spatial correlation
function.
You can install the released version of StempCens from CRAN with:
install.packages("StempCens")
StempCens
package provides five functions:
CovarianceM
: Computes the spatio-temporal covariance matrix for
balanced data.EffectiveRange
: Computes the effective range for an isotropic
spatial correlation function.EstStempCens
: Returns the maximum likelihood estimates of the
unknown parameters.PredStempCens
: Performs spatio-temporal prediction in a set of new
spatial locations for fixed time points.CrossStempCens
: Performs cross-validation, which measure the
performance of the predictive model on new test dataset.DiagStempCens
: Returns measures and graphics for diagnostic
analysis.This is a basic example which shows you how to solve a problem using
functions EstStempCens
(parameter estimation) and PredStempCens
(prediction in new locations):
library(StempCens)
# Initial parameter values
beta <- c(-1,1.50)
phi <- 5; rho <- 0.60
tau2 <- 0.80; sigma2 <- 2
# Simulating data
n1 <- 17 # Number of spatial locations
n2 <- 5 # Number of temporal index
set.seed(12345)
x.co <- round(runif(n1,0,10),9) # X coordinate
y.co <- round(runif(n1,0,10),9) # Y coordinate
coord <- cbind(x.co,y.co) # Cartesian coordinates without repetitions
coord2 <- cbind(rep(x.co,each=n2),rep(y.co,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) # Covariates
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,0.50,"pow.exp")
# Data
require(mvtnorm)
y <- as.vector(rmvnorm(1,mean=as.vector(media),sigma=Cov))
data <- data.frame(coord2,time2,y,x)
names(data) <- c("x.coord","y.coord","time","yObs","x1","x2")
# Splitting the dataset
local.est <- coord[-c(4,13),]
data.est <- data[data$x.coord%in%local.est[,1]&data$y.coord%in%local.est[,2],]
data.valid <- data[data$x.coord%in%coord[c(4,13),1]&data$y.coord%in%coord[c(4,13),2],]
# Censored
perc <- 0.10
y <- data.est$yObs
aa <- sort(y); bb <- aa[1:(perc*nrow(data.est))]; cutof <- bb[perc*nrow(data.est)]
cc <- matrix(1,nrow(data.est),1)*(y<=cutof)
y[cc==1] <- cutof
data.est <- cbind(data.est[,-c(4,5,6)],y,cc,data.est[,c(5,6)])
names(data.est) <- c("x.coord","y.coord","time","yObs","censored","x1","x2")
# Estimation
y <- data.est$yObs
x <- cbind(data.est$x1,data.est$x2)
cc <- data.est$censored
time2 <- matrix(data.est$time)
coord2 <- data.est[,1:2]
LI <- y; LI[cc==1] = -Inf # Left-censored
LS <- y
est_teste <- EstStempCens(y, x, cc, time2, coord2, LI, LS, init.phi=3.5,
init.rho=0.5, init.tau2=1, kappa=0.5, type.S="pow.exp",
IMatrix=TRUE, M=20, perc=0.25, MaxIter=300, pc=0.20)
# Prediction
locPre <- data.valid[,1:2]
timePre <- matrix(data.valid$time)
xPre <- cbind(data.valid$x1,data.valid$x2)
pre_teste <- PredStempCens(est_teste, locPre, timePre, xPre)
library(ggplot2)
Model <- rep(c("y Observed","y Predicted"),each=10)
station <- rep(rep(c("Station 1", "Station 2"),each=5),times=2)
xcoord1 <- rep(seq(1:5),4)
ycoord1 <- c(data.valid$yObs,pre_teste$predValues)
data2 <- data.frame(Model,station,xcoord1,ycoord1)
ggplot(data=data2,aes(x=xcoord1,y=ycoord1)) + geom_line(aes(color=Model)) +
facet_wrap(station~.,nrow=2) + labs(x="",y="") + theme(legend.position="bottom")
For the diagnostic analysis in the EstStempCens
function the parameter
input IMatrix
needs to be TRUE
.
diag <- DiagStempCens(est, type.diag="location", diag.plot = TRUE, ck=1)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.