bliss: Bayesian functional Linear regression with Sparse Step functions

Description

A method for the Bayesian Functional Linear Regression model (functions-on-scalar), including two estimators of the coefficient function and an estimator of its support. A representation of the posterior distribution is also available.

Author(s)

Maintainer: Paul-Marie Grollemund paul_marie.grollemund@uca.fr

Other contributors:

• Isabelle Sanchez isabelle.sanchez@inra.fr [contractor]

• Meili Baragatti meili.baragatti@supagro.fr [contractor]

See Also

Useful links:

• https://github.com/pmgrollemund/bliss

• Report bugs at https://github.com/pmgrollemund/bliss/issues

between

Description

Check if a number belong to a given interval.

Usage

value %between% interval

Arguments

Value

a logical value.

Examples

1 %between% c(0,2)
2 %between% c(0,2)
3 %between% c(0,2)

BIC_model_choice

Description

Model selection with BIC criterion.

Usage

BIC_model_choice(Ks, iter, data, verbose = T)

Arguments

Value

A numerical vector, the BIC values for the Bliss model for different K value.

Examples

param_sim <- list(Q=1,n=100,p=c(50),grids_lim=list(c(0,1)))
data      <- sim(param_sim,verbose=TRUE)
iter = 1e2
Ks <- 1:5

res_BIC <- BIC_model_choice(Ks,iter,data)
plot(res_BIC,xlab="K",ylab="BIC")

Bliss_Gibbs_Sampler

Description

A Gibbs Sampler algorithm to sample the posterior distribution of the Bliss model.

Usage

Bliss_Gibbs_Sampler(data, param, verbose = FALSE)

Arguments

Value

a list containing :

trace

a matrix, the trace of the Gibbs Sampler.

param

a list containing parameters used to run the function.

Examples

param_sim <- list(Q=1,n=25,p=50,grids_lim=list(c(0,1)),iter=2e2,K=2)
data_sim <- sim(param_sim,verbose=FALSE)
res_Bliss_Gibbs_Sampler <- Bliss_Gibbs_Sampler(data_sim,param_sim)
theta_1 <- res_Bliss_Gibbs_Sampler$trace[1,]
theta_1

Bliss_Simulated_Annealing

Description

A Simulated Annealing algorithm to compute the Bliss estimate.

Usage

Bliss_Simulated_Annealing(
  beta_sample,
  posterior_sample,
  param,
  verbose_cpp = FALSE
)

Arguments

Value

a list containing:

Bliss_estimate

a numerical vector, corresponding to the Bliss estimate of the coefficient function.

Smooth_estimate

a numerical vector, which is the posterior expectation of the coefficient function for each time points.

trace

a matrix, the trace of the algorithm.

Examples

data(data1)
data(param1)
data(res_bliss1)
param1$Q <- length(data1$x)
param1$grids <- data1$grids
param1$p <- sapply(data1$grids,length)

posterior_sample <- res_bliss1$posterior_sample
beta_sample <- compute_beta_sample(posterior_sample,param1)

res_sann <- Bliss_Simulated_Annealing(beta_sample,posterior_sample,param1)

build_Fourier_basis

Description

Define a Fourier basis to simulate functional covariate observations.

Usage

build_Fourier_basis(grid, dim, per = 2 * pi)

Arguments

Details

See the sim_x function.

Value

a matrix. Each row is an functional observation evaluated on the grid time points.

Examples

# See the function \code{sim_x}.

change_grid

Description

Compute a function (evaluated on a grid) on a given (finer) grid.

Usage

change_grid(fct, grid, new_grid)

Arguments

Value

a numerical vector, the approximation of the function on the new grid.

Examples

grid <- seq(0,1,l=1e1)
new_grid <- seq(0,1,l=1e2)
fct <- 3*grid^2 + sin(grid*2*pi)
plot(grid,fct,type="o",lwd=2,cex=1.5)
lines(new_grid,change_grid(fct,grid,new_grid),type="o",col="red",cex=0.8)

choose_beta

Description

Compute a coefficient function for the Function Linear Regression model.

Usage

choose_beta(param)

Arguments

Details

Several shapes are available.

