Package 'lgpr'

Description

Interpretable nonparametric modeling of longitudinal data using additive Gaussian process regression. Contains functionality for inferring covariate effects and assessing covariate relevances. Models are specified using a convenient formula syntax, and can include shared, group-specific, non-stationary, heterogeneous and temporally uncertain effects. Bayesian inference for model parameters is performed using 'Stan' (rstan). The modeling approach and methods are described in detail in Timonen et al. (2021).

Core functions

Main functionality of the package consists of creating and fitting an additive GP model:

• lgp: Specify and fit an additive GP model with one command.

• create_model: Define an lgpmodel object.

• sample_model: Fit a model by sampling the posterior distribution of its parameters and create an lgpfit object.

• pred: Computing model predictions and inferred covariate effects after fitting a model.

• relevances: Assessing covariate relevances after fitting a model.

• prior_pred: Prior predictive sampling to check if your prior makes sense.

Visualization

• plot_pred: Plot model predictions.

• plot_components: Visualize inferred model components.

• plot_draws: Visualize parameter draws.

• plot_data: Visualize longitudinal data.

Data

The data that you wish to analyze with 'lgpr' should be in an R data.frame where columns correspond to measured variables and rows correspond to observations. Some functions that can help working with such data frames are:

• new_x: Creating new test points where the posterior distribution of any function component or sum of all components, or the posterior predictive distribution can be computed after model fitting.

• Other functions: add_factor, add_factor_crossing, add_dis_age, adjusted_c_hat.

Vignettes and tutorials

See https://jtimonen.github.io/lgpr-usage/index.html. The tutorials focus on code and use cases, whereas the Mathematical description of lgpr models vignette describes the statistical models and how they can be customized in 'lgpr'.

Citation

Run citation("lgpr") to get citation information.

Feedback

Bug reports, PRs, enhancement ideas or user experiences in general are welcome and appreciated. Create an issue in Github or email the author.

Author(s)

Juho Timonen (first.last at iki.fi)

References

1. Timonen, J. et al. (2021). lgpr: an interpretable non-parametric method for inferring covariate effects from longitudinal data. Bioinformatics, url.

2. Carpenter, B. et al. (2017). Stan: A probabilistic programming language. Journal of Statistical Software 76(1).

See Also

Useful links:

• https://github.com/jtimonen/lgpr

• Report bugs at https://github.com/jtimonen/lgpr/issues

Description

Creates the disease-related age covariate vector based on the disease initiation times and adds it to the data frame

Usage

add_dis_age(data, t_init, id_var = "id", time_var = "age")

Arguments

Value

A data frame with one column added. The new column will be called dis_age. For controls, its value will be NaN.

See Also

Other data frame handling functions: add_factor(), add_factor_crossing(), adjusted_c_hat(), new_x(), split()

Description

Easily add a categorical covariate to a data frame

Usage

add_factor(data, x, id_var = "id")

Arguments

Value

A data frame with one column added. The new column will have same name as the variable passed as input x.

See Also

Other data frame handling functions: add_dis_age(), add_factor_crossing(), adjusted_c_hat(), new_x(), split()

Description

Add a crossing of two factors to a data frame

Usage

add_factor_crossing(data, fac1, fac2, new_name)

Arguments

Value

a data frame

See Also

Other data frame handling functions: add_dis_age(), add_factor(), adjusted_c_hat(), new_x(), split()

Description

Creates the c_hat input for lgp, so that it accounts for normalization between data points in the "poisson" or "nb" observation model

Usage

adjusted_c_hat(y, norm_factors)

Arguments

Value

a vector of length n, which can be used as the c_hat input to the lgp function

See Also

Other data frame handling functions: add_dis_age(), add_factor(), add_factor_crossing(), new_x(), split()

Description

Apply variable scaling

Usage

apply_scaling(scaling, x, inverse = FALSE)

Arguments

Value

a similar object as x

See Also

Other variable scaling functions: create_scaling()

Description

Character representations of different formula objects

Usage

## S4 method for signature 'lgpexpr'
as.character(x)

## S4 method for signature 'lgpterm'
as.character(x)

## S4 method for signature 'lgpformula'
as.character(x)

Arguments

Value

a character representation of the object

Description

See the Mathematical description of lgpr models vignette for more information about the connection between different options and the created statistical model.

Usage

create_model(
  formula,
  data,
  likelihood = "gaussian",
  prior = NULL,
  c_hat = NULL,
  num_trials = NULL,
  options = NULL,
  prior_only = FALSE,
  verbose = FALSE,
  sample_f = !(likelihood == "gaussian")
)

Arguments

Value

An object of class lgpmodel, containing the Stan input created based on parsing the specified formula, prior, and other options.

See Also

Other main functions: draw_pred(), get_draws(), lgp(), pred(), prior_pred(), sample_model()

Description

Parse the covariates and model components from given data and formula

Usage

create_model.covs_and_comps(data, model_formula, x_cont_scl, verbose)

Arguments

Value

parsed input to Stan and covariate scaling, and other info

See Also

Other internal model creation functions: create_model.formula(), create_model.likelihood(), create_model.prior()

Description

Checks if formula is in advanced format and translates if not.

Usage

