R/BayesX.R

Defines functions transformBayesX controlBayesX setupBayesX samplerBayesX process_mfile resultsBayesX sx sx.construct sx.construct.default do.xt sx.construct.pspline.smooth.spec sx.construct.random.smooth.spec sx.construct.rps.smooth.spec sx.construct.kriging.smooth.spec construct.shrw sx.construct.offset.smooth.spec sx.construct.spatial.smooth.spec make_by rmf help.map.name sfoofun splitme resplit

####################################
## (1) BayesX specific functions. ##
####################################
transformBayesX <- function(x, ...)
{
  family <- bayesr.family(attr(x, "family"))
  call <- x$call; x$call <- NULL

  if(family$cat) {
    ylevels <- attr(x, "ylevels")
    reference <- attr(x, "reference")
    grid <- attr(x, "grid")
    rn <- attr(attr(x, "model.frame"), "response.name")
    if(is.null(ylevels)) {
      ylevels <- rn
      reference <- "NA"
      names(x) <- paste(names(x)[1], ylevels, sep = ":")
      for(j in seq_along(x))
        x[[j]]$cat.formula <- x[[j]]$formula
    } else {
      f <- as.formula(paste("~ -1 +", rn))
      rm <- as.data.frame(model.matrix(f, data = attr(x, "model.frame")))
      cn <- colnames(rm) <- rmf(colnames(rm))
      rm <- rm[, !grepl(reference, cn), drop = FALSE]
      mf <- cbind(attr(x, "model.frame"), rm)
      attr(mf, "response.name") <- rn
      x <- x[ylevels != reference]
      attr(x, "model.frame") <- mf
      attr(x, "grid") <- grid
    }
    attr(x, "ylevels") <- ylevels
    attr(x, "reference") <- reference
  }

  attr(attr(x, "model.frame"), "orig.names") <- names(attr(x, "model.frame"))
  names(attr(x, "model.frame")) <- rmf(names(attr(x, "model.frame")))
  attr(attr(x, "model.frame"), "response.name") <- rmf(attr(attr(x, "model.frame"), "response.name"))

  tBayesX <- function(obj, ...) {
    if(!any(c("formula", "fake.formula", "response") %in% names(obj))) {
      nx <- names(obj)
      nx <- nx[nx != "call"]
      if(is.null(nx)) nx <- 1:length(obj)
      if(length(unique(nx)) < length(obj)) nx <- 1:length(obj)
      for(j in nx)
        obj[[j]] <- tBayesX(obj[[j]], ...)
    } else {
      obj <- randomize(obj)
      obj$response <- rmf(obj$response)
      if(!is.null(dim(obj$X)))
        colnames(obj$X) <- rmf(colnames(obj$X))
      if(length(obj$pterms))
        obj$pterms <- rmf(obj$pterms)
      if(length(obj$sterms))
        obj$sterms <- rmf(obj$sterms)

      if(length(obj$smooth)) stop("arbitrary smooths not supported yet, please use sx()!")
      if(length(obj$sx.smooth)) {
        for(j in seq_along(obj$sx.smooth)) {
          obj$sx.smooth[[j]]$term <- rmf(obj$sx.smooth[[j]]$term)
          tmp <- sx.construct(obj$sx.smooth[[j]],
            if(is.null(attr(obj, "model.frame"))) {
              attr(x, "model.frame")
            } else attr(obj, "model.frame"))
          attr(tmp, "specs") <- obj$sx.smooth[[j]]
          obj$sx.smooth[[j]] <- tmp
        }
      }
    }
    obj
  }

  x <- tBayesX(x, ...)

  if(inherits(x, "bayesr.input") & !any(c("formula", "fake.formula", "response") %in% names(x))) {
    n <- length(x)
  } else {
    x <- list(x)
    n <- 1
  }

  names(x) <- rep(family$names, length.out = n)
  if(!is.null(family$bayesx$order)) {
    mf <- attr(x, "model.frame")
    class(x) <- "list"
    x <- x[rev(family$bayesx$order)]
    class(x) <- c("bayesr.input", "list")
    attr(x, "model.frame") <- mf; rm(mf)
  }
  attr(x, "call") <- call
  attr(x, "family") <- family

  x
}


controlBayesX <- function(n.iter = 1200, thin = 1, burnin = 200,
  seed = NULL, predict = "light", model.name = "bayesr", data.name = "d",
  prg.name = NULL, dir = NULL, verbose = FALSE, cores = NULL, ...)
{
  if(is.null(seed))
    seed <- '##seed##'
  stopifnot(burnin < n.iter)
  if(is.null(model.name))
    model.name <- 'bayesr'
  if(is.null(data.name))
    data.name <- 'd'
  if(is.null(prg.name))
    prg.name <- paste(model.name, 'prg', sep = '.')
  if(!grepl(".prg", prg.name))
    prg.name <- paste(prg.name, "prg", sep = ".")
  if(is.null(dir)) {
    dir.create(dir <- tempfile())
    attr(dir, "unlink") <- TRUE
  } else dir <- path.expand(dir)
  if(!file.exists(dir)) dir.create(dir)
  if(is.null(cores)) cores <- 1

  cvals <- list(
    "prg" = list(
      "iterations" = n.iter, "burnin" = burnin, "step" = thin,
      "setseed" = seed, "predict" = predict
    ),
    "setup" = list(
      "main" = c(rep(FALSE, 3), rep(TRUE, 2)), "model.name" = model.name, "data.name" = data.name,
      "prg.name" = prg.name, "dir" = dir, "verbose" = verbose, "cores" = cores
    )
  )

  cvals
}

