R/BAMLSS.R

Defines functions stg BIC.bamlss AIC.bamlss sum_diag2 sum_diag scale_model.frame scale.model.matrix XnotinY XinY match.index matrix_inv blockMatrixDiagonal h_response model.response2 plot.bamlss.residuals residuals.bamlss create.dp kfitted score term.labels2 term.labels coef.bamlss coef.bamlss confint.bamlss continue samples grep2 fitted.bamlss results.bamlss.default get_sterms_labels get_pterms_labels has_pterms has_sterms has_response has_intercept model.terms terms.bamlss.formula has_dot drop.terms.bamlss print.bamlss.formula formula.bamlss.terms logLik.bamlss DIC.bamlss print.bamlss print.summary.bamlss summary.bamlss delete.NULLs delete.args bamlss_factor2d_plot bamlss_random_plot plot.bamlss.effect.default plot.bamlss.effect plot.bamlss.results plot.bamlss.results plot.bamlss traceplot2 smooth.construct.fdl.smooth.spec Predict.matrix.strandom.smooth smooth.construct.str.smooth.spec Predict.matrix.kriging.smooth smooth.construct.kr.smooth.spec krDesign2D krDesign1D Predict.matrix.harmon.smooth smooth.construct.ha.smooth.spec boost.nnet.predict Predict.matrix.nnet.smooth smooth.construct.nnet.smooth.spec n Predict.matrix.lasso.smooth smooth.construct.la.smooth.spec blockstand la rs.plot Predict.matrix.rs.smooth smooth.construct.rs.smooth.spec rs smooth.construct.rsc.smooth.spec rsc s2 fitted.bamlss predict.bamlss predict.bamlss c95 add.partial compute_s.effect fitted_matrix quick_quantiles c.bamlss process.chains randomize formula_hierarchical formula_insert formula_hcheck formula_rm_at formula_at formula_and response.name formula_extend fake.formula all.labels.formula terms.formula2 splitFormula all.vars.formula make_fFormula bamlss.formula.cat get_formula_envir bamlss.formula as.list.Formula complete.bamlss.family bamlss.family rm_infinite bamlss.model.frame samplestats get.all.parnames bamlss.setup bamlss par2list parameters smooth.construct.bamlss.terms smooth.construct model.matrix.bamlss.terms extract.design.construct model.search do.XWX make.prior make.fit.fun sparse.matrix.fit.fun sparse.solve sparse.backsolve sparse.forwardsolve sparse.chol sparse.setup sparse.matrix.ordering sparse.matrix.index design.construct match.index.ff print.bamlss.frame bamlss.frame

Documented in bamlss bamlss.family bamlss.formula bamlss.frame bamlss.model.frame boost.nnet.predict c95 coef.bamlss confint.bamlss continue fitted.bamlss la model.matrix.bamlss.terms n parameters plot.bamlss plot.bamlss.residuals plot.bamlss.results predict.bamlss print.summary.bamlss randomize residuals.bamlss results.bamlss.default s2 samples samplestats smooth.construct smooth.construct.bamlss.terms summary.bamlss terms.bamlss.formula

## Create a 'bamlss.frame'.
bamlss.frame <- function(formula, data = NULL, family = "gaussian",
  weights = NULL, subset = NULL, offset = NULL, na.action = na.omit,
  contrasts = NULL, knots = NULL, specials = NULL, reference = NULL,
  model.matrix = TRUE, smooth.construct = TRUE, ytype = c("matrix", "vector", "integer"),
  scale.x = FALSE, scale.d = FALSE, ...)
{
  ## Parse family object.
  family <- bamlss.family(family, ...)

  ## Parse formula.
  if(!inherits(formula, "bamlss.formula"))
    formula <- bamlss.formula(formula, family, specials, env = parent.frame())
  if(!is.null(attr(formula, "orig.formula")))
    formula <- attr(formula, "orig.formula")

  ## Setup return object.
  bf <- list()
  bf$call <- match.call()

  ## Create the model frame.
  bf$model.frame <- bamlss.model.frame(formula, data, family, weights,
    subset, offset, na.action, specials, contrasts)

  if(!inherits(bf$model.frame, "ffdf")) {
    ## Type of y.
    ytype <- match.arg(ytype)

    ## Process categorical responses and assign 'y'.
    cf <- bamlss.formula.cat(formula, family, bf$model.frame, reference)
    if(!is.null(cf) & (ytype != "integer")) {
      rn <- response.name(formula, hierarchical = FALSE, na.rm = TRUE)
      hrn <- response.name(formula, hierarchical = TRUE, na.rm = TRUE)
      orig.formula <- formula
      formula <- cf$formula
      if(ytype == "matrix") {
        if(is.factor(bf$model.frame[[rn[1]]])) {
          f <- as.formula(paste("~ -1 +", rn[1]), env = NULL)
          bf$y <- bf$model.frame[rn]
          bf$y[rn[1]] <- model.matrix(f, data = bf$model.frame)
          colnames(bf$y[[rn[1]]]) <- rmf(gsub(rn[1], "", colnames(bf$y[[rn[1]]])))
          bf$y[[rn[1]]] <- bf$y[[rn[1]]][, rmf(c(names(formula), cf$reference))]
        } else {
          bf$y <- bf$model.frame[rn]
        }
      } else {
        bf$y <- bf$model.frame[rn]
        if(is.factor(bf$model.frame[[rn[1]]])) {
          if(ytype == "integer")
            bf$y[[1]] <- as.integer(bf$y[[1]]) - if(nlevels(bf$model.frame[[rn[1]]]) < 3) 1L else 0L
        }
      }
      if(length(hrn) > length(rn)) {
        ynot <- hrn[!(hrn %in% rn)]
        bf$y <- cbind(bf$y, bf$model.frame[ynot])
      }
      attr(bf$y, "reference") <- cf$reference
      family$names <- names(formula)
      family$links <- rep(family$links, length.out = length(formula))
      names(family$links) <- names(formula)
      attr(formula, "orig.formula") <- orig.formula
    } else {
      rn <- response.name(formula, hierarchical = FALSE, keep.functions = TRUE)
      rn <- rn[rn %in% names(bf$model.frame)]
      bf$y <- bf$model.frame[rn]
      for(j in rn) {
        if(is.factor(bf$y[[j]]) & (ytype == "matrix")) {
          f <- as.formula(paste("~ -1 +", j), env = NULL)
          bf$y[j] <- model.matrix(f, data = bf$model.frame)
        }
        if(is.factor(bf$y[[j]]) & (ytype == "integer")) {
          bf$y[[j]] <- as.integer(bf$y[[j]]) - if(nlevels(bf$y[[j]]) < 3) 1L else 0L
        }
      }
    }
  } else {
    rn <- response.name(formula, hierarchical = FALSE, keep.functions = TRUE)
    rn <- rn[rn %in% names(bf$model.frame)]
    bf$y <- bf$model.frame[rn]
  }

  bf$formula <- formula
  attr(bf$formula, "response.name") <- rn

  ## Add the terms object.
  bf$terms <- terms.bamlss.formula(formula, data = data, drop = FALSE, specials = specials, ...)

  ## Process possible score and hess functions.
  if(!is.null(score <- family$score)) {
    if(is.function(score))
      score <- list(score)
    family$score <- rep(score, length.out = length(formula))
    names(family$score) <- names(formula)
  }
  if(!is.null(hess <- family$hess)) {
    if(is.function(hess))
      hess <- list(hess)
    family$hess <- rep(hess, length.out = length(formula))
    names(family$hess) <- names(formula)
  }

  ## Add more functions to family object.
  bf$family <- complete.bamlss.family(family)

  if(inherits(bf$model.frame, "data.frame") & scale.d)
    bf$model.frame <- scale_model.frame(bf$model.frame, not = rn)

  ## Assign the 'x' master object.
  bf$x <- design.construct(bf$terms, data = bf$model.frame, knots = knots,
    model.matrix = model.matrix, smooth.construct = smooth.construct, model = NULL,
    scale.x = scale.x, specials = specials, ...)

  bf$knots <- knots

  ## Assign class and return.
  class(bf) <- c("bamlss.frame", "list")

  return(bf)
}


## Simple print method for 'bamlss.frame'
print.bamlss.frame <- function(x, ...)
{
  cat("'bamlss.frame' structure:", "\n")  
  nx <- c("call", "model.frame", "formula", "family", "terms", "x", "y", "knots")
  nx <- c(nx, names(x)[!(names(x) %in% nx)])
  for(i in nx) {
    if(!is.null(x[[i]])) {
      cat("  ..$", i, "\n")
      if(i == "x") {
        for(j in names(x[[i]])) {
          cat("  .. ..$", j, "\n")
          if(!all(c("formula", "fake.formula") %in% names(x[[i]][[j]]))) {
            for(k in names(x[[i]][[j]])) {
              cat("  .. .. ..$", k, "\n")
              for(d in names(x[[i]][[j]][[k]])) {
               cat("  .. .. .. ..$", d, "\n")
              }
            }
          } else {
            for(k in names(x[[i]][[j]]))
              cat("  .. .. ..$", k, "\n")
          }
        }
      }
      if(i == "y") {
        for(j in names(x[[i]])) {
          cat("  .. ..$", j, "\n")
        }
      }
    }
  }
  invisible(NULL)
}


## ff version for indexing.
match.index.ff <- function(x)
{
#  nodups <- ffwhich(x, !duplicated(x))
#  ind <- ffdfmatch(x, x[nodups, , drop = FALSE])
#  ord <- fforder(ind)
#  sindex <- ind[ord]
#  
#  return(list("match.index" = ind, "nodups" = nodups, "order" = ord, "sorted.index" = sindex, "uind" = ind[nodups]))
## FIXME: ff support!
  match.index(x)
}