create_model.formula(formula, data, verbose = FALSE)

Arguments

Value

an object of class lgpformula

See Also

Other internal model creation functions: create_model.covs_and_comps(), create_model.likelihood(), create_model.prior()

Description

Parse the response variable and its likelihood model

Usage

create_model.likelihood(
  data,
  likelihood,
  c_hat,
  num_trials,
  y_name,
  sample_f,
  verbose
)

Arguments

Value

a list of parsed options

See Also

Other internal model creation functions: create_model.covs_and_comps(), create_model.formula(), create_model.prior()

Description

Parse the given modeling options

Usage

create_model.options(options, verbose)

Arguments

Value

a named list of parsed options

Description

Parse given prior

Usage

create_model.prior(prior, stan_input, verbose)

Arguments

Value

a named list of parsed options

See Also

Other internal model creation functions: create_model.covs_and_comps(), create_model.formula(), create_model.likelihood()

Description

Helper function for plots

Usage

create_plot_df(object, x = "age", group_by = "id")

Arguments

Value

a data frame

Description

Create a standardizing transform

Usage

create_scaling(x, name)

Arguments

Value

an object of class lgpscaling

See Also

Other variable scaling functions: apply_scaling()

Description

Using the same parametrization as Stan. More info here.

Usage

dinvgamma_stanlike(x, alpha, beta, log = FALSE)

qinvgamma_stanlike(p, alpha, beta)

Arguments

Value

density/quantile value

See Also

Other functions related to the inverse-gamma distribution: plot_invgamma(), priors

Description

Draw pseudo-observations from predictive distribution. If pred contains draws from the component posterior (prior) distributions, then the output is draws from the posterior (prior) predictive distribution. If pred is not specified, then whether output draws are from prior or posterior predictive distribution depends on whether fit is created using the lgp option prior_only=TRUE or not.

Usage

draw_pred(fit, pred = NULL)

Arguments

Value

An array with shape $S x P$, where $S$ is the number of draws that pred contains and $P$ is the length of each function draw. Each row $s = 1, \ldots, S$ of the output is one vector drawn from the predictive distribution, given parameter draw $s$.

See Also

Other main functions: create_model(), get_draws(), lgp(), pred(), prior_pred(), sample_model()

Description

Quick way to create an example lgpfit, useful for debugging

Usage

example_fit(
  formula = y ~ id + age + age | SEX + age | LOC,
  likelihood = "gaussian",
  chains = 1,
  iter = 30,
  num_indiv = 6,
  num_timepoints = 5,
  ...
)

Arguments

Value

An lgpfit object created by fitting the example model.

Description

Print a fit summary.

Usage

fit_summary(fit, ignore_pars = c("f_latent", "eta", "teff_raw", "lp__"))

Arguments

Value

object invisibly.

Description

An S4 class to represent analytically computed predictive distributions (conditional on hyperparameters) of an additive GP model

Usage

## S4 method for signature 'GaussianPrediction'
show(object)

## S4 method for signature 'GaussianPrediction'
component_names(object)

## S4 method for signature 'GaussianPrediction'
num_components(object)

## S4 method for signature 'GaussianPrediction'
num_paramsets(object)

## S4 method for signature 'GaussianPrediction'
num_evalpoints(object)

Arguments

Methods (by generic)

• show(GaussianPrediction): Print a summary about the object.

• component_names(GaussianPrediction): Get names of components.

• num_components(GaussianPrediction): Get number of components.

• num_paramsets(GaussianPrediction): Get number of parameter combinations (different parameter vectors) using which predictions were computed.

• num_evalpoints(GaussianPrediction): Get number of points where predictions were computed.

Slots

f_comp_mean

component means

f_comp_std

component standard deviations

f_mean

signal mean (on normalized scale)

f_std

signal standard deviation (on normalized scale)

y_mean

predictive mean (on original data scale)

y_std

predictive standard deviation (on original data scale)

x

a data frame of points (covariate values) where the function posteriors or predictive distributions have been evaluated

See Also

Prediction

Description

Uses extract with permuted = FALSE and inc_warmup = FALSE.

Usage

get_draws(object, draws = NULL, reduce = NULL, ...)

Arguments

Value

The return value is always a 2-dimensional array of shape num_param_sets x num_params.

See Also

Other main functions: create_model(), draw_pred(), lgp(), pred(), prior_pred(), sample_model()

Description

NOTE: It is not recommended for users to call this. Use pred instead.

Usage

get_pred(fit, draws = NULL, reduce = NULL, verbose = TRUE)

Arguments

Value

an object of class GaussianPrediction or Prediction

Description

These have STAN_kernel_* counterparts. These R versions are provided for reference and are not optimized for speed. These are used when generating simulated data, and not during model inference.

Usage

kernel_eq(x1, x2, alpha = 1, ell)

kernel_ns(x1, x2, alpha = 1, ell, a)

kernel_zerosum(x1, x2, M)

kernel_bin(x1, x2, pos_class = 0)

kernel_cat(x1, x2)

kernel_varmask(x1, x2, a, vm_params)

kernel_beta(beta, idx1_expand, idx2_expand)

Arguments

Value

A matrix of size $n$ x $m$.

Functions