setupBayesX <- function(x, control = controlBayesX(...), ...)
{
  names(attr(x, "model.frame")) <- rmf(names(attr(x, "model.frame")))

  mf.names <- names(attr(x, "model.frame"))

  family <- attr(x, "family")
  x$call <- x$family <- NULL
  lhs <- family$bayesx$lhs

  ## Handling of weights.
  add.weights <- FALSE
  if(!is.null(family$bayesx$weights)) {
    weights0 <- attr(x, "model.frame")$weights
    rn <- attributes(attr(x, "model.frame"))$response.name
    for(j in names(family$bayesx$weights)) {
      weights1 <- family$bayesx$weights[[j]](attr(x, "model.frame")[[rn[1]]])
      if(!is.null(weights0)) weights1 <- weights1 * weights0
      attr(x, "model.frame")[[paste(j, "weights", sep = "")]] <- weights1
      rm(weights1)
    }
	  rm(weights0)
    mf.names <- names(attr(x, "model.frame"))
    add.weights <- TRUE
  }
  if(length(grep("weights", mf.names)))
    add.weights <- TRUE
  zero <- family$bayesx$zero
  if(is.null(zero)) zero <- FALSE
  if(zero) {
    attr(x, "model.frame")[["ybinom"]] <- 1 * (attr(x, "model.frame")[[rn]] > 0)
  }
  quantile <- family$bayesx$quantile

  dirichlet <- family$family == "dirichlet"
  family$bayesx[c("order", "lhs", "rm.number", "weights", "quantile", "zero")] <- NULL

  args <- list(...)
  model.name <- control$setup$model.name
  data.name <- rmf(control$setup$data.name)
  prg.name <- control$setup$prg.name
  dir <- control$setup$dir
  cores <- control$setup$cores
  if(is.null(cores)) cores <- 1

  BayesX_data <- function(obj, h = FALSE, id = NULL) {
    if(!is.null(obj$cat.formula) & !h) {
      obj$response <- formula_respname(obj$cat.formula)
      obj$response.vec <- attr(x, "model.frame")[[obj$response]]
    }
    if(zero & !is.null(attr(obj$formula, "name"))) {
      if(attr(obj$formula, "name") == "pi" & !h) {
        obj$response <- "ybinom"
        obj$response.vec <- attr(x, "model.frame")[["ybinom"]]
      }
    }
    if(h)
      obj$response.vec <- attr(x, "model.frame")[[obj$response]]
    X <- NULL
    if(!is.null(obj$X)) {
      if(ncol(obj$X) > 0) {
        X <- obj$X
        if(length(i <- grep("Intercept", colnames(X), fixed = TRUE)))
          X <- obj$X[, -i, drop = FALSE]
        X <- if(ncol(X) > 0) as.data.frame(X) else NULL
      }
    }
    if(k2 <- length(obj$sterms)) {
      X <- if(!is.null(X)) {
        cbind(X, attr(x, "model.frame")[, obj$sterms, drop = FALSE])
      } else {
        attr(x, "model.frame")[, obj$sterms, drop = FALSE]
      }
      k1 <- ncol(X)
      colnames(X)[(k1 - k2 + 1):k1] <- obj$sterms
      for(sj in obj$sterms) {
        if(is.factor(X[, sj])) {
          f <- as.formula(paste("~ -1 +", sj))
          fX <- as.data.frame(model.matrix(f, data = X))
          X <- cbind(X, fX)
          X <- X[, unique(colnames(X))]
        }
      }
    }
    if(!is.null(X)) {
      cnX <- colnames(X)
      if(ncol(X) > 0)
        X <- X[, unique(cnX[!grepl("Intercept", cnX, fixed = TRUE)]), drop = FALSE]
      if(!is.null(obj$response)) {
        X[[obj$response]] <- obj$response.vec
      }
    } else {
      X <- list()
      if(!is.null(obj$response))
        X[[obj$response]] <- obj$response.vec
      X <- as.data.frame(X)
    }
    if(ncol(X) < 1) {
      X <- NULL
    } else {
      yf <- if(is.null(obj$response)) FALSE else is.factor(X[[obj$response]])
      X <- as.data.frame(X)
      for(j in 1:ncol(X)) {
        if(is.factor(X[, j])) {
          warn <- getOption("warn")
          options(warn = -1)
          tx <- as.integer(as.character(X[, j]))
          options("warn" = warn)
          X[, j] <- if(!any(is.na(tx))) tx else as.integer(X[, j])
          X <- X[order(X[, j]), , drop = FALSE]
        }
      }
      if(yf) {
        if(!h) {
          if(!is.null(obj$response))
            X[[obj$response]] <- as.integer(as.factor(X[[obj$response]])) - 1
        }
        X <- X[order(X[[obj$response]]), , drop = FALSE]
      }
    }
    if(!h) {
      wi <- paste(id, "weights", sep = "")
      if(!is.null(id) & (wi %in% mf.names)) {
        X[[wi]] <- attr(x, "model.frame")[[wi]]
      } else {
        if("weights" %in% mf.names)
          X$ModelWeights <- attr(x, "model.frame")[["weights"]]
      }
      if("offset" %in% mf.names)
        X$ModelOffset <- attr(x, "model.frame")[["offset"]]
    }
    if(h) {
      X <- unique(X)
      if(!is.null(obj$response))
        X <- X[order(X[[obj$response]]), , drop = FALSE]
    }

    return(X)
  }

  nx <- names(x)
  n <- length(x)

  count <- 1
  d <- prg <- NULL
  dpath0 <- file.path(dir, paste(dname0 <- paste(data.name, 0, sep = ''), "raw", sep = '.'))
  for(j in n:1) {
    if(!"fake.formula" %in% names(x[[j]])) {
      for(i in 1:length(x[[j]])) {
        d2 <- NULL
        if(i < 2) {
          d2 <- BayesX_data(x[[j]][[i]], id = nx[j])
          if(is.null(d)) {
            d <- d2
          } else {
            if(!is.null(d2))
              d <- cbind(d, d2)
          }
          d2 <- NULL
          d <- as.data.frame(d)
          d <- d[, unique(names(d)), drop = FALSE]
          x[[j]][[i]]$dname <- dname0
        } else d2 <- BayesX_data(x[[j]][[i]], i > 1)
        if(is.null(x[[j]][[i]]$hlevel))
          x[[j]][[i]]$hlevel <- i
        if(!is.null(d2)) {
          dpath <- file.path(dir, paste(dname <- paste(data.name, count, sep = ''),
            "raw", sep = '.'))
          x[[j]][[i]]$dname <- dname
          write.table(d2, file = dpath, quote = FALSE, row.names = FALSE, col.names = TRUE)
          prg <- c(prg,
            paste('dataset ', data.name, count, sep = ''),
            paste(data.name, count, '.infile using ', dpath, sep = '')
          )
          count <- count + 1
        }
      }
    } else {
      d2 <- BayesX_data(x[[j]], id = nx[j])
      if(is.null(d)) {
        d <- d2
      } else {
        if(!is.null(d2))
          d <- cbind(d, BayesX_data(x[[j]], id = nx[j]))
      }
      d2 <- NULL
      if(!is.null(d)) {
        d <- as.data.frame(d)
        d <- d[, unique(names(d)), drop = FALSE]
      }
      x[[j]]$dname <- dname0
      x[[j]]$hlevel <- 1
    }
  }
  response.name <- attr(attr(x, "model.frame"), "response.name")
  write.table(d, file = dpath0, quote = FALSE, row.names = FALSE, col.names = TRUE)
  prg <- c(
    paste('dataset', dname0),
    paste(dname0, '.infile using ', dpath0, sep = ''),
	prg
	)

  prg <- c(paste('logopen using', file.path(dir, paste(model.name, 'log', sep = '.'))), "", prg) 
  prg <- c(prg, paste('\nmcmcreg', model.name))
  prg <- c(prg, paste(model.name, '.outfile = ',
    if(cores < 2) file.path(dir, model.name) else '##outfile##',
    sep = ''))

  make_eqn <- function(x, ctr = TRUE, id = NULL) {
    n <- length(x)
    eqn <- NULL
    for(j in n:1) {
      if(!"fake.formula" %in% names(x[[j]])) {
        eqn <- c(eqn, make_eqn(x[[j]], ctr, id = j))
        ctr <- FALSE
      } else {
        teqn <- paste(model.name, '.hregress ',
          if(family$cat | !is.null(lhs)) {
            if(x[[j]]$hlevel > 1) {
              x[[j]]$response
            } else {
              if(!is.null(lhs)) {
                if(nx[if(is.null(id)) j else id] %in% names(lhs)) {
                  lhs[nx[if(is.null(id)) j else id]]
                } else {
                  if(is.null(x[[j]]$response)) response.name[1] else x[[j]]$response
                }
              } else formula_respname(x[[j]]$cat.formula)
            }
          } else {
            if(zero & x[[j]]$hlevel < 2 & !is.null(x[[j]]$formula)) {
              if(attr(x[[j]]$formula, "name") == "pi") {
                "ybinom"
              } else {
                if(is.null(x[[j]]$response)) response.name[1] else x[[j]]$response
              }
            } else {
              if(is.null(x[[j]]$response)) response.name[1] else x[[j]]$response
            }
          }, sep = '')
        et <- x[[j]]$pterms
        fctr <- attr(x[[j]]$formula, "control")
        if(is.null(fctr)) fctr <- ""
        if(length(i <- grep("Intercept", et, fixed = TRUE)))
          et[i] <- "const"
        if(length(x[[j]]$sx.smooth)) {
          et <- c(et, unlist(x[[j]]$sx.smooth))
        }
        if(x[[j]]$hlevel < 2) {
          if("offset" %in% mf.names)
            et <- c(et, "ModelOffset(offset)")
        }
        if(length(et))
          teqn <- paste(teqn, '=', paste(et, collapse = ' + '))
        if(x[[j]]$hlevel < 2 & add.weights) {
          wi <- paste(nx[if(is.null(id)) j else id], "weights", sep = "")
          if(wi %in% mf.names) {
            teqn <- paste(teqn, "weight", wi)
          } else {
            if("weights" %in% mf.names)
              teqn <- paste(teqn, "weight ModelWeights")
          }
        }
        if(ctr) {
          ok <- control$setup$main
          c2 <- control$prg[!ok]
          ok <- sapply(names(c2), function(x) { !any(grepl(x, fctr)) })
          c2 <- c2[ok]
          if(length(c2)) {
            teqn <- paste(teqn, ", ", paste(names(c2), "=", unlist(c2),
              sep = "", collapse = " "), sep = "")
          }
          ctr <- FALSE
        } else teqn <- paste(teqn, ",", sep = "")
        ok <- TRUE
        eqtj <- if(family$cat & (if(is.null(id)) j else id) != 1) 3 else 2
        if(!is.null(x[[j]]$hlevel)) {
          if(!any(grepl("hlevel", fctr))) {
            teqn <- paste(teqn, " hlevel=", if(x[[j]]$hlevel > 1) 2 else 1, sep = "")
            if(x[[j]]$hlevel > 1) {
              if(!any(grepl("family", fctr)))
                teqn <- paste(teqn, " family=gaussian_re", sep = "")
              if(!any(grepl("family", fctr))) {
                teqn <- paste(teqn, " equationtype=", family$bayesx[[nx[if(is.null(id)) j else id]]][eqtj], sep = "")
                ok <- FALSE
              }
            }
            if(x[[j]]$hlevel < 2 & any(grepl("mean", family$bayesx[[nx[if(is.null(id)) j else id]]]))) {
              if(any(grepl("predict", names(control$prg))) & !zero) {
                teqn <- paste(teqn, " predict=", control$prg$predict, " setseed=", control$prg$setseed, sep = "")
              }
            }
          } else ok <- grepl("1", fctr[grepl("hlevel", fctr)])
        }
        if(ok) {
          if(dirichlet) {
            teqn <- paste(teqn, " nrcat=", family$ncat, " family=dirichlet equationtype=",
              if(j < 2) "mean" else paste("alpha", j, sep = ""), sep = "")
          } else {
            if(!any(grepl("family", fctr)))
              teqn <- paste(teqn, " family=", family$bayesx[[nx[if(is.null(id)) j else id]]][1], sep = "")
            if(!any(grepl("equationtype", fctr)))
					    teqn <- paste(teqn, " equationtype=", family$bayesx[[nx[if(is.null(id)) j else id]]][eqtj], sep = "")
            if(!is.null(quantile))
              teqn <- paste(teqn, " quantile=", quantile, sep = "")
          }
        }
        if(length(fctr)) {
          fctr <- paste(fctr, collapse = " ")
          if(nchar(fctr) > 0)
            teqn <- paste(teqn, fctr)
        }
        if(!is.null(x[[j]]$dname)) {
          if(!any(grepl("using", fctr)))
            teqn <- paste(teqn, " using ", x[[j]]$dname, sep = "")
        }
        eqn <- c(eqn, "", teqn)
      }
    }
    eqn
  }

  prg_extras <- function(x) {
    n <- length(x)
    prgex <- NULL
    for(j in n:1) {
      if(!"fake.formula" %in% names(x[[j]])) {
        prgex <- c(prgex, prg_extras(x[[j]]))
      } else {
        if(length(x[[j]]$sx.smooth)) {
          for(k in x[[j]]$sx.smooth) {
            if(!is.null(attr(k, "write"))) {
              if(!is.null(attr(k, "map.name"))) {
                if(!any(grepl(paste("map", attr(k, "map.name")), prgex)))
                  prgex <- c(prgex, attr(k, "write")(dir))
              } else prgex <- c(prgex, attr(k, "write")(dir))
            }
          }
        }
      }
    }
    if(!is.null(prgex)) prgex <- c(prgex, "")
    prgex
  }
  
  prg <- c(prg_extras(x), prg, make_eqn(x))

  if(zero) {
    prg <- c(prg, "", paste(model.name, ".hregress ybinom, family=zero_adjusted setseed=",
      control$prg$setseed," predict=light using ", x[[1]]$dname, sep = "" ))
  }

  prg <- gsub("(random", "(hrandom", prg, fixed = TRUE)
  for(i in 1:5)
    prg <- gsub(paste("(psplinerw", i, sep = ""), "(pspline", prg, fixed = TRUE)
  prg <- c(prg, "", paste(model.name, "getsample", sep = "."))
  prg <- c(
    paste('%% BayesX program created by BayesR: ', as.character(Sys.time()), sep = ''),
    paste('%% usefile ', file.path(dir, prg.name), sep = ''), "",
    prg
  )

  return(list("prg" = prg, "control" = control))
}