## Compute the 'bamlss.frame' 'x' master object.
design.construct <- function(formula, data = NULL, knots = NULL,
  model.matrix = TRUE, smooth.construct = TRUE, binning = FALSE,
  before = TRUE, gam.side = NULL, model = NULL, drop = NULL,
  scale.x = TRUE, absorb.cons = NULL, sparse.cons = 0, specials = NULL, ...)
{
  if(!model.matrix & !smooth.construct)
    return(NULL)

  if(is.null(gam.side))
    gam.side <- if(binning) FALSE else TRUE

  if(inherits(formula, "bamlss.frame")) {
    data <- if(is.null(data)) model.frame(formula) else data
    formula <- formula(formula)
  }
  if(!inherits(formula, "bamlss.terms")) {
    if(!inherits(formula, "bamlss.formula"))
      formula <- bamlss.formula(formula, ...)
    if(inherits(formula, "bamlss.formula"))
      formula <- terms.bamlss.formula(formula, data = data, ...)
  }
  formula <- formula.bamlss.terms(formula)
  if(is.null(data))
    stop("data needs to be supplied!")

  no_ff <- !inherits(data, "ffdf")

  if(!is.character(data) & no_ff) {
    if(!inherits(data, "data.frame"))
      data <- as.data.frame(data)
  }
  if(is.character(data)) {
    ## data <- read.table.ffdf(file = data,
    ##   na.strings = "", header = TRUE, sep = ",")
    ## FIXME: ff data.frames!
    data <- read.table(file = data, header = TRUE, ...)
  }
  if(inherits(data, "ffdf")) {
    before <- TRUE
    gam.side <- FALSE
    if(is.null(binning))
      binning <- TRUE
  }
  if(!is.null(model))
    formula <- model.terms(formula, model)
  if(!binning)
    binning <- NULL

  assign.design <- function(obj, dups = NULL)
  {
    if(!is.null(dups) & no_ff) {
      if(any(dups)) {
        mi <- match.index(data[, all.vars(obj$fake.formula), drop = FALSE])
        obj[names(mi)] <- mi
        data <- subset(data, !dups)
      }
    }
    obj$binning <- binning
    if(!all(c("formula", "fake.formula") %in% names(obj)))
      return(obj)
    if(model.matrix) {
      if(!inherits(data, "ffdf")) {
        obj$model.matrix <- model.matrix(drop.terms.bamlss(obj$terms,
          sterms = FALSE, keep.response = FALSE, data = data, specials = specials), data = data)
        if(ncol(obj$model.matrix) > 0) {
          if(scale.x)
            obj$model.matrix <- scale.model.matrix(obj$model.matrix)
        } else obj$model.matrix <- NULL
      } else {
        mm_terms <- drop.terms.bamlss(obj$terms,
          sterms = FALSE, keep.response = FALSE, data = NULL, specials = specials)
        mm_intercept <- attr(mm_terms, "intercept") > 0
        mm_vars <- all.vars.formula(mm_terms)
        if(mm_intercept) {
          ## FIXME: ff support!
          ## obj$model.matrix <- ffdf("Intercept" = ff(1, length = nrow(data)))
          obj$model.matrix <- cbind("Intercept" = rep(1, length = nrow(data)))
        }
        if(!is.null(mm_vars)) {
          if(!all(mm_vars %in% colnames(data)))
            stop("variables missing in supplied data!")
          if(mm_intercept) {
            for(v in mm_vars)
              obj$model.matrix[[v]] <- data[[v]]
          } else {
            obj$model.matrix <- data[mm_vars]
          }
          if(is.numeric(binning)) {
            for(v in mm_vars) {
              obj$model.matrix[[v]] <- round(obj$model.matrix[[v]], binning)
            }
          }
        }
        if(!is.null(obj$model.matrix)) {
          bind <- match.index.ff(obj$model.matrix)
          obj$model.matrix <- as.matrix(as.data.frame(obj$model.matrix[bind$nodups, , drop = FALSE]))
          if(is.null(colnames(obj$model.matrix)) & ncol(obj$model.matrix) < 2) {
            if(mm_intercept)
              colnames(obj$model.matrix) <- "(Intercept)"
            if(!is.null(mm_vars))
              colnames(obj$model.matrix) <- mm_vars
          }
          attr(obj$model.matrix, "binning") <- bind
        }
      }
    }
    if(smooth.construct) {
      tx <- drop.terms.bamlss(obj$terms,
        pterms = FALSE, keep.response = FALSE, data = data, specials = specials)
      sid <- unlist(attr(tx, "specials"))
      if(!length(sid))
        sid <- NULL
      smt <- NULL
      if(!is.null(sid)) {
        sterms <- sterm_labels <- attr(tx, "term.labels")[sid]
        sterms <- lapply(sterms, function(x) { eval(parse(text = x)) })
        nst <- NULL
        for(j in seq_along(sterms)) {
          sl <- sterms[[j]]$label
          if(is.null(sl))
            sl <- sterm_labels[j]
          nst <- c(nst, sl)
        }
        names(sterms) <- nst
        for(tsm in sterms) {
          if(is.null(tsm$xt))
            tsm$xt <- list()
          if(is.null(tsm$xt$binning))
            tsm$xt$binning <- binning
          if(!is.null(tsm$xt$binning)) {
            if(!is.logical(tsm$xt$binning)) {
              for(tsmt in tsm$term) {
                if(!inherits(data, "ffdf")) {
                  if(!is.factor(data[[tsmt]]))
                    data[[tsmt]] <- round(data[[tsmt]], digits = tsm$xt$binning)
                } else {
                  if(is.numeric(binning))
                    data[[tsmt]] <- round(data[[tsmt]], digits = tsm$xt$binning)
                }
              }
            }
          }
        }
        no.mgcv <- NULL
        smooth <- list()
        for(tsm in sterms) {
          special <- FALSE
          if(!is.null(tsm$special))
            special <- tsm$special
          if(!special) {
            if(is.null(tsm$xt))
              tsm$xt <- list()
            if(is.null(tsm$xt$binning))
              tsm$xt$binning <- binning
            acons <- TRUE
            if(inherits(tsm, "tensor.smooth.spec")) {
              if(!is.null(tsm$margin[[1]]$xt$center))
                acons <- tsm$margin[[1]]$xt$center
            } else {
              if(!is.null(tsm$xt$center))
                acons <- tsm$xt$center
            }
            tsm$xt$center <- acons
            tsm$xt$before <- before
            if(!is.null(tsm$xt$binning)) {
              term.names <- c(tsm$term, if(tsm$by != "NA") tsm$by else NULL)
              if(!inherits(data, "ffdf")) {
                tsm$binning <- match.index(data[, term.names, drop = FALSE])
                tsm$binning$order <- order(tsm$binning$match.index)
                tsm$binning$sorted.index <- tsm$binning$match.index[tsm$binning$order]
              } else {
                tsm$binning <- match.index.ff(data[term.names])
                attr(tsm, "ff") <- TRUE
              }
              if(!inherits(data, "ffdf")) {
                smt <- smoothCon(tsm, if(before) data[tsm$binning$nodups, term.names, drop = FALSE] else data,
                  knots, absorb.cons = if(is.null(absorb.cons)) acons else absorb.cons, sparse.cons = sparse.cons)
                smooth <- c(smooth, smt)
              } else {
                xdata <- as.data.frame(data[tsm$binning$nodups, term.names, drop = FALSE])
                if(!inherits(xdata, "data.frame")) {
                  xdata <- data.frame(xdata)
                  names(xdata) <- term.names
                }
                smt <- smoothCon(tsm, xdata,
                  knots, absorb.cons = if(is.null(absorb.cons)) acons else absorb.cons,
                  sparse.cons = sparse.cons)
                smooth <- c(smooth, smt)
              }
            } else {
              smt <- smoothCon(tsm, data, knots,
                absorb.cons = if(is.null(absorb.cons)) acons else absorb.cons,
                sparse.cons = sparse.cons)
              smooth <- c(smooth, smt)
            }
          } else {
            if(is.null(tsm$by))
              tsm$by <- "NA"
            if(inherits(tsm, "mrf.smooth.spec")) {
              if(!is.null(tsm$xt$map)) {
                vl <- levels(data[[tsm$term]])
                mapn <- names(tsm$xt$map)
                if(!all(mapn %in% vl))
                  levels(data[[tsm$term]]) <- c(vl, mapn[!(mapn %in% vl)])
                tsm$xt$polys <- as.list(tsm$xt$map)
              }
            }
            if((tsm$by != "NA") & is.factor(data[[tsm$by]])) {
              fm <- model.matrix(as.formula(paste("~ -1 +", tsm$by)), data = data)
              tlab <- tsm$label
              byvar <- tsm$by
              for(jj in 1:ncol(fm)) {
                tsm$by <- colnames(fm)[jj]
                tsm$label <- gsub(byvar, colnames(fm)[jj], tlab, fixed = TRUE)
                data[[colnames(fm)[jj]]] <- fm[, jj]
                smt2 <- smooth.construct(tsm, data, knots)
                if(inherits(tsm, "no.mgcv") | inherits(smt2, "no.mgcv")) {
                  no.mgcv <- c(no.mgcv, list(smt2))
                } else {
                  class(smt2) <- c(class(smt2), "mgcv.smooth")
                  smt <- list(smt2)
                  smooth <- c(smooth, smt)
                }
              }
            } else {
              smt2 <- smooth.construct(tsm, data, knots)
              if(inherits(tsm, "no.mgcv") | inherits(smt2, "no.mgcv")) {
                no.mgcv <- c(no.mgcv, if(!inherits(smt2, "smooth.list")) list(smt2) else smt2)
              } else {
                class(smt2) <- c(class(smt2), "mgcv.smooth")
                smt <- if(!inherits(smt2, "smooth.list")) list(smt2) else smt2
                smooth <- c(smooth, smt)
              }
            }
          }
        }
        if(length(smooth) > 0) {
          if(gam.side) {
            if(is.null(obj$model.matrix)) {
              Xp <- model.matrix(drop.terms.bamlss(obj$terms,
                sterms = FALSE, keep.response = FALSE, data = data, specials = specials), data = data)
              smooth <- try(gam.side(smooth, Xp, tol = .Machine$double.eps^.5), silent = TRUE)
            } else {
              smooth <- try(gam.side(smooth, obj$model.matrix, tol = .Machine$double.eps^.5), silent = TRUE)
            }
            if(inherits(smooth, "try-error")) {
              cat("---\n", smooth, "---\n")
              if(binning)
                stop("gam.side() produces an error when binning, try to set before = FALSE or set gam.side = FALSE!")
              else
                stop("gam.side() produces an error, try to set gam.side = FALSE!")
            }
          }
          sme <- NULL
          if(smooth.construct)
            sme <- expand.t2.smooths(smooth)
          if(is.null(sme)) {
            original.smooth <- NULL
          } else {
            original.smooth <- smooth
            smooth <- sme
            rm(sme)
          }
        }
        for(j in seq_along(smooth))
          smooth[[j]][["X.dim"]] <- ncol(smooth[[j]]$X)
        if(!is.null(no.mgcv))
          smooth <- c(smooth, no.mgcv)
        if(length(smooth))
          obj$smooth.construct <- smooth
      }
    }
    if(!is.null(obj$smooth.construct)) {
      sl <- NULL
      for(j in seq_along(obj$smooth.construct)) {
        slj <- obj$smooth.construct[[j]]$label
        if(!is.null(obj$smooth.construct[[j]]$by)) {
          if(obj$smooth.construct[[j]]$by != "NA") {
            if(grepl(pat <- paste("):", obj$smooth.construct[[j]]$by, sep = ""), slj, fixed = TRUE)) {
              slj <- gsub(pat, paste(",by=", obj$smooth.construct[[j]]$by, "):", sep = ""), slj, fixed = TRUE)
              slj <- strsplit(slj, "", fixed = TRUE)[[1]]
              if(slj[length(slj)] == ":")
                slj <- slj[-length(slj)]
              slj <- paste(slj, collapse = "")
              obj$smooth.construct[[j]]$label <- slj
            }
          }
        } else obj$smooth.construct[[j]]$by <- "NA"
        sl <- c(sl, slj)
      }
      if(length(unique(sl)) < length(sl)) {
        sld <- sl[duplicated(sl)]
        for(j in seq_along(sld)) {
          for(jj in which(sl == sld[j])) {
            clj <- class(obj$smooth.construct[[jj]])
            clj <- strsplit(clj, ".", fixed = TRUE)[[1]][1]
            if(clj == "random")
              clj <- "re"
            sl[jj] <- paste(sl[jj], clj, sep = ":")
          }
        }
      }
      names(obj$smooth.construct) <- sl
    }
    if(!is.null(drop)) {
      take <- c("model.matrix", "smooth.construct")[c(model.matrix, smooth.construct)]
      obj[!(names(obj) %in% take)] <- NULL
    }

    obj
  }

  if(!all(c("formula", "fake.formula") %in% names(formula))) {
    for(j in seq_along(formula)) {
      if(!all(c("formula", "fake.formula") %in% names(formula[[j]]))) {
        for(i in seq_along(formula[[j]])) {
          formula[[j]][[i]] <- assign.design(formula[[j]][[i]],
            if(i > 1) duplicated(data[, all.vars(formula[[j]][[i]]$fake.formula), drop = FALSE]) else NULL)
        }
      } else formula[[j]] <- assign.design(formula[[j]])
    }
  } else formula <- assign.design(formula)

  if((!all(c("formula", "fake.formula") %in% names(formula))) & smooth.construct) {
    for(i in seq_along(formula)) {
      if(!all(c("formula", "fake.formula") %in% names(formula[[i]]))) {
        for(j in seq_along(formula[[i]])) {
          if(!is.null(formula[[i]][[j]]$smooth.construct)) {
            for(k in seq_along(formula[[i]][[j]]$smooth.construct)) {
              if(is.null(formula[[i]][[j]]$smooth.construct[[k]]$fit.fun))
                formula[[i]][[j]]$smooth.construct[[k]]$fit.fun <- make.fit.fun(formula[[i]][[j]]$smooth.construct[[k]])
              if(is.null(formula[[i]][[j]]$smooth.construct[[k]]$prior)) {
                priors <- make.prior(formula[[i]][[j]]$smooth.construct[[k]])
                formula[[i]][[j]]$smooth.construct[[k]]$prior <- priors$prior
                formula[[i]][[j]]$smooth.construct[[k]]$grad <- priors$grad
                formula[[i]][[j]]$smooth.construct[[k]]$hesss <- priors$hess
              }
            }
          }
        }
      } else {
        if(!is.null(formula[[i]]$smooth.construct)) {
          for(j in seq_along(formula[[i]]$smooth.construct)) {
            if(is.null(formula[[i]]$smooth.construct[[j]]$fixed))
              formula[[i]]$smooth.construct[[j]]$fixed <- FALSE
            if(length(formula[[i]]$smooth.construct[[j]]$S)) {
              for(sj in seq_along(formula[[i]]$smooth.construct[[j]]$S)) {
                if(!is.list(formula[[i]]$smooth.construct[[j]]$S[[sj]]) & !is.function(formula[[i]]$smooth.construct[[j]]$S[[sj]])) {
                  nc <- ncol(formula[[i]]$smooth.construct[[j]]$S[[sj]])
                  formula[[i]]$smooth.construct[[j]]$S[[sj]] <- formula[[i]]$smooth.construct[[j]]$S[[sj]] + diag(1e-05, nc, nc)
                }
              }
            }
            if(is.null(formula[[i]]$smooth.construct[[j]]$fit.fun))
              formula[[i]]$smooth.construct[[j]]$fit.fun <- make.fit.fun(formula[[i]]$smooth.construct[[j]])
            if(is.null(formula[[i]]$smooth.construct[[j]]$prior)) {
              priors <- make.prior(formula[[i]]$smooth.construct[[j]])
              formula[[i]]$smooth.construct[[j]]$prior <- priors$prior
              formula[[i]]$smooth.construct[[j]]$grad <- priors$grad
              formula[[i]]$smooth.construct[[j]]$hess <- priors$hess
            }
          }
        }
      }
    }
  } else {
    if(!is.null(formula$smooth.construct)) {
      for(j in seq_along(formula$smooth.construct)) {
        if(is.null(formula$smooth.construct[[j]]$fixed))
          formula$smooth.construct[[j]]$fixed <- FALSE
        if(length(formula[[i]]$smooth.construct[[j]]$S)) {
          for(sj in seq_along(formula$smooth.construct[[j]]$S)) {
            if(!is.list(formula$smooth.construct[[j]]$S[[sj]]) & !is.function(formula$smooth.construct[[j]]$S[[sj]])) {
              nc <- ncol(formula$smooth.construct[[j]]$S[[sj]])
              formula$smooth.construct[[j]]$S[[sj]] <- formula$smooth.construct[[j]]$S[[sj]] + diag(1e-05, nc, nc)
            }
          }
        }
        if(is.null(formula$smooth.construct[[j]]$fit.fun))
          formula$smooth.construct[[j]]$fit.fun <- make.fit.fun(formula$smooth.construct[[j]])
        if(is.null(formula$smooth.construct[[j]]$prior)) {
          priors <- make.prior(formula$smooth.construct[[j]])
          formula$smooth.construct[[j]]$prior <- priors$prior
          formula$smooth.construct[[j]]$grad <- priors$grad
          formula$smooth.construct[[j]]$hess <- priors$hess
        }
      }
    }
  }

  attr(formula, "specials") <- NULL
  attr(formula, ".Environment") <- NULL
  class(formula) <- "list"
  if(!is.null(drop)) {
    if(drop & (length(formula) < 2))
      formula <- formula[[1]]
  }

  return(formula)
}