• kernel_eq(): Uses the exponentiated quadratic kernel.

• kernel_ns(): Uses the non-stationary kernel (input warping + squared exponential).

• kernel_zerosum(): Uses the zero-sum kernel. Here, x1 and x2 must be integer vectors (integers denoting different categories). Returns a binary matrix.

• kernel_bin(): Uses the binary (mask) kernel. Here, x1 and x2 must be integer vectors (integers denoting different categories). Returns a binary matrix.

• kernel_cat(): Uses the categorical kernel. Here, x1 and x2 must be integer vectors (integers denoting different categories). Returns a binary matrix.

• kernel_varmask(): Computes variance mask multiplier matrix. NaN's in x1 and x2 will be replaced by 0.

• kernel_beta(): Computes the heterogeneity multiplier matrix. NOTE: idx_expand needs to be given so that idx_expand[j]-1 tells the index of the beta parameter that should be used for the $j$th observation. If observation $j$ doesn't correspond to any beta parameter, then idx_expand[j] should be 1.

Description

An S4 class to represent input for kernel matrix computations

Usage

## S4 method for signature 'KernelComputer'
show(object)

## S4 method for signature 'KernelComputer'
num_components(object)

## S4 method for signature 'KernelComputer'
num_evalpoints(object)

## S4 method for signature 'KernelComputer'
num_paramsets(object)

## S4 method for signature 'KernelComputer'
component_names(object)

Arguments

Methods (by generic)

• show(KernelComputer): Print a summary about the object.

• num_components(KernelComputer): Get number of components.

• num_evalpoints(KernelComputer): Get number of evaluation points.

• num_paramsets(KernelComputer): Get number of parameter sets.

• component_names(KernelComputer): Get component names.

Slots

input

Common input (for example parameter values).

K_input

Input for computing kernel matrices between data points (N x N). A list.

Ks_input

Input for computing kernel matrices between data and output points (P x N). A list.

Kss_input

Input for computing kernel matrices between output points (P x P). A list, empty if full_covariance=FALSE.

comp_names

Component names (character vector).

full_covariance

Boolean value determining if this can compute full predictive covariance matrices (or just marginal variance at each point).

no_separate_output_points

Boolean value determining if Ks_input and Kss_input are the same thing. Using this knowledge can reduce unnecessary computations of kernel matrices.

STREAM

external pointer (for calling 'Stan' functions)

Description

Creates an additive Gaussian process model using create_model and fits it using sample_model. See the Mathematical description of lgpr models vignette for more information about the connection between different options and the created statistical model.

Usage

lgp(
  formula,
  data,
  likelihood = "gaussian",
  prior = NULL,
  c_hat = NULL,
  num_trials = NULL,
  options = NULL,
  prior_only = FALSE,
  verbose = FALSE,
  sample_f = !(likelihood == "gaussian"),
  quiet = FALSE,
  skip_postproc = sample_f,
  ...
)

Arguments

Value

Returns an object of the S4 class lgpfit.

Model formula syntax

There are two ways to define the model formula:

1. Using a common formula-like syntax, like in y ~ age + age|id + sex. Terms can consist of a single variable, such as age, or an interaction of two variables, such as age|id. In single-variable terms, the variable can be either continuous (numeric) or categorical (factor), whereas in interaction terms the variable on the left-hand side of the vertical bar (|) has to be continuous and the one on the right-hand side has to be categorical. Formulae specified using this syntax are translated to the advanced format so that

   • single-variable terms become gp(x) if variable x is numeric and zs(x) if x is a factor

   • interaction terms x|z become gp(x)*zs(z)

2. Using the advanced syntax, like in y ~ gp(age) + gp(age)*zs(id) + het(id)*gp_vm(disAge). This creates lgprhs objects, which consist of lgpterms, which consist of lgpexprs. This approach must be used if creating nonstationary, heterogeneous or temporally uncertain components.

Either one of the approaches should be used and they should not be mixed.

Defining priors

The prior argument must be a named list, like list(alpha=student_t(4), wrp=igam(30,10)). See examples in tutorials. Possible allowed names are

• "alpha" = component magnitude parameters

• "ell" = component lengthscale parameters

• "wrp" = input warping steepness parameters

• "sigma" = noise magnitude (Gaussian obs. model)

• "phi" = inv. overdispersion (negative binomial obs. model)

• "gamma" = overdispersion (beta-binomial obs. model)

• "beta" = heterogeneity parameters

• "effect_time" = uncertain effect time parameters

• "effect_time_info" = additional options for the above

See priors for functions that can be used to define the list elements. If a parameter of a model is not given in this list, a default prior will be used for it.

When to not use default priors

It is not recommended to use default priors blindly. Rather, priors should be specified according to the knowledge about the problem at hand, as in any Bayesian analysis. In lgpr this is especially important when

1. Using a non-Gaussian likelihood or otherwise setting sample_f = TRUE. In this case the response variable is not normalized, so the scale on which the data varies must be taken into account when defining priors of the signal magnitude parameters alpha and possible noise parameters (sigma, phi, gamma). Also it should be checked if c_hat is set in a sensible way.