samplerBayesX <- function(x, ...)
{
  if(any(grepl("##seed##", x$prg, fixed = TRUE)))
    x$prg <- gsub("##seed##", round(runif(1L) * .Machine$integer.max), x$prg, fixed = TRUE)

  dir <- x$control$setup$dir
  if(!is.null(attr(dir, "unlink"))) {
    if(attr(dir, "unlink"))
      on.exit(unlink(dir))
  }
  prg.name <- x$control$setup$prg.name
  model.name <- x$control$setup$model.name
  cores <- x$control$setup$cores

  if(cores > 1) {
    k <- 1
    while(file.exists(cdir <- file.path(dir, paste("Core", k, sep = "")))) {
      k <- k + 1
    }
    dir.create(cdir)
    x$prg <- gsub('##outfile##', prg <- file.path(cdir, model.name), x$prg, fixed = TRUE)
    prg.name <- paste(gsub(".prg", "", prg.name, fixed = TRUE), k, ".prg", sep = "")
    if(length(i <- grep('usefile', x$prg)))
      x$prg[i] <- paste('%% usefile', file.path(dir, prg.name))
    x$prg <- gsub(file.path(dir, paste(model.name, 'log', sep = '.')),
      file.path(cdir, paste(model.name, 'log', sep = '.')), x$prg, fixed = TRUE)
    x$prg <- gsub(paste(model.name, ".", sep = ""), paste(model.name, k, ".", sep = ""),
      x$prg, fixed = TRUE)
    x$prg[grep("mcmcreg", x$prg)] <- paste("mcmcreg ", model.name, k, sep = "")

    dir <- cdir
  }

  prg <- file.path(dir, prg.name)
writeLines(x$prg)
  writeLines(x$prg, file.path(dir, prg.name))

  warn <- getOption("warn")
  options(warn = -1)
  ok <- BayesXsrc::run.bayesx(prg = prg, verbose = x$control$setup$verbose, ...)
  options("warn" = warn)
  if(length(i <- grep("error", ok$log, ignore.case = TRUE))) {
    errl <- gsub("^ +", "", ok$log[i])
    errl <- gsub(" +$", "", errl)
    errl <- encodeString(errl, width = NA, justify = "left")
    errl <- paste(" *", errl)
    errm <- paste("an error occurred running the BayesX binary! The following messages are returned:\n",
      paste(errl, collapse = "\n", sep = ""), sep = "")
    warning(errm, call. = FALSE)
  }
  if(length(i <- grep("-nan", ok$log, ignore.case = TRUE))){
    warning("the BayesX engine returned NA samples, please check your model specification! In some cases it can be helpful to center continuous covariates!", call. = FALSE)
  }
  samples <- NULL
  mfile <- grep(paste(x$control$setup$model.name, "_R.r", sep = ""), dir(dir), value = TRUE, fixed = TRUE)
  if(length(mfile)) {
    mfile <- readLines(file.path(dir, mfile))
    mfile <- process_mfile(mfile)
    mspecs <- list("effects" = list(), "model" = list())
    for(j in 1:nrow(mfile$effects)) {
      ts <- read.table(mfile$effects[j, "samples"], header = TRUE)
      ts$intnr <- NULL
      names(ts) <- paste(mfile$effects[j, "terms"], "[", 1:ncol(ts), "]", sep = "")
      samples <- cbind(samples, as.matrix(ts))
      if(!is.na(mfile$effects[j, "varsamples"])) {
        ts <- read.table(mfile$effects[j, "varsamples"], header = TRUE)
        ts$intnr <- NULL
        names(ts) <- paste(mfile$effects[j, "terms"], ":var[", 1:ncol(ts), "]", sep = "")
        samples <- cbind(samples, as.matrix(ts))
      }
      mspecs$effects[[mfile$effects[j, "terms"]]] <- list(
        "basis" = if(!is.na(mfile$effects[j, "basis"])) {
            eval(parse(file = mfile$effects[j, "basis"]))
         } else function(x) { x },
        "family" = mfile$effects[j, "family"],
        "eqntype" = mfile$effects[j, "eqntype"],
        "filetype" = mfile$effects[j, "filetype"],
        "hlevel" =  as.integer(mfile$effects[j, "hlevel"])
      )
    }
    dic <- if(!is.null(mfile$model$dic)) read.table(mfile$model$dic, header = TRUE) else NULL
    samples <- cbind(samples, "dic" = dic$dic, "pd" = dic$pd)
    samples <- as.mcmc(samples)
    attr(samples, "model.specs") <- mspecs
  }

  samples
}