Value

A numerical vector which corresponds to the coefficient function at given times points (grid).

Examples

### smooth
param <- list(p=100,grid=seq(0,1,length=100),shape="smooth")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="l")
### random_smooth
param <- list(p=100,grid=seq(0,1,length=100),shape="random_smooth")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="l")
### simple
param <- list(p=100,grid=seq(0,1,length=100),shape="simple")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="s")
### simple_bis
param <- list(p=100,grid=seq(0,1,length=100),shape="simple_bis")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="s")
### random_simple
param <- list(p=100,grid=seq(0,1,length=100),shape="random_simple")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="s")
### sinusoid
param <- list(p=100,grid=seq(0,1,length=100),shape="sinusoid")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="l")
### flat_sinusoid
param <- list(p=100,grid=seq(0,1,length=100),shape="flat_sinusoid")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="l")
### sharp
param <- list(p=100,grid=seq(0,1,length=100),shape="sharp")
beta_function <- choose_beta(param)
plot(param$grid,beta_function,type="l")

compute_beta_posterior_density

Description

Compute the posterior density of the coefficient function.

Usage

compute_beta_posterior_density(beta_sample, param)

Arguments

Details

The posterior densities correponds to approximations of the marginal posterior distribitions (of beta(t) for each t). The sample is thinned in order to reduce the correlation and the computational time of the function kde2d.

Value

An approximation of the posterior density on a two-dimensional grid (corresponds to the result of the kde2d function).

Examples

library(RColorBrewer)
data(data1)
data(param1)
data(res_bliss1)
param1$grids <- data1$grids
param1$p <- sapply(data1$grids,length)
param1$Q <- length(data1$x)

density_estimate <- compute_beta_posterior_density(res_bliss1$beta_sample,param1)

compute_beta_sample

Description

Compute the posterior coefficient function from the posterior sample.

Usage

compute_beta_sample(posterior_sample, param)

Arguments

Value

a matrix containing the coefficient function posterior sample.

Examples

data(data1)
data(param1)
data(res_bliss1)
param1$grids <- data1$grids
param1$p <- sapply(data1$grids,length)
param1$Q <- length(data1$x)
beta_sample <- compute_beta_sample(posterior_sample=res_bliss1$posterior_sample,
                                   param=param1)

compute_chains_info

Description

Compute summaries of Gibbs Sampler chains.

Usage

compute_chains_info(chain, param)

Arguments

Value

Return a list containing the estimates of mu and sigma_sq, the Smooth estimate and the chain autocorrelation for mu, sigma_sq and beta.

Examples

a=1

compute_random_walk

Description

Compute a (Gaussian) random walk.

Usage

compute_random_walk(n, p, mu, sigma, start = rep(0, n))

Arguments

Details

See the sim_x function.

Value

a matrix where each row is a random walk.

Examples

# see the sim_x() function.

compute_starting_point_sann

Description

Compute a starting point for the Simulated Annealing algorithm.

Usage

compute_starting_point_sann(beta_expe)

Arguments

Value

a matrix with 3 columns : "m", "l" and "b". The two first columns define the begin and the end of the intervals and the third gives the mean values of each interval.

Examples

data(res_bliss1)
mystart<-compute_starting_point_sann(apply(res_bliss1$beta_sample[[1]],2,mean))

corr_matrix

Description

Compute an autocorrelation matrix.

Usage

corr_matrix(diagonal, ksi)

Arguments

Value

a symmetric matrix.

Examples

### Test 1 : weak autocorrelation
ksi     <- 1
diagVar <- abs(rnorm(100,50,5))
Sigma   <- corr_matrix(diagVar,ksi^2)
persp(Sigma)
### Test 2 : strong autocorrelation
ksi     <- 0.2
diagVar <- abs(rnorm(100,50,5))
Sigma   <- corr_matrix(diagVar,ksi^2)
persp(Sigma)

a list of data

Description

A data object for bliss model

Usage

data1

Format

a list of data

y

y coordinate

x

x coordinate

betas

the coefficient function used to generate the data

grids

the grid of the observation times

determine_intervals

Description

Determine for which intervals a function is nonnull.

Usage