2. Using a model that contains a gp_ns(x) or gp_vm(x) expression in its formula. In this case the corresponding covariate x is not normalized, and the prior for the input warping steepness parameter wrp must be set according to the expected width of the window in which the nonstationary effect of x occurs. By default, the width of this window is about 36, which has been set assuming that the unit of x is months.

See Also

Other main functions: create_model(), draw_pred(), get_draws(), pred(), prior_pred(), sample_model()

Description

An S4 class to represent an lgp expression

Slots

covariate

name of a covariate

fun

function name

See Also

See operations for performing arithmetics on lgprhs, lgpterm and lgpexpr objects.

Description

An S4 class to represent the output of the lgp function

Usage

## S4 method for signature 'lgpfit'
show(object)

## S4 method for signature 'lgpfit'
component_names(object)

## S4 method for signature 'lgpfit'
num_components(object)

## S4 method for signature 'lgpfit'
postproc(object, verbose = TRUE)

## S4 method for signature 'lgpfit'
contains_postproc(object)

## S4 method for signature 'lgpfit'
clear_postproc(object)

## S4 method for signature 'lgpfit'
get_model(object)

## S4 method for signature 'lgpfit'
get_stanfit(object)

## S4 method for signature 'lgpfit'
is_f_sampled(object)

## S4 method for signature 'lgpfit,missing'
plot(x, y)

Arguments

Methods (by generic)

• show(lgpfit): Print information and summary about the fit object.

• component_names(lgpfit): Get names of model components.

• num_components(lgpfit): Get number of model components. Returns a positive integer.

• postproc(lgpfit): Apply postprocessing. Returns an updated lgpfit object (copies data).

• contains_postproc(lgpfit): Check if object contains postprocessing information.

• clear_postproc(lgpfit): Returns an updated (copies data) lgpfit object without any postprocessing information.

• get_model(lgpfit): Get the stored lgpmodel object. Various properties of the returned object can be accessed as explained in the documentation of lgpmodel.

• get_stanfit(lgpfit): Get the stored stanfit object. Various properties of the returned object can be accessed or plotted as explained here or in the documentation of stanfit.

• is_f_sampled(lgpfit): Determine if inference was done by sampling the latent signal f (and its components).

• plot(x = lgpfit, y = missing): Visualize parameter draws using plot_draws.

Slots

stan_fit

An object of class stanfit.

model

An object of class lgpmodel.

num_draws

Total number of parameter draws.

postproc_results

A named list containing possible postprocessing results.

See Also

For extracting parameter draws, see get_draws, or the rstan methods for stanfit objects.

For more detailed plotting functions, see plot_draws, plot_beta, plot_warp, plot_effect_times

Description

An S4 class to represent an lgp formula

Slots

terms

an object of class lgprhs

y_name

name of the response variable

call

original formula call

See Also

See operations for performing arithmetics on lgprhs, lgpterm and lgpexpr objects.

Description

An S4 class to represent an additive GP model

Usage

## S4 method for signature 'lgpmodel'
show(object)

## S4 method for signature 'lgpmodel'
parameter_info(object, digits = 3)

## S4 method for signature 'lgpmodel'
component_info(object)

## S4 method for signature 'lgpmodel'
num_components(object)

## S4 method for signature 'lgpmodel'
covariate_info(object)

## S4 method for signature 'lgpmodel'
component_names(object)

## S4 method for signature 'lgpmodel'
is_f_sampled(object)

Arguments

Methods (by generic)

• show(lgpmodel): Print information and summary about the object. Returns object invisibly.

• parameter_info(lgpmodel): Get a parameter summary (bounds and priors). Returns a data.frame.

• component_info(lgpmodel): Get a data frame with information about each model component.

• num_components(lgpmodel): Get number of model components. Returns a positive integer.

• covariate_info(lgpmodel): Get covariate information.

• component_names(lgpmodel): Get names of model components.

• is_f_sampled(lgpmodel): Determine if inference of the model requires sampling the latent signal f (and its components).

Slots

formula

An object of class lgpformula

data

The original unmodified data.

stan_input

The data to be given as input to rstan::sampling

var_names

List of variable names grouped by type.

var_scalings

A named list with fields

• y - Response variable normalization function and its inverse operation. Must be an lgpscaling object.

• x_cont - Continuous covariate normalization functions and their inverse operations. Must be a named list with each element is an lgpscaling object.

var_info

A named list with fields

• x_cat_levels - Names of the levels of categorical covariates before converting from factor to numeric.

info

Other info in text format.

sample_f

Whether the signal f is sampled or marginalized.

full_prior

Complete prior information.

Description

An S4 class to represent the right-hand side of an lgp formula

Slots

summands

a list of one or more lgpterms

See Also

See operations for performing arithmetics on lgprhs, lgpterm and lgpexpr objects.

Description

An S4 class to represent variable scaling

Slots

loc

original location (mean)

scale

original scale (standard deviation)

var_name

variable name

Description

An S4 class to represent a data set simulated using the additive GP formalism

Usage

## S4 method for signature 'lgpsim'
show(object)

## S4 method for signature 'lgpsim,missing'
plot(x, y, ...)

Arguments

Methods (by generic)

• show(lgpsim): Show summary of object.

• plot(x = lgpsim, y = missing): Plot the data and generating process. For more information see plot_sim.