process_mfile <- function(x)
{
  predict <- dic <- NULL
  if(length(i <- grep("predict.res", x, ignore.case = TRUE))) {
    predict <- gsub(";", "", gsub("\\s", "", strsplit(x[i], "=", fixed = TRUE)[[1]][2]))
    x <- x[-i]
  }
  if(length(i <- grep("DIC.res", x, ignore.case = TRUE))) {
    dic <- gsub(";", "", gsub("\\s", "", strsplit(x[i], "=", fixed = TRUE)[[1]][2]))
    x <- x[-i]
  }
  n <- length(x)
  x <- gsub("equation=", "family=", x, fixed = TRUE)
  x <- gsub("-", "_", x, fixed = TRUE)
  samples <- varsamples <- eqntype <- family <- basis <- terms <- hlevel <- filetype <- rep(NA, length = n)

  c2c <- function(x) {
    if(any(grepl("(", x, fixed = TRUE))) {
      i <- regexpr("\\((.*)\\)", x)
      x2 <- substring(x, i + 1, i + attr(i, "match.length") - 2)
      x3 <- gsub(",", ";", x2, fixed = TRUE)
      x <- gsub(x2, x3, x, fixed = TRUE)
    }
    x
  }

  for(j in 1:n) {
    specs <- strsplit(c2c(x[j]), ",")[[1]]
    if(length(tmp <- grep("pathsamples", specs, value = TRUE))) {
      samples[j] <- gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2])
    }
    if(length(tmp <- grep("pathvarsample", specs, value = TRUE))) {
      varsamples[j] <- gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2])
    }
    if(length(tmp <- grep("equationtype", specs, value = TRUE))) {
      eqntype[j] <- gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2])
    }
    if(length(tmp <- grep("filetype", specs, value = TRUE))) {
      filetype[j] <- gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2])
    }
    if(length(tmp <- grep("family", specs, value = TRUE))) {
      family[j] <- tolower(gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2]))
    }
    if(length(tmp <- grep("pathbasis", specs, value = TRUE))) {
      basis[j] <- gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2])
    }
    if(length(tmp <- grep("hlevel", specs, value = TRUE))) {
      hlevel[j] <- as.integer(gsub("\\s", "", strsplit(tmp, "=", fixed = TRUE)[[1]][2]))
    } else hlevel[j] <- 1
    if(length(tmp <- grep("term", specs, value = TRUE))) {
      tt <- strsplit(tmp, "=", fixed = TRUE)[[1]]
      tt <- paste(tt[2:length(tt)], collapse = "=")
      tt <- gsub(";", ",", tt, fixed = TRUE)
      tt <- gsub("^\\s+|\\s+$", "", tt)
      terms[j] <- gsub("\\s", "+", tt)
      
      ## FIXME: this is just a brute force fix for multinom hierarchical structures!
      ##        Searching for the correct category in the samples path.
      if(grepl("multinom", family[j])) {
        cat <- strsplit(samples[j], "_", fixed = TRUE)[[1]]
        cat <- cat[4]
        terms[j] <- paste(terms[j], cat, sep = ":")
      }

    }
  }

  for(j in 1:n) {
    if(grepl("gaussian_re", family[j], fixed = TRUE) |
      grepl("Gaussian_random_effect", family[j], fixed = TRUE))
    {
      i <- grep(eqntype[j], eqntype)
      i <- i[i != j]
      ok <- TRUE
      for(l in i) {
        if(hlevel[l] < 2 & ok) {
          family[j] <- family[l]
          ok <- FALSE
        }
      }
    }
  }

  ## FIXME: this is just a brute force fix for multinom hierarchical structures!
  ##        Searching for the correct category by matching the equationtype.
  ##        Not possible for more than 2 modeled categories!
  if(any(grepl("multinom", family))) {
    for(j in seq_along(hlevel)) {
      if(hlevel[j] > 1 & !grepl(":", terms[j], fixed = TRUE)) {
        et <- eqntype[j]
        ok <- TRUE
        for(i in seq_along(eqntype)) {
          if(i != j & et == eqntype[i]) {
            if(grepl(":", terms[i], fixed = TRUE) & ok) {
              terms[j] <- paste(terms[j], strsplit(terms[i], ":", fixed = TRUE)[[1]][2], sep = ":")
              ok <- FALSE
            }
          }
        }
      }
    }
  }

  known_paramaters <- c("mu", "sigma", "sigma2", "pi", "delta", "tau",
    "a", "b", "p", "df", "rho", "lambda", "alpha", "nu")

  ft <- sapply(strsplit(family, "_"), function(x) {
    x <- gsub("zip", "lambda", x)
    if(length(x) > 1) {
      if(x[2] %in% known_paramaters) x[2] else x[1]
    } else x
  })

  ft <- paste(ft, eqntype, sep = "")

  terms <- paste(terms, ft, sep = ":")

  if(any(i <- duplicated(terms))) {
    for(j in terms[i]) {
      tt <- grep(j, terms, value = TRUE, fixed = TRUE)
      tt <- paste(tt, 1:length(tt), sep = ":")
      terms[terms == j] <- tt
    }
  }
  terms <- paste(terms, hlevel, sep = ":")

  rval <- list("effects" = cbind(terms, filetype, samples, varsamples, basis, family, eqntype, hlevel),
    "model" = list("predict" = predict, "dic" = dic))
  rval
}


