Package 'bamlss'

Description

The bamlss package is a general tool for complex Bayesian regression modeling with structured additive predictors based on Markov chain Monte Carlo simulation. The design of this package substantially focuses on maximum flexibility and easy integration of new code and/or standalone systems. The package makes heavy use of mgcv infrastructures to build up all necessary model matrices and information from which it is relatively easy for the user to construct estimation algorithms or interfaces to existing software packages. The package can also be seen as an harmonized framework for regression modeling since it does not restrict to any type of problem. The main function in this package is bamlss, which is a wrapper function that calls optimizer and/or sampling functions for fitting Bayesian additive models for location scale and shape (and beyond). These model fitting functions can be exchanged by the user. Moreover, the package contains numerous functions for creating post-estimation results like summary statistics and effect plots etc.

Author(s)

Nikolaus Umlauf, Nadja Klein, Achim Zeileis.

References

Umlauf N, Klein N, Zeileis A (2019). BAMLSS: Bayesian Additive Models for Location, Scale and Shape (and Beyond). Journal of Computational and Graphical Statistics, 27(3), 612–627. doi:10.1080/10618600.2017.1407325

Umlauf N, Klein N, Simon T, Zeileis A (2021). bamlss: A Lego Toolbox for Flexible Bayesian Regression (and Beyond). Journal of Statistical Software, 100(4), 1–53. doi:10.18637/jss.v100.i04

See Also

bamlss, bamlss.frame

Description

This is the main model fitting function of the package. Function bamlss() is a wrapper function that parses the data and the model formula, or extended bamlss.formula, as well as the bamlss.family into a bamlss.frame. The bamlss.frame then holds all model matrices and information that is needed for setting up estimation engines. The model matrices are based on mgcv infrastructures, i.e., smooth terms are constructed using smooth.construct and smoothCon. Therefore, all mgcv model term constructors like s, te, t2 and ti can be used. Identifiability conditions are imposed using function gam.side.

After the bamlss.frame is set up function bamlss() applies optimizer and/or sampling functions. These functions can also be provided by the user. See the details below on how to create new engines to be used with function bamlss().

Finally, the estimated parameters and/or samples are used to create model output results like summary statistics or effect plots. The computation of results may also be controlled by the user.

Usage

bamlss(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, ...)

Arguments

Details

The main idea of this function is to provide infrastructures that make it relatively easy to create estimation engines for new problems, or write interfaces to existing software packages.

The steps that are performed within the function are:

• First, the function parses the data, the formula or the extended bamlss.formula as well as the bamlss.family into a model frame like object, the bamlss.frame. This object holds all necessary model matrices and information that is needed for subsequent model fitting engines. Per default, all package mgcv smooth term constructor functions like s, te, t2 and ti can be used (see also function smooth.construct), however, even special user defined constructors can be included, see the manual of bamlss.frame.

• In a second step, the bamlss.frame can be transformed, e.g., if a mixed model representation of smooth terms is needed, see function randomize.

• Then an optimizer function is started, e.g., a function that finds posterior mode estimates of the parameters. A convention for model fitting engines is that such functions should have the following arguments:

  optimizer(x, y, family, start, weights, offset, ...)

  Internally, function bamlss() will send the x object that holds all model matrices, the response y object, the family object, starting values for the parameters, possible weights and offsets of the created bamlss.frame to the optimizer function (see the manual of bamlss.frame for more details on the x, y and other objects). The job of the optimizer is to return a named numeric vector of optimum parameters. The names of the parameters should be such that they can be uniquely mapped to the corresponding model matrices in x. See function parameters for more details on parameter names. The default optimizer function is opt_bfit. The optimizer can return more information than only the optimum parameters. It is possible to return a list, the convention here is that an element named "parameters" then holds the named vector of estimated parameters. Possible other return values could be fitted values, the Hessian matrix, information criteria or information about convergence of the algorithm, etc. Note that the parameters are simply added to the bamlss.frame in an (list) entry named parameters.

• After the optimization step, a sampler function is started. The arguments of such sampler functions are the same as for the optimizer functions

  sampler(x, y, family, start, weights, offset, ...)

  Sampler functions must return a matrix of samples, each row represents one iteration and the matrix can be coerced to mcmc objects. The function may return a list of samples, e.g., if multiple chains are returned each list entry then holds one sample matrix of one chain. The column names of the sample matrix should be the same as the names of estimated parameters. For a possible naming convention see function parameters, which ensures unique mapping of samples with the model matrices in the x object of the bamlss.frame. The samples are added to the bamlss.frame in an (list) entry named samples.

• Next, the samplestats function is applied. This function can compute any quantity from the samples and the x object, the arguments of such functions are

  samplestats(samples, x, y, family, ...)

  where argument samples are the samples returned from the sampler function, and x, y and family are the same objects as passed to the optimizer and or sampler functions. For example, the default function in bamlss() for this task is also called samplestats and returns the mean of the log-likelihood and the log-posterior computed of all samples, as well as the DIC.

• The last step is to compute more complex information about the model using the results function. The arguments of such results functions are

  results(bamlss.frame, ...)

  here, the full bamlss.frame including possible parameters and samples is passed to the function within bamlss(). The default function for this task is results.bamlss.default which returns an object of class "bamlss.results" for which generic plotting functions are and a summary function is provided. Hence, the user can control the output of the model, the plotting and summary statistics, too.

Note that function transform(), optimizer(), sampler(), samplestats() and results() can be provided from the bamlss.family object, e.g., if a bamlss.family object has an element named "optimizer", which represents a valid optimizer function such as opt_bfit, exactly this optimizer function will be used as a default when using the family.

Value

An object of class "bamlss". The object is in principle only a slight extension of a bamlss.frame, i.e., if an optimizer is applied it will hold the estimated parameters in an additional element named "parameters". If a sampler function is applied it will additionally hold the samples in an element named "samples". The same mechanism is used for results function.

If the optimizer function computes additional output next to the parameters, this will be saved in an element named "model.stats". If a samplestats function is applied, the output will also be saved in the "model.stats" element.

Additionally, all functions that are called are saved as attribute "functions" in the returned object.

Author(s)

Nikolaus Umlauf, Nadja Klein, Achim Zeileis.

References

Umlauf N, Klein N, Zeileis A (2019). BAMLSS: Bayesian Additive Models for Location, Scale and Shape (and Beyond). Journal of Computational and Graphical Statistics, 27(3), 612–627. doi:10.1080/10618600.2017.1407325

Umlauf N, Klein N, Simon T, Zeileis A (2021). bamlss: A Lego Toolbox for Flexible Bayesian Regression (and Beyond). Journal of Statistical Software, 100(4), 1–53. doi:10.18637/jss.v100.i04

See Also

bamlss.frame, family.bamlss, bamlss.formula, randomize, bamlss.engine.setup, opt_bfit, sam_GMCMC, continue, coef.bamlss, parameters, predict.bamlss, plot.bamlss

Examples

## Not run: ## Simulated data example.
d <- GAMart()
f <- num ~ s(x1) + s(x2) + s(x3) + te(lon, lat)
b <- bamlss(f, data = d)
summary(b)
plot(b)
plot(b, which = 3:4)
plot(b, which = "samples")

## Use of optimizer and sampler functions:
## * first run optimizer,
b1 <- bamlss(f, data = d, optimizer = opt_bfit, sampler = FALSE)
print(b1)
summary(b1)