Slots

data

the actual data

response

name of the response variable in the data

components

the drawn function components

kernel_matrices

the covariance matrices for each gp

info

A list with fields

• par_ell the used lengthscale parameters

• par_cont the parameters used to generate the continuous covariates

• p_signal signal proportion

effect_times

A list with fields

• true possible true effect times that generate the disease effect

• observed possible observed effect times

Description

An S4 class to represent one formula term

Slots

factors

a list of at most two lgpexprs

See Also

See operations for performing arithmetics on lgprhs, lgpterm and lgpexpr objects.

Description

Print a model summary.

Usage

model_summary(object, digits = 3)

param_summary(object, digits = 3)

Arguments

Value

object invisibly.

Description

Replaces a continuous variable x in the data frame, and possibly another continuous variable x_ns derived from it, with new values, for each level of a grouping factor (usually id)

Usage

new_x(data, x_values, group_by = "id", x = "age", x_ns = NULL)

Arguments

Value

a data frame containing the following columns

• all factors in the original data

• x

• x_ns (unless it is NULL)

See Also

Other data frame handling functions: add_dis_age(), add_factor(), add_factor_crossing(), adjusted_c_hat(), split()

Description

Operations on formula terms and expressions

Usage

## S4 method for signature 'lgprhs,lgprhs'
e1 + e2

## S4 method for signature 'lgpterm,lgpterm'
e1 + e2

## S4 method for signature 'lgprhs,lgpterm'
e1 + e2

## S4 method for signature 'lgpterm,lgpterm'
e1 * e2

Arguments

Value

The behaviour and return type depend on the types of e1 and e2. You can

• sum (+) two lgprhs's to yield an lgprhs

• sum (+) two lgpterm's to yield an lgprhs

• sum (+) an lgprhs and an lgpterm to yield an lgprhs

• multiply (*) two lgpterm's to yield an lgpterm

Description

Data frames specified in arguments df, and df_err must have a format where

• The first column is the grouping factor (usually id).

• The second column is the x-axis variable (usually age).

• The third column is the coloring factor. If name of the third column is NA, coloring is not done.

• A column named y must contain the y-axis variable (not for df_err).

• A column named lower (upper) must contain the lower (upper) bound of error bar (only for df_err).

• The posterior draw using which the fit has been computed can be specified with a factor named _draw_ (only for df).

Usage

plot_api_c(
  df,
  df_err = NULL,
  alpha = 1,
  alpha_err = 0.2,
  no_err = FALSE,
  no_line = FALSE
)

Arguments

Value

A ggplot object.

See Also

Other internal plot API functions: plot_api_g()

Description

Data frames specified in arguments df_data, df_signal, df_fit, and df_fit_err must have a format where

• the first column is the grouping factor (usually id)

• the second column is the x-axis variable (usually age)

• a column named y must contain the y-axis variable (not for df_fit_err)

• a column named lower (upper) must contain the lower (upper) bound of error bar (only for df_fit_err)

• a column named draw must be a factor that specifies the posterior draw using which the fit has been computed (only for df_fit)

Usage

plot_api_g(
  df_data,
  df_signal = NULL,
  df = NULL,
  df_err = NULL,
  teff_signal = NULL,
  teff_obs = NULL,
  i_test = NULL,
  color_signal = color_palette(2)[1],
  color = color_palette(2)[2],
  color_err = colorset("red", "light_highlight"),
  color_vlines = colorset("gray", "mid_highlight"),
  alpha = 1,
  alpha_err = 0.5,
  nrow = NULL,
  ncol = NULL,
  y_transform = function(x) x
)

Arguments

Value

A ggplot object.

See Also

Other internal plot API functions: plot_api_c()

Description

This calls plot_f for all model components.

Usage

plot_components(
  fit,
  pred = NULL,
  group_by = "id",
  t_name = "age",
  MULT_STD = 2,
  verbose = TRUE,
  draws = NULL,
  reduce = function(x) base::mean(x),
  color_by = NA,
  no_err = FALSE,
  ylim = NULL,
  draw = TRUE,
  nrow = NULL,
  ncol = NULL,
  gg_add = NULL,
  x = NULL,
  ...
)

Arguments

Value

a list of ggplot objects invisibly

See Also

Other main plot functions: plot_draws(), plot_pred()

Description

Vizualizing longitudinal data

Usage

plot_data(
  data,
  x_name = "age",
  y_name = "y",
  group_by = "id",
  facet_by = NULL,
  color_by = NULL,
  highlight = NULL,
  main = NULL,
  sub = NULL
)

Arguments

Value

a ggplot object

Description

Visualize the distribution of parameter draws

Usage

plot_draws(
  fit,
  type = "intervals",
  regex_pars = c("alpha", "ell", "wrp", "sigma", "phi", "gamma"),
  ...
)

plot_beta(fit, type = "dens", verbose = TRUE, ...)

plot_warp(
  fit,
  num_points = 300,
  window_size = 48,
  color = colorset("red", "dark"),
  alpha = 0.5
)

plot_effect_times(fit, type = "areas", verbose = TRUE, ...)

Arguments

Value

a ggplot object or list of them

Functions