resultsBayesX <- function(x, samples, ...)
{
  family <- attr(x, "family")
  grid <- attr(x, "grid")
  response.name <- attr(attr(x, "model.frame"), "response.name")
  if(is.null(grid)) grid <- 100
  if(is.function(family))
    family <- family()
  mspecs <- attr(samples, "model.specs")
  ylevels <- attr(x, "ylevels")
  reference <- attr(x, "reference")
  if(family$cat)
    ylevels <- ylevels[ylevels != reference]

  rename.p <- function(x) {
    if(family$family %in% c("binomial", "multinomial", "quant", "poisson", "zip", "dirichlet")) {
      foo <- switch(family$family,
        "binomial" = function(x) gsub("binomial", "pi", x),
        "multinomial" = function(x) {
          if(any(grepl("):", x, fixed = TRUE))) {
            x <- strsplit(x, "):", fixed = TRUE)
            x <- sapply(x, function(x2) {
              rn <- if(reference == "NA") "" else response.name[1]
              if(length(x2) > 1) {
                x2[2] <- gsub(rn, "", gsub("multinom:", "", x2[2]))
                x2[2] <- gsub("multinomialprobit:", "", x2[2])
                x2 <- paste(x2, collapse = "):")
              } else {
                x2 <- gsub(rn, "", gsub("multinom:", "", x2))
                x2 <- gsub("multinomialprobit:", "", x2)
              }
              x2
            })
            x <- unlist(x)
          } else {
            x <- gsub(response.name[1], "", gsub("multinom:", "", x))
            x <- gsub("multinomialprobit:", "", x)
          }
          x <- gsub("multinommean:", "", x, fixed = TRUE)
          x <- gsub("multinommeanservant:", "", x, fixed = TRUE)
          x
        },
        "quant" = function(x) gsub("quantreg", "mu", x),
        "poisson" = function(x) gsub("poisson", "lambda", x),
        "zip" = function(x) gsub("zip", "lambda", x),
        "dirichlet" = function(x) {
          x <- gsub("dirichletmean", "dirichletalpha1", x)
          x <- gsub("dirichletalpha", "alpha", x)
          x
        }
      )
      x <- foo(x)
    }
    if(grepl("zero-adjusted", family$family)) {
      x <- gsub("binomial", "pi", x)
    }
    x
  }

  if(is.null(mspecs) & is.list(samples))
    mspecs <- attr(samples[[1]], "model.specs")
  names(mspecs$effects) <- rename.p(names(mspecs$effects))

  createBayesXresults <- function(obj, samples, id = NULL, sid = FALSE) {
    if(inherits(samples[[1]], "mcmc.list"))
      samples <- do.call("c", samples)
    chains <- length(samples)
    rval <- vector(mode = "list", length = chains)
    snames <- rename.p(colnames(samples[[1]]))
    if(is.null(snames))
      snames <- attributes(samples[[1]])$dimnames[[2]]
    for(j in 1:chains) {
      if(any(grepl("dic", snames))) {
        DIC <- mean(as.numeric(samples[[j]][, grepl("dic", snames)]))
        pd <- mean(as.numeric(samples[[j]][, grepl("pd", snames)]))
      } else {
        DIC <- pd <- NA
      }

      ## Compute model term effects.
      param.effects <- effects <- effects.hyp <- NULL
      fitted.values <- 0

      ## remove number from id
      id2 <- id
      if(!is.null(family$bayesx$rm.number)) {
        for(i in 1:10)
          id2 <- gsub(i, "", id2)
      }
      if(!is.null(obj$hlevel) & FALSE) {
        if(obj$hlevel > 1)
          id2 <- "gaussian"
      }

      ## Parametric effects.
      if(k <- ncol(obj$X)) {
        nx <- obj$pterms
        nx <- gsub("Intercept", "const", nx, fixed = TRUE)
        pt <- paste(nx, collapse = "+")
        pt <- paste(pt, paste(id2, family$bayesx[[id]][2], sep = ""), if(obj$hlevel > 1) 2 else 1, sep = ":")
        if(any(grepl(pt, snames, fixed = TRUE))) {
          samps <- as.matrix(samples[[j]][, grepl(pt, snames, fixed = TRUE)], ncol = k)
          nx <- gsub("const", "(Intercept)", nx, fixed = TRUE)
          colnames(samps) <- nx
          qu <- t(apply(samps, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE))
          sd <- drop(apply(samps, 2, sd, na.rm = TRUE))
          me <- drop(apply(samps, 2, mean, na.rm = TRUE))
          param.effects <- cbind(me, sd, qu)
          rownames(param.effects) <- nx
          colnames(param.effects) <- c("Mean", "Sd", "2.5%", "50%", "97.5%")
          fitted.values <- as.vector(fitted.values + obj$X %*% param.effects[, 1])
          attr(param.effects, "samples") <- as.mcmc(samps)
        } else warning("please check the BayesX files!")
      }

      ## Smooth terms.
      if(length(i <- grep("sx(", names(mspecs$effects), fixed = TRUE))) {
        sx.smooth <- mspecs$effects[i]
        sx.terms <- unlist(lapply(obj$sx.smooth, function(x) {
          attr(x, "specs")$label
        }))
        sx.extras <- lapply(obj$sx.smooth, function(x) {
          cx <- class(xs <- attr(x, "specs"))
          if(cx %in% "rps.smooth.spec") {
            return(list(paste("sx(", xs$term, ")", sep = ""), paste("sx(", xs$by, ")", sep = "")))
          } else {
            if(cx %in% "re.smooth.spec" & xs$by != "NA") {
              return(list(paste("sx(", xs$by, ")", sep = ""),
                paste("sx(", xs$term, ",by=", xs$by, ")", sep = ""), xs$by))
            } else return(NULL)
          }
        })
        if(length(sx.extras)) {
          sx.terms <- unique(c(sx.terms, unlist(sapply(sx.extras, function(x) { x[[1]] } ))))
        }
        sx.check <- FALSE
        if(!is.null(sx.terms)) {
          nxc <- sapply(strsplit(names(sx.smooth), ":", fixed = TRUE), function(x) { x[1] })
          for(char in c(".", "_")) {
            sx.terms <- gsub(char, "", sx.terms, fixed = TRUE)
            nxc <- gsub(char, "", nxc, fixed = TRUE)
          }
          if(any(take <- nxc %in% sx.terms)) {
            sx.check <- TRUE
            sx.smooth <- sx.smooth[take]
          }
        }
        for(idj in c(id2, "gaussian")) {
          if(sx.check & length(i <- grep(paste(paste(idj, family$bayesx[[id]][2], sep = ""),
              if(obj$hlevel > 1) 2 else 1, sep = ":"), names(sx.smooth), fixed = TRUE))) {
            if(!is.list(effects))
              effects <- list()
            for(i in names(sx.smooth[i])) {
              tn <- grep(i, snames, fixed = TRUE, value = TRUE)
              psamples <- as.matrix(samples[[j]][, snames %in% tn], ncol = length(tn))

              ## Possible variance parameter samples.
              if(length(vs <- grep(":var[", colnames(psamples), fixed = TRUE))) {
                vsamples <- psamples[, vs[1]]
                psamples <- psamples[, -vs, drop = FALSE]
              }

              ## Prediction matrix.
              re.check <- NULL
              if(length(sx.extras)) {
                re.check <- unlist(sapply(sx.extras, function(x) {
                  ii <- strsplit(i, ":")[[1]][1]
                  return(if(ii %in% x[[1]]) {
                      if(length(x) > 2) x[2:3] else x[[2]]
                    } else NULL)
                }))
              }

              is.re.slope <- FALSE
              if(is.null(re.check)) {
                tn0 <- strsplit(i, ":")[[1]][1]
              } else {
                if(length(re.check) < 2) {
                  tn0 <- re.check[1]
                } else {
                  tn0 <- re.check[1]
                  is.re.slope <- TRUE
                }
              }

              tn <- gsub(")", "", gsub("sx(", "", tn0, fixed = TRUE), fixed = TRUE)
              tn <- strsplit(tn, ",", fixed = TRUE)[[1]]
              tn <- tn[!grepl("by", tn)]
              tn1 <- eval(parse(text = tn0))

              basis <- sx.smooth[[i]]$basis

              if(tn1$by != "NA") {
                tn <- c(tn, tn1$by)
                get.X <- function(data) {
                  X <- if(ncol(data) > 1) {
                    if(!is.factor(data[, ncol(data)]) & !is.re.slope) {
                      as.numeric(data[, ncol(data)]) * basis(data[, 1:(ncol(data) - 1), drop = FALSE])
                    } else {
                      basis(data[, 1:(ncol(data) - 1), drop = FALSE])
                    }
                  } else basis(data[, 1, drop = FALSE])
                  return(X)
                }
              } else {
                get.X <- basis
              }

              X <- sx.smooth[[i]]$basis(attr(x, "model.frame")[1, tn, drop = FALSE])
            
              get.mu <- function(X, g) {
                X %*% as.numeric(g)
              }

              ## Compute final smooth term object.
              stype <- attr(X, "type")
              class(tn1) <- paste(stype, "smooth.spec", sep = ".")

              fst <- try(compute_term(tn1, get.X = get.X, get.mu = get.mu,
                psamples = psamples, vsamples = vsamples, asamples = NULL, FUN = NULL,
                snames = snames, effects.hyp = effects.hyp, fitted.values = fitted.values,
                data = attr(x, "model.frame")[, tn, drop = FALSE], grid = grid,
                hlevel = obj$hlevel, sx = TRUE, re.slope = is.re.slope), silent = TRUE)
           
              if(!inherits(fst, "try-error")) {
                attr(fst$term, "specs")$get.mu <- get.mu
                attr(fst$term, "specs")$basis <- sx.smooth[[i]]$basis
                if(sid)
                  attr(fst$term, "specs")$label <- paste(attr(fst$term, "specs")$label, id, sep = ":")
                if(is.re.slope) {
                  attr(fst$term, "specs")$by <- re.check[2]
                  attr(fst$term, "specs")$re.slope <- TRUE
                }

                ## Add term to effects list.
                effects[[paste(tn0, stype, sep = ":")]] <- fst$term
                effects.hyp <- fst$effects.hyp
                fitted.values <- fst$fitted.values
                rm(fst)
              } else warning(paste("problems computing effect object for term ", tn1$label, "!", sep = ""))
            }
          }
        }
      }

      ## Compute partial residuals.
      if(!is.null(effects) & obj$hlevel < 2) {
        if(length(obj$response)) {
          if(obj$response %in% names(attr(x, "model.frame"))) {
            if(!is.na(link <- family$links[id])) {
              link <- make.link2(link)
            } else link <- NULL
            effects <- partial.residuals(effects, attr(x, "model.frame")[[obj$response]],
              fitted.values, family)
          }
        }
      }

      ## Stuff everything together.
      rval[[j]] <- list(
        "model" = list("DIC" = DIC, "pd" = pd, "N" = nrow(attr(x, "model.frame")),
        "formula" = obj$formula, "hlevel" = obj$hlevel), "param.effects" = param.effects, "effects" = effects,
        "effects.hyp" = effects.hyp,
        "fitted.values" = if(length(fitted.values) < 2) NA else fitted.values 
      )

      class(rval[[j]]) <- "bayesr"
    }
    names(rval) <- paste("Chain", 1:chains, sep = "_")
    if(length(rval) < 2) {
      rval <- rval[[1]]
    }
    attr(rval, "fixed.names") <- TRUE
    class(rval) <- "bayesr"
    return(rval)
  }

  nx <- names(x)
  nx <- nx[!(nx %in% c("call", "family"))]
  if(is(samples, "mcmc"))
    samples <- as.mcmc.list(list(samples))
  if(!all(c("formula", "fake.formula") %in% names(x))) {
    rval <- list()
    fn <- family$names
    if(is.null(names(x))) {
      nx <- paste("h", 1:length(x), sep = "")
      names(x) <- nx
    }
    if(length(unique(nx)) != length(nx)) {
      nx <- paste(rep(nx, length.out = length(nx)), 1:length(nx), sep = ":")
      names(x) <- nx
    }
    if(family$cat) {
      fn <- rmf(ylevels)
      names(x) <- nx <- fn
    }
    if(length(fn) != length(nx))
      fn <- nx
    for(j in seq_along(nx)) {
      if(!all(c("formula", "fake.formula") %in% names(x[[nx[j]]]))) {
        rval[[nx[j]]] <- list()
        nh <- NULL
        for(ii in seq_along(x[[nx[j]]])) {
          nh <- c(nh, x[[nx[j]]][[ii]]$hlevel)
          rval[[nx[j]]][[ii]] <- createBayesXresults(x[[nx[j]]][[ii]],
            samples, id = fn[j], sid = length(nx) > 1)
        }
        if(any(dups <- duplicated(nh))) {
          dups <- nh[dups]
          for(dn in dups)
            nh[nh == dn] <- paste(nh[nh == dn], 1:length(nh[nh == dn]), sep = "-")
        }
        names(rval[[nx[j]]]) <- paste("h", nh, sep = "")
        attr(rval[[nx[j]]], "hlevel") <- TRUE
      } else {
        rval[[nx[j]]] <- createBayesXresults(x[[nx[j]]], samples, id = fn[j], sid = length(nx) > 1)
      }
      class(rval[[nx[j]]]) <- "bayesr"
    }
    if(length(nx) > 1)
      names(rval) <- fn
    else
      rval <- rval[[1]]
    attr(rval, "fixed.names") <- TRUE
    class(rval) <- "bayesr"
    return(rval)
  } else {
    return(createBayesXresults(x, samples, id = family$names))
  }
}