## Functions for sparse matrices.
sparse.matrix.index <- function(x, ...)
{
  if(is.null(dim(x)))
    return(NULL)
  index <- apply(x, 1, function(x) {
    which(x != 0)
  })
  if(length(index) < 1)
    return(NULL)
  if(is.list(index)) {
    n <- max(sapply(index, length))
    index <- lapply(index, function(x) {
      if((nx <- length(x)) < n)
        x <- c(x, rep(-1L, length = n - nx))
      x
    })
    index <- do.call("rbind", index)
  } else {
    index <- if(is.null(dim(index))) {
      matrix(index, ncol = 1)
    } else t(index)
  }
  storage.mode(index) <- "integer"
  index
}


## Bandwidth minimization permutation.
sparse.matrix.ordering <- function(x, ...)
{
  x <- as.spam(x)
  i <- try(ordering(chol.spam(x)), silent = TRUE)
  if(inherits(i, "try-error"))
    i <- 1:nrow(x)
  return(i)
}


## Setup sparse indices for various algorithms.
sparse.setup <- function(x, S = NULL, ...)
{
  symmetric <- nrow(x) == ncol(x)
  index.matrix <- sparse.matrix.index(x, ...)
  if(!symmetric)
    x <- crossprod(x)
  if(!is.null(S)) {
    if(!is.list(S))
      S <- list(S)
    for(j in seq_along(S)) {
      x <- x + if(length(S[[j]]) < 1) 0 else { if(is.function(S[[j]])) S[[j]](c("b" = rep(0, attr(S[[j]], "npar")))) else S[[j]] }
    }
  }
  index.crossprod <- if(!symmetric) sparse.matrix.index(x, ...) else NULL
  setup <- list(
    "matrix" = index.matrix,
    "crossprod" = index.crossprod
  )
  if(!is.null(index.crossprod)) {
    # make block.index only if coefficients do not overlap
    tmp <- setup$crossprod[!duplicated(setup$crossprod), , drop = FALSE]
    l <- nrow(tmp)
    if(any(unique(tmp[duplicated(tmp, MARGIN = 0)]) > 0)){
      return(setup)
    } else {
      if((l > 1) & (l <= nrow(setup$crossprod))) {
        setup$block.index <- split(tmp, 1:l)
        setup$block.index <- lapply(1:l, function(i) setup$block.index[[i]][setup$block.index[[i]] > 0])
        setup$is.diagonal <- all(sapply(setup$block.index, length) == 1)
      }
    }
  }
  return(setup)
}


## Sparse cholesky decomposition,
## returns the lower triangle.
sparse.chol <- function(x, index, ...)
{
  if(all(dim(x) < 2))
    return(sqrt(x))

  imat <- index[["matrix"]]
  p <- index[["ordering"]]

  # imat: index matrix of a[p,p]
  # ??? check for positive definiteness?
  n <- nrow(x)
  l <- matrix(0, nrow = n, ncol = n)
  
  # First column simplified (no elements to sum up)
  l[1, 1] <- (x[p[1], p[1]])^0.5
  for(i in imat[1,][imat[1,]>1]) {
    l[i, 1] <- x[p[i], p[1]] / l[1, 1]
  }
  c <- 1
  for(j in p[2:(n-1)]) {
    c <- c + 1
    l[c, c] <- (x[j, j] - sum(l[c, 1:(c - 1)]^2))^0.5
    # use only non-zero entries in lower subdiagonal
    for(i in imat[c,][imat[c,] > c]) {   
      l[i, c] <- (x[p[i], p[c]] -
                    sum(l[i, 1:(c - 1)] * l[c, 1:(c - 1)])) / l[c, c]
    }
  }
  # last column simplified: no subdiagonal - maybe still leave in loop?
  l[n, n] <- (x[p[n], p[n]] - sum(l[n, 1:(n - 1)]^2))^0.5
  j <- c(1:n)[p]
  return(l[j,j])
}


## Sparse forward substitution.
## L %*% x = bn
## with bn = t(P) %*% b
sparse.forwardsolve <- function(l, x, index, ...)
{
  if(all(dim(l) < 2))
    return(x / l)
  imat <- index[["matrix"]]
  p <- index[["ordering"]]
  n <- ncol(l)
  Pt <- diag(n)[p,]
  xn <- Pt %*% x
  y <- matrix(rep(NA, n), ncol=1)
  y[1] <- xn[1]/l[1, 1]
  for(i in 2:n){
    y[i] <- xn[i]/l[i, i]  
    if(max(imat[i,]) > 0){
      y[i] <- y[i] - sum(l[i, imat[i,][imat[i,] > 0] ] * y[imat[i,][imat[i,] > 0]])/l[i, i]
    }
  }
  return(y)
}


## Sparse backward substitution.
## t(L) %*% xn = x,
## with x = P %*% xn.
sparse.backsolve <- function(r, x, index = NULL, ...)
{
  if(all(dim(x) < 2))
    return(r / x)
  imat <- index[["matrix"]]
  p <- index[["ordering"]]
  n <- ncol(r)
  P <- diag(n)[,p]
  xn <- rep(NA, n)
  xn[n] <- x[n]/r[n,n]
  for(i in (n-1):1){
    xn[i] <- x[i]/r[i,i]
    if(max(imat[i,]) > 0){
      xn[i] <- xn[i] - sum(r[imat[i,][imat[i,] > 0],i ] * xn[imat[i,][imat[i,] > 0]])/r[i, i]
    }
  }
  x <- P %*% xn
  return(x)
}


## Sparse matrix solve.
sparse.solve <- function(a, b, index, ...)
{
  if(all(dim(a) < 2))
    return(b / a)
  id <- if(!("crossprod" %in% names(index))) "matrix" else "crossprod"
  L <- sparse.chol(a, index = list("matrix" = index[[id]], "ordering" = index[["ordering"]]), ...)
  y <- sparse.forwardsolve(L, b, index = list("matrix" = index$forward, "ordering" = index[["ordering"]]), ...)
  z <- sparse.backsolve(L, y, list("matrix" = index$backward, "ordering" = index[["ordering"]]), ...)
  return(z)
}


## Computation of fitted values with index matrices.
sparse.matrix.fit.fun <- function(X, b, index = NULL)
{
  if(!is.null(index)) {
    if(nrow(index) != nrow(X))
      return(drop(X %*% b))
  }
  fit <- if(inherits(X, "dgCMatrix") | is.null(index) | inherits(X, "Matrix")) {
    drop(X %*% b)
  } else .Call("sparse_matrix_fit_fun", X, b, index, PACKAGE = "bamlss")
  return(fit)
}


## The model term fitting function.
make.fit.fun <- function(x, type = 1)
{
  if(inherits(x, "nnet.smooth"))
    return(x$fit.fun)

  ff <- function(X, b, expand = TRUE, no.sparse.setup = FALSE) {
    if(!is.null(names(b))) {
      b <- if(!is.null(x$pid)) b[x$pid$b] else get.par(b, "b")
    }
    if(inherits(X, "spam") | inherits(X, "Matrix")) {
      f <- as.matrix(X %*% b)
    } else {
      what <- if(type < 2) "matrix" else "grid.matrix"
      f <- if(is.null(x$sparse.setup[[what]]) | no.sparse.setup) {
        drop(X %*% b)
      } else sparse.matrix.fit.fun(X, b, x$sparse.setup[[what]])
    }
    if(!is.null(x$binning$match.index) & expand) {
      if(inherits(x$binning$match.index, "ff")) {
        ## f <- as.ff(f)
        ## FIXME: ff support!
        return(f[x$binning$match.index])
      }
      f <- f[x$binning$match.index]
    }
    if(!is.null(x$xt$force.center))
      f <- f - mean(f, na.rm = TRUE)
    return(as.numeric(f))
  }
  return(ff)
}


