R/useful_commands.R

In spTimer: Spatio-Temporal Bayesian Modelling

##
## Summary Statistics another built in
##
spT.MCMC.stat<-function(x, nBurn=0)
{
  options(warn=-1)
  if(!inherits(x, "spT")){
    stop("\n# Error: provide valid posterior output \n")
  }
  model<-x$model
  nItr<-x$iterations
  if(!is.null(nBurn)){
    if(nBurn == 0){
      nBurn<-0
      nItr<-nItr-x$nBurn
    }
    if(nBurn != 0){ 
      nBurn<-nBurn
      nItr<-nItr-x$nBurn
    }
  }
  if((nBurn+x$nBurn) >= (nItr+x$nBurn)){
    cat("# Number of Iterations:            ", nItr+x$nBurn, "\n")
    cat("# Number of Burn-in (fitted model):", x$nBurn, "\n")
    cat("# More Burn-in:                    ", nBurn, "\n")
    cat("# Total Number of Burn-in:         ", nBurn+x$nBurn, "\n")
    stop("\n# Error: iteration (",nItr+x$nBurn,") is less than or equal to total burn-in (",nBurn+x$nBurn,") \n")
  }
  cat("\n")
  cat("# Model:", model, "\n")
  if(is.null(model)==TRUE){
    stop("\n# Error: need to define the model")
  }
  else if(model=="AR"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-dim(x$mu_lp)[[1]]
    p<-dim(x$betap)[[1]]
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
    }
    #
    para<-spT.Summary.Stat(para)
    dimnames(para)[[1]][1:(1+p+3)]<-c(dimnames(x$X)[[2]],"rho",
                                      "sig2eps","sig2eta","phi")
    round(para,4)
  }
  else if(model == "GPP"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-x$r
    p<-x$p
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
    }
    #
    para<-spT.Summary.Stat(para)
    dimnames(para)[[1]][1:(1+p+2+1)]<-c(dimnames(x$X)[[2]],"rho",
                                        "sig2eps","sig2eta","phi")
    round(para,4)
  }
  else if(model == "GP"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-x$r
    p<-x$p
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
    }
    #
    para<-spT.Summary.Stat(para)
    dimnames(para)[[1]]<-c(dimnames(x$X)[[2]],
                           "sig2eps","sig2eta","phi")
    round(para,4)
  }
  else{
    
  }
}
##
## MCMC plot another built in 
##
spT.MCMC.plot<-function(x, nBurn=0, ACF="FALSE", PARTIAL.acf="FALSE")
{
  options(warn=-1)
  if(!inherits(x, "spT")){
    stop("\n# Error: provide valid posterior output \n")
  }
  model<-x$model
  nItr<-x$iterations
  if(!is.null(nBurn)){
    if(nBurn == 0){
      nBurn<-0
      nItr<-nItr-x$nBurn
    }
    if(nBurn != 0){ 
      nBurn<-nBurn
      nItr<-nItr-x$nBurn
    }
  }
  if((nBurn+x$nBurn) >= (nItr+x$nBurn)){
    cat("# Number of Iterations:            ", nItr+x$nBurn, "\n")
    cat("# Number of Burn-in (fitted model):", x$nBurn, "\n")
    cat("# More Burn-in:                    ", nBurn, "\n")
    cat("# Total Number of Burn-in:         ", nBurn+x$nBurn, "\n")
    stop("\n# Error: iteration (",nItr+x$nBurn,") is less than or equal to total burn-in (",nBurn+x$nBurn,") \n")
  }
  cat("\n")
  cat("# Model:", model, "\n")
  if(is.null(model)==TRUE){
    stop("\n# Error: need to define the model")
  }
  else if(model=="AR"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-x$r
    p<-x$p
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1+1)]<-c(dimnames(x$X)[[2]],"rho",
                                            "sig2eps","sig2eta","phi","nu")
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1)]<-c(dimnames(x$X)[[2]],"rho",
                                          "sig2eps","sig2eta","phi")
    }
    #
    #
  }
  else if(model == "GPP"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-dim(x$mu_lp)[[1]]
    p<-dim(x$betap)[[1]]
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1+1)]<-c(dimnames(x$X)[[2]],"rho",
                                            "sig2eps","sig2eta","phi","nu")
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$rhop[(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1)]<-c(dimnames(x$X)[[2]],"rho",
                                          "sig2eps","sig2eta","phi")
    }
    #
  }
  else if(model == "GP"){
    cat("\n")
    cat("# Number of Iterations:   ", nItr+x$nBurn, "\n")
    cat("# Total number of Burn-in:", nBurn+x$nBurn, "\n")
    cat("\n")
    if(nItr <= nBurn){
      stop("\n# Error: iteration (",nItr,") is less than or equal to nBurn (",nBurn,") \n")
    }
    r<-x$r
    p<-x$p
    #
    if(x$cov.fnc=="matern"){
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]),t(x$nup[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1+1)]<-c(dimnames(x$X)[[2]],
                                            "sig2eps","sig2eta","phi","nu")
    }
    else{
      para<-rbind((x$betap[,(nBurn+1):nItr]),
                  t(x$sig2ep[(nBurn+1):nItr]),
                  t(x$sig2etap[(nBurn+1):nItr]),
                  t(x$phip[(nBurn+1):nItr]))
      dimnames(para)[[1]][1:(1+p+2+1)]<-c(dimnames(x$X)[[2]],
                                          "sig2eps","sig2eta","phi")
    }
    #
    #
  }
  else{
    
  }
  #
  ##
  dev.new()
  for(i in 1:dim(para)[[1]]){
    MCMC.plot.obj(para[i,],name=c(dimnames(para)[[1]][i]),ACF=ACF,PARTIAL.acf=PARTIAL.acf)
    par(ask=TRUE)
  }
}
##
## segment plot for upper and lower limit
##  
spT.segment.plot<-function(obs, est, up, low, limit=NULL){
  #
  tmp<-cbind(obs,est,up,low)
  tmp<-na.omit(tmp)
  if(is.null(limit)==TRUE){
    plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions", pch="*",
         xlim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)),
         ylim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)))
  }
  #
  else{
    plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions",
         xlim=c(limit[1],limit[2]),ylim=c(limit[1],limit[2]),pch="*")
  }
  #
  segments(tmp[,1],tmp[,2],tmp[,1],tmp[,3])
  segments(tmp[,1],tmp[,2],tmp[,1],tmp[,4])
  #  
}
##
## hit and false alarm function for forecast
## code for hit and false 
spT.validation2<-function(z,zhat,cutoff,names=FALSE){
  #
  obs <- z; fore <- zhat; tol <- cutoff;
  tmp<-cbind(obs,fore)
  tmp<-na.omit(tmp)
  #
  c11<-tmp[tmp[,1]<=tol & tmp[,2]<=tol,]
  c11<-length(c11)/2
  c12<-tmp[tmp[,1]<=tol & tmp[,2]>tol,]
  c12<-length(c12)/2  
  c21<-tmp[tmp[,1]>tol & tmp[,2]<=tol,]
  c21<-length(c21)/2  
  c22<-tmp[tmp[,1]>tol & tmp[,2]>tol,]
  c22<-length(c22)/2
  mat<-matrix(c(c11,c21,c12,c22),2,2)
  dimnames(mat)[[1]]<-c(paste("[Obs:<=",tol,"]"),paste("[Obs:> ",tol,"]")) 
  dimnames(mat)[[2]]<-c(paste("[Forecast:<=",tol,"]"),paste("[Forecast:>",tol,"]")) 
  # Type I error (false positive): mat[1,2]
  # Type II error (false negative): mat[2,1]
  # True positive: mat[1,1]
  # True negative: mat[2,2]
  TPR <- mat[1,1]/sum(mat[,1]) # True Positive Rate, Sensitivity, Hit rate, Recall
  FPR <- mat[1,2]/sum(mat[,2]) # False Positive Rate, False alarm
  FNR <- 1-TPR # False Negative Rate, Miss rate 
  TNR <- 1-FPR # True Negative Rate, Specificity 
  Prevalence <- sum(mat[,1])/sum(c(mat))
  Accuracy <- sum(diag(mat))/sum(c(mat))
  Precision <- mat[1,1]/sum(mat[1,]) # positive predictive value
  FOR <- mat[2,1]/sum(mat[1,]) # false ommission rate
  LRp <- TPR/FPR # positive likelihood ratio 
  LRn <- FNR/TNR # negative likelihood ratio
  FDR <- mat[1,2]/sum(mat[1,]) # false discovery rate 
  NPV <- mat[2,2]/sum(mat[2,]) # negative predictive value
  DOR <- LRp/LRn # diagnostic odds ratio 
  F1score <- 2/(1/TPR+1/Precision)
  top<-2*(mat[1,1]*mat[2,2]-mat[1,2]*mat[2,1])
  bot<-mat[1,2]^2+mat[2,1]^2+2*mat[1,1]*mat[2,2]+(mat[1,2]+mat[2,1])*sum(diag(mat))
  S<-round(top/bot,4) # Heidke Skill
  x<-list(
    TPR=TPR, FPR=FPR, FNR=FNR, TNR=TNR, Prevalence=Prevalence, Accuracy=Accuracy,
    Precision=Precision, FOR=FOR, LRp=LRp, LRn=LRn, FDR=FDR, NPV=NPV, DOR=DOR,
    F1.Score=F1score,Heidke.Skill=S)
  ##
  if(names==TRUE){
    cat("##\n True Positive Rate/ Sensitivity/ Hit rate/ Recall \n False Positive Rate/ False alarm \n False Negative Rate/ Miss rate \n True Negative Rate/ Specificity \n Prevalence \n Accuracy \n Precision/ Positive Predictive Value \n False Ommission Rate \n Positive Likelihood Ratio \n Negative Likelihood Ratio \n False Discovery Rate \n Negative Predictive Value \n Diagnostic Odds Ratio \n F1score \n Heidke Skill \n##\n") 
  }
  ##
  unlist(x)
}
##
##
## For data split
##
spT.subset<-function (data, var.name, s = NULL, reverse = FALSE) 
{
  #
  # This function is to select and deduct the sites used to
  # fit or valid the model
  # Input: 	data
  #		s = the site numbers to be selected/deselected, e.g., c(4,7,10)
  
  if (missing(var.name)) {
    stop("Error: need to define the var.name")
  }
  if(!is.data.frame(data)){
    stop("Error: data should be in data.frame")
  }
  if(!var.name %in%names(data)){
    stop("Error: var.name should be in data")
  }
  data<-data[,c(var.name,names(data)[!names(data)%in%var.name])]
  s.index <- unique(data[, var.name[1]])
  if (reverse == FALSE) {
    data <- data[data[,1] %in% (s.index[s.index %in% s]),]
    data
  }
  else {
    data <- data[data[,1] %in% (s.index[!s.index %in% s]),]
    data
  }
}
##
## grid coordinates
##
spT.grid.coords<-function(Longitude=c(max,min),Latitude=c(max,min),by=c(NA,NA))
{
  max.lon <- Longitude[1]
  min.lon <- Longitude[2]
  max.lat <- Latitude[1]
  min.lat <- Latitude[2]
  if(is.na(by[1]) || is.na(by[2])){
    stop('Error: need to specify grid dimension, n x m')
  }
  knots.lon <- seq(max.lon,min.lon,length=by[1])
  knots.lat <- seq(max.lat,min.lat,length=by[2])
  #knots.coords <- cbind(rep(knots.lon,by),c(outer(rep(1,by),knots.lat)))
  knots.coords <- cbind(c(outer(rep(1,by[2]),knots.lon)),rep(knots.lat,by[1]))
  as.matrix(knots.coords)
}
##
## Percentage Coverage
##
spT.pCOVER<-function(z=NULL,zup=NULL,zlow=NULL,zsample=NULL,level=95)
{
  if(is.null(z)){
    stop('\n # Provide observed values\n')
  }
  if(!is.null(zsample)){
    z<-c(z)
    low<-(1-level/100)/2
    up<-1-low
    x<-apply(zsample,1,quantile,prob=c(low,up))
    x<-rbind(x,z)
    x<-t(x)
    x<-na.exclude(x)
    y<-x[x[, 1] <= x[, 3] & x[, 3] <= x[, 2], ]
    round(length(y[,1])/length(x[,1])*100,2)
  }
  else if(is.null(zsample)){
    z<-as.matrix(z)
    zlow<-as.matrix(zlow)
    zup<-as.matrix(zup)
    x<-cbind(zlow,z,zup)
    x<-na.exclude(x)
    y<- x[x[,1] <= x[,2] & x[,2] <= x[,3],]
    round(length(y[,1])/length(x[,1])*100,2)
  }
}
##
## Geodetic distance in K.M. or Miles
##
spT.geodist<-function(Lon, Lat, KM=TRUE){
  #
  # This function is for geodistance using C/C++
  #
  n<-length(Lat)
  if(KM==TRUE){	
    ds<-.C("GeoDist_km", as.integer(n),as.double(Lat),
           as.double(Lon),out=matrix(double(n*n),n,n))$out
  }
  else {
    ds<-.C("GeoDist_miles", as.integer(n),as.double(Lat),
           as.double(Lon),out=matrix(double(n*n),n,n))$out
  }
  #dis<-matrix(ds,length(Lat),length(Lat))
  #return(dis)
  ds
}
##
## Geodetic distance: another approach: two points
##
spT.geo.dist <- function(point1, point2){
  #
  # The following program computes the distance on the surface 
  # of the earth between two points point1 and point2. 
  # Both the points are of the form (Longitude, Latitude)
  #
  R <- 6371
  p1rad <- point1 * pi/180
  p2rad <- point2 * pi/180
  d <- sin(p1rad[2])*sin(p2rad[2])+cos(p1rad[2])*cos(p2rad[2])*cos(abs(p1rad[1]-p2rad[1]))	
  d <- acos(d)
  R*d
}
##
## Geodetic distance: using points 1:(lon, lat) 2:(lon, lat)
##
spT.geo_dist <- function(points)
{
  point1 <- c(points[1], points[2])
  point2 <- c(points[3], points[4])
  #The following program computes the distance on the surface 
  #The argument points is of the form (Long, Lat, Long, Lat)
  dist <- 0
  if(sum((point1 - point2)^2) > 1e-05) {
    R <- 6371
    p1rad <- (point1 * pi)/180
    p2rad <- (point2 * pi)/180
    d <- sin(p1rad[2]) * sin(p2rad[2]) + cos(p1rad[2]) * cos(p2rad[
      2]) * cos(abs(p1rad[1] - p2rad[1]))
    d <- acos(d)
    dist <- R * d
  }
  dist
}
##
## To keep the long lat positions based on tol.dist
##
spT.keep.morethan.dist <- function(coords, tol.dist=100)  
{
  #
  # coords must have two columns named long and lat 
  #
  a <- as.data.frame(coords)
  names(a) <- c("long","lat")
  n <- nrow(coords)
  c1 <- rep(1:n, each=n)
  c2 <- rep(1:n, n)
  b <- matrix(NA, nrow=n, ncol=n)
  w <- as.vector(upper.tri(b))
  bigmat <- matrix(0, nrow=n*n, ncol=7)
  bigmat[, 1] <- c1
  bigmat[, 2] <- c2
  bigmat[, 3] <- a$long[c1]
  bigmat[, 4] <- a$lat[c1]
  bigmat[, 5] <- a$long[c2]
  bigmat[, 6] <- a$lat[c2]
  ubmat <- bigmat[w, ]
  ubmat[,7] <- as.vector(apply(ubmat[,3:6], 1, spT.geo_dist)) 
  
  v <- ubmat[,7]
  w <- ubmat[v<tol.dist, ]
  
  z <- unique(w[,1])
  a <- coords[-z,  ]
  a
}
##
## This function extracts the data using formula
##
Formula.matrix <- function(formula, data, na.rm=FALSE)
{
  if(na.rm == FALSE){
    mt <- terms.formula(formula, data=data, specials=c("sp","tp"))
    if(missing(data)) data <- sys.frame(sys.parent())
    spcheck<-attr(mt, "specials")$sp
    tpcheck<-attr(mt, "specials")$tp
    mf <- model.frame(mt, data=data, na.action=na.pass)
  }
  else{
    mt <- terms.formula(formula, data=data, specials=c("sp","tp"))
    if(missing(data)) data <- sys.frame(sys.parent())
    spcheck<-attr(mt, "specials")$sp
    tpcheck<-attr(mt, "specials")$tp
    mf <- model.frame(mt, data=data)
  }
  # null model support
  "%w/o%" <- function(x, y) x[!x %in% y]
  x.names <- attr(mt, "term.labels") # X variable names
  X <- if (!is.empty.model(mt)) model.matrix(mt, mf)
  if((!is.null(spcheck)) & (is.null(tpcheck))){
    # only spatially varying  
    x.names.sp <- x.names[attr(mt, "specials")$sp-1]
    X.sp <- X[,c(x.names.sp)]
    x.names <- dimnames(X)[[2]]
    x.names <- x.names %w/o% x.names.sp
    X <- X[,c(x.names)]
    X <- as.matrix(X)
    dimnames(X)[[2]]<-as.list(x.names) 
    X.sp <- as.matrix(X.sp) 
    dimnames(X.sp)[[2]]<-as.list(x.names.sp) 
    X.tp <- NULL; x.names.tp<-NULL
    if(dim(X)[[2]]==0){
      X<-as.matrix(rep(0, dim(mf)[[1]])); 
      x.names <- "(Intercept)"
    }  
  }
  else if((!is.null(tpcheck)) & (is.null(spcheck))) {
    # only temporally varying
    x.names.tp <- x.names[attr(mt, "specials")$tp-1]
    X.tp <- X[,c(x.names.tp)]
    x.names <- dimnames(X)[[2]]
    x.names <- x.names %w/o% x.names.tp
    X <- X[,c(x.names)]
    X <- as.matrix(X)
    dimnames(X)[[2]]<-as.list(x.names) 
    X.tp <- as.matrix(X.tp) 
    dimnames(X.tp)[[2]]<-as.list(x.names.tp) 
    X.sp <- NULL; x.names.sp<-NULL 
    if(dim(X)[[2]]==0){
      X<-as.matrix(rep(0, dim(mf)[[1]])); 
      x.names <- "(Intercept)"
    }  
  }
  else if((!is.null(spcheck)) & (!is.null(tpcheck))){
    # both spatially and temporally varying
    x.names.sp <- x.names[attr(mt, "specials")$sp-1]
    x.names.tp <- x.names[attr(mt, "specials")$tp-1]
    X.sp <- X[,c(x.names.sp)]
    X.tp <- X[,c(x.names.tp)]
    x.names <- dimnames(X)[[2]]
    x.names <- x.names %w/o% x.names.sp
    x.names <- x.names %w/o% x.names.tp
    X <- X[,c(x.names)]
    X <- as.matrix(X)
    dimnames(X)[[2]]<-as.list(x.names) 
    X.sp <- as.matrix(X.sp) 
    dimnames(X.sp)[[2]]<-as.list(x.names.sp) 
    X.tp <- as.matrix(X.tp) 
    dimnames(X.tp)[[2]]<-as.list(x.names.tp) 
    if(dim(X)[[2]]==0){
      X<-as.matrix(rep(0, dim(mf)[[1]])); 
      x.names <- "(Intercept)"
    }  
  }
  else {
    x.names <- dimnames(X)[[2]] 
    X.sp <- NULL; x.names.sp<-NULL; X.tp <- NULL; x.names.tp<-NULL 
  }
  Y <- as.matrix(model.response(mf, "numeric")) # Y matrix
  #
  return(list(Y=Y, X=X, x.names=x.names, X.sp=X.sp, x.names.sp=x.names.sp, X.tp=X.tp, x.names.tp=x.names.tp))
}
##
## This function extracts the coords using formula
##
Formula.coords <- function(formula, data, na.rm=TRUE)
{
  mt <- terms(formula, data=data)
  if(missing(data)) data <- sys.frame(sys.parent())
  mf <- model.frame(mt, data=data)
  # null model support
  X <- if (!is.empty.model(mt)) model.matrix(mt, mf)
  X <- as.matrix(X[,2:3])         # X matrix
  X <- unique(X); dimnames(X)<-NULL
  X
}
##
## Summary Statistics
##
spT.Summary.Stat <- function(y) 
{
  
  ## define Upper and lower limit function
  up.low.limit<-function(y,limit)
  {
    #
    y<-sort(y)
    y<-y[limit]
    y
    #
  }
  #
  if(is.vector(y)==TRUE){
    y <- matrix(y[!is.na(y)])
  }
  else{
    y <- t(y)
  } 
  N <- length(y[1,  ])
  nItr <- length(y[, 1])
  z <- matrix(nrow = N, ncol = 5)
  dimnames(z) <- list(dimnames(y)[[2]], c("Mean","Median","SD","Low2.5p","Up97.5p"))
  #
  if (nItr < 40) {
    stop("\n##\n# Error: number of samples must be >= 40\n##")
  }
  z[, 1] <- apply(y,2,mean)
  z[, 2] <- apply(y,2,median)
  z[, 3] <- apply(y,2,sd)
  #z[, 4] <- apply(y,2,quantile,0.025)
  #z[, 5] <- apply(y,2,quantile,0.975)
  nl <- as.integer(nItr * 0.025)
  nu <- as.integer(nItr * 0.975)
  z[, 4] <- apply(y,2,up.low.limit,limit=nl)
  z[, 5] <- apply(y,2,up.low.limit,limit=nu)
  
  as.data.frame(z)
}
##
## Predictive Model Choice Criteria
##
PMCC<-function(z=NULL, z.mean=NULL, z.sd=NULL, z.samples=NULL)
{
  #
  # Predictive Model Choice Criteria
  #
  if(is.null(z)){
    stop("Error: need to provide z values.")
  }
  #
  if(is.null(z.mean) | is.null(z.sd)){
    if(!is.null(z.mean)){
      stop("Error: need to provide z.sd value.")
    }   
    if(!is.null(z.sd)){
      stop("Error: need to provide z.mean value.")
    }   
    if(is.null(z.samples)){
      stop("Error: need to provide z.samples value.")
    }
  }
  #
  if(!is.null(z.samples)){
    if ( !is.matrix(z.samples)) {
      stop("Error: z.samples must be a (N x nItr) matrix")
    }
    if (dim(z.samples)[1] != length(z)) {
      stop("Error: observations in z.samples in each iteration must be equal to length of z")
    }
    if ( dim(z.samples)[2] < 40) {
      stop("Error: samples are too small to obtain summary statistics")
    }
    sum.stat = matrix(NA,length(c(z)),6)
    sum.stat[,1:5] = as.matrix(spT.Summary.Stat(z.samples))
    sum.stat[,6] = c(z)
    sum.stat = sum.stat[!is.na(sum.stat[,6]),]
    goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,6])^2),2)
    penalty = round(sum(sum.stat[,3]^2),2)
    pmcc =  round(goodness.of.fit + penalty,2)
    out = NULL
    out$pmcc = pmcc; 
    out$goodness.of.fit = goodness.of.fit
    out$penalty = penalty
    #class(out) <- "PMCC"
    out
  }
  else{
    if(is.null(z.mean) | is.null(z.sd)){
      stop("Error: need to provide z.mean and/or z.sd values.")
    }
    if(length(c(z)) != length(c(z.mean))){
      stop("Error: z and z.mean should be in same length.")
    }
    if(length(c(z)) != length(c(z.sd))){
      stop("Error: z and z.sd should be in same length.")
    }
    sum.stat = matrix(NA,length(c(z)),3)
    sum.stat[,1] = c(z)
    sum.stat[,2] = c(z.mean)
    sum.stat[,3] = c(z.sd)
    sum.stat = sum.stat[!is.na(sum.stat[,1]),]
    goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,2])^2),2)
    penalty = round(sum(sum.stat[,3]^2),2)
    pmcc =  round(goodness.of.fit + penalty,2)
    out = NULL
    out$pmcc = pmcc; 
    out$goodness.of.fit = goodness.of.fit
    out$penalty = penalty
    #class(out) <- "PMCC"
    out
  }
}