########################################
## (2) BayesX model term construction ##
########################################
sx <- function(x, z = NULL, bs = "ps", by = NA, ...)
{
  by <- deparse(substitute(by), backtick = TRUE, width.cutoff = 500)
  call <- match.call()

  available.terms <- c(
    "rw1", "rw2",
    "season",
    "ps", "psplinerw1", "psplinerw2", "pspline",
    "te", "pspline2dimrw2", "te1", "pspline2dimrw1",
    "kr", "kriging",
    "gk", "geokriging",
    "gs", "geospline",
    "mrf", "spatial",
    "bl", "baseline",
    "factor",    
    "ridge", "lasso", "nigmix",
    "re", "ra", "random",
    "cs", "catspecific",
    "offset",
    "generic",
    "rps", "hrandom_pspline"
  )
  if(!bs %in% available.terms) stop(paste("basis type", sQuote(bs), "not supported by BayesX"))

  if(bs %in% c("rsps", "hrandom_pspline")) {
    bs <- "rsps"
    x <- deparse(substitute(x), backtick = TRUE, width.cutoff = 500)
    rcall <- paste("R2BayesX:::r(x = ", by, ", bs = ", sQuote(bs), ", by = ", x, ", ...)", sep = "")
    rval <- eval(parse(text = rcall))
  } else {
    if(length(grep("~", term <- deparse(call$x))) && bs %in% c("re", "ra", "random")) {
      x <- deparse(substitute(x), backtick = TRUE, width.cutoff = 500)
      rcall <- paste("R2BayesX:::r(x = ", x, ", by = ", by, ", ...)", sep = "")
      rval <- eval(parse(text = rcall))
    } else {
      k <- -1
      m <- NA
      xt <- list(...)
      if(is.null(xt$lambda))
        xt$lambda <- 100
      if("m" %in% names(xt))
        stop("argument m is not allowed, please see function s() using this specification!")
      if("k" %in% names(xt))
        stop("argument k is not allowed, please see function s() using this specification!")
      if(!is.null(xt$xt))
        xt <- xt$xt
      warn <- getOption("warn")
      options("warn" = -1)
      if(by != "NA" && is.vector(by) && length(by) < 2L && !is.na(as.numeric(by))) {
        xt["b"] <- by
        by <- "NA"
      }
      options("warn" = warn)
      if(bs %in% c("pspline2dimrw1", "pspline2dimrw2", "te",
        "gs", "geospline", "kr", "gk", "kriging", "geokriging")) {
        if(by != "NA")
          stop(paste("by variables are not allowed for smooths of type bs = '", bs, "'!", sep = ""))
      }
      if(bs == "te1") bs <- "pspline2dimrw1"
      if(!is.null(xt$knots))
        xt$nrknots <- xt$knots
      if(bs %in% c("ps", "te", "psplinerw1", "psplinerw2", "pspline",
        "pspline2dimrw2", "pspline2dimrw1", "gs", "geospline")) {
        if(!is.null(xt$degree))
          m <- xt$degree
        if(!is.null(xt$order)) {
          if(is.na(m))
            m <- c(3L, xt$order)
          else
            m <- c(m[1L], xt$order)
        }
        if(length(m) < 2L && (bs %in% c("gs", "geospline")))
          m <- c(3L, 1L)
        if(length(m) < 2L && is.na(m))
          m <- c(3L, 2L)   
        if(is.null(xt$order) && length(m) < 2L)
          m <- c(m, 2L)
        if(is.null(xt$order) && length(m) < 2L)
          m <- c(m, 1L)
        m[1L] <- m[1L] - 1L
        if(!is.null(xt$nrknots))
          k <- xt$nrknots + m[1L]
        else {
          if(bs %in% c("ps", "psplinerw1", "psplinerw2", "pspline"))
            k <- 20L + m[1L]
          else
            k <- 10L + m[1L]
        }
      }
      if(bs %in% c("kr", "gk", "kriging", "geokriging")) {
        m <- c(1L, 1L)
        if(!is.null(xt$nrknots)) {
          k <- xt$nrknots
        } else k <- -1L
      }
      if(!is.null(xt$map))
        xt$map.name <- rmf(as.character(call$map))
      xt[c("degree", "order", "knots", "nrknots")] <- NULL
      if(!length(xt))
        xt <- NULL
      if(!is.null(call$z)) 
        term <- c(term, deparse(call$z))
      rval <- mgcv::s(x, z, k = k, bs = bs, m = m, xt = xt)
      rval$term <- term
      rval$dim <- length(term)
      rval$by <- by
      rval$label <- paste("sx(", paste(term, collapse = ",", sep = ""),
        if(by != "NA") paste(",by=", by, sep = "") else NULL, ")", sep = "")
    }
  }

  return(rval)
}