## * afterwards, start sampler with staring values,
b2 <- bamlss(f, data = d, start = coef(b1), optimizer = FALSE, sampler = sam_GMCMC)
print(b2)
summary(b2)

## Continue sampling.
b3 <- continue(b2, n.iter = 12000, burnin = 0, thin = 10)
plot(b3, which = "samples")
plot(b3, which = "max-acf")
plot(b3, which = "max-acf", burnin = 500, thin = 4)

## End(Not run)

Description

A single class and corresponding methods encompassing all bamlss.family distributions (from the bamlss package) using the workflow from the distributions3 package.

Usage

BAMLSS(family, ...)

Arguments

Details

The constructor function BAMLSS sets up a distribution object, representing a distribution from the BAMLSS (Bayesian additive model of location, scale, and shape) framework by the corresponding parameters plus a family attribute, e.g., gaussian_bamlss for the normal distribution or binomial_bamlss for the binomial distribution. The parameters employed by the family vary across the families but typically capture different distributional properties (like location, scale, shape, etc.).

All parameters can also be vectors, so that it is possible to define a vector of BAMLSS distributions from the same family with potentially different parameters. All parameters need to have the same length or must be scalars (i.e., of length 1) which are then recycled to the length of the other parameters.

For the BAMLSS distribution objects there is a wide range of standard methods available to the generics provided in the distributions3 package: pdf and log_pdf for the (log-)density (PDF), cdf for the probability from the cumulative distribution function (CDF), quantile for quantiles, random for simulating random variables, and support for the support interval (minimum and maximum). Internally, these methods rely on the usual d/p/q/r functions provided in bamlss, see the manual pages of the individual families. The methods is_discrete and is_continuous can be used to query whether the distributions are discrete on the entire support or continuous on the entire support, respectively.

See the examples below for an illustration of the workflow for the class and methods.

Value

A BAMLSS distribution object.

See Also

bamlss.family

Examples

## package and random seed
library("distributions3")
set.seed(6020)

## three Weibull distributions
X <- BAMLSS("weibull", lambda = c(1, 1, 2), alpha = c(1, 2, 2))
X

## moments (FIXME: mean and variance not provided by weibull_bamlss)
## mean(X)
## variance(X)

## support interval (minimum and maximum)
support(X)
is_discrete(X)
is_continuous(X)

## simulate random variables
random(X, 5)

## histograms of 1,000 simulated observations
x <- random(X, 1000)
hist(x[1, ], main = "Weibull(1,1)")
hist(x[2, ], main = "Weibull(1,2)")
hist(x[3, ], main = "Weibull(2,2)")

## probability density function (PDF) and log-density (or log-likelihood)
x <- c(2, 2, 1)
pdf(X, x)
pdf(X, x, log = TRUE)
log_pdf(X, x)

## cumulative distribution function (CDF)
cdf(X, x)

## quantiles
quantile(X, 0.5)

## cdf() and quantile() are inverses
cdf(X, quantile(X, 0.5))
quantile(X, cdf(X, 1))

## all methods above can either be applied elementwise or for
## all combinations of X and x, if length(X) = length(x),
## also the result can be assured to be a matrix via drop = FALSE
p <- c(0.05, 0.5, 0.95)
quantile(X, p, elementwise = FALSE)
quantile(X, p, elementwise = TRUE)
quantile(X, p, elementwise = TRUE, drop = FALSE)

## compare theoretical and empirical mean from 1,000 simulated observations
## (FIXME: mean not provided by weibull_bamlss)
## cbind(
##   "theoretical" = mean(X),
##   "empirical" = rowMeans(random(X, 1000))
## )

Description

These functions can be useful when setting up new model fitting engines that are based on the setup function bamlss.engine.setup. See the examples.

Usage

## Functions to extract parameter states.
get.par(x, what = NULL)
get.state(x, what = NULL)
set.par(x, replacement, what)

## Function for setting starting values.
set.starting.values(x, start)

Arguments

See Also

bamlss.engine.setup

Examples

## Create a bamlss.frame.
d <- GAMart()
bf <- bamlss.frame(num ~ s(x1) + s(x2) + te(lon,lat), data = d, family = "gaussian")
names(bf$x$mu$smooth.construct)

## Use the setup function for
## adding state elements.
bf$x <- bamlss.engine.setup(bf$x, df = c("s(x1)" = 1, "s(x2)" = 3))
names(bf$x$mu$smooth.construct)

## Extract regression coefficients.
get.state(bf$x$mu$smooth.construct[["te(lon,lat)"]], "b")

## Extract smoothing variances.
get.state(bf$x$mu$smooth.construct[["te(lon,lat)"]], "tau2")

## More examples.
state <- bf$x$mu$smooth.construct[["te(lon,lat)"]]$state
get.par(state$parameters, "b")
get.par(state$parameters, "tau2")

state$parameters <- set.par(state$parameters, c(0.1, 0.5), "tau2")
get.par(state$parameters, "tau2")

## Setting starting values.
start <- c("mu.s.s(x1).b" = 1:9, "mu.s.s(x1).tau2" = 0.1)
bf$x <- set.starting.values(bf$x, start = start)
get.state(bf$x$mu$smooth.construct[["s(x1)"]], "b")
get.state(bf$x$mu$smooth.construct[["s(x1)"]], "tau2")

Description

This function takes the x object of a bamlss.frame and adds additional objects that are useful for model fitting engines. This is applied only for 'regular' terms, e.g., as created by s and te. For special model terms the corresponding smooth.construct method is in charge of this (see also the examples for function bfit).

Usage

bamlss.engine.setup(x, update = "iwls", propose = "iwlsC_gp",
  do.optim = NULL, df = NULL, parametric2smooth = TRUE, ...)

Arguments

Details

For each model term in the "smooth.construct" object of the x list (as returned from bamlss.frame), this function adds a named list called "state" with the following entries:

• "parameters": A numeric vector. Regression coefficients are named with "b", smooth variances are named with "tau2".

• "fitted.values": Given the "parameters", the actual fitted values of the model term.

• "edf": Given the smoothing variances, the actual equivalent degrees of freedom (edf) of the model term.

• "do.optim": Should an optimizer function try to find optimum smoothing variances?

The state will be changed in each iteration and can be passed outside an updating function.

Additionally, if missing in the xt argument of a model term (see, e.g., function s for xt) the function adds the corresponding log-prior and its first and second order derivatives w.r.t. regression coefficients in functions grad() and hess().

Also, objects named "lower" and "upper" are added to each model term. These indicate the lower and upper boundaries of the parameter space.

Value

A transformed x list, as returned from function bamlss.frame.

See Also

bamlss.frame, bfit, GMCMC, get.par, set.par, get.state

Examples

d <- GAMart()
bf <- bamlss.frame(num ~ s(x1) + s(x2), data = d, family = "gaussian")
names(bf$x$mu$smooth.construct)
bf$x <- bamlss.engine.setup(bf$x, df = c("s(x1)" = 1, "s(x2)" = 3))
names(bf$x$mu$smooth.construct)
names(bf$x$mu$smooth.construct[["s(x1)"]])
names(bf$x$mu$smooth.construct[["s(x1)"]]$state)
sapply(bf$x$mu$smooth.construct, function(x) {
  c(x$state$edf, get.state(x, "tau2"))
})

Description

This function creates an extended BAMLSS formula. In combination with a bamlss.family object, each parameter of the response distribution is linked to a single formula. If no formula is supplied for a parameter, a simple intercept only model is created. Moreover, the function identifies hierarchical structures, see the examples. This function is useful for creating complex model.frames for (hierarchical) multi-parameter models and is used by function bamlss.frame.