## The prior function.
make.prior <- function(x, sigma = 0.1)
{
  prior <- NULL
  if(!is.null(x$xt$prior)) {
    prior <- x$xt$prior
    if(is.character(x$xt$prior)) {
      prior <- tolower(prior)
      if(!(prior %in% c("ig", "hc", "sd", "hn", "hn.lasso")))
        stop(paste('smoothing variance prior "', prior, '" not supported!', sep = ''))
    }
  } else {
    prior <- "ig"
  }

  if(!is.null(x$margin)) {
    if(!is.null(x$margin[[1]]$xt)) {
      xt <- x$margin[[1]]$xt
      if(is.null(names(xt)) & (length(xt) == 1)) {
        if(length(xt[[1]]) > 0)
          xt <- xt[[1]]
      }
      x$xt <- c(x$xt, xt)
    }
  }

  if(!is.function(prior)) {
    rval <- list()

    a <- if(is.null(x$xt[["a"]])) {
      if(is.null(x[["a"]])) 1e-04 else x[["a"]]
    } else x$xt[["a"]]
    b <- if(is.null(x$xt[["b"]])) {
      if(is.null(x[["b"]])) 1e-04 else x[["b"]]
    } else x$xt[["b"]]
    theta <- if(is.null(x$xt[["theta"]])) {
      x[["theta"]]
    } else x$xt[["theta"]]

    if(is.null(theta)) {
      theta <- switch(prior,
        "sd" = 0.00877812,
        "hc" = 0.01034553,
        "hn" = 0.1457644
      )
    }

    fixed <- if(is.null(x$fixed)) FALSE else x$fixed

    igs <- log((b^a)) - log(gamma(a))
    var_prior_fun <- switch(prior,
      "ig" = function(tau2) { igs + (-a - 1) * log(tau2) - b / tau2 },
      "hc" = function(tau2) { -log(1 + tau2 / (theta^2)) - 0.5 * log(tau2) - log(theta^2) },
      "sd" = function(tau2) { -0.5 * log(tau2) + 0.5 * log(theta) - (tau2 / theta)^(0.5) },
      "hn" = function(tau2) { -0.5 * log(tau2) - tau2 / (2 * theta^2) },
      "hn.lasso0" = function(tau2) { -0.2257913 - log(sigma) - tau2^2/(2 * sigma^2) },
      "hn.lasso" = function(tau2) {
        theta <- sqrt(pi) / (sigma * sqrt(2))
        log(2 * theta / pi) - (tau2^2 * theta^2) / pi
      }
    )
    
    rval$prior <- function(parameters) {
      if(is.null(x$pid)) {
        if(!is.null(names(parameters))) {
          gamma <- get.par(parameters, "b")
          tau2 <- get.par(parameters, "tau2")
        } else {
          gamma <- parameters
          tau2 <- numeric(0)
        }
      } else {
        gamma <- parameters[x$pid$b]
        tau2 <-  parameters[x$pid$tau2]
      }
      if(fixed | !length(tau2)) {
        lp <- sum(dnorm(gamma, sd = 1000, log = TRUE))
      } else {
        if(length(tau2) < 2) {
          K <- if(is.function(x$S[[1]])) x$S[[1]](c(parameters, x$fixed.hyper)) else x$S[[1]]
          if(is.null(x$rank))
            x$rank <- qr(K)$rank
          lp <- -log(tau2) * x$rank / 2 + drop(-0.5 / tau2 * t(gamma) %*% K %*% gamma) + var_prior_fun(tau2)
        } else {
          ld <- 0
          P <- if(inherits(x$X, "Matrix")) Matrix(0, ncol(x$X), ncol(x$X)) else 0
          for(j in seq_along(tau2)) {
            P <- P + 1 / tau2[j] * if(is.function(x$S[[j]])) x$S[[j]](c(parameters, x$fixed.hyper)) else x$S[[j]]
            ld <- ld + var_prior_fun(tau2[j])
          }
          ##lp <- dmvnorm(gamma, sigma = matrix_inv(P), log = TRUE) + ld
          dP <- determinant(P, logarithm = TRUE)
          dP <- as.numeric(dP$modulus) * as.numeric(dP$sign)
          lp <- 0.5 * dP - 0.5 * (t(gamma) %*% P %*% gamma) + ld
        }
      }
      return(as.numeric(lp))
    }

    attr(rval$prior, "var_prior") <- prior

    rval$grad <- function(score = NULL, parameters, full = TRUE) {
      gamma <- get.par(parameters, "b")
      tau2 <-  get.par(parameters, "tau2")
      grad2 <- NULL
      if(x$fixed | !length(tau2)) {
        grad <- rep(0, length(gamma))
      } else {
        grad <- 0; grad2 <- NULL
        for(j in seq_along(tau2)) {
          tauS <- -1 / tau2[j] * if(is.function(x$S[[j]])) x$S[[j]](c(parameters, x$fixed.hyper)) else x$S[[j]]
          grad <- grad + tauS %*% gamma
          if(full & !is.null(tau2[j])) {
            grad2 <- c(grad2, drop(-x$rank[j] / (2 * tau2[j]) - 1 / (2 * tau2[j]^2) * (if(is.function(x$S[[j]])) x$S[[j]](c(parameters, x$fixed.hyper)) else x$S[[j]]) %*% gamma + (-x$a - 1) / tau2[j] + x$b / (tau2[j]^2)))
            x$X <- cbind(x$X, 0)
          }
          grad <- drop(grad)
        }
      }
      if(!is.null(score)) {
        grad <- if(!is.null(x$xbin.ind)) {
          drop(crossprod(x$X[x$xbin.ind, , drop = FALSE], score)) + c(grad, grad2)
        } else drop(crossprod(x$X, score)) + c(grad, grad2)
      } else grad <- c(grad, grad2)
      return(grad)
    }

    rval$hess <- function(score = NULL, parameters, full = FALSE) {
      tau2 <- get.par(parameters, "tau2")
      if(x$fixed | !length(tau2)) {
        k <- length(get.par(parameters, "b"))
        hx <- matrix(0, k, k)
      } else {
        hx <- 0
        for(j in seq_along(tau2)) {
          hx <- hx + (1 / tau2[j]) * if(is.function(x$S[[j]])) x$S[[j]](c(parameters, x$fixed.hyper)) else x$S[[j]]
        }
      }
      return(hx)
    }

    return(rval)
  } else {
    return(prior(x))
  }
}


## Fast block diagonal crossproduct with weights.
do.XWX <- function(x, w, index = NULL)
{
  if(is.null(index) | inherits(x, "dgCMatrix")) {
    rval <- crossprod(x / w, x)
  } else {
    if(is.null(dim(index)))
      index <- matrix(index, ncol = 1)
    rval <- .Call("do_XWX", x, w, index, PACKAGE = "bamlss")
  }
  rval
}


## Get the model.frame.
model.frame.bamlss <- model.frame.bamlss.frame <- function(formula, ...) 
{
  dots <- list(...)
  nargs <- dots[match(c("data", "na.action", "subset"), names(dots), 0L)]
  mf <- if(length(nargs) || is.null(formula$model.frame)) {
    fcall <- formula$call
    fcall[[1L]] <- quote(bamlss.model.frame)
    fcall[names(nargs)] <- nargs
    env <- environment(formula$formula)
    if(is.null(env))
      env <- parent.frame()
    ft <- eval(fcall[["formula"]], env)
    if(!is.null(attr(ft, "orig.formula"))) {
      fcall["formula"] <- parse(text = paste("attr(", fcall["formula"], ", 'orig.formula')", sep = ""))
    }
    nf <- names(fcall)
    nf <- nf[!(nf %in% names(formals(bamlss.model.frame)))]
    nf <- nf[nchar(nf) > 0]
    fcall[nf] <- NULL
    fcall["drop.unused.levels"] <- FALSE
    if(is.null(fcall["family"]))
      fcall["family"] <- parse(text = "gaussian_bamlss()")
    eval(fcall, env)
  } else formula$model.frame
  mf
}


## Search for parts in models, optionally extract.
model.search <- function(x, what, model = NULL, part = c("x", "formula", "terms"),
  extract = FALSE, drop = FALSE)
{
  if(!inherits(x, "bamlss.formula") & !inherits(x, "bamlss.frame"))
    stop("x must be a 'bamlss.formula' or 'bamlss.frame' object!")
  part <- match.arg(part)
  if(is.null(x[[part]]))
    return(FALSE)
  x <- model.terms(x, model = model, part = part)
  elmts <- c("formula", "fake.formula")
  nx <- names(x)
  rval <- list()
  for(i in nx) {
    if(!all(elmts %in% names(x[[i]]))) {
      rval[[i]] <- list()
      for(j in names(x[[i]])) {
        rval[[i]][[j]] <- if(is.null(x[[i]][[j]][[what]])) FALSE else TRUE
        if(extract & rval[[i]][[j]])
          rval[[i]][[j]] <- x[[i]][[j]][[what]]
      }
    } else {
      rval[[i]] <- if(is.null(x[[i]][[what]])) FALSE else TRUE
      if(extract & rval[[i]])
        rval[[i]] <- x[[i]][[what]]
    }
  }
  if(!extract) {
    rval <- unlist(rval)
  } else {
    if(drop & (length(rval) < 2))
      rval <- rval[[1]]
  }
  rval
}


## Wrapper for design construct extraction.
extract.design.construct <- function(object, data = NULL,
  knots = NULL, model = NULL, drop = TRUE, what = c("model.matrix", "smooth.construct"),
  specials = NULL, ...)
{
  if(!inherits(object, "bamlss.frame") & !inherits(object, "bamlss.formula") & !inherits(object, "bamlss.terms"))
    stop("object must be a 'bamlss.frame', 'bamlss.formula' or 'bamlss.terms' object!")
  what <- match.arg(what)
  model.matrix <- what == "model.matrix"
  smooth.construct <- what == "smooth.construct"
  if(inherits(object, "bamlss.frame")) {
    if(!is.null(data)) {
      object$model.frame <- NULL
      object <- design.construct(object, data = data, knots = knots,
        model.matrix = model.matrix, smooth.construct = smooth.construct,
        model = model, drop = drop, specials = specials, ...)
    } else {
      if(!all(model.search(object, what, model, part = "x"))) {
        object <- design.construct(object, model.matrix = model.matrix,
          smooth.construct = smooth.construct, model = model, drop = TRUE,
          specials = specials, ...)
      } else {
        object <- model.search(object, what, model, extract = TRUE, drop = drop, part = "x")
      }
    }
  } else {
    if(is.null(data))
      stop("argument data is missing!")
    object <- design.construct(object, data = data, knots = knots,
      model.matrix = model.matrix, smooth.construct = smooth.construct, model = model,
      drop = drop, specials = specials, ...)
  }
  if(!is.null(drop)) {
    if(length(object) & drop & (length(object) < 2))
      object <- object[[1]]
  }
  if(!length(object))
    return(NULL)
  mostattributes(object) <- NULL
  attr(object, "orig.formula") <- NULL
  if(what == "model.matrix") {
    if(is.list(object)) {
      for(j in seq_along(object)) {
        if(is.list(object[[j]])) {
          if((length(object[[j]]) < 2) & (names(object[[j]]) == "model.matrix")) {
            object[[j]] <- object[[j]]$model.matrix
          }
        }
      }
    }
  }
  return(object)
}


## Model matrix extractor.
model.matrix.bamlss.frame <- model.matrix.bamlss.formula <- model.matrix.bamlss.terms <- function(object, data = NULL, model = NULL, drop = TRUE, scale.x = FALSE, ...)
{
  extract.design.construct(object, data = data,
    knots = NULL, model = model, drop = drop, what = "model.matrix",
    scale.x = scale.x)
}


## Extract smooth constructs.
smooth.construct <- function(object, data, knots, ...)
{
  UseMethod("smooth.construct")
}

smooth.construct.bamlss.frame <- smooth.construct.bamlss.formula <- smooth.construct.bamlss.terms <- function(object, data = NULL, knots = NULL, model = NULL, drop = TRUE, ...)
{
  extract.design.construct(object, data = data,
    knots = knots, model = model, drop = drop, what = "smooth.construct")
}