• plot_draws(): visualizes the distribution of any set of model parameters (defaults to kernel hyperparameters and possible observation model parameters)

• plot_beta(): visualizes the distribution of the individual-specific disease effect magnitude parameter draws

• plot_warp(): visualizes the input warping function for different draws of the warping steepness parameter

• plot_effect_times(): visualizes the input warping function for different parameter draws

See Also

Other main plot functions: plot_components(), plot_pred()

Description

Visualize input warping function with several steepness parameter values

Usage

plot_inputwarp(wrp, x, color = colorset("red", "dark"), alpha = 0.5)

Arguments

Value

a ggplot object

Description

Plot the inverse gamma-distribution pdf

Usage

plot_invgamma(
  alpha,
  beta,
  by = 0.01,
  log = FALSE,
  IQR = 0.95,
  return_quantiles = FALSE,
  linecolor = colorset("red", "dark"),
  fillcolor = colorset("red", "mid")
)

Arguments

Value

a ggplot object

See Also

Other functions related to the inverse-gamma distribution: dinvgamma_stanlike(), priors

Description

• Function draws at data points can be visualized using plot_pred. If the pred argument is NULL, it is computed using the pred function with x=NULL.

• The total signal f or any of its additive components can be plotted using plot_f.

Usage

plot_pred(
  fit,
  pred = NULL,
  group_by = "id",
  t_name = "age",
  MULT_STD = 2,
  verbose = TRUE,
  draws = NULL,
  reduce = function(x) base::mean(x),
  x = NULL,
  ...
)

plot_f(
  fit,
  pred = NULL,
  group_by = "id",
  t_name = "age",
  MULT_STD = 2,
  verbose = TRUE,
  draws = NULL,
  reduce = function(x) base::mean(x),
  comp_idx = NULL,
  color_by = NA,
  x = NULL,
  ...
)

Arguments

Value

a ggplot object

See Also

Other main plot functions: plot_components(), plot_draws()

Description

Visualize an lgpsim object (simulated data)

Usage

plot_sim(
  simdata,
  group_by = "id",
  x_name = "age",
  h_name = "h",
  y_name = "y",
  comp_idx = NULL,
  color_by = NA,
  verbose = TRUE,
  ...
)

Arguments

Value

a ggplot object

Description

Graphical posterior predictive checks

Usage

ppc(fit, data = NULL, fun = default_ppc_fun(fit), verbose = TRUE, ...)

Arguments

Value

a ggplot object

See Also

Introduction to graphical posterior predictive checks: here. Prior predictive check can be done by calling prior_pred and then bayesplot::pp_check().

Description

• If fit is for a model that marginalizes the latent signal f (i.e. is_f_sampled(fit) is FALSE), this computes the analytic conditional posterior distributions of each model component, their sum, and the conditional predictive distribution. All these are computed for each (hyper)parameter draw (defined by draws), or other parameter set (obtained by a reduction defined by reduce). Results are stored in a GaussianPrediction object which is then returned.

• If fit is for a model that samples the latent signal f (i.e. is_f_sampled(fit) is TRUE), this will extract these function samples, compute their sum, and a version of the sum f that is transformed through the inverse link function. If x is not NULL, the function draws are extrapolated to the points specified by x using kernel regression. Results are stored in a Prediction object which is then returned.

Usage

pred(
  fit,
  x = NULL,
  reduce = function(x) base::mean(x),
  draws = NULL,
  verbose = TRUE,
  STREAM = get_stream(),
  c_hat_pred = NULL,
  force = FALSE,
  debug_kc = FALSE
)

Arguments

Value

An object of class GaussianPrediction or Prediction.

See Also

Other main functions: create_model(), draw_pred(), get_draws(), lgp(), prior_pred(), sample_model()

Description

An S4 class to represent prior or posterior draws from an additive function distribution.

Usage

## S4 method for signature 'Prediction'
show(object)

## S4 method for signature 'Prediction'
component_names(object)

## S4 method for signature 'Prediction'
num_components(object)

## S4 method for signature 'Prediction'
num_paramsets(object)

## S4 method for signature 'Prediction'
num_evalpoints(object)

Arguments

Methods (by generic)

• show(Prediction): Print a summary about the object.

• component_names(Prediction): Get names of components.

• num_components(Prediction): Get number of components.

• num_paramsets(Prediction): Get number of parameter combinations (different parameter vectors) using which predictions were computed.

• num_evalpoints(Prediction): Get number of points where predictions were computed.

Slots

f_comp

component draws

f

signal draws

h

predictions (signal draws + scaling factor c_hat, transformed through inverse link function)

x

a data frame of points (covariate values) where the functions/predictions have been evaluated/sampled

extrapolated

Boolean value telling if the function draws are original MCMC draws or if they have been created by extrapolating such draws.

See Also

GaussianPrediction

Description

These functions take an lgpmodel object, and

• prior_pred samples from the prior predictive distribution of the model

• sample_param_prior samples only its parameter prior using sampling

Usage

prior_pred(
  model,
  verbose = TRUE,
  quiet = FALSE,
  refresh = 0,
  STREAM = get_stream(),
  ...
)

sample_param_prior(model, verbose = TRUE, quiet = FALSE, ...)

Arguments

Value