Usage

bamlss.formula(formula, family = NULL, specials = NULL, env = NULL, ...)

Arguments

Value

A named list of class "bamlss.formula". Each list entry specifies a model, e.g., for one parameter of a provided bamlss.family object. Each entry (parameter model) then holds:

See Also

bamlss, bamlss.frame, bamlss.family

Examples

## Simple formula without family object.
f <- bamlss.formula(y ~ x1 + s(x2))
print(f)
print(str(f))

## Complex formula with list of formulae.
f <- list(
  y1 ~ x1 + s(x2),
  y2 ~ x3 + te(lon, lat),
  u ~ x4 + x1
)

f <- bamlss.formula(f)
print(f)
print(names(f))

## Same formula but using extended formulae
## of package Formula.
f <- y1|y2|u ~ x1 + s(x2)|x3 + te(lon,lat)|x4 + x1
f <- bamlss.formula(f)
print(f)
print(names(f))

## Using a bamlss family object, e.g., gaussian_bamlss().
## The family has two parameters, mu and sigma, for
## each parameter one formula is returned. If no
## formula is specified an intercept only model is
## generated for the respective parameter.
f <- bamlss.formula(y ~ x1 + s(x2), family = gaussian_bamlss)

## Note, same as:
f <- bamlss.formula(y ~ x1 + s(x2), family = "gaussian")
print(f)

## Specify model for parameter sigma, too
f <- list(
  y ~ x1 + s(x2),
  sigma ~ x2 + te(lon, lat)
)
f <- bamlss.formula(f, family = "gaussian")
print(f)

## With complex hierarchical structures,
## each parameter is another list of formulae,
## indicated by the h1,...,hk, names.
f <- list(
  y ~ x1 + s(x2) + id1,
  sigma ~ x2 + te(lon, lat) + id2,
  id1 ~ s(x3) + x4 + s(id3),
  id3 ~ x5 + s(x5, x6),
  id2 ~ x7
)
f <- bamlss.formula(f, family = "gaussian")
print(f)

Description

This function parses the data and the model formula, or extended bamlss.formula, as well as the bamlss.family into a bamlss.frame object. The bamlss.frame then holds all model matrices and information that is needed for setting up estimation engines.

Usage

bamlss.frame(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, ...)

Arguments

Details

The function parses the data, the formula or the extended bamlss.formula as well as the bamlss.family into a model frame like object, the bamlss.frame. This object holds all necessary model matrices and information that is needed for model fitting engines. Per default, all package mgcv smooth term constructor functions like s, te, t2 and ti can be used (see also function smooth.construct), however, even special user defined constructors can be included, see the examples below.

Function bamlss.frame() uses function model.matrix.bamlss.frame to compute all design matrices for simple linear parts, all smooth terms are parsed with function smooth.construct.bamlss.frame.

It is also possible to create a "bamlss.frame" using hierarchical formulae, see the example below.

Value

An list of class "bamlss.frame" with the following elements:

See Also

bamlss, bamlss.formula, bamlss.family, smooth.construct.bamlss.frame, model.matrix.bamlss.frame

Examples

## Create a 'bamlss.frame'.
d <- GAMart()
f <- list(
  num ~ fac + s(x1) + s(x2) + te(lon, lat),
  sigma ~ id + s(x2) + s(x3)
)
bf <- bamlss.frame(f, data = d, family = "gaussian")

## Show parts of the 'bamlss.frame'.
print(bf)

## Categorical responses.
f <- list(
  cat ~ fac + s(x1) + s(x2)
)

bf <- bamlss.frame(f, data = d, family = "multinomial", reference = "low")
print(bf)

## The response is a matrix per default.
head(bf$y)

## 0/1 responses.
d <- cbind(d, model.matrix(~ -1 + cat, data = d))

f <- list(
  catnone ~ fac + s(x1),
  catlow ~ s(x2),
  catmedium ~ s(x3)
)

bf <- bamlss.frame(f, data = d, family = "multinomial")
print(bf)

## Hierarchical structures.
f <- list(
  num ~ s(x1) + s(x2) + id,
  id ~ te(lon, lat),
  sigma ~ s(x1) + fac
)

bf <- bamlss.frame(f, data = d, family = "gaussian")
print(bf)

## Special model term constructors,
## set up "new" constructor function and eval
## with bamlss.frame().
s77 <- function(...) {
  sm <- s(...)
  sm$label <- paste("s77(", paste(sm$term, collapse = ","), ")", sep = "")
  sm
}

f <- list(
  num ~ s77(x1) + s(x2) + id,
  sigma ~ s77(x1)
)

bf <- bamlss.frame(f, data = d, family = "gaussian", specials = "s77")
print(bf)
names(bf$x$mu$smooth.construct)

Description

Wrapper function for applying bootstrap estimation using gradient boosting.

Usage

## Bootstrap boosting.
bboost(..., data, type = 1, cores = 1,
  n = 2, prob = 0.623, fmstop = NULL,
  trace = TRUE, drop = FALSE, replace = FALSE)

## Plotting function.
bboost_plot(object, col = NULL)

## Predict method.
## S3 method for class 'bboost'
predict(object, newdata, ..., cores = 1, pfun = NULL)

Arguments

Value

A list of bamlss objects.

See Also

bamlss, boost, lasso, BayesX

Examples

## Not run: ## Simulate data.
set.seed(123)
d <- GAMart()

## Estimate model.
f <- num ~ s(x1) + s(x2) + s(x3) + s(lon,lat)

## Function for evaluation of hold out sample
## criterion to find the optimum mstop.
fmstop <- function(model, data) {
  p <- predict(model, newdata = data, model = "mu")
  mse <- NULL
  for(i in 1:nrow(model$parameters))
    mse <- c(mse, mean((data$num - p[, i])^2))
  list("MSE" = mse, "mstop" = which.min(mse))
}

## Bootstrap boosted models.
b <- bboost(f, data = d, n = 50, cores = 3, fmstop = fmstop)

## Plot hold out sample MSE.
bboost_plot(b)

## Predict for each bootstrap sample.
nd <- data.frame("x2" = seq(0, 1, length = 100))
p <- predict(b, newdata = nd, model = "mu", term = "x2")
plot2d(p ~ x2, data = nd)

## End(Not run)

Description

Some simple shortcuts to model fitting engines.

Usage

## BayesX.
bayesx2(...)

## Gradient boosting.
boost2(...)

## Lasso.
lasso2(...)

Arguments

Value

A bamlss object.

See Also

bamlss, boost, lasso, BayesX

Examples

## Not run: ## Simulate data.
set.seed(123)
d <- GAMart()

## Estimate model.
f <- num ~ s(x1) + s(x2) + s(x3) + s(lon,lat)

## Boosted model.
b <- boost2(f, data = d)

## Plot estimated effects.
plot(b)

## End(Not run)

Description

Small helper function that computes the 2.5% and 97.5% quantiles and the mean of a vector. Useful for example when using function predict.bamlss.

Usage

c95(x)

Arguments

See Also

predict.bamlss, coef.bamlss

Examples

x <- rnorm(100)
c95(x)

## Not run: ## Example computing predictions.
set.seed(123)
d <- data.frame("x" = seq(-3, 3, length = 30))
d$y <- sin(d$x) + rnorm(30, sd = 0.3)