## Extract/initialize parameters.
parameters <- function(x, model = NULL, start = NULL, fill = c(0, 0.0001),
  list = FALSE, simple.list = FALSE, extract = FALSE, ...)
{
  if(inherits(x, "bamlss") | extract) {
    if(!is.null(x$parameters)) {
      if(is.null(model)) {
        if(list) {
          return(x$parameters)
        } else {
          if(inherits(x$parameters, "data.frame") | inherits(x$parameters, "matrix")) {
            args <- list(...)
            mstop <- if(is.null(args$mstop)) nrow(x$parameters) else args$mstop
            return(x$parameters[mstop, ])
          } else return(unlist(x$parameters))
        }
      } else {
        if(is.list(x$parameters)) {
          if(list) return(x$parameters[model]) else return(unlist(x$parameters[model]))
        } else {
          if(!is.character(model))
            model <- names(x$terms)[model]
          rp <- grep(paste(model, ".", sep = ""), names(x$parameters), fixed = TRUE, value = TRUE)
          if(inherits(x$parameters, "data.frame") | inherits(x$parameters, "matrix")) {
            rp <- grep(paste(model, ".", sep = ""), colnames(x$parameters), fixed = TRUE, value = TRUE)
            args <- list(...)
            mstop <- if(is.null(args$mstop)) nrow(x$parameters) else args$mstop
            return(x$parameters[mstop, rp])
          } else return(x$parameters[rp])
        }
      }
    }
  }
  if(inherits(x, "bamlss.frame")) {
    if(is.null(x$x)) {
      x <- design.construct(x, data = x$model.frame,
        knots = x$knots, model.matrix = TRUE, smooth.construct = TRUE, model = NULL)
    } else x <- x$x
  }
  fill <- rep(fill, length.out = 2)
  if(!is.null(start)) {
    if(is.list(start))
      start <- unlist(start)
  }
  par <- list()
  for(i in names(x)) {
    par[[i]] <- list()
    if(!all(c("formula", "fake.formula") %in% names(x[[i]]))) {
      for(j in names(x[[i]])) {
        par[[i]][[j]] <- list()
        if(!is.null(x[[i]][[j]]$model.matrix)) {
          nc <- ncol(x[[i]][[j]]$model.matrix)
          if(simple.list) {
            par[[i]][[j]]$p <- fill[1]
          } else {
            par[[i]][[j]]$p <- rep(fill[1], length = nc)
            if(is.null(cn <- colnames(x[[i]][[j]]$model.matrix)))
              cn <- paste("b", 1:nc, sep = "")
            names(par[[i]][[j]]$p) <- cn
            if(!is.null(start)) {
              if(length(ii <- grep(paste(i, j, "p", sep = "."), names(start), fixed = TRUE))) {
                spar <- start[ii]
                spn <- names(spar)
                cn2 <- paste(i, j, "p", cn, sep = ".")
                take <- which(spn %in% cn2)
                if(length(take)) {
                  par[[i]][[j]]$p[which(cn2 %in% spn)] <- spar[take]
                }
              }
            }
          }
        }
        if(!is.null(x[[i]][[j]]$smooth.construct)) {
          par[[i]][[j]]$s <- list()
          for(k in names(x[[i]][[j]]$smooth.construct)) {
            if(simple.list) {
              par[[i]][[j]]$s[[k]] <- fill[1]
            } else {
              if(!is.null(x[[i]][[j]]$smooth.construct[[k]]$rand)) {
                tpar1 <- rep(fill[1], ncol(x[[i]][[j]]$smooth.construct[[k]]$rand$Xr))
                tpar2 <- rep(fill[1], ncol(x[[i]][[j]]$smooth.construct[[k]]$Xf))
                cn1 <- colnames(x[[i]][[j]]$smooth.construct[[k]]$rand$Xr)
                cn2 <- colnames(x[[i]][[j]]$smooth.construct[[k]]$Xf)
                if(is.null(cn1))
                  cn1 <- paste("b", 1:length(tpar1), ".re", sep = "")
                if(is.null(cn2))
                  cn2 <- paste("b", 1:length(tpar2), ".fx", sep = "")
                names(tpar1) <- cn1
                names(tpar2) <- cn2
                tpar <- c(tpar1, tpar2)
              } else {
                nfill <- if(is.null(x[[i]][[j]]$smooth.construct[[k]]$special.npar)) {
                  ncol(x[[i]][[j]]$smooth.construct[[k]]$X)
                } else x[[i]][[j]]$smooth.construct[[k]]$special.npar
                tpar <- rep(fill[1], nfill)
                cn <- colnames(x[[i]][[j]]$smooth.construct[[k]]$X)
                if(is.null(cn))
                  cn <- paste("b", 1:length(tpar), sep = "")
                names(tpar) <- cn
              }
              if(length(x[[i]][[j]]$smooth.construct[[k]]$S)) {
                tpar3 <- NULL
                for(kk in seq_along(x[[i]][[j]]$smooth.construct[[k]]$S)) {
                  tpar3 <- c(tpar3, fill[2])
                }
                names(tpar3) <- paste("tau2", 1:length(tpar3), sep = "")
                tpar <- c(tpar, tpar3)
              }
              par[[i]][[j]]$s[[k]] <- tpar
              if(!is.null(start)) {
                if(length(ii <- grep(paste(i, j, "s", k, sep = "."), names(start), fixed = TRUE))) {
                  spar <- start[ii]
                  cn <- names(par[[i]][[j]]$s[[k]])
                  if(length(tau2 <- grep("tau2", names(spar)))) {
                    tau2 <- spar[tau2]
                    if(length(jj <- grep("tau2", cn, fixed = TRUE))) {
                      tau2 <- rep(tau2, length.out = length(jj))
                      par[[i]][[j]]$s[[k]][jj] <- tau2
                    }
                  }
                  if(any(b <- !grepl("tau2", names(spar)))) {
                    b <- spar[b]
                    if(any(jj <- !grepl("tau2", cn, fixed = TRUE))) {
                      b <- rep(b, length.out = sum(jj))
                      par[[i]][[j]]$s[[k]][jj] <- b
                    }
                  }
                }
              }
            }
          }
        }
      }
    } else {
      if(!is.null(x[[i]]$model.matrix)) {
        if(ncol(x[[i]]$model.matrix) > 0) {
          if(simple.list) {
            par[[i]]$p <- fill[1]
          } else {
            nc <- ncol(x[[i]]$model.matrix)
            par[[i]]$p <- rep(fill[1], length = nc)
            if(is.null(cn <- colnames(x[[i]]$model.matrix)))
              cn <- paste("b", 1:nc, sep = "")
            names(par[[i]]$p) <- cn
            if(!is.null(start)) {
              if(length(ii <- grep(paste(i, "p", sep = "."), names(start), fixed = TRUE))) {
                spar <- start[ii]
                spn <- names(spar)
                cn2 <- paste(i, "p", cn, sep = ".")
                take <- which(spn %in% cn2)
                if(length(take)) {
                  par[[i]]$p[which(cn2 %in% spn)] <- spar[take]
                }
              }
            }
          }
        }
      }
      if(!is.null(x[[i]]$smooth.construct)) {
        par[[i]]$s <- list()
        for(k in names(x[[i]]$smooth.construct)) {
          re.effect <- !is.null(x[[i]]$smooth.construct[[k]]$rand)
          if(re.effect) {
            if(is.logical(x[[i]]$smooth.construct[[k]]$rand))
              re.effect <- FALSE
          }
          if(re.effect) {
            tpar1 <- rep(fill[1], ncol(x[[i]]$smooth.construct[[k]]$rand$Xr))
            tpar2 <- rep(fill[1], ncol(x[[i]]$smooth.construct[[k]]$Xf))
            cn1 <- colnames(x[[i]]$smooth.construct[[k]]$rand$Xr)
            cn2 <- colnames(x[[i]]$smooth.construct[[k]]$Xf)
            if(is.null(cn1))
              cn1 <- paste("b", 1:length(tpar1), ".re", sep = "")
            if(is.null(cn2))
              cn2 <- paste("b", 1:length(tpar2), ".fx", sep = "")
            names(tpar1) <- cn1
            names(tpar2) <- cn2
            tpar <- c(tpar1, tpar2)
          } else {
            nfill <- if(is.null(x[[i]]$smooth.construct[[k]]$special.npar)) {
              ncol(x[[i]]$smooth.construct[[k]]$X)
            } else x[[i]]$smooth.construct[[k]]$special.npar
            tpar <- rep(fill[1], nfill)
            cn <- colnames(x[[i]]$smooth.construct[[k]]$X)
            if(is.null(cn))
              cn <- paste("b", 1:length(tpar), sep = "")
            if(!simple.list)
              names(tpar) <- cn
          }
          if(length(x[[i]]$smooth.construct[[k]]$S)) {
            tpar3 <- NULL
            for(kk in seq_along(x[[i]]$smooth.construct[[k]]$S)) {
              tpar3 <- c(tpar3, fill[2])
            }
            if(!simple.list)
              names(tpar3) <- paste("tau2", 1:length(tpar3), sep = "")
            tpar <- c(tpar, tpar3)
          }
          if(!is.null(x[[i]]$smooth.construct[[k]]$special.mpar)) {
            tpar <- x[[i]]$smooth.construct[[k]]$special.mpar()
          }
          par[[i]]$s[[k]] <- tpar
          if(!is.null(start)) {
            if(length(ii <- grep(paste(i, "s", k, sep = "."), names(start), fixed = TRUE))) {
              spar <- start[ii]
              cn <- names(par[[i]]$s[[k]])
              if(length(tau2 <- grep("tau2", names(spar)))) {
                tau2 <- spar[tau2]
                if(length(jj <- grep("tau2", cn, fixed = TRUE))) {
                  tau2 <- rep(tau2, length.out = length(jj))
                  par[[i]]$s[[k]][jj] <- tau2
                }
              }
              if(any(b <- !grepl("tau2", names(spar)))) {
                b <- spar[b]
                if(any(jj <- !grepl("tau2", cn, fixed = TRUE))) {
                  b <- rep(b, length.out = sum(jj))
                  par[[i]]$s[[k]][jj] <- b
                }
              }
            }
          }
        }
      }
    }
  }
  if(!is.null(model))
    par <- par[model]
  if(!list)
    par <- unlist(par)
  return(par)
}


par2list <- function(x)
{
  xl <- list()
  nx <- names(x)
  npl <- strsplit(nx, ".", fixed = TRUE)
  np <- unique(sapply(npl, function(x) { x[1] }))
  for(j in np) {
    xl[[j]] <- list()
    tmp <- grep(paste(j, ".", sep = ""), nx, value = TRUE)
    tmp2 <- unique(sapply(strsplit(tmp, ".", fixed = TRUE), function(x) { x[2] }))
    for(jj in tmp2) {
      tmp3 <- grep(paste(j, jj, sep = "."), nx, fixed = TRUE, value = TRUE)
      if(jj == "p") {
        xl[[j]][[jj]] <- x[tmp3]
        names(xl[[j]][[jj]]) <- gsub(paste(j, ".", jj, ".", sep = ""), "", names(xl[[j]][[jj]]), fixed = TRUE)
      } else {
        tmp4 <- grep(paste(j, jj, sep = "."), nx, fixed = TRUE, value = TRUE)
        tmp4 <- unique(sapply(strsplit(tmp4, ".", fixed = TRUE), function(x) { x[3] }))
        xl[[j]][[jj]] <- list()
        for(jjj in tmp4) {
          xl[[j]][[jj]][[jjj]] <- x[grep(paste(j, ".", jj, ".", jjj, ".", sep = ""), nx, fixed = TRUE)]
          names(xl[[j]][[jj]][[jjj]]) <- gsub(paste(j, ".", jj, ".", jjj, ".", sep = ""), "", names(xl[[j]][[jj]][[jjj]]), fixed = TRUE)
        }
      }
    }
  }
  xl
}


## Main bamlss().
bamlss <- function(formula, family = "gaussian", data = NULL, start = NULL, knots = NULL,
  weights = NULL, subset = NULL, offset = NULL, na.action = na.omit, contrasts = NULL,
  reference = NULL, transform = NULL, optimizer = NULL, sampler = NULL, samplestats = NULL, results = NULL,
  cores = NULL, sleep = NULL, combine = TRUE, model = TRUE, x = TRUE, light = FALSE, ...)
{
  ## Search for functions in family object.
  family <- bamlss.family(family, ...)
  if(!is.null(family$transform) & is.null(transform))
    transform <- family$transform
  if(!is.null(family$optimizer) & is.null(optimizer))
    optimizer <- family$optimizer
  if(!is.null(family$sampler) & is.null(sampler))
    sampler <- family$sampler
  if(!is.null(family$samplestats) & is.null(samplestats))
    samplestats <- family$samplestats
  if(!is.null(family$results) & is.null(results))
    results <- family$results

  ## Switch for light variant.
  if(light) {
    results <- FALSE
    samplestats <- FALSE
  }

  ## Setup all processing functions.
  foo <- list("transform" = transform, "optimizer" = optimizer,
    "sampler" = sampler, "samplestats" = samplestats, "results" = results)
  nf <- names(foo)
  default_fun <- c("no.transform", "bfit", "GMCMC", "samplestats", "results.bamlss.default")
  functions <- list()
  for(j in 1:length(foo)) {
    if(is.null(foo[[nf[j]]])) {
      foo[[nf[j]]] <- if(default_fun[j] != "no.transform") {
        get(default_fun[j], envir = asNamespace("bamlss"))
      } else FALSE
    }
    if(is.list(foo[[nf[j]]])) {
      args <- foo[[nf[j]]]
      fun <- default_fun[j]
      functions[[nf[j]]] <- function(x, ...) {
        args <- c(args, list(...))
        args$x <- x
        do.call(fun, args)
      }
    } else functions[[nf[j]]] <- foo[[nf[j]]]
    if(!is.function(functions[[nf[j]]])) {
      if(!is.logical(functions[[nf[j]]])) {
        stop(paste("argument", nf[j], "is not a function!"))
      } else {
        if(functions[[nf[j]]]) {
          functions[[nf[j]]] <- get(default_fun[j])
        }
      }
    }
  }

  ## Create the 'bamlss.frame'.
  bf <- match.call(expand.dots = TRUE)
  bf[c("transform", "optimizer", "sampler", "samplestats",
    "results", "cores", "sleep", "combine", "model", "x")] <- NULL
  bf[[1]] <- as.name("bamlss.frame")
  bf <- eval(bf, envir = parent.frame())

  ## Transform.
  if(is.function(functions$transform)) {
    tbf <- functions$transform(bf, ...)
    bf[names(tbf)] <- tbf
    rm(tbf)
  }

  ## Start optimizer.
  if(is.function(functions$optimizer)) {
    opt <- functions$optimizer(x = bf$x, y = bf$y, family = bf$family,
      start = start, weights = model.weights(bf$model.frame),
      offset = model.offset(bf$model.frame), ...)
    if(!is.list(opt)) {
      if(inherits(opt, "numeric")) {
        opt <- list("parameters" = drop(opt))
      } else stop("the optimizer should return the parameters as named numeric vector!")
    }
    if(is.null(opt$parameters))
      stop("the optimizer must return an element $parameters!")
    if(inherits(opt$parameters, "data.frame") | inherits(opt$parameters, "matrix")) {
      if(is.null(colnames(opt$parameters)))
        stop("the returned parameters must be a named numeric data.frame or matrix!")
    } else {
      if(is.null(names(opt$parameters)))
        stop("the returned parameters must be a named numeric vector!")
    }
    bf$parameters <- opt$parameters
    if(!is.null(opt$fitted.values))
      bf$fitted.values <- opt$fitted.values
    if(!is.null(opt$hessian))
      bf$hessian <- opt$hessian
    ne <- names(opt)
    bf$model.stats <- opt[ne[!(ne %in% c("parameters", "fitted.values", "hessian"))]]
    rm(opt)
  }

  ## Start sampling.
  if(is.function(functions$sampler)) {
    if(is.null(cores)) {
      bf$samples <- functions$sampler(x = bf$x, y = bf$y, family = bf$family,
        weights = model.weights(bf$model.frame),
        offset = model.offset(bf$model.frame),
        start = if(is.null(bf$parameters)) start else unlist(bf$parameters),
        hessian = bf$hessian, ...)
    } else {
      parallel_fun <- function(j) {
        if(j > 1 & !is.null(sleep)) Sys.sleep(sleep)
        functions$sampler(x = bf$x, y = bf$y, family = bf$family,
          weights = model.weights(bf$model.frame),
          offset = model.offset(bf$model.frame),
          start = if(is.null(bf$parameters)) start else unlist(bf$parameters),
          hessian = bf$hessian, ...)
      }
      bf$samples <- parallel::mclapply(1:cores, parallel_fun, mc.cores = cores)
    }
    if(!inherits(bf$samples, "mcmc")) {
      if(!is.null(bf$samples)) {
        if(is.list(bf$samples)) {
          bf$samples <- as.mcmc.list(lapply(bf$samples, as.mcmc))
        } else {
          bf$samples <- as.mcmc(bf$samples)
        }
      }
    }

    ## Process samples.
    bf$samples <- process.chains(bf$samples, combine)

    ## Optionally, compute more model stats from samples.
    if(is.function(functions$samplestats)) {
      ms <- functions$samplestats(samples = bf$samples, x = bf$x, y = bf$y, family = bf$family, ...)
      if(is.null(bf$model.stats)) {
        bf$model.stats <- list("sampler" = list())
        bf$model.stats$sampler[names(ms)] <- ms
      } else {
        bf$model.stats <- list("optimizer" = bf$model.stats, "sampler" = list())
        bf$model.stats$sampler[names(ms)] <- ms
      }
    } else {
      if(!is.null(bf$model.stats))
        bf$model.stats <- list("optimizer" = bf$model.stats)
    }
  } else {
    if(!is.null(bf$model.stats))
      bf$model.stats <- list("optimizer" = bf$model.stats)
  }

  ## Compute results.
  if(is.function(functions$results))
    bf$results <- try(functions$results(bf, bamlss = TRUE,  ...))

  ## Save the model frame?
  if(!model | light)
    bf$model.frame <- NULL

  ## Save 'x' master object?
  if(!x)
    bf$x <- NULL
  if(light) {
    if(!is.null(bf$x)) {
      for(j in names(bf$x)) {
        if(length(bf$x[[j]]$smooth.construct)) {
          for(i in seq_along(bf$x[[j]]$smooth.construct)) {
            bf$x[[j]]$smooth.construct[[i]][c("X", "S", "Xr", "Xf", "binning", "prior", "grad", "hess")] <- NULL
            bf$x[[j]]$smooth.construct[[i]][["xt"]][["binning"]] <- NULL
            if(!is.null(bf$x[[j]]$smooth.construct[[i]][["margin"]])) {
              for(jj in seq_along(bf$x[[j]]$smooth.construct[[i]][["margin"]])) {
                bf$x[[j]]$smooth.construct[[i]][["margin"]][[jj]][c("X", "S", "Xr", "Xf", "binning", "prior", "grad", "hess")] <- NULL
                bf$x[[j]]$smooth.construct[[i]][["margin"]][[jj]][["xt"]][["binning"]] <- NULL
              }
            }
          }
        }
      }
    }
    bf$y <- NULL
    bf$fitted.values <- NULL
    attr(bf$formula, ".Environment") <- NULL
    attr(bf$terms, ".Environment") <- NULL
  }

  bf$call <- match.call()
  class(bf) <- c("bamlss", "bamlss.frame", "list")
  attr(bf, "functions") <- functions

  bf
}


