R/predict.sptmRes.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' Prediction for SPTMres objects
#'
#' Update the seasons for LU
#' @param SPTMresobj result of SpTMixture function.
#' @keywords updating season
#' @export
#' @examples
#' data("wq_analysis_week2")
#' SPTMData(wq.raw.obs, frequency = "quarter")
predict.sptmRes = function(SPTMresobj, transform = "none"){

  if(transform == "none"){
    fT = function(x){x}
  }
  if(transform == "log"){
    fT = function(x){exp(x)}
  }

  theta  = SPTMresobj$Result$theta
  p.sites = SPTMresobj$Result$latent

  n.tempFun = 1 * as.numeric(length(dim(SPTMresobj$data$data$season))==2) +
    (1-as.numeric(length(dim(SPTMresobj$data$data$season))==2)) * dim(SPTMresobj$data$data$season)[2]

  if(length(theta) < (ncol(p.sites) * (2 + n.tempFun))){theta = c(theta, rep(theta[length(theta)], ncol(p.sites)-1))}



  season.sites = SPTMresobj$seasonSites


  if(!is.null(SPTMresobj$data$data$basis)){
    y.p = meanFuncSplines(theta, n.mixt=ncol(p.sites), basis = SPTMresobj$data$data$basis)
  }else{
    y.p = meanFunc(theta, n.mixt=ncol(p.sites), p.sites=p.sites, season.sites = season.sites)
  }


  idx.sigma = (length(theta) - ncol(p.sites) + 1):length(theta)

  if(!is.null(SPTMresobj$GibbsOut)){
    theta[idx.sigma] = log(sqrt(theta[idx.sigma]))
  }

  df = SPTMresobj$data$data

  SiteID = unique(df$obs.data$ID)

  n.sites = length(SiteID)

  n.mixt = ncol(p.sites)

  dfPred = NULL

  ListCl = apply(p.sites, 1,which.max)

  rangeCl = matrix(NA, nrow = ncol(p.sites), ncol = 2)
  for(i in sort(unique(ListCl))){
    Y = df$obs.data[which(df$obs.data$ID %in% SiteID[which(ListCl == i)]),]
    rangeCl[i,] = as.character(range(Y$date[!is.na(Y$obs)]))
  }

  for(i in 1:length(SiteID)){

    obs = df$obs.data$obs[df$obs.data$ID %in% SiteID[i]]
    idx.k = which.max(p.sites[i,])
    mat = matrix(p.sites[i,], ncol=1)
    idx.na = which(is.na(colSums(y.p)))
    pred= y.p %*% mat
    if(length(idx.na)>0){pred= y.p[,-idx.na] %*% mat[-idx.na]}
    dates = df$obs.data$date[df$obs.data$ID %in% SiteID[i]]
    ci = cbind(pred - 2*exp(theta[idx.sigma][idx.k]),pred + 2*exp(theta[idx.sigma][idx.k]))

    dfTemp = data.frame(obs = fT(obs))
    dfTemp$date = dates
    dfTemp$ID = SiteID[i]
    dfTemp$pred = fT(pred)
    dfTemp$CIlow = fT(ci[,1])
    dfTemp$CIup = fT(ci[,2])
    dfTemp$Cluster = which.max(p.sites[i,])

    rangeDates = which(dfTemp$date < rangeCl[idx.k,1] | dfTemp$date > rangeCl[idx.k,2])
    if(length(rangeDates)>0){
      dfTemp[rangeDates,c("pred","CIlow", "CIup")] <- NA
    }


    dfPred = rbind(dfPred, dfTemp)

  }

  dfPred$Cluster = as.factor(dfPred$Cluster)

  return(dfPred)

}