## Estimate model and compute predictions.
## with c95().
b <- bamlss(y ~ s(x), data = d)
p <- predict(b, model = "mu", FUN = c95)
plot(d)
matplot(d$x, p, type = "l", lty = c(2, 1, 2),
  col = "black", add = TRUE)

## Example extracting coefficients.
coef(b, FUN = c95)

## End(Not run)

Description

Methods to extract coefficients of fitted bamlss objects, either coefficients returned from optimizer functions, or samples from a sampler functions.

Method confint.bamlss() produces credible intervals or parameter samples using quantiles.

Usage

## S3 method for class 'bamlss'
coef(object, model = NULL, term = NULL,
  FUN = NULL, parameters = NULL,
  pterms = TRUE, sterms = TRUE,
  hyper.parameters = TRUE, list = FALSE,
  full.names = TRUE, rescale = FALSE, ...)

## S3 method for class 'bamlss'
confint(object, parm, level = 0.95,
  model = NULL, pterms = TRUE, sterms = FALSE,
  full.names = FALSE, hyper.parameters = FALSE, ...)

Arguments

Value

Depending on argument list and the number of distributional parameters, either a list or vector/matrix of model coefficients.

See Also

bamlss.

Examples

## Not run: ## Simulate data.
d <- GAMart()

## Model formula.
f <- list(
  num ~ s(x1) + s(x2) + s(x3),
  sigma ~ s(x1) + s(x2) + s(x3)
)

## Estimate model.
b <- bamlss(f, data = d)

## Extract coefficients based on MCMC samples.
coef(b)

## Now only the mean.
coef(b, FUN = mean)

## As list without the full names.
coef(b, FUN = mean, list = TRUE, full.names = FALSE)

## Coefficients only for "mu".
coef(b, model = "mu")

## And "s(x2)".
coef(b, model = "mu", term = "s(x2)")

## With optimizer parameters.
coef(b, model = "mu", term = "s(x2)", parameters = TRUE)

## Only parameteric part.
coef(b, sterms = FALSE, hyper.parameters = FALSE)

## For sigma.
coef(b, model = "sigma", sterms = FALSE,
  hyper.parameters = FALSE)

## 95 perc. credible interval based on samples.
confint(b)

## End(Not run)

Description

Function to generate a color legend, the legend may be added to an existing plot or drawn in a separate plotting window.

Usage

colorlegend(color = NULL, ncol = NULL, x = NULL, breaks = NULL, 
  pos = "center", shift = 0.02, side.legend = 1L, side.ticks = 1L,
  range = NULL, lrange = NULL, width = 0.25, height = 0.05,
  scale = TRUE, xlim = NULL, ylim = NULL, plot = NULL, full = FALSE,
  add = FALSE, col.border = "black", lty.border = 1L, lwd.border = 1L,
  ticks = TRUE, at = NULL, col.ticks = "black", lwd.ticks = 1L,
  lty.ticks = 1L, length.ticks = 0.3, labels = NULL,
  distance.labels = 0, col.labels = "black", cex.labels = 1L,
  digits = 2L, swap = FALSE, symmetric = TRUE, xpd = NULL,
  title = NULL, side.title = 2, shift.title = c(0, 0),
  cex.title = 1, ...)

Arguments

Value

A named list with the colors generated, the breaks and the function map, which may be used for mapping of x values to the colors specified in argument colors, please see the examples below.

Examples

## Play with colorlegend.
colorlegend()
colorlegend(side.legend = 2)
colorlegend(side.legend = 2, side.ticks = 2)
colorlegend(height = 2)
colorlegend(width = 1, height = 0.8, scale = FALSE, 
  pos = c(0, 0.2), length.ticks = 0.5)
colorlegend(color = heat.colors, ncol = 9)
colorlegend(color = heat.colors, ncol = 9, swap = TRUE)
colorlegend(pos = "bottomleft")
colorlegend(pos = "topleft")
colorlegend(pos = "topright")
colorlegend(pos = "bottomright")