• prior_pred returns a list with components

  • y_draws: A matrix containing the prior predictive draws as rows. Can be passed to bayesplot::pp_check() for graphical prior predictive checking.

  • pred_draws: an object of class Prediction, containing prior draws of each model component and their sum

  • param_draws: a stanfit object of prior parameter draws (obtained by calling sample_param_prior internally)

• sample_param_prior returns an object of class stanfit

See Also

Other main functions: create_model(), draw_pred(), get_draws(), lgp(), pred(), sample_model()

Description

Convert given prior to numeric format

Usage

prior_to_num(desc)

Arguments

Value

a named list of parsed options

Description

These use the same parametrizations as defined in the 'Stan' documentation. See the docs for gamma and inverse gamma distributions.

Usage

uniform(square = FALSE)

normal(mu, sigma, square = FALSE)

student_t(nu, square = FALSE)

gam(shape, inv_scale, square = FALSE)

igam(shape, scale, square = FALSE)

log_normal(mu, sigma, square = FALSE)

bet(a, b)

Arguments

Value

a named list

See Also

Other functions related to the inverse-gamma distribution: dinvgamma_stanlike(), plot_invgamma()

Examples

# Log-normal prior
log_normal(mu = 1, sigma = 1)

# Cauchy prior
student_t(nu = 1)

# Exponential prior with rate = 0.1
gam(shape = 1, inv_scale = 0.1)

# Create a similar priors as in LonGP (Cheng et al., 2019)
# Not recommended, because a lengthscale close to 0 is possible.
a <- log(1) - log(0.1)
log_normal(mu = 0, sigma = a / 2) # for continuous lengthscale
student_t(nu = 4) # for interaction lengthscale
igam(shape = 0.5, scale = 0.005, square = TRUE) # for sigma

Description

Function for reading the built-in proteomics data

Usage

read_proteomics_data(parentDir = NULL, protein = NULL, verbose = TRUE)

Arguments

Value

a data.frame

Description

Assess component relevances

Usage

relevances(fit, reduce = function(x) base::mean(x), verbose = TRUE, ...)

Arguments

Value

a named vector with length equal to num_comps + 1

Description

S4 generics for lgpfit, lgpmodel, and other objects

Usage

parameter_info(object, digits)

component_info(object)

covariate_info(object)

component_names(object)

get_model(object)

is_f_sampled(object)

get_stanfit(object)

postproc(object, ...)

contains_postproc(object)

clear_postproc(object)

num_paramsets(object)

num_evalpoints(object)

num_components(object)

Arguments

Value

• parameter_info returns a data frame with one row for each parameter and columns for parameter name, parameter bounds, and the assigned prior

• component_info returns a data frame with one row for each model component, and columns encoding information about model components

• covariate_info returns a list with names continuous and categorical, with information about both continuous and categorical covariates

• component_names returns a character vector with component names

• get_model for lgpfit objects returns an lgpmodel

• is_f_sampled returns a logical value

• get_stanfit returns a stanfit (rstan)

• postproc applies postprocessing and returns an updated lgpfit

• clear_postproc removes postprocessing information and returns an updated lgpfit

• num_paramsets, num_evalpoints and num_components return an integer

Functions

• parameter_info(): Get parameter information (priors etc.).

• component_info(): Get component information.

• covariate_info(): Get covariate information.

• component_names(): Get component names.

• get_model(): Get lgpmodel object.

• is_f_sampled(): Determine if signal f is sampled or marginalized.

• get_stanfit(): Extract stanfit object.

• postproc(): Perform postprocessing.

• contains_postproc(): Determine if object contains postprocessing information.

• clear_postproc(): Clear postprocessing information (to reduce size of object).

• num_paramsets(): Get number of parameter sets.

• num_evalpoints(): Get number of points where posterior is evaluated.

• num_components(): Get number of model components.

See Also

To find out which methods have been implemented for which classes, see lgpfit, lgpmodel, Prediction and GaussianPrediction.

Description

• sample_model takes an lgpmodel object and fits it using sampling.

• optimize_model takes an lgpmodel object and fits it using optimizing.

Usage

sample_model(
  model,
  verbose = TRUE,
  quiet = FALSE,
  skip_postproc = is_f_sampled(model),
  ...
)

optimize_model(model, ...)

Arguments

Value

• sample_model returns an object of class lgpfit containing the parameter draws, the original model object, and possible postprocessing results. See documentation of lgpfit for more information.

• optimize_model directly returns the list returned by optimizing. See its documentation for more information.

See Also

Other main functions: create_model(), draw_pred(), get_draws(), lgp(), pred(), prior_pred()

Description

• select performs strict selection, returning either TRUE or FALSE for each component.

• select.integrate is like select, but instead of a fixed threshold, computes probabilistic selection by integrating over a threshold density.

• select_freq performs the selection separately using each parameter draw and returns the frequency at which each component was selected.

• select_freq.integrate is like select_freq, but instead of a fixed threshold, computes probabilistic selection frequencies by integrating over a threshold density.

Usage

select(fit, reduce = function(x) base::mean(x), threshold = 0.95, ...)

select_freq(fit, threshold = 0.95, ...)