##
## To identify the Integer Number
##

#is.wholenumber<-function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

##
## To check sites that has <tol km of distances
## 
spT.check.locations<-function(fit.locations, pred.locations, method, tol){
  #
  #
  if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
                    "maximum", "manhattan", "canberra")){
    stop("\n# Error: correctly define distance.method \n")
  }
  #
  coords.all <- rbind(fit.locations,pred.locations)
  tn.fitsites <- length(fit.locations[, 1])
  nfit.sites <- 1:tn.fitsites
  tn.predsites <- length(coords.all[, 1]) - tn.fitsites
  npred.sites <- (tn.fitsites + 1):(length(coords.all[, 1]))
  #
  if(method=="geodetic:km"){
    coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
  }
  else if(method=="geodetic:mile"){
    coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
  }
  else{
    coords.D <- as.matrix(dist(coords.all, method=method, diag = T, upper = T))
  }  
  coords.D[is.na(coords.D)]<-0
  #
  diag(coords.D)<-NA
  # 
  fdmis<-coords.D[nfit.sites, npred.sites]
  if(is.matrix(fdmis)==TRUE){
    fdmis<-cbind(c(t(fdmis)),1:dim(fdmis)[[2]],sort(rep(1:dim(fdmis)[[1]],dim(fdmis)[[2]]))) # add pred sites and fitted sites
    fdmis<-fdmis[fdmis[,1] < tol,]
    #
    if(!is.na(fdmis[1])==TRUE){
      cat("#\n# Tolerance Limit:", paste(tol))
      cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
      fdmis<-matrix(fdmis,(length(fdmis)/3),3) 
      for(i in 1:dim(fdmis)[[1]]){
        print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",fdmis[i,2]," Fitted location:", fdmis[i,3],""))
      }
      cat("#\n# Romove the locations and run again. \n#\n")
      dimnames(fdmis)[[2]]<-c('distance','pred_location','fit_location')
      fdmis
    }
    else{
      cat("#\n# Tolerance Limit (unit):", paste(tol))
      #cat("\n# Fitted and Predicted location distances are alright \n#\n")  
      cat("\n# Location distances are alright \n#\n")  
    }
  }
  else{
    fdmis<-cbind(c(fdmis),1:length(fdmis)) # 
    fdmis<-fdmis[fdmis[,1] < tol,]
    if(!is.na(fdmis[1])==TRUE){
      cat("#\n# Tolerance Limit:", paste(tol))
      cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
      fdmis<-as.matrix(fdmis) 
      for(i in 1:dim(fdmis)[[1]]){
        print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",1," Fitted location:", fdmis[i,2],""))
      }
      cat("#\n# Romove the locations and run again. \n#\n")
      dimnames(fdmis)[[2]]<-c('distance','fit_location')
      fdmis
    }
    else{
      cat("#\n# Tolerance Limit (unit):", paste(tol))
      #cat("\n# Fitted and Predicted location distances are alright \n#\n")  
      cat("\n# Location distances are alright \n#\n")  
    }
  }
}
##
## To check sites inside codes
## 
spT.check.sites.inside<-function(coords, method, tol=0.1){
  #
  #
  if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
                    "maximum", "manhattan", "canberra")){
    stop("\n# Error: correctly define distance.method \n")
  }
  #
  #
  if(method=="geodetic:km"){
    fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=TRUE))
  }
  else if(method=="geodetic:mile"){
    fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=FALSE))
  }
  else{
    fdm<- as.matrix(dist(coords, method, diag = TRUE, upper = TRUE))
  } 
  #
  diag(fdm)<-NA
  # 
  fdm<-cbind(c(fdm),1:dim(fdm)[[2]],sort(rep(1:dim(fdm)[[1]],dim(fdm)[[2]]))) 
  #
  fdm<-fdm[!is.na(fdm[,1]),]
  #
  #tol <- 0.01
  fdmis<-fdm[fdm[,1] < tol,]
  #
  if(!is.na(fdmis[1])==TRUE){
    cat("#\n# There are some locations very close:")
    print(paste("( < ",tol," unit) to each other."))
    fdmis<-matrix(fdmis,(length(fdmis)/3),3) 
    for(i in 1:dim(fdmis)[[1]]){
      print(paste("Distance (unit):", round(fdmis[i,1],2)," site:",fdmis[i,2]," site:", fdmis[i,3],""))
    }
    dimnames(fdmis)[[2]]<-c('dist_km','pred_site','fit_site')
    fdmis
    cat("#\n# Romove the sites and run again. \n#\n")
    stop("Error: Termination.")
  }
  #
  #tol <- 0.5
  #fdmis<-fdm[fdm[,1] < tol,]
  #
  #if(!is.na(fdmis[1])==TRUE){
  # cat("#\n# Warnings: There are some locations very close ( < 0.5 unit) to each other.\n#\n")
  #}
  #
}
##
##
#check.sites.inside.pred<-spT.check.locations
##
## convert seconds into min. hour. and day
##
fnc.time<-function(t)
{
  #
  if(t < 60){
    t <- round(t,2)
    tt <- paste(t," - Sec.")
    cat(paste("##\n# Elapsed time:",t,"Sec.\n##\n"))
  } 
  #
  if(t < (60*60) && t >= 60){
    t1 <- as.integer(t/60)
    t <- round(t-t1*60,2) 
    tt <- paste(t1," - Mins.",t," - Sec.")
    cat(paste("##\n# Elapsed time:",t1,"Min.",t,"Sec.\n##\n"))
  }
  #
  if(t < (60*60*24) && t >= (60*60)){
    t2 <- as.integer(t/(60*60))
    t <- t-t2*60*60
    t1 <- as.integer(t/60)
    t <- round(t-t1*60,2) 
    tt <- paste(t2," - Hour/s.",t1," - Mins.",t," - Sec.")
    cat(paste("##\n# Elapsed time:",t2,"Hour/s.",t1,"Min.",t,"Sec.\n##\n"))
  }
  #
  if(t >= (60*60*24)){
    t3 <- as.integer(t/(60*60*24))
    t <- t-t3*60*60*24
    t2 <- as.integer(t/(60*60))
    t <- t-t2*60*60
    t1 <- as.integer(t/60)
    t <- round(t-t1*60,2)
    tt <- paste(t3," - Day/s.",t2," - Hour/s.",t1," - Mins.",t," - Sec.")
    cat(paste("##\n# Elapsed time:",t3,"Day/s.",t2,"Hour/s.",t1,"Mins.",t,"Sec.\n##\n"))
  }
  #
  tt
}
##
## validation criteria
##
spT.validation <- function(z, zhat, names=FALSE)
{
  ##
  ## Validation Mean Squared Error (VMSE) 
  ##
  VMSE <- function(z, zhat)
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- c(z-zhat)
    u <- x[!is.na(x)]
    round(sum(u^2)/length(u), 4)
  }
  ##
  ## Root Mean Squared Error (RMSE) 
  ##
  RMSE<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- c(z-zhat)
    u <- x[!is.na(x)]
    round(sqrt(sum(u^2)/length(u)), 4)
  }
  ##
  ## Mean Absolute Error (MAE) 
  ##
  MAE<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- abs(c(zhat-z))
    u <- x[!is.na(x)]
    round(sum(u)/length(u), 4)
  }
  ##
  ## Mean Absolute Percentage Error (MAPE) 
  ##
  MAPE<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- abs(c(zhat-z))/abs(z)
    u <- x[!is.na(x)]
    u <- u[!is.infinite(u)]
    round(sum(u)/length(u)*100, 4)
  }
  ##
  ## Bias (BIAS) 
  ##
  BIAS<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- c(zhat-z)
    u <- x[!is.na(x)]
    round(sum(u)/length(u), 4)
  }
  ##
  ## Relative Bias (rBIAS) 
  ##
  rBIAS<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- c(zhat-z)
    u <- x[!is.na(x)]
    round(sum(u)/(length(u)*mean(z,na.rm=TRUE)), 4)
  }
  ##
  ## Relative Mean Separation (rMSEP) 
  ##
  rMSEP<-function (z,zhat) 
  {
    z<-as.matrix(z)
    zhat<-as.matrix(zhat)
    x <- c(zhat-z)
    u <- x[!is.na(x)]
    y <- c(mean(zhat,na.rm=TRUE)-z)
    v <- y[!is.na(y)]
    round(sum(u^2)/sum(v^2), 4)
  }
  ##
  if(names==TRUE){
    cat("##\n Mean Squared Error (MSE) \n Root Mean Squared Error (RMSE) \n Mean Absolute Error (MAE) \n Mean Absolute Percentage Error (MAPE) \n Bias (BIAS) \n Relative Bias (rBIAS) \n Relative Mean Separation (rMSEP)\n##\n") 
  }
  ##
  out<-NULL
  out$MSE<-VMSE(c(z), c(zhat))
  out$RMSE<-RMSE(c(z), c(zhat))
  out$MAE<-MAE(c(z), c(zhat))
  out$MAPE<-MAPE(c(z), c(zhat))
  out$BIAS<-BIAS(c(z), c(zhat))
  out$rBIAS<-rBIAS(c(z), c(zhat))
  out$rMSEP<-rMSEP(c(z), c(zhat))
  unlist(out)
  
}
##
## MCMC plots for individual values (trace, density, acf)
##
MCMC.plot.obj<-function(post_val, nItr, nBurn=0, 
                        name=c('....'), ACF="FALSE", PARTIAL.acf="FALSE") 
{
  nBurn=nBurn+1; x<-post_val;
  #if(missing(initial_val)){ 
  #   y<-0
  #}
  #else{
  #   y<-initial_val
  #} 
  if(missing(nItr)){ 
    nItr<-length(post_val)
  }
  if(ACF==FALSE & PARTIAL.acf==FALSE){
    #windows()
    #dev.new()
    par(mfrow=c(1,2))
    #plot(nBurn:nItr, x[nBurn:nItr], xlab = "Iterations", 
    #   ylab = paste("Values of  (", name, ")", sep=''), type = "l",
    #   ylim=c(min(y,x[nBurn:nItr]),max(y,x[nBurn:nItr])),main='Trace plot')
    #abline(h = y, lty = 2,col="blue")
    #plot(density(x[nBurn:nItr]),main='Density plot',xlim=c(min(y,x[nBurn:nItr]),max(y,x[nBurn:nItr])))
    #abline(v = y, lty = 2,col="blue")
    plot(nBurn:nItr, x[nBurn:nItr], xlab = "Iterations", 
         ylab = paste("Values of  (", name, ")", sep=''), type = "l",
         ylim=c(min(x[nBurn:nItr]),max(x[nBurn:nItr])),main='Trace plot')
    plot(density(x[nBurn:nItr]),main='Density plot',xlim=c(min(x[nBurn:nItr]),max(x[nBurn:nItr])))
  }	
  else if(ACF==TRUE & PARTIAL.acf==TRUE){
    #windows()
    #dev.new() 
    par(mfrow=c(1,4))
    plot(nBurn:nItr, x[nBurn:nItr], xlab = "Iterations", 
         ylab = paste("Values of  (", name, ")", sep=''), type = "l", 
         ylim=c(min(x[nBurn:nItr]), max(x[nBurn:nItr])),main='Trace plot')
    #abline(h = y, lty = 2,col="blue")
    plot(density(x[nBurn:nItr]),main='Density plot',xlim=c(min(x[nBurn:nItr]), max(x[nBurn:nItr])))
    #abline(v = y, lty = 2,col="blue")
    acf(x[nBurn:nItr],main='ACF plot',col="red")
    pacf(x[nBurn:nItr],main='Partial ACF plot',col="red")
  }
  else if(ACF=="TRUE" & PARTIAL.acf=="FALSE"){
    #dev.new() 
    par(mfrow=c(1,3))
    plot(nBurn:nItr, x[nBurn:nItr], xlab = "Iterations", 
         ylab = paste("Values of  (", name, ")", sep=''), type = "l", 
         ylim=c(min(x[nBurn:nItr]), max(x[nBurn:nItr])),main='Trace plot')
    #abline(h = y, lty = 2, col="blue")
    plot(density(x[nBurn:nItr]),main='Density plot',xlim=c(min(x[nBurn:nItr]), max(x[nBurn:nItr])))
    #abline(v = y, lty = 2, col="blue")
    acf(x[nBurn:nItr],main='ACF plot',col="red")
  }
  else {
    
  }
}
##
## Multiple imputation for initial values using Amelia with m=1
##
#imputation.z<-function(x){
#
# x is the r*T by n matrix
#
#  x<-amelia(x,m=1)$imputation[[1]]
#  x
#}
##
## To use in coda package
##
as.mcmc.spT<-function(x, ...){
  
  if (x$combined.fit.pred == TRUE) {
    stop("\n# Error: not available for combined fit-prediction option. Please read the MCMC samples from the *.txt file manually.")
  }
  model <- x$model
  if (is.null(model) == TRUE) {
    stop("\n# Error: need to define the model")
  }
  else if (model == "AR" | model == "truncatedAR") {
    r <- x$r
    p <- x$p
    #          
    if(x$cov.fnc=="matern"){
      para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$phip), t(x$nup))
      dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "phi", "nu")
    }
    else {
      para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$phip))
      dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "phi")
    }
    #
    para<-t(para)
    para<-mcmc(para)
    para
  }
  else if (model == "GPP" | model == "truncatedGPP") {
    r <- x$r
    p <- x$p
    #          
    if(x$cov.fnc=="matern"){
      if((!is.null(x$sp.covariate.names))){  
        para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "sig2beta", "phi", "nu")
      }
      else{
        para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "phi", "nu")
      }
    }
    else {
      if((!is.null(x$sp.covariate.names))){  
        para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "sig2beta", "phi")
      }
      else{
        para <- rbind((x$betap), t(x$rhop), t(x$sig2ep), t(x$sig2etap), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]],"rho", "sig2eps", "sig2eta", "phi")
      }
    }
    #
    para<-t(para)
    para<-mcmc(para)
    para
  }
  else if (model == "GP" | model == "truncatedGP") {
    r <- x$r
    p <- x$p
    #          
    if(x$cov.fnc=="matern"){
      if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi", "nu")
      }
      else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi", "nu")
      }
      else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi", "nu")
      } 
      else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi", "nu")
      }
      else {
        stop("Error")
      } 
    }
    else {
      if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi")
      }
      else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi")
      }
      else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi")
      } 
      else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){  
        para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
        dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi")
      }
      else {
        stop("Error")
      } 
    }
    #
    para<-t(para)
    para<-mcmc(para)
    para
  }
  else {
  }
}
##
## use of model.frame
##
model.frame.spT<-function(formula, ...){
  #   tmp<-cbind(formula$Y,formula$X[,-1])
  #   dimnames(tmp)[[2]]<-dimnames(attr(terms(formula$call),"factors"))[[1]]
  #   tmp
  if(formula$combined.fit.pred==TRUE){
    stop("\n# Error: not useful for output with combined fit and predict \n")
  }
  else{
    model.frame(formula$call,formula$data,na.action=na.pass)
  }
}
##
## use of model.matrix
##
model.matrix.spT<-function(object, ...){
  if(object$combined.fit.pred==TRUE){
    stop("\n# Error: not useful for output with combined fit and predict \n")
  }
  else{
    Formula.matrix(object$call,object$data)[[2]]
  }
}
##
## to include dimnames: function used in function summary statistics for spatially varying coef
##
sp.dimname.fnc<-function(x,names,n,q){
  dimnames(x)[[2]][1:(n*q)]<-1:(n*q)
  for(i in 1:q){
    dimnames(x)[[2]][(1+(i-1)*n):(n*i)]<-rep(paste(names[i],"site",1:n,sep=""))
  }
  x
}
##
## to include dimnames: function used in function summary statistics for temporally varying coef
##
tp.dimname.fnc<-function(x,names,u,T){
  dimnames(x)[[2]][1:(u*T)]<-1:(u*T)
  for(i in 1:u){
    dimnames(x)[[2]][(1+(i-1)*T):(T*i)]<-rep(paste(names[i],"time",1:T,sep=""))
  }
  x
}
##
## use of summary
##
summary.spT<-function(object, digits=4, package="spTimer", coefficient=NULL, ...){
  #
  if(package=="coda"){
    if(object$combined.fit.pred==TRUE){
      stop("\n# Error: coda package is not useful for output with combined fit and predict \n")
    }
    else{
      if(is.null(coefficient)){
        cat("\n#### MCMC summary statistics using coda package ####\n")
        if(!is.null(object$sp.covariate.names)) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
        tmp<-as.mcmc(object)
        summary(tmp, ...)
      }
      else if(coefficient=="spatial"){
        if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
        cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
        tmp<-as.mcmc(t(object$betasp))
        if(object$model=="GPP"){n<-object$knots} 
        else{n<-object$n}
        tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q) 
        summary(tmp, ...)
      }
      else if(coefficient=="temporal"){
        if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
        cat("\n#### MCMC summary statistics for the temporally varying dynamic coefficients ####\n")
        tmp<-as.mcmc(t(object$betatp))
        tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T) 
        summary(tmp, ...)
      }
      else if(coefficient=="rho"){
        if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
        cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
        tmp<-as.mcmc(t(object$rhotp))
        dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="") 
        summary(tmp, ...)
      }
      else{
        stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
      }
    }
  }
  else{
    if(package=="spTDyn"){
      coefficient <- coefficient
    }
    #
    else if("Xsp"%in%names(object) | "Xtp"%in%names(object)){
      coefficient <- coefficient   
    }
    else{
      coefficient <- NULL
    }
    #
    if(is.null(coefficient)){
      if((!is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
      else if((!is.null(object$sp.covariate.names)) & (is.null(object$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
      else if((is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
      else { }
      #ans<-NULL
      #ans$Model=object$model;ans$Call=object$call;ans$Iterations=object$iterations;ans$nBurn=object$nBurn;ans$PMCC=object$pmcc;ans$Parameters=round(object$parameter,digits=digits)
      #ans
      print(object)
      cat("-----------------------------------------------------"); cat('\n');
      cat("Parameters:\n")
      print(round(object$parameter,digits=digits)); #cat("\n");
      cat("-----------------------------------------------------"); cat('\n');
    }
    else if(coefficient=="spatial"){
      if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
      cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
      tmp<-as.mcmc(t(object$betasp))
      if(object$model=="GPP"){n<-object$knots} 
      else{n<-object$n}
      tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q) 
      summary(tmp, ...)
    }
    else if(coefficient=="temporal"){
      if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
      cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
      tmp<-as.mcmc(t(object$betatp))
      tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T) 
      summary(tmp, ...)
    }
    else if(coefficient=="rho"){
      if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
      cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
      tmp<-as.mcmc(t(object$rhotp))
      dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="") 
      summary(tmp, ...)
    }
    else{
      stop("Error: the argument coefficient only takes character 'spatial' and 'temporal'.")
    }
  }
}
##
## use of plot
##
plot.spT<-function(x, residuals=FALSE, coefficient=NULL, ...){
  if(as.logical(residuals)==FALSE){
    if(x$combined.fit.pred==TRUE){
      if(!is.null(x$sp.covariate.names)) {cat("## \n# Spatially varying parameters are not included.\n##\n ")}
      cat("\n## Only residual plots are available for output with combined fit and predict option \n\n")
      plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
      par(ask=TRUE)
      qqnorm(residuals(x));qqline(residuals(x),lty=2)
    }
    else{
      if(is.null(coefficient)){
        if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
        else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
        else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
        else {  }
        tmp<-as.mcmc(x)
        plot(tmp, ...)
      }
      else if(coefficient=="spatial"){
        if(is.null(x$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
        tmp<-as.mcmc(t(x$betasp))
        if(x$model=="GPP"){n<-x$knots} 
        else{n<-x$n}
        tmp<-sp.dimname.fnc(x=tmp,names=x$sp.covariate.names,n=n,q=x$q) 
        plot(tmp, ...)
      }
      else if(coefficient=="temporal"){
        if(is.null(x$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
        tmp<-as.mcmc(t(x$betatp))
        tmp<-tp.dimname.fnc(x=tmp,names=x$tp.covariate.names,u=x$u,T=x$T) 
        plot(tmp, ...)
      }
      else if(coefficient=="rho"){
        if(is.null(x$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
        tmp<-as.mcmc(t(x$rhotp))
        dimnames(tmp)[[2]]<-paste("rho",1:x$u,sep="") 
        plot(tmp, ...)
      }
      else{
        stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
      }
    } 
  }
  else{
    if(!is.null(x$sp.covariate.names)) {cat("## \n# Models fitted with spatially varying parameters.\n## \n ")}
    plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
    par(ask=TRUE)
    qqnorm(residuals(x));qqline(residuals(x),lty=2)
  } # end of 2nd if 
}
##
## use of contour
##
#contour.spT<-function(x, surface="Mean", time=c(1), ...){
#	z<-array(fitted(x)[,paste(surface)],dim=c(x$T*x$r,x$n))
#	xyz<-cbind(x$coords,c(z[time,]))
#	xyz<-interp(x=xyz[,1], y=xyz[,2], z=xyz[,3], xo=seq(min(xyz[,1]), max(xyz[,1]), length = 150),
#		 yo=seq(min(xyz[,2]), max(xyz[,2]), length = 150),linear = TRUE, extrap=FALSE, 
#		 duplicate = "error", dupfun = NULL, ncp = NULL)
#	contour(xyz, ...) 
#}
##
## use of coefficients
##
coef.spT<-function(object, digits=4, ...){
  round(t(object$parameter)[1,],digits=digits)
}
##
## use of residuals
##
residuals.spT<-function(object, ...){
  #if(object$combined.fit.pred==TRUE){
  #   stop("\n# Error: not useful for output with combined fit and predict")
  #}
  #else{
  if(object$scale.transform=="NONE"){
    tmp<-object$Y-object$fitted[,1]
    tmp
  }
  else if(object$scale.transform=="SQRT"){
    tmp<-sqrt(object$Y)-object$fitted[,1]
    tmp
  }
  else if(object$scale.transform=="LOG"){
    tmp<-log(object$Y)-object$fitted[,1]
    tmp
  }
  else{
  }
  #}
}
##
## use of formula
##
formula.spT<-function(x, ...){
  x$call
}
##
## use of terms
##
terms.spT<-function(x, ...){
  terms(x$call)
}
##
## For spatially varying beta
##
sp<-function(x){
  class(x)<-"spBeta"
  x
}
##
## For temporally varying beta
##
tp<-function(x){
  class(x)<-"tpBeta"
  x
}
##
## find close locations
##
spT.ObsGridLoc<-function(obsLoc, gridLoc, distance.method="geodetic:km",
                         plot=FALSE)
{
  #
  #
  if(!distance.method %in% c("geodetic:km", "geodetic:mile", "euclidean",
                             "maximum", "manhattan", "canberra")){
    stop("\n# Error: correctly define distance.method \n")
  }
  #
  coords.all <- rbind(obsLoc, gridLoc)
  n1.sites <- 1:dim(obsLoc)[[1]]
  n2.sites <- (dim(obsLoc)[[1]] + 1):(dim(coords.all)[[1]])
  #
  if(distance.method=="geodetic:km"){
    coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
  }
  else if(distance.method=="geodetic:mile"){
    coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
  }
  else{
    coords.D <- as.matrix(dist(coords.all, distance.method, diag = T, upper = T))
  }  
  #
  coords.D[is.na(coords.D)]<-0
  diag(coords.D)<-NA
  # 
  fdmis<-coords.D[n1.sites, n2.sites]
  min.fnc<-function(x){
    x<-cbind(x,1:length(x))
    x<-x[order(x[,1]),] 
    x<-x[1,]
    x  
  }
  fdmis<-t(apply(fdmis,1,min.fnc))
  fdmis<-cbind(obsLoc,gridLoc[fdmis[,2],],fdmis[,2],fdmis[,1])
  dimnames(fdmis)[[2]]<-c("obsLon","obsLat","gridLon","gridLat","gridNum","dist")
  if(plot==TRUE){
    plot(fdmis[,1:2],pch="*")
    points(fdmis[,3:4],pch=19,col=2)
    legend("bottomleft",pch=c(8,19),col=c(1,2),legend=c("Observation locations", "Grid locations"),cex=0.9)
  }
  fdmis
}
##
## combine grid data and coords for spTimer input, gridData should be in array with dim = lon, lat, day, year 
##
spT.gridTodata<-function(gridData, gridLoc=NULL, gridLon=NULL, gridLat=NULL)
{
  # gridData should be in array with dim = lon, lat, day, year
  # for morethan one array use list argument
  # gridLoc should be obtain using expand.grid() argument
  
  options(warn=-1)
  if(!is.null(gridLoc)){
    lonlat<-gridLoc
    dimnames(lonlat)[[2]]<-c("lon","lat")
  }
  else if((!is.null(gridLon)) & (!is.null(gridLat))){ 
    lonlat<-expand.grid(gridLon, gridLat)
    dimnames(lonlat)[[2]]<-c("lon","lat")
  }
  else{
    stop("Error: check grid location input in the function")
  }
  if(is.list(gridData)){
    n<-length(gridData)
    if(n>1){
      dat<-NULL
      for(i in 1:n){
        grid<-c(gridData[[i]])
        dat<-cbind(dat,grid)
      }
      for(i in 1:n){ dimnames(dat)[[2]][i]<-paste("grid",i,sep="") }
      dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
                 lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),dat)
    }
    else{
      dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
                 lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),grid=c(gridData))
    }
  }
  else{
    dat<-cbind(lon=rep(lonlat[,1],length(c(gridData))/dim(lonlat)[[1]]),
               lat=rep(lonlat[,2],length(c(gridData))/dim(lonlat)[[1]]),grid=c(gridData))
  }
  dat<-dat[order(dat[,1],dat[,2]),]
  dat
}
##
## combine observation and grid data for model fitting 
##
spT.ObsGridData<-function(obsData, gridData, obsLoc, gridLoc, distance.method="geodetic:km")
{
  # obsData should be in data form
  # gridData should be in array form with dim = lon, lat, day, year
  
  obsLoc<-as.matrix(obsLoc)
  gridLoc<-as.matrix(gridLoc)
  dimnames(obsLoc)<-NULL
  dimnames(gridLoc)<-NULL
  
  if((is.list(gridData)) & (is.list(gridLoc))){
    n1<-length(gridData)
    n2<-length(gridLoc)
    if(n1 != n2){ stop("\n Error: list of gridData and gridLoc should be same") }
    datt<-NULL
    for(i in 1:n1){
      gLoc<-gridLoc[[i]]
      gLoc<-gLoc[order(gLoc[,1],gLoc[,2]),]
      comb<-spT.ObsGridLoc(obsLoc=obsLoc, gridLoc=gLoc, distance.method=distance.method)
      grid<-spT.gridTodata(gridData=gridData[[i]], gridLoc=gLoc)
      grid<-cbind(gridNum=sort(rep(1:dim(gLoc)[[1]],dim(grid)[[1]]/dim(gLoc)[[1]])),grid)
      dat<-NULL
      for(j in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==j,]) }
      dimnames(dat)[[2]]<-c(paste("gridNum",i,sep=""),paste("grid.lon",i,sep=""),paste("grid.lat",i,sep=""),paste("grid",i,sep=""))
      if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n obsData and gridData time points mismatch.\n")}
      datt<-cbind(datt,dat)
    }
    dat<-cbind(obsData,datt); rm(datt);
    dat
  }
  else{
    comb<-spT.ObsGridLoc(obsLoc=obsLoc, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),], distance.method=distance.method)
    grid<-spT.gridTodata(gridData=gridData, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),])
    grid<-cbind(gridNum=sort(rep(1:dim(gridLoc)[[1]],dim(grid)[[1]]/dim(gridLoc)[[1]])),grid)
    dat<-NULL
    for(i in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==i,]) }
    dimnames(dat)[[2]][2:3]<-c("grid.lon","grid.lat")
    if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n   obsData and gridData time points mismatch.\n")}
    dat<-cbind(obsData,dat)
    dat
  }
}
##
## Print function
##
print.spT<-function(x, ...) {
  cat("-----------------------------------------------------"); cat('\n');
  cat("Model: "); cat(x$model); cat('\n');
  cat("Call: "); print(x$call); #cat('\n')
  cat("Iterations: "); cat(x$iterations); cat("\n")
  cat("nBurn: "); cat(x$nBurn); cat("\n")
  cat("Acceptance rate for phi (%): "); cat(x$accept); cat("\n")
  cat("-----------------------------------------------------"); cat('\n');
  #cat("PMCC: "); cat("\n");
  print(x$PMCC); 
  cat("-----------------------------------------------------"); cat('\n');
  #cat("Parameters:\n")
  #print(x$parameter); #cat("\n");
  #cat("-----------------------------------------------------"); cat('\n');
  cat("Computation time: "); cat(x$computation.time); cat("\n")
}
##
## fitted function
##
fitted.spT<-function(object, ...){
  x<-data.frame(object$fitted)
  x
}
##
## use of package forecast
##
as.forecast.object<-function(object, site=1, level=c(80,95), ...){
  # object is the output from the predict.spT for forecast
  x<-NULL
  x$model<-list(Name = object$model)
  x$method<-paste(object$model,"spatio-temporal model, site:",site)
  x$x<-ts(object$obsData)[,site]
  x$fitted<-ts(object$fittedData)[,site]
  x$residuals<-ts(object$residuals)[,site]
  x$mean<-ts(object$Mean,start=c(dim(object$obsData)[[1]]+1,1))[,site]
  x$level<-level
  x$upper<-array(apply(object$fore.samples,1,quantile,probs=c(level/100)),dim=c(length(level),dim(object$Mean)[[1]],dim(object$Mean)[[2]]))
  x$upper<-t(x$upper[,,site])
  x$lower<-array(apply(object$fore.samples,1,quantile,probs=c(1-level/100)),dim=c(length(level),dim(object$Mean)[[1]],dim(object$Mean)[[2]]))
  x$lower<-t(x$lower[,,site])
  cat(paste("Forecast for site:",site),"\n")
  class(x)<-"forecast"
  x
}
##
## use of confint()
##
confint.spT<-function(object, parm, level=0.95, ...){
  #
  x<-as.mcmc(object)
  up<-level+(1-level)/2
  low<-(1-level)/2
  FUN <- function(x){quantile(x,probs=c(low,up))}
  out<-apply(x,2,FUN=FUN)
  out<-t(out)
  if(missing(parm)){
    out
  }
  else{
    if(length(parm)>1){ 
      out<-as.matrix(out[dimnames(out)[[1]] %in% parm,])
      out
    }
    else{
      out<-t(as.matrix(out[dimnames(out)[[1]] %in% parm,]))
      dimnames(out)[[1]]<-parm
      out
    }
  }
}
##
## Gamma prior
##
Gamm<-function(a=NA,b=NA){
  out<-matrix(c(a,b),1,2)
  class(out)<-"Gamma"
  out
}
##
## Normal prior
##
Norm<-function(mu=NA,sig=NA){
  out<-matrix(c(mu,sig),1,2)
  class(out)<-"Normal"
  out
}
##
## Uniform prior
##
Unif<-function(low=NA,up=NA){
  out<-matrix(c(low,up),1,2)
  class(out)<-"Uniform"
  out
}
##
## Truncated Gaussian code
##
truncated.fnc<-function(Y, at=0, lambda=NULL, both=FALSE){
  #
  #
  # at is left-tailed
  #
  if(is.null(lambda)){
    stop("Error: define truncation parameter lambda properly using list ")
  }
  if(is.null(at)){
    stop("Error: define truncation point properly using list ")
  }
  if(at < 0){
    stop("Error: currently truncation point only can take value >= zero ")
  }
  zm <- cbind(Y,Y-at,1)
  zm[zm[, 2] <= 0, 3] <- 0
  zm[zm[, 3] == 0, 2] <- - rhnorm(nrow(zm[zm[,3]==0,]),sd(zm[zm[, 3] != 0, 1],na.rm=TRUE))
  ck <- min(zm[,2],na.rm=TRUE)
  zm[,2] <- zm[,2]-ck
  zm[,2] <- zm[,2]^(1/lambda)
  #zm[, 2] <- zm[, 2]^(1/lambda)
  #zm[zm[, 3] == 0, 2] <- - rhnorm(nrow(zm[zm[,3]==0,]),sd(zm[zm[, 3] != 0, 1]^(1/lambda)))
  #
  if (both == TRUE) {
    list(zm=zm,at2=ck)
  }
  else {
    list(zm=c(zm[,2]),at2=ck)
  }
  #
}
##
## for truncated model
##
reverse.truncated.fnc<-function(Y, at=0, lambda=NULL, at2=0){
  #
  # at is left-tailed
  #
  if(is.null(lambda)){
    stop("Error: define truncation parameter lambda properly using list ")
  }
  #zm <- Y
  #zm[zm <= 0]<- 0
  #zm <- zm^(lambda)
  #zm <- zm + at
  zm <- Y
  zm <- zm^(lambda)
  zm <- zm + at2
  zm[zm <= 0]<- 0
  zm <- zm + at
  zm
  #
}
##
## probability below threshold (for truncated)
##
prob.below.threshold <- function(out, at){
  # out is the N x nItr MCMC samples
  fnc<-function(x, at){
    length(x[x <= at])/length(x)
  }   
  apply(out,1,fnc,at=at)
}
##
##
##