determine_intervals(beta_fct)

Arguments

Value

a matrix with 3 columns : "begin", "end" and "value". The two first columns define the begin and the end of the intervals and the third gives the mean values of each interval.

Examples

data(data1)
data(param1)
# result of res_bliss1<-fit_Bliss(data=data1,param=param1)
data(res_bliss1)
intervals <- determine_intervals(res_bliss1$Bliss_estimate[[1]])
plot(data1$grids[[1]],res_bliss1$Bliss_estimate[[1]],type="s")
for(k in 1:nrow(intervals)){
   segments(data1$grids[[1]][intervals[k,1]],intervals[k,3],
           data1$grids[[1]][intervals[k,2]],intervals[k,3],col=2,lwd=4)
}

do_need_to_reduce

Description

Determine if it is required to reduce the size of the grid time points for each functional covariate.

Usage

do_need_to_reduce(param)

Arguments

Value

a boolean value.

Examples

data(param1)
param1$p <- sapply(data1$grids,length)

do_need_to_reduce(param1)

dposterior

Description

Compute (non-normalized) posterior densities for a given parameter set.

Usage

dposterior(posterior_sample, data, theta = NULL)

Arguments

Details

If the theta is NULL, the posterior density is computed from the MCMC sample given in the posterior_sample.

Value

Return the (log) posterior density, the (log) likelihood and the (log) prior density for the given parameter set.

Examples

data(data1)
data(param1)
# result of res_bliss1<-fit_Bliss(data=data1,param=param1)
data(res_bliss1)
# Compute the posterior density of the MCMC sample :
res_poste <- dposterior(res_bliss1$posterior_sample,data1)

fit_Bliss

Description

Fit the Bayesian Functional Linear Regression model (with Q functional covariates).

Usage

fit_Bliss(
  data,
  param,
  sann = TRUE,
  compute_density = TRUE,
  support_estimate = TRUE,
  sann_trace = FALSE,
  verbose = TRUE
)

Arguments

Value

return a list containing:

alpha

a list of Q numerical vector. Each vector is the function alpha(t) associated to a functional covariate. For each t, alpha(t) is the posterior probabilities of the event "the support covers t".

beta_posterior_density

a list of Q items. Each item contains a list containing information to plot the posterior density of the coefficient function with the image function.

grid_t

a numerical vector: the x-axis.

grid_beta_t

a numerical vector: the y-axis.

density

a matrix: the z values.

new_beta_sample

a matrix: beta sample used to compute the posterior densities.

beta_sample

a list of Q matrices. The qth matrix is a posterior sample of the qth functional covariates.

Bliss_estimate

a list of numerical vectors corresponding to the Bliss estimates of each functional covariates.

data

a list containing the data.

posterior_sample

a list of information about the posterior sample: the trace matrix of the Gibbs sampler, a list of Gibbs sampler parameters and the posterior densities.

support_estimate

a list of support estimates of each functional covariate.

support_estimate_fct

another version of the support estimates.

trace_sann

a list of Q matrices which are the trace of the Simulated Annealing algorithm.

Examples

# see the vignette BlissIntro.

image_Bliss

Description

Plot an approximation of the posterior density.

Usage

image_Bliss(beta_posterior_density, param = list(), q = 1, to_print = TRUE)

Arguments

Examples

data(data1)
data(param1)
data(res_bliss1)

image_Bliss(res_bliss1$beta_posterior_density,param1,q=1)

integrate_trapeze

Description

Trapezoidal rule to approximate an integral.

Usage

integrate_trapeze(x, y)

Arguments

Value

a numerical value, the approximation.

Examples

x <- seq(0,1,le=1e2)
integrate_trapeze(x,x^2)

integrate_trapeze(data1$grids[[1]],t(data1$x[[1]]))

interpretation_plot

Description

Provide a graphical representation of the functional data with a focus on the detected periods with the Bliss method.

Usage

interpretation_plot(data, Bliss_estimate, q = 1, centered = FALSE, cols = NULL)

Arguments

Examples