## Take x values for the color legend.
x <- runif(100, -2, 2)
colorlegend(color = diverge_hcl, x = x)
colorlegend(color = diverge_hcl, x = x, at = c(-1.5, 0, 1.5))
colorlegend(color = diverge_hcl, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"))
colorlegend(color = rainbow_hcl, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5)
colorlegend(color = heat_hcl, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5,
  lwd.border = 2, lwd.ticks = 2, cex.labels = 1.5, font = 2)
colorlegend(color = topo.colors, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5,
  lwd.border = 2, lwd.ticks = 2, cex.labels = 1.5, font = 2,
  col.border = "green3", col.ticks = c(2, 5, 2), 
  col.labels = c(6, 4, 3))
colorlegend(color = diverge_hsv, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5,
  lwd.border = 2, lwd.ticks = 2, cex.labels = 1.5, font = 2,
  col.border = "green3", col.ticks = c(2, 5, 2), 
  col.labels = c(6, 4, 3), lty.border = 2, lty.ticks = c(2, 3, 2))
colorlegend(color = diverge_hsv, x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5,
  lwd.border = 2, lwd.ticks = 2, cex.labels = 1.5, font = 2,
  col.border = "green3", col.ticks = c(2, 5, 2), 
  col.labels = c(6, 4, 3), lty.border = 2, lty.ticks = c(2, 3, 2),
  ncol = 3)
colorlegend(color = c("red", "white", "red"), x = x, at = c(-1.5, 0, 1.5),
  labels = c("low", "middle", "high"), length.ticks = 1.5,
  lwd.border = 2, lwd.ticks = 2, cex.labels = 1.5, font = 2,
  col.border = "green3", col.ticks = c(2, 5, 2), 
  col.labels = c(6, 4, 3), lty.border = 2, lty.ticks = c(2, 3, 2),
  ncol = 3, breaks = c(-2, -1, 1, 2))
colorlegend(color = diverge_hcl, x = x, range = c(-3, 3))
colorlegend(color = diverge_hcl, x = x, range = c(-3, 3), lrange = c(-6, 6))


## Combine plot with color legend.
n <- 100
x <- y <- seq(-3, 3, length.out = n)
z <- outer(sin(x), cos(x)) 
pal <- colorlegend(color = diverge_hcl, x = z, plot = FALSE)
op <- par(no.readonly = TRUE)
par(mar = c(4.1, 4.1, 1.1, 1.1))
layout(matrix(c(1, 2), nrow = 1), widths = c(1, 0.3))
image(x = x, y = y, z = z, col = pal$colors, breaks = pal$breaks)
par(mar = c(4.1, 0.1, 1.1, 3.1))
colorlegend(color = diverge_hcl, x = z, plot = TRUE, full = TRUE,
  side.legend = 2, side.ticks = 2)
par(op)


## Another example with different plot.
n <- 50
x <- sin(seq(-3, 3, length.out = n)) 
pal <- colorlegend(color = diverge_hcl, x = x, plot = FALSE)
op <- par(no.readonly = TRUE)
par(mar = c(7.1, 4.1, 1.1, 1.1))
barplot(x, border = "transparent", col = pal$map(x))
colorlegend(color = diverge_hcl, x = x, plot = FALSE, add = TRUE,
  xlim = c(0, 60), ylim = c(-1, 1), pos = c(0, -0.15), xpd = TRUE,
  scale = FALSE, width = 60, height = 0.15,
  at = seq(min(x), max(x), length.out = 9))
par(op)

Description

This function takes a bamlss object which was created using a sampler function and continues sampling from the last state of the MCMC chain.

Usage

continue(object, cores = NULL, combine = TRUE,
  sleep = NULL, results = TRUE, ...)

Arguments

Value

A bamlss object.

See Also

bamlss

Examples

## Not run: ## Simulate and run model with MCMC.
set.seed(123)
d <- GAMart()
b <- bamlss(num ~ s(x1) + s(x2) + s(x3) + te(lon,lat), data = d)

## Continue sampling.
a <- continue(b)

## Plot all samples.
## plot(a, which = "samples")

## End(Not run)

Description

This function takes a fitted Cox model, i.e., a model estimated by opt_Cox or sam_Cox and computes predictions given a new data set or the original values. Survival probabilities are computed using numerical integration, therefore, computation may take some time. To avoid problems with computer memory, the prediction of survival probabilities can be split into chunks and computed parallel on different cores.

Usage

cox_predict(object, newdata,
  type = c("link", "parameter", "probabilities"),
  FUN = function(x) { mean(x, na.rm = TRUE) },
  time = NULL, subdivisions = 100, cores = NULL,
  chunks = 1, verbose = FALSE, ...)

Arguments

Value

Depending on the type of function provided in argument FUN, a numeric vector or matrix.

See Also

sam_Cox, cox_bamlss, surv_transform, simSurv, bamlss, predict.bamlss

Examples

## Not run: library("survival")
set.seed(123)

## Simulate survival data.
d <- simSurv(n = 500)

## Formula of the survival model, note
## that the baseline is given in the first formula by s(time).
f <- list(
  Surv(time, event) ~ s(time) + s(time, by = x3),
  gamma ~ s(x1) + s(x2)
)

## Cox model with continuous time.
## Note the the family object cox_bamlss() sets
## the default optimizer and sampler function!
## First, posterior mode estimates are computed
## using function opt_Cox(), afterwards the
## sampler sam_Cox() is started.
b <- bamlss(f, family = "cox", data = d)

## Predict survival probabilities P(T > t).
p <- predict(b, type = "probabilities",
  time = 3, subdivisions = 100, FUN = c95)

## End(Not run)

Description

This function creates simulated data based on a crazy function.

Usage

Crazy(n = 1000)

Arguments

Value

A data frame with response y and covariate x.

See Also

GAMart

Examples

d <- Crazy(1000)
head(d)
plot(d)

Description

The function computes the continuous rank probability score (CRPS). Note that the function uses numerical integration, for highly efficient computation please see the scoringRules package.

Usage

CRPS(object, newdata = NULL,
  interval = c(-Inf, Inf), FUN = mean,
  term = NULL, ...)

Arguments

References

Gneiting T, Raftery AE (2007). Strictly Proper Scoring Rules, Prediction, and Estimation." Journal of the American Statistical Association, 102(477), 359–378. doi:10.1198/016214506000001437cd ..

Gneiting T, Balabdaoui F, Raftery AE (2007). Probabilistic Forecasts, Calibration and Sharpness. Journal of the Royal Statistical Society B, 69(2), 243–268. doi:10.1111/j.1467-9868.2007.00587.x

Examples

## Not run: ## Simulate data.
d <- GAMart()

## Model only including covariate x1.
b1 <- bamlss(num ~ s(x1), data = d)

## Now, also including x2 and x2.
b2 <- bamlss(num ~ s(x1) + s(x2) + s(x3), data = d)

## Compare using the CRPS score.
CRPS(b1)
CRPS(b2)

## End(Not run)

Description

This function interfaces keras infrastructures for high-level neural networks. The function can be used as a standalone model fitting engine such as bamlss or as an on top model engine to capture special features in the data that could not be captures by other model fitting engines.

Usage

## Deep distributional neural net.
ddnn(object, optimizer = "adam",
  learning_rate = 0.01,
  epochs = 100, batch_size = NULL,
  nlayers = 2, units = 100, activation = "relu",
  l1 = NULL, l2 = NULL,
  validation_split = 0.2, early_stopping = TRUE, patience = 50,
  verbose = TRUE, ...)

## Predict method.
## S3 method for class 'ddnn'
predict(object, newdata,
  model = NULL, type = c("link", "parameter"),
  drop = TRUE, ...)

## CV method for optimizing
## the number of epochs using
## the CRPS.
cv_ddnn(formula, data, folds = 10,
  min_epochs = 300, max_epochs = 400,
  interval = c(-Inf, Inf), ...)

Arguments

Details

The default keras model is a sequential model with two hidden layers with "relu" activation function and 100 units in each layer. Between each layer is a dropout layer with 0.1 dropout rate.

Value

For function ddnn() an object of class "ddnn". Note that extractor functions fitted and residuals.bamlss can be applied. For function predict.ddnn() a list or vector of predicted values.

WARNINGS

The deep learning infrastructure is experimental!

See Also

bamlss.frame, bamlss

Examples

## Not run: ## Simulate data.
set.seed(123)
n <- 300
x <- runif(n, -3, 3)
fsigma <- -2 + cos(x)
y <- sin(x) + rnorm(n, sd = exp(fsigma))

## Setup model formula.
f <- list(
  y ~ x,
  sigma ~ x
)

## Fit neural network.
library("keras")
b <- ddnn(f, epochs = 2000)

## Plot estimated functions.
par(mfrow = c(1, 2))
plot(x, y)
plot2d(fitted(b)$mu ~ x, add = TRUE)
plot2d(fitted(b)$sigma ~ x,
  ylim = range(c(fitted(b)$sigma, fsigma)))
plot2d(fsigma ~ x, add = TRUE, col.lines = "red")

## Predict with newdata.
nd <- data.frame(x = seq(-6, 6, length = 100))
nd$p <- predict(b, newdata = nd, type = "link")

par(mfrow = c(1, 2))
plot(x, y, xlim = c(-6, 6), ylim = range(c(nd$p$mu, y)))
plot2d(p$mu ~ x, data = nd, add = TRUE)
plot2d(p$sigma ~ x, data = nd,
  ylim = range(c(nd$p$sigma, fsigma)))
plot2d(fsigma ~ x, add = TRUE, col.lines = "red")

## Plot quantile residuals.
e <- residuals(b)
plot(e)

## End(Not run)

Description

Generic function returning the deviance information criterion (DIC) of a fitted model object.

Usage

DIC(object, ...)

Arguments

Examples

## Not run: d <- GAMart()
b1 <- bamlss(num ~ s(x1), data = d)
b2 <- bamlss(num ~ s(x1) + s(x2), data = d)
DIC(b1, b2)

## End(Not run)

Description

disttree Families for MVN with Cholesky Parameterization

Usage

dist_mvnchol(k, r = k - 1L, type = c("basic", "modified", "chol"), ...)

Arguments

Details

NOTE: These functions are under development!! disttree families that models a multivariate Normal (Gaussian) distribution by (modified) Cholesky decomposition of the covariance matrix.

Value

a bamlss family.

Description

The function shows available optimizer and sampling engines for a given family object.

Usage

engines(family, ...)

Arguments

Examples

engines(gaussian_bamlss, "gamma", cox_bamlss)

Description

Family objects in bamlss specify the information that is needed for using (different) model fitting engines, e.g., the parameter names and corresponding link functions, the density function, derivatives of the log-likelihood w.r.t. the predictors, and so forth. The optimizer or sampler functions that are called by bamlss must know how much information is needed to interpret the model since the family objects are simply passed through. Family objects are also used for computing post-modeling statistics, e.g., for residual diagnostics or random number generation. See the details and examples.

Usage

## Family objects in bamlss:
ALD_bamlss(..., tau = 0.5, eps = 0.01)
beta_bamlss(...)
binomial_bamlss(link = "logit", ...)
cnorm_bamlss(...)
cox_bamlss(...)
dw_bamlss(...)
DGP_bamlss(...)
dirichlet_bamlss(...)
ELF_bamlss(..., tau = 0.5)
gaussian_bamlss(...)
gaussian2_bamlss(...)
Gaussian_bamlss(...)
gamma_bamlss(...)
logNN_bamlss(...)
multinomial_bamlss(...)
mvnorm_bamlss(k = 2, ...)
mvnormAR1_bamlss(k = 2, ...)
poisson_bamlss(...)
gpareto_bamlss(...)
glogis_bamlss(...)
AR1_bamlss(...)
beta1_bamlss(ar.start, ...)
nbinom_bamlss(...)
ztnbinom_bamlss(...)
lognormal_bamlss(...)
weibull_bamlss(...)
Sichel_bamlss(...)
GEV_bamlss(...)
gumbel_bamlss(...)
mix_bamlss(f1, f2, ...)
ZANBI_bamlss(...)

## Extractor functions:
## S3 method for class 'bamlss'
family(object, ...)
## S3 method for class 'bamlss.frame'
family(object, ...)

Arguments

Details

The following lists the minimum requirements on a bamlss family object to be used with bamlss and bamlss.frame:

• The family object must return a list of class "family.bamlss".

• The object must contain the family name as a character string.

• The object must contain the names of the parameters as a character string, as well as the corresponding link functions as character string.

For most optimizer and sampling functions at least the density function, including a log argument, should be provided. When using generic model fitting engines like opt_bfit or sam_GMCMC, as well as for computing post-modeling statistics with function samplestats, and others, it is assumed that the density function in a family object has the following arguments:

d(y, par, log = FALSE, ...)

where argument y is the response (possibly a matrix) and par is a named list holding the evaluated parameters of the distribution, e.g., using a normal distribution par has two elements, one for the mean par$mu and one for the standard deviation par$sigma. The dots argument is for passing special internally used objects, depending on the type of model this feature is usually not needed.

Similarly, for derivative based algorithms, e.g. using iteratively weighted least squares (IWLS, see function opt_bfit, the family object holds derivative functions evaluating derivatives of the log-likelihood w.r.t. the predictors (or expectations of derivatives). For each parameter, these functions also hold the following arguments:

score(y, par, ...)

for computing the first derivative of the log-likelihood w.r.t. a predictor and

hess(y, par, ...)

for computing the negative second derivatives. Within the family object these functions are organized in a named list, see the examples below.

In addition, for the cumulative distribution function (p(y, par, ...)), for the quantile function (q(y, par, ...)) or for creating random numbers (r(n, par, ...)) the same structure is assumed. See, e.g., the code of function gaussian.bamlss().

Some model fitting engines can initialize the distributional parameters which oftentimes leads to much faster convergence. The initialize functions are again organized within a named list, one entry for each parameter, similar to the score and hess functions, e.g., see the code of family object gaussian.bamlss.

Using function bamlss, residuals.bamlss and predict.bamlss the family objects may also specify the transform()er, optimizer(), sampler(), samplestats(), results(), residuals() and predict() function that should be used with this family. See for example the setup of cox_bamlss.

For using specialized estimation engines like sam_JAGS it is recommended to supply any extra arguments needed by those engines with an additional list entry within the family object, e.g., using gaussian_bamlss() with sam_JAGS the family objects holds special details in an element named "bugs".

The examples below illustrate this setup. See also the code of the bamlss family functions.

See Also

bamlss, bamlss.frame

Examples

## New family object for the normal distribution,
## can be used by function opt_bfit() and sam_GMCMC().
normal_bamlss <- function(...) {
  f <- list(
    "family" = "normal",
    "names" = c("mu", "sigma"),
    "links" = c("identity", "log"),
    "d" = function(y, par, log = FALSE) {
      dnorm(y, mean = par$mu, sd = par$sigma, log = log)
    },
    "score" = list(
      "mu" = function(y, par, ...) {
        drop((y - par$mu) / (par$sigma^2))
      },
      "sigma" = function(y, par, ...) {
        drop(-1 + (y - par$mu)^2 / (par$sigma^2))
      }
    ),
    "hess" = list(
      "mu" = function(y, par, ...) {
        drop(1 / (par$sigma^2))
      },
      "sigma" = function(y, par, ...) { 
        rep(2, length(y))
      }
    )
  )
  class(f) <- "family.bamlss"
  return(f)
}

## Not run: ## Test on simulated data.
d <- GAMart()
b <- bamlss(num ~ s(x1) + s(x2) + s(x3),
  data = d, family = "normal")
plot(b)

## Compute the log-likelihood using the family object.
f <- family(b)
sum(f$d(y = d$num, par = f$map2par(fitted(b)), log = TRUE))

## For using JAGS() more details are needed.
norm4JAGS_bamlss <- function(...) {
  f <- normal_bamlss()
  f$bugs <- list(
    "dist" = "dnorm",
    "eta" = BUGSeta,
    "model" = BUGSmodel,
    "reparam" = c(sigma = "1 / sqrt(sigma)")
  )
  return(f)
}

## Now with opt_bfit() and sam_JAGS().
b <- bamlss(num ~ s(x1) + s(x2) + s(x3), data = d,
  optimizer = opt_bfit, sampler = sam_JAGS, family = "norm4JAGS")
plot(b)

## End(Not run)

Description

This data set includes weekly fatalities in Austria from 2000 to 46 weeks in 2020. The data is taken from the Eurostat data base.

Usage

data("fatalities")

Format

The fatalities data contains the following variables:

num:

Integer, the number of fatalities.

year:

Integer, the corresponding year fatalities are recorded.

week:

Integer, the corresponding week fatalities are recorded..

References

Eurostat Database (2020). Population and social conditions, demography and migration, mortality, weekly deaths, deaths by week and NUTS 3 region, Austria https://ec.europa.eu/eurostat/

Examples

data("fatalities")
plot(num ~ week, data = fatalities)

Description

Function to compute fitted values for bamlss models. The function calls predict.bamlss to compute fitted values from samples.

Usage

## S3 method for class 'bamlss'
fitted(object, model = NULL, term = NULL,
  type = c("link", "parameter"), samples = TRUE,
  FUN = c95, nsamps = NULL, ...)

Arguments

Value

Depending on arguments model, FUN and the structure of the bamlss model, a list of fitted values or simple vectors or matrices of fitted values.

See Also

bamlss, predict.bamlss.

Examples

## Not run: ## Generate some data.
d <- GAMart()

## Model formula.
f <- list(
  num ~ s(x1) + s(x2) + s(x3) + te(lon,lat),
  sigma ~ s(x1) + s(x2) + s(x3) + te(lon,lat)
)

## Estimate model.
b <- bamlss(f, data = d)

## Fitted values returned from optimizer.
f1 <- fitted(b, model = "mu", samples = FALSE)

## Fitted values returned from sampler.
f2 <- fitted(b, model = "mu", samples = TRUE, FUN = mean)

plot(f1, f2)

## End(Not run)

Description

This function creates artificial GAM-type data.frames. The function is mainly used for testing purposes.

Usage

GAMart(n = 500, sd = 0.1, seed = FALSE,
  ti = c("none", "vcm", "main", "both"))

Arguments

Examples

d <- GAMart()
head(d)

## Not run: b <- bamlss(num ~ s(x1) + s(x2) + s(x3) + te(lon,lat), data = d)
plot(b)

## End(Not run)

Description

The functions searches in the gamlss.dist namespace for available distributions. It returns a named list of faily functions which can be used with bamlss.

Usage

gamlss_distributions(type = c("continuous", "discrete"))

Arguments

See Also

bamlss

Examples

## Not run: 
dists <- gamlss_distributions(type = "continuous")
print(dists)

## End(Not run)

Description

Function to get a family.bamlss object to be used for fitting. The main purpose of this function is to ease the handling of extra arguments to the family object.

Usage

gF(x, ...)

Arguments

Value

A family.bamlss object.

See Also

family.bamlss.

Examples

f <- gF(gaussian)
print(f)

Description

This dataset is taken from the Regression Book and is about prices of used VW Golf cars.

Usage

data("Golf")

Format

The Golf data contains the following variables:

price:

Numeric, sale price in 1000 Euro.

age:

Numeric, age of the car in month.

kilometer:

Numeric, kilometer reading in 1000 kilometers.

tia:

Numeric, month until the next TIA appointment (German TUEV).

abs:

Factor, does the car have abs?

sunroof:

Factor, does the car have a sunroof?

References

Fahrmeir, L., Kneib, T., Lang, S. and Marx, B. (2013). Regression - Models, Methods and Applications, Springer. https://www.uni-goettingen.de/de/551357.html.

Examples

data("Golf")
plot(price ~ age, data = Golf)
plot(price ~ kilometer, data = Golf)

Description

This function downloads and compiles the HOMSTART-project data set. The data is downloaded from the Zentralanstalt fuer Meteorologie und Geodynamik (ZAMG, http://www.zamg.ac.at) and funded by the Austrian Climate Research Programme (ACRP) and is free for research purposes.

Usage

homstart_data(dir = NULL, load = TRUE, tdir = NULL)

Arguments

Value

A data frame containing the following variables:

References

Nemec J, Gruber C, Chimani B, Auer I (2012). Trends in extreme temperature indices in Austria based on a new homogenised dataset. International Journal of Climatology. DOI 10.1002/joc.3532.

Nemec J, Chimani B, Gruber C, Auer I (2011). Ein neuer Datensatz homogenisierter Tagesdaten. OEGM Bulletin, 1, 19–20. https://www.meteorologie.at/docs/OEGM_bulletin_2011_1.pdf

Umlauf N, Mayr G, Messner J, Zeileis A (2012). Why does it always rain on me? A spatio-temporal analysis of precipitation in Austria. Austrian Journal of Statistics, 41(1), 81–92. doi:10.17713/ajs.v41i1.190

Examples

## Not run: homstart_data(load = TRUE)
head(homstart)

## End(Not run)

Description

Family object to fit a flexible additive joint model for longitudinal and survival data under a Bayesian approach as presented in Koehler et al. (2017a, b). All parts of the joint model can be specified as structured additive predictors. See the details and examples.

Usage

## JM family object.
jm_bamlss(...)

## "bamlss.frame" transformer function
## to set up joint models.
jm_transform(x, y, data, terms, knots, formula, family, subdivisions = 25,
  timedependent = c("lambda", "mu", "alpha", "dalpha"), timevar = NULL, 
  idvar = NULL, alpha = .Machine$double.eps, mu = NULL, sigma = NULL, 
  sparse = TRUE, nonlinear = FALSE, edf_alt = FALSE, start_mu = NULL, 
  k_mu = 6, ...)
                    
## Posterior mode optimizing engine.     
opt_JM(x, y, start = NULL, weights = NULL, offset = NULL, 
  criterion = c("AICc", "BIC", "AIC"), maxit = c(100, 1), 
  nu = c("lambda" = 0.1, "gamma" = 0.1, "mu" = 1, "sigma" = 1,
    "alpha" = 1, "dalpha" = 1),
  update.nu = FALSE, eps = 0.0001, alpha.eps = 0.001, ic.eps = 1e-08, 
  nback = 40, verbose = TRUE, digits = 4, ...)

jm_mode(x, y, start = NULL, weights = NULL, offset = NULL, 
  criterion = c("AICc", "BIC", "AIC"), maxit = c(100, 1), 
  nu = c("lambda" = 0.1, "gamma" = 0.1, "mu" = 1, "sigma" = 1,
    "alpha" = 1, "dalpha" = 1),
  update.nu = FALSE, eps = 0.0001, alpha.eps = 0.001, ic.eps = 1e-08, 
  nback = 40, verbose = TRUE, digits = 4, ...)

## Sampler function.
sam_JM(x, y, family, start = NULL, weights = NULL, offset = NULL, 
  n.iter = 1200, burnin = 200, thin = 1, verbose = TRUE, digits = 4, 
  step = 20, ...)

jm_mcmc(x, y, family, start = NULL, weights = NULL, offset = NULL, 
  n.iter = 1200, burnin = 200, thin = 1, verbose = TRUE, digits = 4, 
  step = 20, ...)
  
## Predict function, set to default in jm_bamlss().
jm_predict(object, newdata,
  type = c("link", "parameter", "probabilities", "cumhaz", "loglik"),
  dt, steps, id, FUN = function(x) { mean(x, na.rm = TRUE) }, 
  subdivisions = 100, cores = NULL, chunks = 1,
  verbose = FALSE,  ...)
  
## Survival plot.
jm_survplot(object, id = 1, dt = NULL, steps = 10, 
  points = TRUE, rug = !points)

Arguments

Details

We refer to the papers of Koehler et al. (2017a, b) for details on the flexible additive joint model. In short, we model the hazard of subject $i$ an event at time $t$ as

$h_{i}(t)= \exp [\eta_{\lambda i}(t)+ \eta_{\gamma i}+\eta_{\alpha i}(\eta_{\mu i}(t), t) ]$

with predictor $\eta_{\lambda}$ for all survival covariates that are time-varying or have a time-varying coefficient (including the log baseline hazard), predictor $\eta_{\gamma}$ for baseline survival covariates, predictor $\eta_{\alpha}$ representing the potentially time-varying or nonlinear association between the longitudinal marker $\eta_{\mu}$ and the hazard. The longitudinal response $y_{ij}$ at time points $t_{ij}$ is modeled as

$y_{ij}=\eta_{\mu i}(t_{ij})+e_{ij}$

with independent normal errors $N(0, \exp[\eta_{\sigma i}(t_{ij})]^2)$.

Each predictor $\eta_{ki}$ is a structured additive predictor, i.e. a sum of functions of covariates $\eta_{ki} = \sum_{m=1}^{M_k} f_{km}(\bm{x}_{ki})$. Each of these functions can be modeled parametrically or using basis function evaluations from the smooth constructors in mgcv such as s, te and ti and can include smooth time-varying, random or spatial effects. For the Bayesian estimation of these effects we specify corresponding priors: For linear or parametric terms we use vague normal priors, smooth and random effect terms are regularized by placing generic multivariate normal priors on the coefficients and for anisotropic smooths, when multiple smoothing variance parameters are involved, more complex prior are in place (cf. Koehler et al., 2017a). We use inverse Gamma hyper-priors, i.e. IG(0.001, 0.001) to obtain an inverse Gamma full conditional for the variance parameters. We estimate the posterior mode by maximizing the log-posterior of the model using a Newton-Raphson procedure, the posterior mean is obtained via derivative-based Metropolis-Hastings sampling. We recommend to use posterior mode estimates for a quick model assessment. In order to draw correct inferences from the model, posterior mean estimates should be computed. We approximate integration in the survival part of the likelihood using trapezoidal rule. For posterior mode estimation.

A variety specifications of the association $\eta_{\alpha i}(\eta_{\mu i}(t), t)$ are possible with an important distinction between associations which are nonlinear in $\eta_{\mu}$ for nonlinear = TRUE (Koehler et al. 2017b) or linear where $\eta_{\alpha i}(\eta_{\mu i}(t), t) = \eta_{\alpha i}(t)\eta_{\mu i}(t)$ for nonlinear = FALSE (Koehler et al. 2017a).

Note

The indicator nonlinear for associations with are linear or nonlinear in $\eta_{\mu}$ was named interaction] in earlier versions stages of the development.

References

Koehler M, Umlauf N, Beyerlein, A., Winkler, C. Ziegler, A.-G., Greven S (2017). Flexible Bayesian Additive Joint Models with an Application to Type 1 Diabetes Research. Biometrical Journal. doi:10.1002/bimj.201600224

Meike Koehler, Nikolaus Umlauf, and Sonja Greven (2018). Nonlinear association structures in flexible Bayesian additive joint models. Statistics in Medicine. doi:10.1002/sim.7967

See Also

bamlss, bamlss.frame.

Examples

## Not run: 

set.seed(123)
## Simulate survival data
## with random intercepts/slopes and a linear effect of time,
## constant association alpha and no effect of the derivative
d <- simJM(nsub = 200, long_setting = "linear", 
  alpha_setting = "constant",
  dalpha_setting = "zero", full = FALSE)

## Formula of the according joint model
f <- list(
  Surv2(survtime, event, obs = y) ~ s(survtime, bs = "ps"),
  gamma ~ s(x1, bs = "ps"),
  mu ~ obstime + s(id, bs = "re") + 
    s(id, obstime, bs = "re"),
  sigma ~ 1,
  alpha ~ 1,
  dalpha ~ -1
)

## Joint model estimation
## jm_bamlss() sets the default optimizer and sampler function.
## First, posterior mode estimates are computed using function
## opt_JM(), afterwards the sampler sam_JM() is started.
b <- bamlss(f, data = d, family = "jm",
           timevar = "obstime", idvar = "id")

## Plot estimated effects.
plot(b)

## Predict event probabilities for two individuals
## at 12 time units after their last longitudinal measurement.
## The event probability is conditional on their survival
## up to their last observed measurement.
p <- predict(b, type = "probabilities", id = c(1, 2), dt = 12, FUN = c95)
print(p)

## Plot of survival probabilities and
## corresponding longitudinal effects
## for individual id.
jm_survplot(b, id = 3)
jm_survplot(b, id = 30)


## Simulate survival data
## with functional random intercepts and a nonlinear effect 
## of time, time-varying association alpha and no effect 
## of the derivative.
## Note: This specification is the simJM default.
d <- simJM(nsub = 200, full = FALSE)

## Formula of the according joint model
## number of knots for the smooth nonlinear effect of time
klong <- 8
f <- list(
  Surv2(survtime, event, obs = y) ~ s(survtime, bs = "ps"),
  gamma ~ s(x1, bs = "ps"),
  mu ~ ti(id, bs = "re") + 
    ti(obstime, bs = "ps", k = klong) + 
    ti(id, obstime, bs = c("re", "ps"), 
      k = c(nlevels(d$id), klong)) + 
    s(x2, bs = "ps"),
  sigma ~ 1,
  alpha ~ s(survtime, bs = "ps"),
  dalpha ~ -1
)

## Estimating posterior mode only using opt_JM()
b_mode <- bamlss(f, data = d, family = "jm",
                 timevar = "obstime", idvar = "id", 
                 sampler = FALSE)

## Estimating posterior means using sam_JM() 
##  with starting values generated from posterior mode 
b_mean <- bamlss(f, data = d, family = "jm",
                 timevar = "obstime", idvar = "id", optimizer = FALSE,
                 start = parameters(b_mode), results = FALSE)

## Plot effects.
plot(b_mean, model = "alpha")



## Simulate survival data
## with functional random intercepts and an association nonlinear in mu
set.seed(234)
d <- simJM(nsub = 300, long_setting = "functional", alpha_setting = "nonlinear", 
           nonlinear = TRUE, full = FALSE, probmiss = 0.9)

## Calculate longitudinal model to obtain starting values for mu
long_df <- 7
f_start <- y ~ ti(id, bs = "re") + ti(obstime, bs = "ps", k = long_df) + 
               ti(id, obstime, bs = c("re", "ps"), k = c(nlevels(d$id), long_df)) + 
               s(x2, bs = "ps")
b_start <- bamlss(f_start, data = d, sampler = FALSE)
mu <- predict(b_start)$mu

## Fit joint model with nonlinear association (nonlinear = TRUE)
f <- list(
  Surv2(survtime, event, obs = y) ~ s(survtime, bs = "ps"),
  gamma ~ x1,
  mu ~ ti(id, bs = "re") + ti(obstime, bs = "ps", k = long_df) + 
       ti(id, obstime, bs = c("re", "ps"), k = c(nlevels(d$id), long_df)) + 
       s(x2, bs = "ps"),
  sigma ~ 1,
  alpha ~ 1,
  dalpha ~ -1
)
b <- bamlss(f, data = d, family = "jm", timevar = "obstime", idvar = "id", 
            nonlinear = TRUE, start_mu = mu, 
            n.iter = 6000, burnin = 2000, thin = 2)
plot(b)
samplestats(b$samples)

## End(Not run)

Description

Smooth constructors and optimizer for Lasso penalization with bamlss. The penalization is based on a Taylor series approximation of the Lasso penalty.

Usage

## Smooth constructor function.
la(formula, type = c("single", "multiple"), ...)

## Single Lasso smoothing parameter optimizer.
opt_lasso(x, y, start = NULL, adaptive = TRUE, lower = 0.001, upper = 1000,
  nlambda = 100, lambda = NULL,  multiple = FALSE, verbose = TRUE,
  digits = 4, flush = TRUE, nu = NULL, stop.nu = NULL,
  ridge = .Machine$double.eps^0.5, zeromodel = NULL, ...)

lasso(x, y, start = NULL, adaptive = TRUE, lower = 0.001, upper = 1000,
  nlambda = 100, lambda = NULL,  multiple = FALSE, verbose = TRUE,
  digits = 4, flush = TRUE, nu = NULL, stop.nu = NULL,
  ridge = .Machine$double.eps^0.5, zeromodel = NULL, ...)

## Lasso transformation function to set
## adaptive weights from an unpenalized model.
lasso_transform(x, zeromodel, nobs = NULL, ...)

## Plotting function for opt_lasso() optimizer.
lasso_plot(x, which = c("criterion", "parameters"),
  spar = TRUE, model = NULL, name = NULL, mstop = NULL,
  retrans = FALSE, color = NULL, show.lambda = TRUE,
  labels = NULL, digits = 2, ...)

## Extract optimum stopping iteration for opt_lasso() optimizer.
## Based on the minimum of the information criterion.
lasso_stop(x)

## Extract retransformed Lasso coefficients.
lasso_coef(x, ...)