BayesSummaryStatLM: An R package for Bayesian Linear Models for Big Data and Data Science

Recent developments in data science and big data research have produced an abundance of large data sets that are too big to be analyzed in their entirety, due to limits on either computer memory or storage capacity. Here, we introduce our R package 'BayesSummaryStatLM' for Bayesian linear regression models with Markov chain Monte Carlo implementation that overcomes these limitations. Our Bayesian models use only summary statistics of data as input; these summary statistics can be calculated from subsets of big data and combined over subsets. Thus, complete data sets do not need to be read into memory in full, which removes any physical memory limitations of a user. Our package incorporates the R package 'ff' and its functions for reading in big data sets in chunks while simultaneously calculating summary statistics. We describe our Bayesian linear regression models, including several choices of prior distributions for unknown model parameters, and illustrate capabilities and features of our R package using both simulated and real data sets.


Introduction
Due to the developments in data science and the exponential growth of big data in recent years, statisticians are facing new challenges in the analysis of large data sets. Here, big data is defined as data sets that are too large to be analyzed in their entirety, due to limits on either computer memory or storage size. The application areas affected by big data are extensive, and include genomics, sustainability, public health, finance, energy and climatology, among many others. Here, we introduce Markov chain Monte Carlo (MCMC) procedures for Bayesian linear regression models with normally distributed errors that use only summary statistics as input. The summary statistics can be combined over subsets of data, which allows for analyzing large data sets in chunks. Since complete data sets do not need to be read into R at once, this overcomes any physical memory limits of a user. Our methods are implemented through the R [1] software package BayesSummaryStatLM (available from the Comprehensive R Archive Network at http://CRAN.R-project.org/), which integrates functions from the R package ff [2]. The ff functions keep only a portion of a data file in resident memory in R, while simultaneously calculating summary statistics on each portion. The summary statistics are then summed over portions to create full data summary statistics.
Similar methods that use summary statistics for Bayesian linear regression models were developed by Ordonez et al. [3] in a big data setting. In this work, the authors introduce a database management system (DBMS) procedure for obtaining summary statistics for full data sets; they then perform Bayesian linear regression on the summary statistics. In addition, Ghosh and Reiter [4] develop methods for Bayesian linear regression based on summary statistics from data subsets, but not in a big data context; instead, the application areas are privacy and secure sharing of information. The methods of Ordonez et al. [3] and Ghosh and Reiter [4] include only one choice of prior distribution for the unknown model parameters, and do not provide software tools to implement their methods. In contrast, we provide an R package for Bayesian linear regression models for big data that includes several choices of prior distributions for the unknown model parameters. In addition, our package can accommodate data that is located in separate files or that is distributed via networks. It can also analyze data that is updated over time; for example, streamed data. We describe each of these user-defined options of our package and provide examples in the following sections.
Our paper is organized as follows. In the Methods section, we introduce the Bayesian linear regression models and summary statistics. In the following section, we describe the main features and functions of our R package; we also illustrate our package using both simulated and real data sets. Our work is summarized in the Discussion section.

Bayesian linear regression model
Bayesian linear regression models are used extensively across a wide range of scientific subjects including economics, finance, and the biological, behavioral and social sciences, among many others; the purpose of linear regression is to model a response variable Y using a set of predictor variables 1 ( ,..., ) σ with known shape and scale parameters.
2) Inverse sigma squared prior for 2 σ , i.e. the Jeffreys prior for 2 σ (see [5]  first.rows The number of rows to read in the first chunk of data. Default = 100,000.

Description of function arguments
next.rows The number of rows to read in the remaining chunks of data. Default = 100,000.

update.summaries
The name of the R object containing previously-calculated summary statistics (if applicable), to be updated with new data. Default = NULL.

Description of returned values
The returned value for the read.regress.data.ff() function is a list containing the summary statistics named xtx, xty and yty (for , and ,  [7] and its function rmvn() for sampling from multivariate normal distributions. The mvnfast package generally speeds sampling compared to other packages due to the use of C++ code through the Rcpp [8,9] and RcppArmadillo [10] packages and the use of efficient numerical algorithms. The arguments used in a call to the function bayes.regress() are the following: bayes.regress(data.values = list(xtx, xty, yty, numsamp.data), beta.prior, sigmasq.prior, Tsamp.out,

zero.intercept)
The descriptions of the available function arguments are given here; the descriptions of the list arguments are given in Tables 1-4. data.values A list that contains input values xtx, xty, yty, numsamp.data; one optional argument and further details in Table 1.
beta.prior A list that contains the arguments for the prior for β ; details in Tables 2   and 3.
sigmasq.prior A list that contains the arguments for the prior for 2 σ ; details in Table 4.

Tsamp.out
The number of MCMC posterior samples to produce. Default = 1,000.
zero.intercept An optional logical parameter with default = FALSE. If zero.intercept = TRUE is specified, the linear regression model sets the y-intercept term 0 β to zero (i.e. removes the 0 β term from the model); by default, the y-intercept term 0 β is included in the model.

Description of returned values
The Further analysis, including plotting and creating posterior summary statistics, can be carried out using the CODA R package (Plummer et al. [11]); examples are given below.

Description of the simulation procedure
Here, we demonstrate our R package by simulating data from the linear regression model (1) with ten predictor variables, with data sample size 100,000,000, and file size 10.43 GB. Note that this file cannot be read into memory in R using standard functions such as read.table(), due to physical memory limits. The predictor data matrix 1 1 0 ( ,..., ) = X X X was simulated from a multivariate normal distribution by the following, where each column vector represents a predictor variable: The variance-covariance matrix Σ is defined as follows: 1, 1,...,10; , We use ρ = 0.2, so that the predictors h X , h = 1,…,10, have moderate correlation levels.
The model parameters β were simulated from a standard normal distribution; these are the parameters that are estimated in our analysis (see Table 5). The dependent values , i y 1,...,100, 000, 000, i = were then simulated from the following linear regression model: The data file was named "sim.multreg". Here, we assign 2 σ = 1.

Illustration of the R functions for simulation data
We read into the function read.regress.data.ff()the data file that contains the simulated X and Y values, named "sim.multreg". At the R prompt, the user enters the following command: The returned value sim.regress.data is a list containing the summary statistics xtx, xty, yty, and the total number of data values numsamp.data = 100,000,000.
These values are then used as input to the function bayes.regress() in the data.values argument, described below.
Next, we illustrate the function bayes.regress() for linear regression model (1) with k = 10. We first specify the priors for β and 2 σ as separate commands in R; these are read into the R function bayes.regress(). Here, we use prior 2) for β and prior 1) for 2 σ , with the default arguments. At the R prompt, the user enters the following commands: Further analysis of the MCMC posterior samples can be carried out using the R package CODA [11]; we illustrate this package in the following commands. Here, the output created above named sim.beta.sigmasq.out is first converted to class "mcmc" using the mcmc command. The history and density plots for the marginal posterior distribution for 1 β are generated using the plot command; results are shown in Figure 1.
The summary statistics for the marginal posterior distribution for 1 β are produced using the summary command, with results shown below. We discard the first 1,000 MCMC posterior samples for burnin. At the R prompt, the user enters the following commands: The returned value for the summary command in R is the following: The returned value for the 95% posterior equal-tail credible interval limits of (-0.  Table 5; we also show the history and density plots for the marginal posterior distribution of 2 σ in Figure 1.

Saving and updating summary statistics
Here, we provide an example of saving and updating summary statistics using the argument update.summaries in our read.regress.data.ff() function. We also present examples of the built-in functions of the R environment for saving and loading objects in the workspace; these features are used when the R session needs to be terminated and restarted at a later time.
We first repeat the command to read in the simulation data file named "sim.multreg" and create summary statistics: The R object named "sim.regress.data" can be saved onto a local directory using the built-in command in R named save.image(), as follows: > save.image(file="c:\simregressdata.RData") This R object can be loaded in a later R session using the built-in command load(), as follows: > load(file="c:\simregressdata.RData") When new data becomes available in a file named "sim2.multreg" (for example), the summary statistics included in the R object "sim.regress.data" are updated using the

Illustration of the R package for real data
Here, we read into the function read.regress.data.ff()the data file that contains the X and Y values, named "air.multreg". At the R prompt, the user enters the following command: commands to those shown above were used to produce posterior summary statistics, posterior history and density plots for β and 2 σ . We show plots for 1 β and 2 σ in Figure   2, and summary statistics for all model parameters in Table 5.

Software computational time
Here, we summarize our computational times for reading in data sets, calculating summary statistics and producing the MCMC posterior samples, using a computer with operating system Windows 7 64-bit and an Intel Core i7-4600U CPU 2.1 GHz Processor; complete results are provided in Table 6. For this, we created simulation data using the procedure similar to that described in Section 4.1 above, and removing the y-intercept

Discussion
We introduced and demonstrated the R package BayesSummaryStatLM for Bayesian linear regression models with MCMC implementation that uses only summary statistics as input. Since the summary statistics can be calculated for data subsets and combined, the complete data set does not need to be read into physical memory at once. As a result,

Prior distributions and full conditional posterior distributions for β
In this section, we introduce three choices of prior distributions for β that are included in the R package.

2) Multivariate Normal prior for β with known mean vector and known covariance matrix
Here, the prior distribution for β is specified by the following: is a known ( 1) 1 vector, specified by the user, is a known ( 1) ( 1) matrix, specified by the user; must be symmetric and positive definite.