data(data1)
data(param1)
# result of res_bliss1 <- fit_Bliss(data=data1,param=param1,verbose=TRUE)
data(res_bliss1)
interpretation_plot(data=data1,Bliss_estimate=res_bliss1$Bliss_estimate,q=1)
interpretation_plot(data=data1,Bliss_estimate=res_bliss1$Bliss_estimate,q=1,centered=TRUE)

lines_bliss

Description

Add a line to a plot obtained with image_Bliss.

Usage

lines_bliss(x, y, col = "black", lty = "solid")

Arguments

Examples

data(data1)
data(param1)
data(res_bliss1)

image_Bliss(res_bliss1$beta_posterior_density,param1,q=1) +
lines_bliss(res_bliss1$data$grids[[1]],res_bliss1$smooth_estimate[[1]])+
lines_bliss(res_bliss1$data$grids[[1]],res_bliss1$Bliss_estimate[[1]],col="purple")

A list of param for bliss model

Description

A list of param for bliss model

Usage

param1

Format

a list of param for bliss model

Q

the number of functional covariates

n

the sample size

p

the number of observation times

beta_shapes

the shapes of the coefficient functions

grids_lim

the range of the observation times

grids

the grids of the observation times

K

the number of intervals for the coefficient function

pdexp

Description

Probability function of a discretized Exponentiel distribution.

Usage

pdexp(a, l_values)

Arguments

Value

a numerical vector, which is the prability function on l_values.

Examples

pdexp(10,seq(0,1,1))

x <- seq(0,10,le=1e3)
plot(x,dexp(x,0.5),lty=2,type="l")
lines(pdexp(0.5,1:10),type="p")

post_treatment_bliss

Description

Compute the post treatment values.

Usage

post_treatment_bliss(posterior_sample, param, data)

Arguments

Value

A list of important post treatment value: BIC, the maximum of the log likelihood and the numbre of parameters.

Examples

data(data1)
data(param1)
data(res_bliss1)

post_treatment_bliss(res_bliss1$posterior_sample,param1,data1)

predict_bliss

Description

Compute predictions.

Usage

predict_bliss(x, grids, burnin, posterior_sample, Smooth_estimate)

Arguments

Value

A vector of predictions for each individual data x.

Examples

data(data1)
data(param1)
data(res_bliss1)

predict_bliss(data1$x,data1$grids,50,res_bliss1$posterior_sample,res_bliss1$smooth_estimate)

predict_bliss_distribution

Description

Compute the distribution of the predictions.

Usage

predict_bliss_distribution(x, grids, burnin, posterior_sample, beta_sample)

Arguments

Value

A matrix containing predictions for each individual data x.

Examples

data(data1)
data(param1)
data(res_bliss1)

predict_bliss_distribution(data1$x,data1$grids,50,res_bliss1$posterior_sample,
   res_bliss1$beta_sample)

Print a bliss Object

Description

Print a bliss Object

Usage

printbliss(x, ...)

Arguments

Examples

# See fit_Bliss() function

reduce_x

Description

Reduce the number of time points.

Usage

reduce_x(data, param)

Arguments

Value

a numerical value, the approximation.

Examples

param <- list(Q=1,n=10,p=c(150),grids_lim=list(c(0,1)))
data <- sim(param)

data(param1)
param1$n <- nrow(data$x[[1]])
param1$p <- sapply(data$grids,length)
param1$Q <- length(data$x)

data <- reduce_x(data,param1)

A result of the BliSS method

Description

A result of the BliSS method

Usage

res_bliss1

Format

a Bliss object (list)

alpha

a list of Q numerical vector. Each vector is the function alpha(t) associated to a functional covariate. For each t, alpha(t) is the posterior probabilities of the event "the support covers t".

beta_posterior_density

a list of Q items. Each item contains a list containing information to plot the posterior density of the coefficient function with the image function.

grid_t

a numerical vector: the x-axis.

grid_beta_t

a numerical vector: the y-axis.

density

a matrix: the z values.

new_beta_sample

a matrix: beta sample used to compute the posterior densities.

beta_sample

a list of Q matrices. The qth matrix is a posterior sample of the qth functional covariates.

Bliss_estimate

a list of numerical vectors corresponding to the Bliss estimates of each functional covariates.

data