## Basic engine setup transformer.
bamlss.setup <- function(x, ...)
{
  list("x" = bamlss.engine.setup(x$x, ...))
}

## family extractor.
family.bamlss <- family.bamlss.frame <- function(object, ...)
{
  return(object$family)
}


## Extract all parameter names.
get.all.parnames <- function(x, rename.p = TRUE)
{
  pn <- names(parameters(if(inherits(x, "bamlss.frame")) x$x else x, list = FALSE))
  if(rename.p)
    pn <- gsub("s.model.matrix", "p.model.matrix", pn, fixed = TRUE)
  pn
}


## Model stats based on samples.
samplestats <- function(samples, x = NULL, y = NULL, family = NULL, logLik = FALSE, ...)
{
  if(inherits(samples, "bamlss")) {
    if(is.null(samples$samples))
      stop("no samples in 'bamlss' object!")
    x <- if(is.null(samples$x)) smooth.construct(samples) else samples$x
    y <- samples$y
    family <- samples$family
    samples <- samples$samples
  }
  what <- c("logLik", "logPost", "DIC", "pd")
  if(inherits(samples, "mcmc.list"))
    samples <- process.chains(samples)
  if(is.null(samples)) return(NULL)
  samples <- as.matrix(samples)
  sn <- colnames(samples)
  stats <- NULL
  taken <- what[what %in% sn]
  if(length(taken)) {
    what <- what[!(what %in% taken)]
    stats <- samples[, taken, drop = FALSE]
    if(logLik) {
      if("loglik" %in% tolower(taken))
        return(stats[, tolower(taken) == "loglik"])
    }
    stats <- as.list(apply(stats, 2, mean, na.rm = TRUE))
  }
  if(is.data.frame(y)) {
    if(ncol(y) < 2)
      y <- y[[1]]
  }
  if(length(what) | logLik) {
    pn <- get.all.parnames(x, rename.p = TRUE)
    pn <- gsub(".p.model.matrix.", ".p.", pn, fixed = TRUE)
    pn <- pn[pn %in% colnames(samples)]
    if(length(pn) & ("DIC" %in% what)) {
      samples <- samples[, pn, drop = FALSE]
      if(is.null(family$p2d) & is.null(family$p2logLik)) {
        nx <- names(x)
        par <- rep(list(0), length = length(x))
        names(par) <- nx
        mpar <- par
        for(i in nx) {
          par[[i]] <- .fitted.bamlss(i, x[[i]], samples)
          par[[i]] <- make.link2(family$links[i])$linkinv(par[[i]])
        }
        msamples <- matrix(apply(samples, 2, mean, na.rm = TRUE), nrow = 1)
        colnames(msamples) <- pn
        for(i in nx)
          mpar[[i]] <- make.link2(family$links[i])$linkinv(.fitted.bamlss(i, x[[i]], msamples))
        tpar <- mpar
        dev <- ll <- rep(NA, ncol(par[[1]]))
        for(j in 1:ncol(par[[1]])) {
          for(i in nx)
            tpar[[i]] <- par[[i]][, j]
          llt <- try(family$loglik(y, tpar), silent = TRUE)
          if(!inherits(llt, "try-error")) {
            ll[j] <- llt
            dev[j] <- -2 * ll[j]
          }
        }
        if(logLik)
          return(ll)
        ll <- try(family$loglik(y, mpar), silent = TRUE)
      } else {
        mpar <- apply(samples, 2, mean, na.rm = TRUE)
        ll <- try(apply(samples, 1, function(x) {
          names(x) <- colnames(samples)
          if(is.null(family$p2logLik)) {
            sum(family$p2d(x, log = TRUE), na.rm = TRUE)
          } else family$p2logLik(x)
        }), silent = TRUE)
        if(inherits(ll, "try-error")) {
          warning("no DIC, cannot evaluate the $p2d() function in 'bamlss' family object!")
        } else {
          if(logLik)
            return(ll)
          dev <- -2 * ll
          ll <- try(if(is.null(family$p2logLik)) {
              sum(family$p2d(mpar, log = TRUE), na.rm = TRUE)
            } else family$p2logLik(mpar), silent = TRUE)
        }
      }
      if(!inherits(ll, "try-error") & !all(!is.finite(dev))) {
        mdev <- -2 * ll
        pd <- mean(dev, na.rm = TRUE) - mdev
        DIC <- mdev + 2 * pd
        if(is.null(stats))
          stats <- list()
        stats$DIC <- DIC
        stats$pd <- pd
      } else {
        warning("no DIC, cannot evaluate the $loglik() function in 'bamlss' family object!")
      }
    }
  }
  
  return(stats)
}


#### -----------------------------------------------------------------------------------------------
#### -----------------------------------------------------------------------------------------------
#### -----------------------------------------------------------------------------------------------
#### -----------------------------------------------------------------------------------------------
#### -----------------------------------------------------------------------------------------------
#### -----------------------------------------------------------------------------------------------
## Could be interesting: http://people.duke.edu/~neelo003/r/
##                       http://www.life.illinois.edu/dietze/Lectures2012/
#########################
## (2) Engine stacker. ##
#########################
#stacker <- function(x, optimizer = bfit0, sampler = samplerJAGS, ...)
#{
#  if(is.function(optimizer) | is.character(optimizer))
#    optimizer <- list(optimizer)
#  if(is.integer(sampler) | is.numeric(sampler)) {
#    n.samples <- as.integer(sampler)
#    sampler <- function(x, ...) { null.sampler(x, n.samples = n.samples) }
#  }
#  if(is.null(sampler))
#    sampler <- null.sampler
#  if(is.function(sampler) | is.character(sampler))
#    sampler <- list(sampler)
#  if(length(optimizer)) {
#    for(j in optimizer) {
#      if(is.character(j)) j <- eval(parse(text = j))
#      if(!is.function(j)) stop("the optimizer must be a function!")
#      x <- j(x, ...)
#    }
#  }
#  if(length(sampler)) {
#    for(j in sampler) {
#      if(is.character(j)) j <- eval(parse(text = j))
#      if(!is.function(j)) stop("the sampler must be a function!")
#      x <- j(x, ...)
#    }
#  }

#  x
#}


"[.bamlss" <- function(x, ...) {
  rval <- NextMethod("[", ...)
  mostattributes(rval) <- attributes(x)
  rval
}


## Create the model.frame.
bamlss.model.frame <- function(formula, data, family = gaussian_bamlss(),
  weights = NULL, subset = NULL, offset = NULL, na.action = na.omit,
  specials = NULL, contrasts.arg = NULL, drop.unused.levels = TRUE, ...)
{
  if(!missing(data)) {
    if(is.character(data)) {
      if(!file.exists(data))
        stop("data path is not existing!")
      ## data <- read.table.ffdf(file = data,
      ##   na.strings = "", header = TRUE, sep = ",")
      ## return(data)
      ## FIXME: ff data.frames!
      data <- read.table(file = data, header = TRUE, ...)
    }
  } else data <- NULL
  if(inherits(formula, "bamlss.frame") | inherits(formula, "bamlss")) {
    if(!is.null(formula$model.frame))
      return(formula$model.frame)
    fcall <- formula$call
    fcall[[1L]] <- quote(bamlss.model.frame)
    env <- environment(formula$formula)
    if(is.null(env))
      env <- parent.frame()
    return(eval(fcall, env))
  } else {
    family <- bamlss.family(family, ...)
    formula <- bamlss.formula(formula, family, specials, env = parent.frame())
  }

  if(is.null(na.action))
    na.action <- get(getOption("na.action"))
  if(missing(data))
    data <- environment(formula)

  if(!is.data.frame(data))
    data <- as.data.frame(data)

  ## Make fake "Formula" object.
  fF <- make_fFormula(formula)
  attr(fF, ".Environment") <- environment(formula)

  ## Resulting terms object.
  mterms <- terms(formula(fF), data = data)

  ## Set up the model.frame.
  data <- list(formula = fF, data = data, subset = subset,
    na.action = na.action, drop.unused.levels = drop.unused.levels, ...)
  data <- do.call("model.frame", data)
  rownames(data) <- NULL

  ## Code from stats model.matrix()
  contr.funs <- as.character(getOption("contrasts"))
  namD <- names(data)
  for(i in namD) {
    if(is.character(data[[i]])) 
      data[[i]] <- factor(data[[i]])
  }
  isF <- vapply(data, function(x) is.factor(x) || is.logical(x), NA)
  isF[attr(mterms, "response")] <- FALSE
  isOF <- vapply(data, is.ordered, NA)
  for(nn in namD[isF]) {
    if(is.null(attr(data[[nn]], "contrasts"))) {
      contrasts(data[[nn]]) <- contr.funs[1 + isOF[nn]]
    }
  }
  if(!is.null(contrasts.arg) && is.list(contrasts.arg)) {
    if(is.null(namC <- names(contrasts.arg))) 
      stop("invalid 'contrasts' argument")
    for(nn in namC) {
      if(is.na(ni <- match(nn, namD))) 
        warning(gettextf("variable '%s' is absent, its contrast will be ignored", nn), domain = NA)
      else {
        ca <- contrasts.arg[[nn]]
        if(is.matrix(ca)) {
          contrasts(data[[ni]], ncol(ca)) <- ca
        } else {
          contrasts(data[[ni]]) <- contrasts.arg[[nn]]
        }
      }
    }
  }

  ## Process weights and offset.
  if(!is.null(weights)) {
    if(!is.list(weights)) {
      weights <- rep(list(weights), length = length(family$names))
      names(weights) <- family$names
    }
    weights <- do.call("cbind", weights)
    colnames(weights) <- names(formula)[1:ncol(weights)]
    if(!is.null(subset)) {
      weights <- if(!is.logical(subset)) {
        weights[subset, , drop = FALSE]
      } else subset(weights, subset)
    }
    if(nrow(weights) < 2)
      weights <- do.call(rbind, replicate(nrow(data), weights, simplify = FALSE))
    for(j in 1:ncol(weights))
      weights[weights[, j] == 0, j] <- .Machine$double.eps
    data[["(weights)"]] <- weights
  }
  if(!is.null(offset)) {
    if(!is.list(offset)) {
      offset <- rep(list(offset), length = length(family$names))
      names(offset) <- family$names
    }
    offset <- do.call("cbind", offset)
    colnames(offset) <- names(formula)[1:ncol(offset)]
    if(!is.null(subset)) {
      offset <- if(!is.logical(subset)) {
        offset[subset, , drop = FALSE]
      } else subset(offset, subset)
    }
    if(nrow(offset) < 2)
      offset <- do.call(rbind, replicate(nrow(data), offset, simplify = FALSE))
    data[["(offset)"]] <- offset
  }

  ## Remove inf values.
  data <- rm_infinite(data)

  ## Assign terms object.
  attr(data, "terms") <- mterms

  ## Check response.
  if(!is.null(family$valid.response)) {
    family$valid.response(model.response(data))
  }

  data
}