sx.construct <- function(object, data)
{
  UseMethod("sx.construct")
}

sx.construct.default <- function(object, data) 
{
  cl <- class(object)
  bs <- gsub(".smooth.spec", "", cl, fixed = TRUE)
  info <- if(cl == paste(bs, "smooth.spec", sep = ".")) {
    paste("with basis", sQuote(bs))
  } else {
    paste("of class", sQuote(cl))
  }
  stop(paste("BayesX does not support smooth terms ", info,
    ", it is recommended to use sx() for specifying smooth terms",
    sep = ""))
}

do.xt <- function(term, object, not = NULL, noco = FALSE)
{
  if(!is.null(object$xt)) {
    names.xt <- names(object$xt)
    if(is.null(not))
      not <- "not"
    count <- 1
    co <- ","
    if(noco)
      co <- NULL
    for(name in names.xt) {
      if(count > 1)
        co <- ","
      if(!name %in% not) {
        if(name %in% c("full", "catspecific", "center", "derivative", "nofixed") || 
          is.logical(object$xt[[name]])) {
          if(is.logical(object$xt[[name]])) {
            if(object$xt[[name]])
              term <- paste(term, co, name, sep = "")
          } else term <- paste(term, co, name, sep = "")
          count <- count + 1
        } else {
          term <- paste(term, co, name, "=", object$xt[name], sep = "")
          count <- count + 1
        }
      }
    }
  }

  return(term)
}

sx.construct.ps.smooth.spec <- sx.construct.psplinerw1.smooth.spec <-
sx.construct.psplinerw2.smooth.spec <- sx.construct.pspline.smooth.spec <-
function(object, data)
{
  if(length(object$p.order) == 1L) 
    m <- rep(object$p.order, 2L)
  else 
    m <- object$p.order
  m[is.na(m)] <- 2L
  object$p.order <- m
  object$p.order[1L] <- object$p.order[1L] + 1L
  if(class(object) == "psplinerw1.smooth.spec")
    object$p.order[2L] <- 1L
  if(class(object) == "psplinerw2.smooth.spec")
    object$p.order[2L] <- 2L
  if(object$bs.dim < 0L)
    object$bs.dim <- 10L
  if(length(object$p.order) > 1L) {
    if(object$p.order[2L] > 2L) {
      warning("order of the difference penalty not supported by BayesX, set to 2!")
      object$p.order <- c(object$p.order[1L], 2L)
    }
  }
  nrknots <- object$bs.dim - object$p.order[1L] + 1L
  if(nrknots < 5L) {
    warning("number of inner knots smaller than 5 not supported by BayesX, set to 5!",
      call. = FALSE)
    nrknots <- 5L
  }
  termo <- object$term
  term <- paste(termo, "(psplinerw", object$p.order[2L], ",nrknots=",
    nrknots, ",degree=", object$p.order[1L], sep = "")
  object$xt[c("knots", "nrknots", "degree")] <- NULL
  term <- paste(do.xt(term, object, NULL), ")", sep = "")
  if(object$by != "NA")
    term <- make_by(term, object, data)
  
  return(term)
}

sx.construct.ra.smooth.spec <- sx.construct.re.smooth.spec <-
sx.construct.random.smooth.spec <- function(object, data)
{
  term <- object$term
  if(is.null(object$ins))
    term <- paste(term, "(random", sep = "")
  else
    term <- paste(term, "(hrandom", sep = "")
  term <- paste(do.xt(term, object, NULL), ")", sep = "")
  if(object$by != "NA")
    term <- make_by(term, object, data)

  return(term)
}

sx.construct.rps.smooth.spec <- function(object, data)
{
  term <- paste(object$term, "(hrandom_pspline,centermethod=meansum2", sep = "")
  term <- paste(do.xt(term, object, NULL), ")", sep = "")
  term <- paste(object$by, term , sep = "*")

  return(term)
}

sx.construct.kr.smooth.spec <- sx.construct.kriging.smooth.spec <- function(object, data)
{
  termo <- object$term
  if(length(termo) < 2L)
    stop("kriging method needs two terms!")
  if(object$bs.dim < 0L)
    object$bs.dim <- 10L
	nrknots <- object$bs.dim
  xt <- object$xt
  if(is.null(xt$full))
    term <- paste(termo[1L], "*", termo[2L], "(kriging,nrknots=", nrknots, sep = "")
  else {
    term <- paste(termo[1L], "*", termo[2L], "(kriging,full", sep = "")    
    object$xt$full <- NULL
  }
  term <- paste(do.xt(term, object, NULL), ")", sep = "")
  if(object$by != "NA")
    term <- make_by(term, object, data)

  return(term)
}