see the description of the object data1.

posterior_sample

a list containing (for each chain) the result of the Bliss_Gibbs_Sampler function.

Smooth_estimate

a list containing the Smooth estimates of the coefficient functions.

support_estimate

a list containing the estimations of the support.

support_estimate_fct

a list containing the estimation of the support.

trace_sann

a list containing (for each chain) the trace of the Simulated Annealing algorithm.

sigmoid

Description

Compute a sigmoid function.

Usage

sigmoid(x, asym = 1, v = 1)

Arguments

Details

see the function sim_x.

Value

a numerical vector.

Examples

## Test 1 :
x <- seq(-7,7,0.1)
y <- sigmoid(x)
plot(x,y,type="l",main="Sigmoid function")
## Test 2 :
x  <- seq(-7,7,0.1)
y  <- sigmoid(x)
y2 <- sigmoid(x,asym=0.5)
y3 <- sigmoid(x,v   =  5)
plot(x,y,type="l",main="Other sigmoid functions")
lines(x,y2,col=2)
lines(x,y3,col=3)

sigmoid_sharp

Description

Compute a sharp sigmoid function.

Usage

sigmoid_sharp(x, loc = 0, ...)

Arguments

Details

see the function sim_x.

Value

a numerical vector.

Examples

## Test 1 :
x <- seq(-7,7,0.1)
y <- sigmoid_sharp(x)
plot(x,y,type="l",main="Sharp sigmoid")
## Test 2 :
x  <- seq(-7,7,0.1)
y  <- sigmoid_sharp(x,loc=3)
y2 <- sigmoid_sharp(x,loc=3,asym=0.5)
y3 <- sigmoid_sharp(x,loc=3,v   =  5)
plot(x,y,type="l",main="Other sharp sigmoids")
lines(x,y2,col=2)
lines(x,y3,col=3)

sim

Description

Simulate a dataset for the Function Linear Regression model.

Usage

sim(param, verbose = FALSE)

Arguments

Value

a list containing:

Q

an integer, the number of functional covariates.

y

a numerical vector, the outcome observations.

x

a list of matrices, the qth matrix contains the observations of the qth functional covariate at time points given by grids.

grids

a list of numerical vectors, the qth vector is the grid of time points for the qth functional covariate.

betas

a list of numerical vectors, the qth vector is the 'true' coefficient function associated to the qth covariate on a grid of time points given with grids.

Examples

library(RColorBrewer)
param <- list(Q=2,n=25,p=c(50,50),grids_lim=list(c(0,1),c(-1,2)))
data <- sim(param)
data$y
cols <- colorRampPalette(brewer.pal(9,"YlOrRd"))(10)
q=2
matplot(data$grids[[q]],t(data$x[[q]]),type="l",lty=1,col=cols)
plot(data$grids[[q]],data$betas[[q]],type="l")
abline(h=0,lty=2,col="gray")

sim_x

Description

Simulate functional covariate observations.

Usage

sim_x(param)

Arguments

Details

Several shape are available for the observations: "Fourier", "Fourier2", "random_walk", "random_sharp", "uniform", "gaussian", "mvgauss", "mvgauss_different_scale", "mvgauss_different_scale2", "mvgauss_different_scale3" and "mvgauss_different_scale4".

Value

a matrix which contains the functional covariate observations at time points given by grid.

Examples

library(RColorBrewer)
### uniform
param <- list(n=15,p=100,grid=seq(0,1,length=100),x_type="uniform")
x <- sim_x(param)
cols <- colorRampPalette(brewer.pal(9,"YlOrRd"))(15)
matplot(param$grid,t(x),type="l",lty=1,col=cols)

support_estimation

Description

Compute the support estimate.

Usage

support_estimation(beta_sample, param)

Arguments

Value

a list containing:

alpha

a numerical vector. The approximated posterior probabilities that the coefficient function support covers t for each time points t.

estimate

a numerical vector, the support estimate.

estimate_fct

a numerical vector, another version of the support estimate.

Examples

data(data1)
data(param1)
data(res_bliss1)
param1$Q <- length(data1$x)

res_support <- support_estimation(res_bliss1$beta_sample,param1)