select.integrate(
  fit,
  reduce = function(x) base::mean(x),
  p = function(x) stats::dbeta(x, 100, 5),
  h = 0.01,
  verbose = TRUE,
  ...
)

select_freq.integrate(
  fit,
  p = function(x) stats::dbeta(x, 100, 5),
  h = 0.01,
  verbose = TRUE,
  ...
)

Arguments

Value

See description.

Description

Printing formula object info using the show generic

Usage

## S4 method for signature 'lgpformula'
show(object)

## S4 method for signature 'lgprhs'
show(object)

## S4 method for signature 'lgpterm'
show(object)

Arguments

Value

the object invisibly

Description

Simulate latent function components for longitudinal data analysis

Usage

sim.create_f(
  X,
  covariates,
  relevances,
  lengthscales,
  X_affected,
  dis_fun,
  bin_kernel,
  steepness,
  vm_params,
  force_zeromean
)

Arguments

Value

a data frame FFF where one column corresponds to one additive component

Description

Create an input data frame X for simulated data

Usage

sim.create_x(
  N,
  covariates,
  names,
  n_categs,
  t_data,
  t_jitter,
  t_effect_range,
  continuous_info
)

Arguments

Value

a list

Description

Simulate noisy observations

Usage

sim.create_y(noise_type, f, snr, phi, gamma, N_trials)

Arguments

Value

A list out, where

• out$h is f mapped through an inverse link function (times N_trials if noise_type is binomial or beta-binomial)

• out$y is the noisy response variable.

Description

Compute all kernel matrices when simulating data

Usage

sim.kernels(
  X,
  types,
  lengthscales,
  X_affected,
  bin_kernel,
  useMaskedVarianceKernel,
  steepness,
  vm_params
)

Arguments

Value

a 3D array

Description

Generate an artificial longitudinal data set.

Usage

simulate_data(
  N,
  t_data,
  covariates = c(),
  names = NULL,
  relevances = c(1, 1, rep(1, length(covariates))),
  n_categs = rep(2, sum(covariates %in% c(2, 3))),
  t_jitter = 0,
  lengthscales = rep(12, 2 + sum(covariates %in% c(0, 1, 2))),
  f_var = 1,
  noise_type = "gaussian",
  snr = 3,
  phi = 1,
  gamma = 0.2,
  N_affected = round(N/2),
  t_effect_range = "auto",
  t_observed = "after_0",
  c_hat = 0,
  dis_fun = "gp_warp_vm",
  bin_kernel = FALSE,
  steepness = 0.5,
  vm_params = c(0.025, 1),
  continuous_info = list(mu = c(pi/8, pi, -0.5), lambda = c(pi/8, pi, 1)),
  N_trials = 1,
  force_zeromean = TRUE
)

Arguments

Value

An object of class lgpsim.

Examples

# Generate Gaussian data
dat <- simulate_data(N = 4, t_data = c(6, 12, 24, 36, 48), snr = 3)

# Generate negative binomially (NB) distributed count data
dat <- simulate_data(
  N = 6, t_data = seq(2, 10, by = 2), noise_type = "nb",
  phi = 2
)

Description

• split_by_factor splits according to given factor

• split_within_factor splits according to given data point indices within the same level of a factor

• split_within_factor_random selects k points from each level of a factor uniformly at random as test data

• split_random splits uniformly at random

• split_data splits according to given data rows

Usage

split_by_factor(data, test, var_name = "id")

split_within_factor(data, idx_test, var_name = "id")

split_within_factor_random(data, k_test = 1, var_name = "id")

split_random(data, p_test = 0.2, n_test = NULL)

split_data(data, i_test, sort_ids = TRUE)

Arguments

Value

a named list with names train, test, i_train and i_test

See Also

Other data frame handling functions: add_dis_age(), add_factor(), add_factor_crossing(), adjusted_c_hat(), new_x()

Description

A very small artificial test data, used mostly for unit tests

Usage

testdata_001

Format

A data frame with 24 rows and 6 variables:

id

individual id, a factor with levels: 1, 2, 3, 4

age

age

dis_age

disease-related age

blood

a continuous variable

sex

a factor with 2 levels: Male, Female

y

a continuous variable

See Also

Other built-in datasets: testdata_002

Description

Medium-size artificial test data, used mostly for tutorials

Usage

testdata_002

Format

A data frame with 96 rows and 6 variables:

id

individual id, a factor with levels: 01-12

age

age

diseaseAge

disease-related age

sex

a factor with 2 levels: Male, Female

group

a factor with 2 levels: Case, Control

y

a continuous variable

See Also

read_proteomics_data

Other built-in datasets: testdata_001

Description

Validate S4 class objects

Usage

validate_lgpexpr(object)

validate_lgpformula(object)

validate_lgpscaling(object)

validate_lgpfit(object)

validate_GaussianPrediction(object)

validate_Prediction(object)

Arguments

Value

TRUE if valid, otherwise reasons for invalidity

Description

Variance masking function

Usage

var_mask(x, stp)

Arguments

Value

a vector of length $n$

See Also

Other kernel utility functions: warp_input()

Description

Input warping function

Usage

warp_input(x, a)

Arguments

Value

a vector of warped inputs $w(x)$, length $n$

See Also

Other kernel utility functions: var_mask()