3) Multivariate Normal prior for β with unknown mean vector and unknown covariance matrix
Here, the following prior distribution is specified for β : Normal is an unknown ( 1) 1 vector, is an unknown ( 1) ( 1) matrix; must be symmetric and positive definite.

1) Inverse Gamma prior for 2 σ with known shape and scale parameters
Here, the prior distribution for 2 σ is: The resulting full conditional posterior distribution for 2 σ is (see [5,6,13,14]): Note that the expression ( ) ( ) Thus, the full conditional posterior distribution for 2 σ can be restated as: 2) Inverse sigma squared prior for 2 σ Here, the following prior distribution is specified for 2 σ , which is the Jeffreys prior for 2 σ (see [5]): The corresponding full conditional posterior distribution for 2 σ is (see [5,6,13,14]):

The Gibbs sampler algorithm
Once the full conditional posterior distributions are determined, the Gibbs sampler proceeds as follows in the R package [5,6].      Tables   Table 1. Summary of available list arguments for the function argument data.values of the function bayes.regress()

List Argument Description xtx
The square matrix T X X of dimension ( 1) ( 1) where k = number of predictors; based on the summaries of the subset data, as shown in Formula (5) above. Default = 1.

xty
The square matrix T X Y of dimension ( 1) where k = number of predictors; based on the summaries of the subset data, as shown in Formula (6) above.
where k = number of predictors; based on the summaries of the subset data, as shown in Formula (7) above. Default = 1.

xtx.inv
Optional argument. This will be used for the Uniform prior distribution, to speed computations. If not provided, the argument xtx will be used to compute the inverse.

numsamp.data
The total number of data points in the full data set.

prec.Cinv
The fixed known prior precision matrix for the Multivariate Normal prior of β . This is automatically calculated by R as the inverse of C and does not need to be specified by the user. If provided, the preference will be given to this matrix. The default is an identity matrix with dimension equal to the length of β .

mu.init
The initial value for μ for the MCMC chain.
The default is a vector of 1's with length equal to the length of β .

Cinv.init
The initial value for 1 − C for the MCMC chain.
The default is an identity matrix with dimension equal to the length of β .  Table 6. Time (in seconds) to read in data and calculate summary statistics using read.regress.data.ff() function, and for producing 11,000 MCMC posterior samples (including burnin) using bayes.regress() function for specific choices of prior distributions. Note that β prior 3) contains three sets of unknown model parameters, including an unknown precision matrix, which adds considerably to the computational time for producing MCMC posterior samples. The results are based on a computer with operating system Windows 7 64-bit and an Intel Core i7-4600U CPU 2.1 GHz Processor.