## Remove Inf values from data.
rm_infinite <- function(x) {
  if(is.null(dim(x))) return(x)
  if(ncol(x) > 0) {
    for(j in 1:ncol(x)) {
      if(any(class(x[, j]) %in% c("numeric", "integer"))) {
        if(any(!is.finite(x[, j]))) {
          warning("infinite values in data, removing these observations in model.frame!")
          x <- x[is.finite(x[, j]), ]
        }
      }
    }
  }
  x
}


## Parse families and get correct family object, depending on type.
bamlss.family <- function(family, type = "bamlss", ...)
{
  family <- if(is.function(family)) family() else {
    if(is.character(family)) {
      if(!is.null(type)) {
        if(!grepl("gF(", family, fixed = TRUE) & !grepl("gF2(", family, fixed = TRUE))
          if(!grepl(type, family))
            family <- paste(family, type, sep = "_")
      }
      family <- eval(parse(text = family[1]))
      if(is.function(family))
        family()
      else family
    } else family
  }
  if(inherits(family, "gamlss.family"))
    family <- tF(family, ...)
  if(!inherits(family, "family.bamlss")) {
    if(!is.character(family)) {
      if(is.null(family$family)) stop("family is specified wrong, no family name available!")
      family <- family$family
    }
    txt <- paste(family, type, sep = if(!is.null(type)) "_" else "")
    txt <- gsub("bamlss.bamlss", "bamlss", txt, fixed = TRUE)
    family <- eval(parse(text = txt[1]))
    family <- family()
  }
  if(is.null(family)) family <- list()

  family
}


complete.bamlss.family <- function(family)
{
  if(is.null(names(family$links)))
    names(family$links) <- family$names

  linkinv <- linkfun <- list()
  for(j in family$names) {
    link <- make.link2(family$links[j])
    linkinv[[j]] <- link$linkinv
    linkfun[[j]] <- link$linkfun
  }

  if(is.null(family$map2par)) {
    family$map2par <- function(eta) {
      if(inherits(eta[[1L]], "ff")) {
        for(j in family$names)
          eta[[j]] <- ff_eval(eta[[j]], FUN = function(x) { linkinv[[j]](x) },
            lower = c(-Inf, -10), upper = c(Inf, 10))
      } else {
        for(j in family$names) {        
          eta[[j]] <- linkinv[[j]](eta[[j]])
          eta[[j]][is.na(eta[[j]])] <- 0
          if(any(jj <- eta[[j]] == Inf))
            eta[[j]][jj] <- 10
          if(any(jj <- eta[[j]] == -Inf))
            eta[[j]][jj] <- -10
        }
      }
      return(eta)
    }
  }

  if(is.null(family$mu)) {
    family$mu <- function(par) { make.link2(family$links[1])$linkinv(par[[1]]) }
  }

  if(is.null(family$loglik)) {
    if(!is.null(family$d))
      family$loglik <- function(y, par, ...) { sum(family$d(y, par, log = TRUE), na.rm = TRUE) }
  }

  err01 <- .Machine$double.eps^(1/3)
  err02 <- err01 * 2
  err11 <- .Machine$double.eps^(1/4)
  err12 <- err11 * 2

  if(is.null(family$score))
    family$score <- list()
  for(i in family$names) {
    if(is.null(family$score[[i]])) {
      fun <- c(
        "function(y, par, ...) {",
        paste("  eta <- linkfun[['", i, "']](par[['", i, "']]);", sep = ""),
        paste("  par[['", i, "']] <- linkinv[['", i, "']](eta + err01);", sep = ""),
        "  d1 <- family$d(y, par, log = TRUE);",
        paste("  par[['", i, "']] <- linkinv[['", i, "']](eta - err01);", sep = ""),
        "  d2 <- family$d(y, par, log = TRUE);",
        "  return((d1 - d2) / err02)",
        "}"
      )
      family$score[[i]] <- eval(parse(text = paste(fun, collapse = "")))
      attr(family$score[[i]], "dnum") <- TRUE
    }
  }

  if(is.null(family$hess))
    family$hess <- list()
  for(i in family$names) {
    if(is.null(family$hess[[i]])) {
      fun <- if(!is.null(attr(family$score[[i]], "dnum"))) {
        c(
          "function(y, par, ...) {",
          paste("  eta <- linkfun[['", i, "']](par[['", i, "']]);", sep = ""),
          paste("  par[['", i, "']] <- linkinv[['", i, "']](eta + err11);", sep = ""),
          paste("  d1 <- family$score[['", i, "']](y, par, ...);", sep = ""),
          paste("  par[['", i, "']] <- linkinv[['", i, "']](eta - err11);", sep = ""),
          paste("  d2 <- family$score[['", i, "']](y, par, ...);", sep = ""),
          "  return(-1 * (d1 - d2) / err12)",
          "}"
        )
      } else {
        c(
          "function(y, par, ...) {",
          paste("  eta <- linkfun[['", i, "']](par[['", i, "']]);", sep = ""),
          paste("  par[['", i, "']] <- linkinv[['", i, "']](eta + err01);", sep = ""),
          paste("  d1 <- family$score[['", i, "']](y, par, ...);", sep = ""),
          paste("  par[['", i, "']] <- linkinv[['", i, "']](eta - err01);", sep = ""),
          paste("  d2 <- family$score[['", i, "']](y, par, ...);", sep = ""),
          "  return(-1 * (d1 - d2) / err02)",
          "}"
        )
      }
      family$hess[[i]] <- eval(parse(text = paste(fun, collapse = "")))
    }
  }

  return(family)
}


## Formula to list().
as.list.Formula <- function(x)
{
  if(!inherits(x, "Formula"))
    x <- as.Formula(x)
  env <- environment(x)
  lhs <- attr(x, "lhs")
  rhs <- attr(x, "rhs")
  nl <- length(lhs)
  nr <- length(rhs)
  if(nl < nr)
    lhs <- c(lhs, rep(list(NA), length = nr - nl))
  if(nr < nl)
    rhs <- c(rhs, rep(list(1), length = nl - nr))
  x <- mapply(c, lhs, rhs, SIMPLIFY = FALSE)
  formula <- list()
  for(i in seq_along(x)) {
    check <- inherits(x[[i]][[1]], "call") | inherits(x[[i]][[1]], "name")
    f <- if(check) {
      as.call(c(as.symbol("~"), x[[i]]))
    } else as.call(c(as.symbol("~"), x[[i]][[2]]))
    formula[[i]] <- eval(f, envir = env)
    attr(formula[[i]], ".Environment") <- NULL
  }
  environment(formula) <- env
  formula
}


## Special formula parser, can deal with multi parameter models
## and hierarchical structures.
bamlss.formula <- function(formula, family = NULL, specials = NULL, env = NULL, ...)
{
  if(is.null(specials))
    specials <- c("s", "te", "t2", "sx", "s2", "rs", "ti", "tx", "tx2", "tx3", "la", "n")
  if(inherits(formula, "bamlss.formula"))
    return(formula)
  if(!is.list(formula)) {
    if(!inherits(formula, "Formula"))
      formula <- as.Formula(formula)
    if(inherits(formula, "Formula"))
      formula <- as.list.Formula(formula)
  }
  if(!is.null(family))
    family <- bamlss.family(family, ...)
  if(!is.list(formula)) formula <- list(formula)
  if(!length(formula)) stop("formula is specified wrong!")
  
  if(is.null(env))
    env <- get_formula_envir(formula)

  complete_formula <- function(formula) {
    if(!is.null(family)) {
      if(length(formula) < length(family$names))
        formula <- c(formula, rep(list(), length = length(family$names) - length(formula)))
    }
    fn <- NULL
    for(j in seq_along(formula)) {
      ft <- if(!inherits(formula[[j]], "formula")) formula[[j]][[1]] else formula[[j]]
      if(!is.null(ft)) {
        yok <- attr(terms(formula(as.Formula(ft), rhs = FALSE)), "response") > 0
        fn <- c(fn, if(yok) all.vars(ft)[1] else NULL)
      }
    }
    fn[fn %in% c("1", "-1")] <- NA

    nas <- which(is.na(fn))
    if(!is.null(family)) {
      if(length(nas))
        fn[nas] <- family$names[nas]
      if(is.null(family$names))
        family$names <- NA
      if(!all(is.na(family$names[1:length(fn)])))
        fn <- family$names[1:length(fn)]
      else
        family$names[1:length(fn)] <- fn
    } else fn[nas] <- paste("par", 1:length(fn[nas]), sep = ".")

    if(is.null(family)) {
      if(length(fn) < length(formula)) {
        k <- length(formula) - length(fn)
        if(k > 1)
          fn <- c(fn, paste("?par", 1:k, sep = ""))
        else
          fn <- c(fn, "?par")
      }
    }
    names(formula) <- fn
    if(!is.null(family)) {
      if(any(i <- is.na(names(formula))))
        names(formula)[i] <- family$names[i]
      for(j in family$names) {
        if(is.null(formula[[j]])) {
          formula[[j]] <- as.formula(paste(j, "1", sep = " ~ "), env = NULL)
        }
      }
    }
    for(j in seq_along(formula)) {
      if(!inherits(formula[[j]], "formula")) {
        if(is.null(names(formula[[j]])))
          names(formula[[j]]) <- paste("h", 1:length(formula[[j]]), sep = "")
      } else {
        attr(formula[[j]], ".Environment") <- NULL
      }
    }
    formula
  }
  formula <- formula_and(formula)
  formula <- formula_at(formula)
  formula <- complete_formula(formula_hierarchical(formula))
  formula <- formula_extend(formula, family, specials)

  environment(formula) <- env
  class(formula) <- c("bamlss.formula", "list")

  formula
}


## Formula environment.
get_formula_envir <- function(formula)
{
  env <- environment(formula)
  if(is.null(env)) {
    get_env <- function(x) {
      if(inherits(x, "list")) {
        env <- NULL
        for(j in x)
          env <- c(env, get_env(j))
        return(env)
      } else return(environment(x))
    }
    env <- get_env(formula)
  }
  if(is.null(env)) env <- .GlobalEnv
  if(is.list(env))
    env <- env[[1]]
  return(env)
}


