R/estimGibbs.R

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

#' Wrapper function for Gibbs sampler
#'
#' Wrapper function for Gibbs sampler
#' @param df.sptmod model data
#' @param method method to be used
#' @param theta.init prior mean of the parameters
#' @param mixture boolean. Is the model a mixture model.
#' @param N.run Number of runs of the sampler.
#' @param beta.range Range of beta for Potts.
#' @param print.res logical. printing or not the steps values.
#' @keywords estimation
#' @export
#' @examples
#' data("wq_analysis_week2")
#' SPTMData(wq.raw.obs, frequency = "quarter")
estimGibbs = function(df.sptmod, method = 'simple', priors = NULL, mixture = TRUE, N.run = 10000, beta.range = c(0.5, 3), print.res = FALSE, debug=FALSE, meanFunc = NULL, Xcov = NULL){

  if(debug){browser()}


  p.comp = df.sptmod$lu

  season = df.sptmod$data$season

  basis = df.sptmod$data$basis

  season.sites=season

  season.lu = updateSeason(season.sites,p.sites= p.comp)

  y = df.sptmod$data$obs.data

  if(method == 'simple'){

    theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
                       sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
    theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
    #theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)

    if(is.null(priors)){
      priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = abs(theta.init[1:(ncol(season.lu))]) / 10),
                    sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
    }else{priors = priors}

    res = gibbSampler(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
                      ID = levels(y$ID), N.run = N.run, mixture = mixture, meanFunc = meanFunc, debug, print.res = print.res)
    p.comp = apply(res$z, 1:2, mean)
    RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
               seasonSites = season.sites,
               GibbsOut = res,
               data = df.sptmod)

  }


  if(method == 'spat'){

    dd = deldir(df.sptmod$coord[,1], df.sptmod$coord[,2],suppressMsge=TRUE)
    tile.dd = tile.list(dd)

    Bds = comPts(tile.dd)

    theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
                       sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
    theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
    theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)

    if(is.null(priors)){
      priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = theta.init[1:(ncol(season.lu))] / 10),
                    sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
    }else{priors = priors}

    res = gibbSPTM(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
                   Coords=NULL ,Bds=Bds, ID = levels(y$ID), N.run = N.run, beta.range = beta.range,
                   meanFunc = meanFunc, debug, print.res = print.res)
    p.comp = apply(res$z, 1:2, mean)
    RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
               seasonSites = season.sites,
               GibbsOut = res,
               data = df.sptmod)
  }

  if(method == 'mixedeffect'){

    dd = deldir(df.sptmod$coord[,1], df.sptmod$coord[,2],suppressMsge=TRUE)
    tile.dd = tile.list(dd)

    Bds = comPts(tile.dd)

    theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
                       sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
    theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
    theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)

    if(is.null(priors)){
      priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = theta.init[1:(ncol(season.lu))] / 10),
                    sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
    }else{priors = priors}

    attr(basis, "indexSeason") <- df.sptmod$data$list.idx

    res = gibbSPTMixed(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
                   Coords=NULL ,Bds=Bds, ID = levels(y$ID), N.run = N.run, beta.range = beta.range,
                   meanFunc = meanFunc, debug, print.res = print.res, Xcov = df.sptmod$covariates)
    p.comp = apply(res$z, 1:2, mean)
    RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
               seasonSites = season.sites,
               GibbsOut = res,
               data = df.sptmod)
  }


  class(RET) <- 'sptmRes'
  return(RET)
}