construct.shrw <- function(object, data, what)
{
  term <- object$term
  term <- paste(term, "(", what, sep = "")
  term <- paste(do.xt(term, object, NULL), ")", sep = "")
  if(object$by != "NA")
    term <- make_by(term, object, data)

  return(term)
}

sx.construct.offset.smooth.spec <- function(object, data)
{
  return(construct.shrw(object, data, "offset"))
}

sx.construct.mrf.smooth.spec <- sx.construct.spatial.smooth.spec <- function(object, data)
{
  if(is.null(object$xt))
    stop("need to supply a map object in argument xt!")  
  map.name <- help.map.name(deparse(substitute(object, env = .GlobalEnv), 
    backtick = TRUE, width.cutoff = 500L))
  if(!is.null(object$xt$map.name))
    map.name <- object$xt$map.name
  if(!is.list(object$xt))
    object$xt <- list(object$xt)
  map.name <- rmf(gsub("\\s", "", paste(map.name, sep = "", collapse = "")))

  map <- object$xt$map
  if(is.null(map)) {
    if(!is.null(object$xt$polys))
      map <- object$xt$polys
    if(!is.null(object$xt$penalty))
      map <- object$xt$penalty
  }
  if(is.null(map))
    map <- object$xt$gra
  if(is.null(map)) {
    if(!is.list(object$xt[[1L]])) {
      if(inherits(object$xt[[1L]], "gra"))
        map <- object$xt[[1L]]
      else
        map <- object$xt
    } else map <- object$xt[[1L]]
    if(is.null(map)) {
      map <- object$xt
      if(is(map, "SpatialPolygonsDataFrame"))
        map <- SPDF2bnd(map)
      if(is.null(map) || (!is.list(map) && !inherits(map, "bnd") || !inherits(map, "gra")))
        stop("need to supply a bnd or graph file object in argument xt!")
    }
  }
  if(is(map, "nb"))
    map <- nb2gra(map)
  if(inherits(map, "SpatialPolygons"))
    map <- sp2bnd(map)
  if(!inherits(map, "bnd") && !inherits(map, "gra")) {
    if(is.list(map))
      class(map) <- "bnd"
    else
      class(map) <- "gra"
  }

  write <- function(dir = NULL) {
    if(is.null(dir)) {
      dir.create(dir <- tempfile())
      on.exit(unlink(dir))
    } else dir <- path.expand(dir)
    if(!file.exists(dir)) dir.create(dir)

    counter <- NULL
    ok <- TRUE
    files <- list.files(dir)
    while(ok) {
      classm <- class(map)
      if(length(classm) > 1L)
        if("list" %in% classm)
          class(map) <- classm[classm != "list"]
      mapfile <- paste(map.name, counter, ".", class(map), sep = "")[1]
      if(any(grepl(mapfile, files))) {
        if(is.null(counter))
          counter <- 0L
        counter <- counter + 1L
      } else ok <- FALSE
    }
    mapfile <- file.path(dir, mapfile)

    prg <- paste("map", map.name)
    if(inherits(map, "bnd")) {
      if(!any(is.na(poly.names <- as.integer(names(map))))) {
        poly.names <- sort(poly.names)
        poly.names <- as.character(poly.names)
      } else poly.names <- sort(names(map))
      map <- map[poly.names]
      class(map) <- "bnd"
      write.bnd(map = map, file = mapfile, replace = TRUE)
      prg <- c(prg, paste(map.name, ".infile using ", mapfile, sep = ""))
    } else {
      if(!is.character(map)) {
        write.gra(map = map, file = mapfile, replace = TRUE)
        prg <- c(prg, paste(map.name, ".infile, graph using ", mapfile, sep = ""))
      } else {
        stopifnot(is.character(map))
        pos <- regexpr("\\.([[:alnum:]]+)$", map)
        fext <- ifelse(pos > -1L, substring(map, pos + 1L), "")
        if(fext == "gra")
          prg <- c(prg, paste(map.name, ".infile, graph using ", path.expand(map), sep = ""))
        else
          prg <- c(prg, paste(map.name, ".infile using ", path.expand(map), sep = ""))
      }
    }
    prg
  }

  term <- object$term
  term <- paste(term, "(spatial,map=", map.name, sep = "")
  term <- paste(do.xt(term, object, c("map", "polys", "penalty", "map.name")), ")", sep = "")
  if(object$by != "NA")
    term <- make_by(term, object, data)

  attr(term, "write") <- write
  attr(term, "map.name") <- map.name

  return(term)
}


make_by <- function(term, object, data)
{
  if(!missing(data) && !is.character(data)) {
    by <- data[[object$by]]
    if(is.factor(by) && nlevels(by) > 1) {
      nocenter <- paste(c("lasso", "nigmix", "ridge", "ra"), ".smooth.spec", sep = "")
      term <- paste(paste(rmf(object$by), rmf(levels(by)), sep = ""), "*", term, sep = "")
      if((k <- length(term)) > 1)
        for(j in 1:k) 
          if(!grepl("center", term[j]) && is.null(object$xt$center))
            if(!class(object) %in% nocenter)
              term[j] <- gsub(")", ",center)", term[j])
    } else term <- paste(rmf(object$by), "*", term, sep = "")
  } else term <- paste(rmf(object$by), "*", term, sep = "")

  return(term)
}

rmf <- function(x) 
{
  for(i in 1L:length(x)) {
    for(char in c("+", "-", "*", ":", "^", "/", " ", "(", ")", "]", "[",
      ",", ".", "<", ">", "?", "!", "'", "#", "~", "`", ";", "=", "&", "$", "@")) {
      x[i] <- gsub(char, "", x[i], fixed = TRUE)
    }
  }

  return(x)
}

help.map.name <- function(x)
{
  if(is.null(x))
    return("")
  x <- splitme(x)
  if(resplit(x[1L:2L]) == "s(") {
    x <- resplit(c("sfoofun", x[2L:length(x)]))
    x <- eval(parse(text = x), envir = parent.frame())
    if(is.null(x))
      x <- "map"
  } else  x <- "map"

  return(x)
}

sfoofun <- function(x, xt = NULL, ...)
{
  if(is.null(x) || is.null(xt))
    return(NULL)
  x <- rval <- deparse(substitute(xt), backtick = TRUE, width.cutoff = 500L)
  if(is.list(xt) && length(xt)>1)
    for(i in 1L:length(xt))
      if(inherits(xt[[i]], "bnd") || inherits(xt[[i]], "gra") || inherits(xt[[i]], "list")) {
        rval <- strsplit(x, ",", " ")[[1L]]
        if(length(rval) > 1L)
          rval <- rval[i]
      }
  rval <- splitme(rval)
  if(length(rval) > 5L)
    if(resplit(rval[1L:5L]) == "list(")
      rval <- rval[6L:length(rval)]
  if(rval[length(rval)] == ")")
    rval <- rval[1L:(length(rval) - 1)]
  if(any(grepl("=", rval)))
    rval <- rval[(grep("=", rval) + 2L):length(rval)]
  rval <- resplit(rval)
   
  return(rval)
}

splitme <- function(x) {
  return(strsplit(x, "")[[1L]])
}

resplit <- function(x) {
  if(!is.null(x))
    x <- paste(x, sep = "", collapse = "")
	
  return(x)
}

Try the BayesR package in your browser

Any scripts or data that you put into this service are public.

BayesR documentation built on May 2, 2019, 6:16 p.m.