## Process categorical responses.
bamlss.formula.cat <- function(formula, family, data, reference)
{
  env <- environment(formula)

  rn <- y <- NULL
  for(j in seq_along(formula)) {
    ft <- if(!inherits(formula[[j]]$formula, "formula")) {
      formula[[j]][[1]]$formula
    } else formula[[j]]$formula
    yok <- attr(terms(formula(as.Formula(ft), rhs = FALSE)), "response") > 0
    if(yok)
      rn <- c(rn, all.vars(ft)[1])
  }
  rn2 <- rn[rn %in% names(data)]
  cat <- !is.null(family$cat) & (length(rn2) > 1)
  if(is.factor(data[[rn2[1]]]) | cat) {
    if((nlevels(data[[rn2[1]]]) > 2) | cat) {
      if(!cat | is.factor(data[[rn2[1]]])) {
        ft <- as.formula(paste("~ -1 +", rn2[1]), env = NULL)
        y <- model.matrix(ft, data = data)
        colnames(y) <- rmf(gsub(rn2[1], "", colnames(y), fixed = TRUE))
        if(is.null(reference)) {
          ty <- table(data[[rn2[1]]])
          reference <- c(names(ty)[ty == max(ty)])[1]
        } else {
          ld <- levels(data[[rn2[1]]])
          reference <- ld[match(gsub(rn2[1], "", reference), ld)]
        }
        if(is.na(reference))
          stop(paste("cannot find reference category within response levels!"))
        reference <- rmf(reference)
        ylevels <- rmf(levels(data[[rn2[1]]]))
        ylevels <- ylevels[ylevels != reference]
        y <- y[, colnames(y) %in% ylevels, drop = FALSE]
      } else {
        y <- data[, rn2]
        ylevels <- colnames(y)
        reference <- ""
      }
      if(length(formula) < ncol(y)) {
        formula <- c(formula, rep(formula, length = ncol(y) - length(formula)))
      }
      if(!(names(formula)[[1]] %in% colnames(y))) {
        names(formula)[[1]] <- colnames(y)[1]
        ft <- if(!inherits(formula[[1]]$formula, "formula")) {
          formula[[1]][[1]]$formula
        } else formula[[1]]$formula
        env <- environment(ft)
        ft <- update(ft, as.formula(paste(colnames(y)[1], ".", sep = "~"), env = NULL))
        environment(ft) <- env
        if(!inherits(formula[[1]]$formula, "formula")) {
          formula[[1]][[1]]$formula <- ft
          formula[[1]][[1]]$response <- colnames(y)[1]
        } else {
          formula[[1]]$formula <- ft
          formula[[1]]$response <- colnames(y)[1]
        }
        ft <- if(!inherits(formula[[1]]$formula, "formula")) {
          formula[[1]][[1]]$fake.formula
        } else formula[[1]]$fake.formula
        ft <- update(ft, as.formula(paste(colnames(y)[1], ".", sep = "~"), env = NULL))
        if(!inherits(formula[[1]]$formula, "formula")) {
          formula[[1]][[1]]$fake.formula <- ft
        } else {
          formula[[1]]$fake.formula <- ft
        }
      }
      if(length(i <- !(names(formula) %in% ylevels))) {
        k <- 1
        ynot <- ylevels[!(ylevels %in% names(formula))]
        for(j in which(i)) {
          names(formula)[[j]] <- ynot[k]
          ft <- if(!all(c("formula", "fake.formula") %in% names(formula[[j]]))) {
            formula[[j]][[1]]$formula
          } else formula[[j]]$formula
          env <- environment(ft)
          ft <- update(ft, as.formula(paste(ynot[k], "~ ."), env = NULL))
          environment(ft) <- env
          if(!inherits(formula[[j]]$formula, "formula")) {
            formula[[j]][[1]]$formula <- ft
            formula[[j]][[1]]$response <- ynot[k]
          } else {
            formula[[j]]$formula <- ft
            formula[[j]]$response <- ynot[k]
          }
          attr(formula[[j]], "name") <- ynot[k]
          k <- k + 1
        }
      }
    }
  }

  rval <- if(!is.null(y)) {
    class(formula) <- "bamlss.formula"
    environment(formula) <- env
    list("formula" = formula, "ylevels" = ylevels, "reference" = reference)
  } else NULL
  rval
}


## Make "Formula" object from fake.formulas.
make_fFormula <- function(formula)
{
  fF <- NULL
  for(j in seq_along(formula)) {
    if(!all(c("formula", "fake.formula") %in% names(formula[[j]]))) {
      for(i in seq_along(formula[[j]]))
        fF <- c(fF, formula[[j]][[i]]$fake.formula)
    } else {
      fF <- c(fF, formula[[j]]$fake.formula)
    }
  }
  fF <- do.call("as.Formula", fF)
  fF
}


all.vars.formula <- function(formula, lhs = TRUE, rhs = TRUE, specials = NULL, intercept = FALSE, type = 1)
{
  env <- environment(formula)
  specials <- unique(c(specials, "s", "te", "t2", "sx", "s2", "rs", "ti", "tx", "tx2", "tx3", "la", "n"))
  tf <- terms(formula, specials = specials, keep.order = TRUE)
  ## sid <- unlist(attr(tf, "specials")) - attr(tf, "response")
  tl <- attr(tf, "term.labels")
  sid <- NULL
  for(j in specials)
    sid <- c(sid, grep2(paste(j, "(", sep = ""), tl, fixed = TRUE))
  if(length(sid))
    sid <- sort(unique(sid))
  vars <- NULL
  if(rhs) {
    if(length(sid)) {
      vars <- tl[-sid]
      if(!length(vars))
        vars <- NULL
      for(j in tl[sid]) {
        tcall <- parse(text = j)[[1]]
        tcall[c("k","fx","bs","m","xt","id","sp","pc","d","mp","np")] <- NULL
        tcall <- eval(tcall)
        vars <- c(vars, tcall$term)
        if(!is.null(tcall$by)) {
          if(tcall$by != "NA")
            vars <- c(vars, tcall$by)
        }
      }
    } else {
      vars <- tl
    }
  }
  if(lhs & (attr(tf, "response") > 0))
    vars <- c(vars, response.name(formula, keep.functions = TRUE))
  if(intercept & (attr(tf, "intercept") > 0))
    vars <- c("1", vars)
  vars <- vars[vars != "."]
  if(length(vars) < 1)
    vars <- NULL
  if(any(i <- grep(":", vars, fixed = TRUE))) {
    dv <- unlist(strsplit(vars[i], ":", fixed = TRUE))
    vars <- c(vars[-i], dv)
  }
  vars <- unique(vars)
  if(type == 1)
    vars <- all.vars(as.formula(paste("~", paste(vars, collapse = "+")), env = NULL))
  unique(vars)
}


## From nlme.
splitFormula <- function(form, sep = "+") 
{
  if(inherits(form, "formula") || mode(form) == "call" && 
    form[[1]] == as.name("~")) 
    return(splitFormula(form[[length(form)]], sep = sep))
  if(mode(form) == "call" && form[[1]] == as.name(sep)) 
    return(do.call("c", lapply(as.list(form[-1]), splitFormula, 
      sep = sep)))
  if(mode(form) == "(") 
    return(splitFormula(form[[2]], sep = sep))
  if(length(form) < 1) 
    return(NULL)
  list(stats::asOneSidedFormula(form))
}


terms.formula2 <- function(formula, specials, keep.order = TRUE, ...)
{
  fs <- splitFormula(formula, sep = c("+", "-"))
  tl <- rep("", length = length(fs))
  adds <- NULL
  for(j in seq_along(fs)) {
    tlj <- attr(terms(fs[[j]], specials = specials), "term.labels")
    if(length(tlj))
      tl[j] <- tlj[1]
    if(length(tlj) > 1)
      adds <- c(adds, tlj[-1])
  }
  tl <- c(tl, adds)
  tl <- tl[tl != ""]
  if(!length(tl))
    tl <- ""
  attr(formula, "term.labels") <- tl
  formula
}


all.labels.formula <- function(formula, specials = NULL, full.names = FALSE)
{
  env <- environment(formula)
  specials <- unique(c("s", "te", "t2", "sx", "s2", "rs", "ti", "tx", "tx2", "tx3", "la", "n", specials))
  tf <- terms.formula2(formula, specials = specials, keep.order = FALSE)
  ## sid <- unlist(attr(tf, "specials")) - attr(tf, "response")
  tl <- attr(tf, "term.labels")
  sid <- NULL
  for(j in specials)
    sid <- c(sid, grep2(paste(j, "(", sep = ""), tl, fixed = TRUE))
  if(length(sid))
    sid <- sort(unique(sid))
  labs <- NULL

  if(length(sid)) {
    labs <- tl[-sid]
    if(full.names & length(labs))
      labs <- paste("p", labs, sep = ".")
    if(!length(labs))
      labs <- NULL
    else
      tl[-sid] <- labs
    tl[sid] <- gsub(" ", "", tl[sid])
    for(j in sid) {
      tcall <- parse(text = tl[j])[[1]]
      tcall[c("k","fx","bs","m","xt","sp","pc","d","mp","np")] <- NULL
      tcall <- eval(tcall)
      tl[j] <- gsub(" ", "", tcall$label)
      if(!is.null(tcall$by)) {
        if(tcall$by != "NA") {
          if(!grepl("by=", tl[j], fixed = TRUE)) {
            tlt <- strsplit(tl[j], "")[[1]]
            tlt <- paste(tlt[1:(length(tlt) - 1)], collapse = "")
            tl[j] <- paste(tlt, ",by=", tcall$by, ")", sep = "")
          }
        }
      }
    }
    if(full.names)
      tl[sid] <- paste("s", tl[sid], sep = ".")
    labs <- tl
  } else labs <- if(full.names) paste("p", tl, sep = ".") else tl
  unique(labs)
}


fake.formula <- function(formula, lhs = TRUE, rhs = TRUE, specials = NULL)
{
  if(all(!lhs & !rhs))
    return(0 ~ 0)
  if(all(rhs)) {
    f <- paste(all.vars.formula(formula, lhs = FALSE, rhs = TRUE, specials, intercept = TRUE, type = 2), collapse = "+")
    if(f == "") f <- "-1"
  }
  if(all(lhs))
    f <- paste(all.vars.formula(formula, lhs = TRUE, rhs = FALSE, type = 2), "~", if(!is.null(f)) f else 0)
  else
    f <- paste("~", if(!is.null(f)) f else 0)
  f <- as.formula(f, env = NULL)
  f
}


## Extend formula by a fake formula with all variables
## to compute a model.frame, create smooth objects.
formula_extend <- function(formula, family, specials = NULL)
{
  if(is.list(formula)) {
    for(j in seq_along(formula))
      formula[[j]] <- formula_extend(formula[[j]], family, specials)
    return(formula)
  } else {
    rn <- response.name(formula)
    ff <- fake.formula(formula, lhs = !(rn %in% family$names), specials = specials)
    return(list("formula" = formula, "fake.formula" = ff))
  }
}


## Get response name.
response.name <- function(formula, hierarchical = TRUE, keep.functions = FALSE, na.rm = FALSE)
{
  rn <- NA
  if(inherits(formula, "bamlss.frame")) {
    if(!is.null(formula$formula)) {
      if(!is.null(attr(formula$formula, "response.name")))
        return(attr(formula$formula, "response.name"))
    }
    formula <- terms(model.frame(formula))
  }
  if(!is.null(attr(formula, "terms")))
    formula <- attr(formula, "terms")
  if(inherits(formula, "formula")) {
    f <- as.Formula(formula)
    f <- formula(f, lhs = TRUE, rhs = FALSE)
    if(keep.functions) {
      cf <- as.character(formula)
      rn <- if(length(cf) < 3) character(0) else cf[2]
    } else {
      rn <- all.vars(f)
    }
  } else {
    if(inherits(formula, "list")) {
      rn <- NULL
      for(i in seq_along(formula)) {
        if(is.null(formula[[i]]$formula) & inherits(formula[[i]], "list")) {
          for(j in seq_along(formula[[i]])) {
            if(!hierarchical & (j > 1)) {
              next
            } else {
              tf <- if(is.null(formula[[i]][[j]]$formula)) {
                formula[[i]][[j]]
              } else formula[[i]][[j]]$formula
              rn <- c(rn, response.name(tf, keep.functions = keep.functions))
            }
          }
        } else {
          rn <- c(rn , response.name(formula[[i]]$formula, keep.functions = keep.functions))
        }
      }
    }
  }
  if(!length(rn))
    rn <- NA
  if(na.rm)
    rn <- rn[!is.na(rn)]
  rn
}

## Search and process "&"
formula_and <- function(formula)
{
  if(nol <- !is.list(formula))
    formula <- list(formula)
  for(j in seq_along(formula)) {
    if(!inherits(formula[[j]], "formula")) {
      formula[[j]] <- formula_and(formula[[j]])
    } else {
      ft <- deparse(formula[[j]])
      if(any(grep("&", ft, fixed = TRUE))) {
        formula[[j]] <- as.list(strsplit(ft, "&", fixed = TRUE)[[1]])
        for(i in seq_along(formula[[j]])) {
          if(!any(grepl("~", formula[[j]][[i]], fixed = TRUE)))
            formula[[j]][[i]] <- paste("~", formula[[j]][[i]])
          formula[[j]][[i]] <- as.formula(formula[[j]][[i]], env = NULL)
        }
      }
    }
  }
  if(nol) {
    names(formula) <- response.name(formula[[1]][[1]])
    if(inherits(formula[[1]], "formula"))
      formula <- formula[[1]]
  }
  formula
}

## Search and process "@"
formula_at <- function(formula)
{
  if(nol <- !is.list(formula))
    formula <- list(formula)
  for(j in seq_along(formula)) {
    if(!