---
abstract: 'Big data are data on a massive scale in terms of volume, intensity, and complexity that exceed the capacity of standard analytic tools. They present opportunities as well as challenges to statisticians. The role of computational statisticians in scientific discovery from big data analyses has been under-recognized even by peer statisticians. This article summarizes recent methodological and software developments in statistics that address the big data challenges. Methodologies are grouped into three classes: subsampling-based, divide and conquer, and online updating for stream data. As a new contribution, the online updating approach is extended to variable selection with commonly used criteria, and their performances are assessed in a simulation study with stream data. Software packages are summarized with focuses on the open source and packages, covering recent tools that help break the barriers of computer memory and computing power. Some of the tools are illustrated in a case study with a logistic regression for the chance of airline delay.'
author:
- 'Chun Wang, Ming-Hui Chen, Elizabeth Schifano, Jing Wu, and Jun Yan'
bibliography:
- 'softrevTechRep.bib'
title: Statistical Methods and Computing for Big Data
---

\#1

[Key words:]{} bootstrap; divide and conquer; external memory algorithm; high performance computing; online update; sampling; software;

Introduction {#sec:intro}
============

A 2011 McKinsey report predicted shortage of talent necessary for organizations to take advantage of big data [@Many:etal:big:2011]. Data now stream from daily life thanks to technological advances, and big data has indeed become a big deal [e.g., @Shaw:why:2014]. In the President’s Corner of the June 2013 issue of AMStat News, the three presidents (elect, current, and past) of the American Statistical Association (ASA) wrote an article titled “The ASA and Big Data” [@Sche:Davi:Rodr:asa:2013]. This article echos the June 2012 column of @Rodr:big:2012 on the recent media focus on big data, and discusses on what the statistics profession needs to do in response to the fact that statistics and statisticians are missing from big data discussions. In the followup July 2013 column, president Marie Davidian further raised the issues of statistics not being recognized as data science and mainstream academic statisticians being left behind by the rise of big data [@Davi:aren:2013]. A white paper prepared by a working group of the ASA called for more ambitious efforts from statisticians to work together with researchers in other fields on national research priorities in order to achieve better science more quickly [@Rudi:etal:disc:2014]. The same concern was expressed in a 2014 president’s address of the Institute of Mathematical Statistics (IMS) [@Yu:let:2014]. President Bin Yu of the IMS called for statisticians to own Data Science by working on real problems such as those from genomics, neuroscience, astronomy, nanoscience, computational social science, personalized medicine/healthcare, finance, and government; relevant methodology/theory will follow naturally.

Big data in the media or the business world may mean differently than what are familiar to academic statisticians [@Jord:Lin:stat:2014]. Big data are data on a massive scale in terms of volume, intensity, and complexity that exceed the ability of standard software tools to manage and analyze [e.g., @Snij:Matz:Reip:big:2012]. The origin of the term “big data” as it is understood today has been traced back in a recent study [@Dieb:pers:2012] to lunch-table conversations at Silicon Graphics in the mid-1990s, in which John Mashey figured prominently [@Mash:big:1998]. Big data are generated by countless online interactions among people, transactions between people and systems, and sensor-enabled machinery. Internet search engines (e.g., Google and YouTube) and social network tools (e.g., Facebook and Twitter) generate billions of activity data per day. Rather than Gigabytes and Terabytes, nowadays, the data produced are estimated by zettabytes, and are growing 40% every day [@Fan:Bife:mini:2013]. In the big data analytics world, a 3V definition by @Lane:3-D:2001 is widely accepted: volume (amount of data), velocity (speed of data in and out), and variety (range of data types and sources). High variety brings nontraditional or even unstructured data types, such as social network sentiments and internet map usage, which calls for new, creative ways to understand the structure of data and even to ask intelligent research questions [e.g., @Jord:Lin:stat:2014]. High volume and high velocity may bring noise accumulation, spurious correlation and incidental homogeneity, creating issues in computational feasibility and algorithmic stability [@Fan:Han:Liu:chal:2014].

Notwithstanding that new statistical thinking and methods are needed for the high variety aspect of big data, our focus is on fitting standard statistical models to big data whose size exceeds the capacity of a single computer from its high volume and high velocity. There are two computational barriers for big data analysis: 1) the data can be too big to hold in a computer’s memory; and 2) the computing task can take too long to wait for the results. These barriers can be approached either with newly developed statistical methodologies and/or computational methodologies. Despite the impression that statisticians are left behind in media discussions or governmental summits on big data, some statisticians have made important contributions and are pushing the frontier. Sound statistical procedures that are scalable computationally to massive datasets have been proposed [@Jord:on:2013]. Examples are subsampling-based approaches [@Klei:Talw:Sark:Jord:scal:2014; @Ma:Maho:Yu:stat:2013; @Lian:etal:resa:2013; @Macl:Adams:fire:2014], divide and conquer approaches [@Lin:Xi:aggr:2011; @Chen:Xie:spli:2014; @Song:Lian:spli:2014; @Neis:Wang:Xing:asym:2013], and online updating approaches [@Schifano2015]. From a computational perspective, much effort has been put into the most active, open source statistical environment, [@R]. Statistician developers are relentless in their drive to extend the reach of into big data [@Rick:stat:2013]. Recent UseR! conferences had many presentations that directly addressed big data, including a 2014 keynote lecture by John Chambers, the inventor of the language [@Cham:inte:2014]. Most cutting edge methods are first and easily implemented in . Given the open source nature of and the active recent development, our focus on software for big data will be on and packages. Revolution Enterprise () is a commercialized version of , but it offers free academic use, so it is also included in our case study and benchmarked. Other commercial software such as , , and will be briefly touched upon for completeness.

The rest of the article is organized as follows. Recent methodological developments in statistics on big data are summarized in Section \[sec:meth\]. Updating formulas for commonly used variable selection criteria in the online setting are developed and their performances studied in a simulation study in Section \[sec:varsel\]. Resources from open source software for analyzing big data with classical models are summarized in Section \[sec:R\]. Commercial software products are presented in Section \[sec:comm\]. A case study on a logistic model for the chance of airline delay is presented in Section \[sec:case\]. A discussion concludes in Section \[sec:disc\].

Statistical Methods {#sec:meth}
===================

The recent methodologies for big data can be loosely grouped into three categories: resampling-based, divide and conquer, and online updating. To put the different methods in a context, consider a dataset with $n$ independent and identically distributed observations, where $n$ is too big for standard statistical routines such as logistic regression.

Subsampling-Based Methods {#sect:stat:resamp}
-------------------------

### Bags of Little Bootstrap

@Klei:Talw:Sark:Jord:scal:2014 proposed the bags of little bootstrap (BLB) approach that provides both point estimates and quality measures such as variance or confidence intervals. It is a combination of subsampling [@Poli:Roma:Wolf:subs:1999], the $m$-out-of-$n$ bootstrap [@Bick:Gotz:van:resa:1997], and the bootstrap [@Efro:boot:1979] to achieve computational efficiency. BLB consists of the following steps. First, draw $s$ subsamples of size $m$ from the original data of size $n$. For each of the $s$ subsets, draw $r$ bootstrap samples of size $n$ instead of $m$, and obtain the point estimates and their quality measures (e.g., confidence interval) from the $r$ bootstrap sample. Then, the $s$ bootstrap point estimates and quality measures are combined (e.g., by average) to yield the overall point estimates and quality measures. In summary, BLB has two nested procedures: the inner procedure applies the bootstrap to a subsample, and the outer procedure combines these multiple bootstrap estimates. The subsample size $m$ was suggested to be $n^{\gamma}$ with $\gamma \in [0.5, 1]$ [@Klei:Talw:Sark:Jord:scal:2014], a much smaller number than $n$. Although the inner bootstrap procedure conceptually generates multiple resampled data of size $n$, what is really needed in the storage and computation is a sample of size $m$ with a weight vector. In contrast to subsampling and the $m$-out-of-$n$ bootstrap, there is no need for an analytic correction (e.g., $\sqrt{m/n}$) to rescale the confidence intervals from the final result. The BLB procedure facilitates distributed computing by letting each subsample of size $m$ be processed by a separate processor. @Klei:Talw:Sark:Jord:scal:2014 proved the consistency of BLB and provided high order correctness. Their simulation study showed good accuracy, convergence rate and the remarkable computational efficiency.

### Leveraging

@Ma:Sun:leve:2015 proposed to use leveraging to facilitate scientific discoveries from big data using limited computing resources. In a leveraging method, one samples a small proportion of the data with certain weights (subsample) from the full sample, and then performs intended computations for the full sample using the small subsample as a surrogate. The key to success of the leveraging methods is to construct the weights, the nonuniform sampling probabilities, so that influential data points are sampled with high probabilities [@Ma:Maho:Yu:stat:2013]. Leveraging methods are different from the traditional subsampling or $m$-out-of-$n$ bootstrap in that 1) they are used to achieve feasible computation even if the simple analytic results are available; 2) they enable visualization of the data when visualization of the full sample is impossible; and 3) they usually use unequal sampling probabilities for subsampling data. This approach is quite unique in allowing pervasive access to extract information from big data without resorting to high performance computing.

### Mean Log-likelihood

@Lian:etal:resa:2013 proposed a resampling-based stochastic approximation approach with an application to big geostatistical data. The method uses Monte Carlo averages calculated from subsamples to approximate the quantities needed for the full data. Motivated from minimizing the Kullback–Leibler (KL) divergence, they approximate the KL divergence by averages calculated from subsamples. This leads to a maximum mean log-likelihood estimation method. The solution to the mean score equation is obtained from a stochastic approximation procedure, where at each iteration, the current estimate is updated based on a subsample of size $m$ drawn from the full data. As $m$ is much smaller than $n$, the method is scalable to big data. @Lian:etal:resa:2013 established the consistency and asymptotic normality of the resulting estimator under mild conditions. In a simulation study, the convergence rate of the method was almost independent of $n$, the sample size of the full data.

### Subsampling-Based MCMC

As a popular general purpose tool for Bayesian inference, Markov chain Monte Carlo (MCMC) for big data is challenging because of the prohibitive cost of likelihood evaluation of every datum at every iteration. @Lian:Kim:boot:2013 extended the mean log-likelihood method to a bootstrap Metropolis–Hastings (MH) algorithm in MCMC. The likelihood ratio of the proposal and current estimate in the MH ratio is replaced with an approximation from the mean log-likelihood based on $k$ bootstrap samples of size $m$. The algorithm can be implemented exploiting the embarrassingly parallel structure and avoids repeated scans of the full dataset in iterations. @Macl:Adams:fire:2014 proposed an auxiliary variable MCMC algorithm called Firefly Monte Carlo (FlyMC) that only queries the likelihoods of a potentially small subset of the data at each iteration yet simulates from the exact posterior distribution. For each data point, a binary auxiliary variable and a strictly positive lower bound of the likelihood contribution are introduced. The binary variable for each datum effectively turn on and off data points in the posterior, hence the “firefly" name. The probability of turning on each datum depends on the ratio of its likelihood contribution and the introduced lower bound. The computational gain depends on that the lower bound is tight enough and that simulation of the auxiliary variables is cheap enough. Because of the need to hold the whole data in computer memory, the size of the data this method can handle is limited.

The pseudo-marginal Metropolis–Hasting algorithm replaces the intractable target (posterior) density in the MH algorithm with an unbiased estimator [@Andr:Robe:pseu:2009]. The the log-likelihood is estimated by an unbiased subsampled version, and an unbiased estimator of the likelihood is obtained by correcting the bias of the exponentiation of this estimator. @Quir:Vill:Kohn:spee:2014 proposed subsampling the data using probability proportional-to-size (PPS) sampling to obtain an approximately unbiased estimate of the likelihood which is used in the M-H acceptance step. The subsampling approach was further improved in @Quir:Vill:Kohn:scal:2015 using the efficient and robust difference estimator form the survey sampling literature.

Divide and Conquer {#sect:stat:dc}
------------------

A divide and conquer algorithm (which may appear under other names such as divide and recombine, split and conquer, or split and merge) generally has three steps: 1) partitions a big dataset into $K$ blocks; 2) processes each block separately (possibly in parallel); and 3) aggregates the solutions from each block to form a final solution to the full data.

### Aggregated Estimating Equations

For a linear regression model, the least squares estimator for the regression coefficient $\beta$ for the full data can be expressed as a weighted average of the least squares estimator for each block with weight being the inverse of the estimated variance matrix. The success of this method for linear regression depends on the linearity of the estimating equations in $\beta$ and that the estimating equation for the full data is a simple summation of that for all the blocks. For general nonlinear estimating equations, @Lin:Xi:aggr:2011 proposed a linear approximation of the estimating equations with the Taylor expansion at the solution in each block, and, hence, reduce the nonlinear estimating equation to the linear case so that the solutions to all the blocks are combined by a weighted average. The weight of each block is the slope matrix of the estimating function at the solution in that block, which is the Fisher information or inverse of the variance matrix if the equations are score equations. @Lin:Xi:aggr:2011 showed that, under certain technical conditions including $K = O(n^{\gamma})$ for some $\gamma \in (0, 1)$, the aggregated estimator has the same limit as the estimator from the full data.

### Majority Voting

@Chen:Xie:spli:2014 consider a divide and conquer approach for generalized linear models (GLM) where both the sample size $n$ and the number of covariates $p$ are large, by incorporating variable selection via penalized regression into a subset processing step. More specifically, for $p$ bounded or increasing to infinity slowly, ($p_n$ not faster than $o(e^{n_k})$, while model size may increase at a rate of $o(n_k)$), they propose to first randomly split the data of size $n$ into $K$ blocks (size $n_k=O(n/K)$). In step 2, penalized regression is applied to each block separately with a sparsity-inducing penalty function satisfying certain regularity conditions. This approach can lead to differential variable selection among the blocks, as different blocks of data may result in penalized estimates with different non-zero regression coefficients. Thus, in step 3, the results from the $K$ blocks are combined by majority vote to create a combined estimator. The implicit assumption is that real effects should be found persistently and therefore should be present even under perturbation by subsampling [e.g. @meinsh.buhl:2010]. The derivation of the combined estimator in step 3 stems from ideas for combining confidence distributions in meta-analysis [@singh.xie.straw:2005; @xie.singh.straw:2011], where one can think of the $K$ blocks as $K$ independent and separate analyses to be combined in a meta-analysis. The authors show under certain regularity conditions that their combined estimator in step 3 is model selection consistent, asymptotically equivalent to the penalized estimator that would result from using all of the data simultaneously, and achieves the oracle property when it is attainable for the penalized estimator from each block [see e.g., @fan.lv:2011]. They additionally establish an upper bound for the expected number of incorrectly selected variables and a lower bound for the expected number of correctly selected variables.

### Screening with Ultrahigh Dimension

Instead of dividing the data into blocks of observations in step 1, @Song:Lian:spli:2014 proposed a split-and-merge (SAM) method that divides the data into subsets of covariates for variable selection in ultrahigh dimensional regression from the Bayesian perspective. This method is particularly suited for big data where the number of covariates $P_n$ is much larger than the sample size $n$, $P_n \gg n$, and possibly increasing with $n$. In step 2, Bayesian variable selection is separately performed on each lower dimensional subset, which facilitates parallel processing. In step 3, the selected variables from each subset are aggregated, and Bayesian variable selection is applied on the aggregated data. The embarrassingly parallel structure in step 2 makes the SAM method applicable to big data problems with millions or more predictors. Posterior consistency is established for correctly specified models and for misspecified models, the latter of which is necessary because it is quite likely that some true predictors are missing. With correct model specification, true covariates will be identified as the sample size becomes large; under misspecified models, all predictors correlated with the response variable will be identified. Compared with the sure independence screening (SIS) approach [@Fan:Lv:sure:2008], the method uses the joint information of multiple predictors in predictor screening while SIS only uses the marginal information of each predictor. Their numerical results show that the SAM approach outperforms competing methods for ultrahigh dimensional regression.

### Parallel MCMC

In the Bayesian framework, it is natural to partition the data into $K$ subsets and run parallel MCMC on each one of them. The prior distribution for each subset is often obtained by taking a power $1/K$ of the prior distribution for whole data in order to preserve the total amount of prior information (which may change the impropriety of the prior). MCMC is run independently on each subset with no communications between subsets (and, thus, embarrassingly parallel), and the resulting samples are combined to approximate samples from the full data posterior distribution. @Neis:Wang:Xing:asym:2013 proposed to use kernel density estimators of the posterior density for each data subset, and estimate the full data posterior by multiplying the subset posterior densities together. This method is asymptotically exact in the sense of being converging in the number of MCMC iterations. @Wang:etal:para:2015 replaced the kernel estimator of @Neis:Wang:Xing:asym:2013 with a random partition tree histogram, which uses the same block partition across all terms in the product representation of the posterior to control the number of terms in the approximation such that it does not explode with $m$. @Scot:etal:baye:2013 proposed a consensus Monte Carlo algorithm, which produces the approximated full data posterior using weighted averages over the subset MCMC samples. The weight used (for Gaussian models) for each subset is the inverse of the variance-covariance matrix of the MCMC samples. The method is effective when the posterior is close to Gaussian but may cause bias when the distribution is skewed or has multi-modes. The consensus Monte Carlo principal is approached from a variational perspective by @Robi:Ange:Jord:vari:2015. The embarrassingly parallel feature of these methods facilitates their implementation in the MapReduce framework that exploits the division and recombination strategy [@Dean:Ghem:mapr:2008]. The final recombination step is implemented in package `parallelMCMCcombine` [@Miro:Conl:para:2014].

Going beyond embarrassingly parallel MCMC remains challenging because of storage issues and communication overheads. General strategies for parallel MCMC such as multiple-proposal MH algorithm [@Cald:gene:2014] and population MCMC [@Song:Wu:Lian:weak:2014] mostly require full data at each node.

Online Updating for Stream Data
-------------------------------

In some applications, data come in streams or large chunks, and a sequentially updated analysis is desirable without storing the data. Motivated from a Bayesian inference perspective, @Schifano2015 extends the work of @Lin:Xi:aggr:2011 in a few important ways. First, they introduce divide-and-conquer-type variance estimates of regression parameters in the linear model and estimating equation settings. These estimates of variability allow for users to make inferences about the true regression parameters based upon the previously developed divide-and-conquer point estimates of the regression parameters. Second, they develop iterative estimating algorithms and statistical inferences for linear models and estimating equations that update as new data arrive. Thus, while the divide-and-conquer setting is quite amenable to parallel processing for each subset, the online-updating approach for data streams is inherently sequential in nature. Their algorithms were designed to be computationally efficient and minimally storage-intensive, as they assume no access/storage of the historical data. Third, the authors address the issue of possible rank deficiencies when dealing with blocks of data, and the uniqueness properties of the combined and cumulative estimators when using a generalized inverse. The authors also provide methods for assessing goodness of fit in the linear model setting, as standard residual-based diagnostics cannot be performed with the cumulative data without access to historical data. Instead, they propose outlier tests relying on predictive residuals, which are based on the predictive values computed from the cumulative estimate of the regression coefficients attained at the previous accumulation point. Additionally, they introduce a new online-updated estimator of the regression coefficients and corresponding estimator of the standard error in the estimating equation setting that takes advantage of information from the previous data. They show theoretically that this new estimator, the cumulative updated estimating equation (CUEE) estimator, is asymptotically consistent, and show empirically that the CUEE estimator is less biased in their finite sample simulations than the cumulatively estimated version of the estimator of @Lin:Xi:aggr:2011.

Criterion-Based Variable Selection with Online Updating {#sec:varsel}
=======================================================

To the best of our knowledge, criterion-based variable selection has not yet been considered in the online updating context. This problem is well worth investigating especially when access/storage of the historical data is limited. Suppose that we have $K$ blocks of data in a sequence with ${{\bf{Y}}}_k$, ${{\bf{X}}}_k$, and $n_k$ being the $n_k$-dimensional vector of responses, the $n_k \times (p + 1)$ matrix of covariates, and the sample size, respectively, for the $k_{th}$ block, $k = 1, \ldots, K$, such that ${{\bf{Y}}}= (Y'_1, Y'_2, \dots, Y'_K )'$ and ${{\bf{X}}}= ({{\bf{X}}}'_1, \dots, {{\bf{X}}}'_k)'$. Consider the standard linear regression model for the whole data with sample size $n = \sum_{i=1}^k n_k$, $${{\bf{Y}}}= {{\bf{X}}}{{\mbox{\boldmath $\beta$}}}+ {{\mbox{\boldmath $\epsilon$}}},$$ where ${{\mbox{\boldmath $\beta$}}}$ is the regression coefficient vector, and ${{\mbox{\boldmath $\epsilon$}}}$ is a normal random vector with mean zero and variance $\theta I_n$. Let $\mathcal{M}$ denote the model space. We enumerate the models in $\mathcal{M}$ by $m=1, 2, ..., 2^p$, where $2^p$ is the dimension of $\mathcal{M}$. For the full model, the least squares estimate of ${{\mbox{\boldmath $\beta$}}}$ and the sum of squared errors based on the $k$th subset is given by $\hat{{{\mbox{\boldmath $\beta$}}}}_{n_k,k}=({{\bf{X}}}'_k{{\bf{X}}}_k)^{-}{{\bf{X}}}'_k{{\bf{Y}}}_k$ and ${\mathsf{SSE}}_{n_k, k}$. In the sequential setting, we only need to store and update the cumulative estimates at each $k$ [see, e.g. @Schifano2015].

Let ${{\mbox{\boldmath $\beta$}}}^{(m)}_{k}=(\beta^{(m)}_{0}, \beta^{(m)}_{1},  \ldots, \beta^{(m)}_{{p_m}})'$ and ${\mathsf{SSE}}^{(m)}_{k}$ denote the cumulative estimates based on all data through subset $k$ for model $m$, where $p_m$ is the number of covariates for model $m$. We further introduce the $(p + 1) \times (p_m +1)$ selection matrix $P^{(m)} = (e_{m_0}, \quad e_{m_1}, \quad \dots \quad e_{m_{p_m}})$, where $e_{m_0}$ is a vector with length $(p+1)$ and the first element as 1, and $e_{m_j}$ denotes a vector of length $(p+1)$ with 1 in the $m_{j}$th position and 0 in every other position for all $j > 0$. Here $(m_1, ..., m_{p_m})$ are not necessarily in sequence, but represents the index of selected variables in the full design matrix $X_k$. Define ${{\bf{X}}}^{(m)}_{k} = {{\bf{X}}}_kP^{(m)}$. Update a $(p_m + 1) \times (p_m + 1)$ matrix $$V^{(m)}_{k}= {{\bf{X}}}^{(m)'}_{k} {{\bf{X}}}^{(m)}_{k} + V^{(m)}_{k-1},$$ where $V^{(m)}_0=0$, and a $(p_m + 1) \times 1$ vector $$A^{(m)}_k = {{\bf{X}}}^{(m)'}_{k}{{\bf{X}}}_k\hat{{{\mbox{\boldmath $\beta$}}}}_{n_k,k}+V^{(m)}_{k-1}\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)}_{k-1},$$ where $\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)}_{0} = 0$. After some algebra, we have $$\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)}_{k} = (V^{(m)}_k)^{-1} A^{(m)}_k,$$ and $$\begin{aligned}
{\mathsf{SSE}}^{(m)}_{k} &= {\mathsf{SSE}}_{n_k k}+\hat{{{\mbox{\boldmath $\beta$}}}}_{n_k k}'{{\bf{X}}}'_k{{\bf{X}}}_k\hat{{{\mbox{\boldmath $\beta$}}}}_{n_k k}+\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)'}_{k-1}V^{(m)}_{k-1}\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)}_{k-1} \\
& \quad - \hat{{{\mbox{\boldmath $\beta$}}}}^{(m)'}_{k}V^{(m)}_{k}\hat{{{\mbox{\boldmath $\beta$}}}}^{(m)}_{k} + {\mathsf{SSE}}^{(m)}_{k-1}.\end{aligned}$$

With $\sigma$ unknown, letting $$\begin{aligned}
B_k^{(m)} & = n \log \frac{2\pi {\mathsf{SSE}}^{(m)}_{k}}{n - p_m - 1},\end{aligned}$$ the Akaike information criterion (AIC) and Bayesian information criterion (BIC) are updated by $$\begin{aligned}
{\mathsf{AIC}}_k^{(m)} &=  B_k^{(m)}  + n + p_m + 1,\\
{\mathsf{BIC}}_k^{(m)} &=  B_k^{(m)}  + n - p_m - 1 + (p_{m}+1)\log n.\end{aligned}$$

To study the Bayesian variable selection criteria, assume a joint conjugate prior for $({{\mbox{\boldmath $\beta$}}}^{(m)},{\theta}^{(m)})$ as follows: ${{\mbox{\boldmath $\beta$}}}^{(m)} | {\theta}^{(m)}$ follows normal distribution with mean ${{\mbox{\boldmath $\mu$}}}_o$, and precision matrix ${{\bf{V}}}_0$, ${\theta}^{(m)}$ follows Inverse Gamma distribution with shape parameter $\nu_0/2$ and scale parameter $\tau_0/2$, e.g, $$\begin{aligned}
  \pi( {{\mbox{\boldmath $\beta$}}}^{(m)},{\theta}^{(m)}| & {{\mbox{\boldmath $\mu$}}}_0, {{\bf{V}}}_0, \nu_0, \tau_0) \\
  & = \pi({{\mbox{\boldmath $\beta$}}}^{(m)}|{\theta}^{(m)},{{\mbox{\boldmath $\mu$}}}_0, {{\bf{V}}}_0) \pi({\theta}^{(m)}|\nu_0,\tau_0),\end{aligned}$$ where ${{\mbox{\boldmath $\mu$}}}_0$ is a prespecified $(p_m+1)$-dimensional vector, ${{\bf{V}}}_0$ is a $(p_m+1)\times (p_m+1)$ positive definite matrix, $\nu_0 > 0$, $\tau_0 > 0$. It can be shown that the deviance information criterion (${\mathsf{DIC}}$) [@Spie:etal:baye:2002] is updated by $${\mathsf{DIC}}_k^{(m)} = n \log\frac{\pi(n-2) {\mathsf{SSE}}_k^{(m)}}{2} + 2n\psi(\frac{n}{2}) + 2p_m +n + 4,$$ where $\psi(x)={\mathrm{d}\!}\log \Gamma(x) / {\mathrm{d}\!}x$ is the digamma function.

We examined the performance of AIC, BIC and DIC under the online updating scenario in a simulation study. Each dataset was generated from linear model $y_i= {{\bf{x}}}_i'{{\mbox{\boldmath $\beta$}}}+ \epsilon_i,$ where $\epsilon_i$’s were independently generated from $N(0, 100)$, $x_i = (1, x_{i1}, x_{i2}, x_{i3}, x_{i4})$ were identically distributed random vectors from a multivariate normal distribution with mean $(1, 0, 0, 0, 0)$ and marginal variances $(0, 16, 9, 0.3, 3)$. Two correlation structures of $(x_{i1}, x_{i2}, x_{i3}, x_{i4})$ were considered: 1) independent and 2) AR(1) with correlation coefficient 0.9. Four different models as determined by the nonzeroness of ${{\mbox{\boldmath $\beta$}}}$ were considered: $(-1, 3, 0, 0, 0)$, $(-1, 3, 0, -1.5, 0)$, $(-1, 3, 2, -1.5, 0)$, and $(-1, 3, 2, -1.5, 1)$. The corresponding signal-to-noise ratios were 1.44, 1.45, 1.81, and 1.83 in the independent case and 1.44, 1.29, 2.85, and 3.33 under the dependent case. The sample size of each block was set as $n_k=100$. The performance of the criteria was investigated with the cumulative estimates at block $k \in \{2, 25, 100\}$. For each scenario, 10,000 independent datasets were generated.

[max width=0.85]{}

[l ccc ccc ccc ccc ccc ccc]{} (lr)[2-10]{} (lr)[11-19]{} & &\
(lr)[2-10]{} (lr)[11-19]{} True & & & & & &\
(lr)[2-4]{}(lr)[5-7]{}(lr)[8-10]{} (lr)[11-13]{}(lr)[14-16]{}(lr)[17-19]{} Model & AIC &BIC &DIC& AIC &BIC &DIC& AIC &BIC &DIC &AIC &BIC &DIC& AIC &BIC &DIC& AIC &BIC &DIC\
\
none &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
$\bf(x_1)$ &**59 &**93 &**59 &**60 &**98 &**60 &**59 &**99 &**59 &**63 &**94 &**62 &**64 &**99 &**64 &**64 &**99 &**64\
($x_2$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2$) &11 &2 &11 &11 &1 &11 &12 &0 &12 &10 &2 &10 &9 &1 &9 &10 &0 &10\
($x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3$) &11 &2 &11 &11 &1 &11 &11 &0 &11 &8 &2 &8 &8 &0 &8 &8 &0 &8\
($x_2, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_3$) &2 &0 &3 &2 &0 &2 &2 &0 &2 &4 &0 &4 &3 &0 &3 &3 &0 &3\
($x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_4$) &11 &2 &11 &11 &0 &11 &11 &0 &11 &9 &2 &9 &8 &0 &9 &8 &0 &8\
($x_2, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_4$) &2 &0 &2 &2 &0 &2 &2 &0 &2 &3 &0 &3 &3 &0 &3 &3 &0 &3\
($x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3, x_4$) &2 &0 &2 &2 &0 &2 &2 &0 &2 &4 &0 &4 &4 &0 &4 &4 &0 &4\
($x_2, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_3, x_4$)&1 &0 &1 &0 &0 &0 &0 &0 &0 &1 &0 &1 &1 &0 &1 &1 &0 &1\
\[5pt\]\
none &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1$) &42 &83 &42 &0 &9 &0 &0 &0 &0 &55 &89 &55 &10 &60 &10 &0 &3 &0\
($x_2$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2$) &8 &2 &8 &0 &0 &0 &0 &0 &0 &11 &3 &11 &10 &4 &10 &1 &2 &1\
($x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
$\bf(x_1, x_3)$ &**28 &**12 &**27 &**71 &**90 &**71 &**70 &**100 &**70 &**13 &**4 &**13 &**50 &**30 &**50 &**69 &**90 &**69\
($x_2, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_3$) &6 &0 &6 &13 &0 &13 &14 &0 &14 &4 &0 &4 &6 &0 &6 &12 &0 &12\
($x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_4$) &8 &2 &8 &0 &0 &0 &0 &0 &0 &10 &3 &10 &14 &6 &14 &3 &5 &3\
($x_2, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_4$) &2 &0 &2 &0 &0 &0 &0 &0 &0 &3 &0 &3 &2 &0 &2 &2 &0 &2\
($x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3, x_4$) &6 &0 &6 &13 &0 &13 &13 &0 &13 &4 &0 &5 &6 &0 &6 &11 &0 &11\
($x_2, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_3, x_4$)&1 &0 &1 &2 &0 &3 &3 &0 &3 &1 &0 &1 &1 &0 &1 &2 &0 &2\
\[5pt\]\
none &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &2 &17 &2 &0 &0 &0 &0 &0 &0\
($x_2$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2$) &50 &85 &50 &0 &9 &0 &0 &0 &0 &64 &74 &64 &28 &83 &28 &1 &29 &1\
($x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &3 &2 &3 &0 &0 &0 &0 &0 &0\
($x_2, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
$\bf(x_1, x_2, x_3)$ &**33 &**13 &**33 &**84 &**90 &**84 &**84 &**100 &**84 &**14 &**3 &**14 &**50 &**14 &**50 &**81 &**67 &**81\
($x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &1 &1 &1 &0 &0 &0 &0 &0 &0\
($x_2, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_4$) &10 &2 &10 &0 &0 &0 &0 &0 &0 &11 &2 &11 &15 &3 &15 &6 &4 &6\
($x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &1 &0 &1 &0 &0 &0 &0 &0 &0\
($x_2, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
$(x_1, x_2, x_3, x_4)$ &7 &0 &7 &15 &0 &15 &16 &0 &16 &4 &0 &5 &7 &0 &7 &13 &0 &13\
\[5pt\]\
none &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 & 0\
($x_1$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &3 &0 &0 &0 &0 &0 &0 &0\
($x_2$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2$) &9 &40 &9 &0 &0 &0 &0 &0 &0 &51 &75 &51 &0 &13 &0 &0 &0 &0\
($x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &4 &6 &4 &0 &0 &0 &0 &0 &0\
($x_2, x_3$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_3$) &6 &6 &6 &0 &0 &0 &0 &0 &0 &7 &1 &7 &0 &0 &0 &0 &0 &0\
($x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &4 &10 &4 &0 &0 &0 &0 &0 &0\
($x_2, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_2, x_4$) &50 &47 &50 &0 &9 &0 &0 &0 &0 &24 &4 &25 &51 &80 &51 &11 &65 &11\
($x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_1, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
($x_2, x_3, x_4$) &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0 &0\
$\bf(x_1, x_2, x_3, x_4)$ &**34 &**7 &**34 &**100 &**91 &**100 &**100 &**100 &**100 &**10 &**1 &**10 &**48 &**7 &**48 &**89 &**35 &**89\
************************************************************************************************************************************************

The percentages of models selected among the $2^4$ models by each of the three criteria are summarized in Table \[tab:sim\]. The entire row in bold represents the true model. Based on the simulation results, BIC performs extremely well when the number of blocks ($k$) is large, which is consistent with known results that the probability of selecting the true model by BIC approaches 1 as $n \rightarrow \infty$ [e.g., @Schw:esti:1978; @Nish:asym:1984]. The BIC also performs better than AIC and DIC when the covariates are independent, even for small sample sizes. When covariates are highly dependent, AIC and DIC provide more reliable results when sample size is small. The performance of AIC and DIC is always very similar. The simulation results also confirm the existing theorem that AIC is not consistent [e.g., @Wood:on:1982]. In the big data setting with large sample size, BIC is generally preferable, especially when the covariates are not highly correlated.

Open Source R and R Packages {#sec:R}
============================

Handling big data is one of the topics of high performance computing. As the most popular open source statistical software, and its adds-on packages provide a wide range of high performance computing; see Comprehensive Archive Network (CRAN) task view on “High-Performance and Parallel Computing with R” [@Edde:cran:2014]. The focus of this section is on how to break the computer memory barrier and the computing power barrier in the context of big data.

Breaking the Memory Barrier {#sec:R:mem}
---------------------------

The size of big data is relative to the available computing resources. The theoretical limit of random access memory (RAM) is determined by the width of memory addresses: 4 gigabyte (GB) ($2^{32}$ bytes) for a 32-bit computer and 16.8 million terabyte ($2^{64}$ bytes) for a 64-bit computer. In practice, however, the latter is limited by the physical space of a computer case, the operating system, and specific software. Individual objects in have limits in size too; an user can hardly work with any object of size close to that limit. @Emer:Kane:dont:2012 suggested that a data set would be considered *large* if it exceeds 20% of RAM on a given machine and *massive* if it exceeds 50%, in which case, even the simplest calculation would consume all the remaining RAM.

Memory boundary can be broken with an external memory algorithms (EMA) [e.g., @Vitt:exte:2001], which conceptually works by storing the data on a disk storage (which has a much greater limit than RAM), and processing one chunk of it at a time in RAM [e.g., @Rpkg:biglm]. The results from each chunk will be saved or updated and the process continues until the entire dataset is exhausted; then, if needed as in an iterative algorithm, the process is reset from the beginning of the data. To implement an EMA for each statistical function, one need to address 1) data management and 2) numerical calculation.

### Data Management

Earlier solutions to oversize data resorted to relational databases. This method depends on an external database management system (DBMS) such as , , , , , , and others. Interfaces to are provided through many packages such as [[sqldf]{}]{} [@Rpkg:sqldf], [[DBI]{}]{} [@Rpkg:dbi], [[RSQLite]{}]{} [@Rpkg:rsqlite], and others. The database approach requires a DBMS to be installed and maintained, and knowledge of structured query language (); an exception for simpler applications is package [[filehash]{}]{} [@Peng:inte:2006], which comes with a simple key-value database implementation itself. The numerical functionality of SQL is quite limited, and calculations for most statistical analyses require copying subsets of the data into objects in facilitated by the interfaces. Extracting chunks from an external DBMS is computationally much less efficient than the more recent approaches discussed below [@Kane:Emer:West:scal:2013].

Two packages, [[bigmemory]{}]{} [@Kane:Emer:West:scal:2013] and [[ff]{}]{} [@Rpkg:ff] provide data structures for massive data while retaining a look and feel of objects. Package [[bigmemory]{}]{} defines a data structure for numeric matrices which uses memory-mapped files to allow matrices to exceed the RAM size on computers with 64-bit operating systems. The underling technology is memory mapping on modern operating systems that associates a segment of virtual memory in a one-to-one correspondence with contents of a file. These files are accessed at a much faster speed than in the database approaches because operations are handled at the operating-system level. The structure has several advantages such as support of shared memory for efficiency in parallel computing, reference behavior that avoids unnecessary temporary copies of massive objects, and column-major format that is compatible with legacy linear algebra packages (e.g., , ) [@Kane:Emer:West:scal:2013]. The package provides utility to read in a csv file to form a object, but it only allows one type of data, numeric; it is a numeric matrix after all.

Package [[ff]{}]{} provides data structures that are stored in binary flat files but behave (almost) as if they were in RAM by transparently mapping only a section (pagesize) of meta data in main memory. Unlike [[bigmemory]{}]{}, it supports ’s standard atomic data types (e.g., double or logical) as well as nonstandard, storage efficient atomic types (e.g., the 2-bit unsigned type allows efficient storage of genomic data as a factor with levels A, T, G, and, C). It also provides class which is like in , and import/export filters for csv files. A binary flat file can be shared by multiple objects in the same or multiple processes for parallel access. Utility functions allow interactive process of selections of big data.

### Numerical Calculation

The data management systems in packages [[bigmemory]{}]{} or [[ff]{}]{} do not mean that one can apply existing functions yet. Even a simple statistical analysis such as linear model or survival analysis will need to be implemented for the new data structures. Chunks of big data will be processed in RAM one at a time, and often, the process needs to be iterated over the whole data. A special case is the linear model fitting, where one pass of the data is sufficient and no resetting from the beginning is needed. Consider a regression model $E[Y] = X\beta$ with $n \times 1$ response $Y$, $n\times p$ model matrix $X$ and $p\times 1$ coefficient $\beta$. The base implementation takes $O(np + p^2)$ memory, which can be reduced dramatically by processing in chunks. The first option is to compute $X'X$ and $X'y$ in increment, and get the least squares estimate of $\beta$, $\hat\beta = (X'X)^{-1}X'Y$. This method is adopted in package [[speedglm]{}]{} [@Rpkg:speedglm]. A slower but more accurate option is to compute the incremental QR decomposition [@Mill:algo:1992] of $X = QR$ to get $R$ and $Q'Y$, and then solve $\beta$ from $R\beta = Q'Y$. This option is implemented in package [[biglm]{}]{} [@Rpkg:biglm]. Function uses only $p^2$ memory of $p$ variables and the fitted object can be updated with more data using . The package also provides an incremental computation of sandwich variance estimator by accumulating a $(p + 1)^2\times (p + 1)^2$ matrix of products of $X$ and $Y$ without a second pass of the data.

In general, a numerical calculation needs an iterative algorithm in computation and, hence, multiple passes of the data are necessary. For example, a GLM fitting is often obtained through the iterated reweighted least squares (IRLS) algorithm. The function in package [[biglm]{}]{} implements the generic IRLS algorithm that can be applied to any specific data management system such as DBMS, [[bigmemory]{}]{}, or [[ff]{}]{}, provided that a function supplies the next chunk of data or zero-row data if there is no more, and resets to the beginning of the data for the next iteration. Specific implementation of the function for object of class and are provided in package [[biganalytics]{}]{} [@Rpkg:biganalytics] and [[ffbase]{}]{} [@Rpkg:ffbase], respectively.

For any statistical analysis on big data making use of the data management system, one would need to implement the necessary numerical calculations like what package [[biglm]{}]{} does for GLM. The family of [[bigmemory]{}]{} provides a collection of functions for objects: [[biganalytics]{}]{} for basic analytic and statistical functions, [[bigtabulate]{}]{} for tabulation operations [@Rpkg:bigtabulate], and [[bigalgebra]{}]{} for matrix operation with the BLAS and LAPACK libraries [@Rpkg:bigalgebra]. Some additional functions for objects are available from other contributed packages, such as [[bigpca]{}]{} for principal component analysis and single-value decomposition [@Rpkg:bigpca], and [[bigrf]{}]{} for random forest [@Rpkg:bigrf]. For objects, package [[ffbase]{}]{} provides basic statistical functions [@Rpkg:ffbase]. Additional functions for objects are provided in other packages, with examples including [[biglars]{}]{} for least angle regression and LASSO [@Rpkg:biglars] and [[PopGenome]{}]{} for population genetic and genomic analysis [@Pfei:popg:2014].

If some statistical analysis, such as generalized estimating equations or Cox proportional hazards model, has not been implemented for big data, then one will need to modify the existing algorithm to implement it. As pointed out by @Kane:Emer:West:scal:2013 [p.5], this would open Pandora’s box of recoding which is not a long-term solution for scalable statistical analyses; this calls for redesign of the next-generation statistical programming environment which could provide seamless scalability through file-backed memory-mapping for big data, help avoid the need for specialized tools for big data management, and allow statisticians and developers to focus on new methods and algorithms.

Breaking the Computing Power Barrier {#sec:R:comp}
------------------------------------

### Speeding Up

As a high level interpreted language, for which most of instructions are executed directly, is infamously slow with loops. Some loops can be avoided by taking advantage of the vectorized functions in or by clever vectorizing with some effort. When vectorization is not straightforward or loops are unavoidable, as in the case of MCMC, acceleration is much desired, especially for big data. The least expensive tool in a programmer’s effort to speed up code is to compile them to byte code with the [[compiler]{}]{} package, which was developed by Luke Tierney and is now part of base . The byte code compiler translates the high-level into a very simple language that can be interpreted by a very fast byte code interpreter, or virtual machine. Starting with  2.14.0 in 2011, the base and recommended packages were pre-compiled into byte-code by default. Users’ functions, expressions, scripts, and packages can be compiled for an immediate boost in speed by a factor of 2 to 5.

Computing bottlenecks can be implemented in a compiled language such as or and interfaced to through ’s foreign language interfaces [@Rext ch.5]. Typical bottlenecks are loops, recursions, and complex data structures. Recent developments have made the interfacing with much easier than it used to be [@Edde:seam:2013]. Package [[inline]{}]{} [@Rpkg:inline] provides functions that wrap (or ) code as strings in and takes care of compiling, linking, and loading by placing the resulting dynamically-loadable object code in the per-session temporary directory used by . For more general usage, package [[Rcpp]{}]{} [@Edde:etal:rcpp:2011] provides classes for many basic data types, which allow straightforward passing of data in both directions. Package [[RcpEigen]{}]{} [@Rpkg:rcppeigen] provides access to the high-performance linear algebra library for a wide variety of matrix methods, various decompositions and support of sparse matrices. Package [[RcppArmadillo]{}]{} [@Edde:Sand:rcpp:2014] connects with , a powerful templated linear algebra library which provides a good balance between speed and ease of use. Package [[RInside]{}]{} [@Rpkg:rinside] gives easy access of objects from by wrapping the existing embedding application programming interface (API) in classes. The [[Rcpp]{}]{} project has revolutionized the integration of with ; it is now used by hundreds of packages.

Diagnostic tools can help identify the bottlenecks in code. Package [[microbenchmark]{}]{} [@Rpkg:microbenchmark] provides very precise timings for small pieces of source code, making it possible to compare operations that only take a tiny amount of time. For a collection of code, run-time of each individual operation can be measured with realistic inputs; the process is known as profiling. Function in does the profiling, but the outputs are not intuitive to understand for many users. Packages [[proftools]{}]{} [@Rpkg:proftools] and [[aprof]{}]{} [@Rpkg:aprof] provide tools to analyze profiling outputs. Packages [[profr]{}]{} [@Rpkg:profr], [[lineprof]{}]{} [@Rpkg:lineprof], and [[GUIProfiler]{}]{} [@Rpkg:GUIProfiler] provide visualization of profiling results.

### Scaling Up

The package system has long embraced integration of parallel computing of various technologies to address the big data challenges. For embarrassingly parallelizable jobs such as bootstrap or simulation, where there is no dependency or communication between parallel tasks, many options are available with computer clusters or multicores. @schm:etal:stat:2009 reviewed the then state-of-the-art parallel computing with , highlighting two packages for cluster use: [[Rmpi]{}]{} [@Yu:rmpi:2002] which provides an interface to the Message Passing Interface (MPI) in parallel computing; [[snow]{}]{} [@Ross:Tier:Li:simp:2007] which provides an abstract layer with the communication details hidden from the end users. Since then, some packages have been developed and some discontinued. Packages [[snowFT]{}]{} [@Rpkg:snowFT] and [[snowfall]{}]{} [@Rpkg:snowfall] extend [[snow]{}]{} with fault tolerance and wrappers for easier development of parallel programs. Package [[multicore]{}]{} [@Rpkg:multicore] provides parallel processing of code on machines with multiple cores or CPUs. Its work and some of [[snow]{}]{} have been incorporated into the base package [[parallel]{}]{}, which was first included in  2.14.0 in 2011. Package [[foreach]{}]{} [@Rpkg:foreach] allows general iteration over elements in a collection without any explicit loop counter. Using loop without side effects facilitates executing the loop in parallel with different parallel mechanisms, including those provided by [[parallel]{}]{}, [[Rmpi]{}]{}, and [[snow]{}]{}. For massive data that exceed the computer memory, a combination of [[foreach]{}]{} and [[bigmemory]{}]{}, with shared-memory data structure referenced by multiple processes, provides a framework with ease of development and efficiency of execution (both in speed and memory) as illustrated by @Kane:Emer:West:scal:2013. Package [[Rdsm]{}]{} provides facilities for distributed shared memory parallelism at the level, and combined with [[bigmemory]{}]{}, it enables parallel processing on massive, out-of-core matrices.

The “Programming with Big Data in ” project (pbdR) enables high-level distributed data parallelism in with easy utilization of large clusters with thousands of cores [@pbdR2012]. Big data are interpreted quite literally to mean that a dataset requires parallel processing either because it does not fit in the memory of a single machine or because its processing time needs to be made tolerable. The project focuses on distributed memory systems where data are distributed across processors and communications between processors are based on MPI. It consists of a collection of packages in a hierarchy. Package [[pbdMPI]{}]{} provides classes to directly interface with MPI to support the Single Program Multiple Data (SPMD) parallelism. Package [[pbdSLAP]{}]{} serves as a mechanism to utilize a subset of functions of scalable dense linear algebra in [@Blac:etal:scal:1997], a subset of routines redesigned with the SPMD style. Package [[pbdBASE]{}]{} contains a set of wrappers of low level functions in , upon which package [[pbdMAT]{}]{} builds to provide distributed dense matrix computing while preserving the friendly and familiar syntax for these computations. Demonstrations on how to use these and other packages from the pbdR are available in package [[pbdDEMO]{}]{}.

A recent, widely adopted open source framework for massive data storage and distributed computing is [@Hadoop]. Its heart is an implementation of the MapReduce programming model first developed at Google [@Dean:Ghem:mapr:2008], which divides the data to distributed systems and computes for each group (the map step), and then recombines the results (the reduce step). It provides fault tolerant and scalable storage of massive datasets across machines in a cluster [@Whit:hado:2011]. The model suits perfectly the embarrassingly parallelizable jobs and the distributed file system helps break the memory boundary. @McCa:West:para:2011 [ch.5–8] demonstrated three ways to combine and . The first is to submit scripts directly to a cluster, which gives the user the most control and the most power, but comes at the cost of a learning curve. The second is a pure solution via package [[Rhipe]{}]{}, which hides the communications to from users. The package (not on CRAN) is from the project, which stands for and Integrated Programming Environment [@Guha:etal:larg:2012]. With [[Rhipe]{}]{}, data analysts only need to write code for the map step and the reduce step [@Guha:etal:larg:2012], and get the power of without leaving . The third approach targets specifically the Elastic MapReduce (EMR) at Amazon by a CRAN package [[segue]{}]{} [@Rpkg:segue], which makes EMR as easy to use as a parallel backend for -style operations. An alternative open source project that connects and is the RHadoop project, which is actively being developed by Revolution Analytics [@RHadoop]. This project is a collection of packages that allow users to manage and analyze data with : [[rhbase]{}]{} provides functions for database management for the HBase distributed database, [[rhdfs]{}]{} provides functions for distributed file system (HDFS), [[rmr]{}]{} provides functions to MapReduce functionality, [[plymr]{}]{} provides higher level data processing for structured data, and the most recent addition [[ravro]{}]{} provides reading and writing functions for files in format, an efficient data serialization system developed at Apache [@Avro].

is a more recent, cousin project of that supports tools for big data related tasks [@Spark]. The functions of and are neither the exactly same nor mutually exclusive, and they often work together. has its own distributed storage system, which is fundamental for any big data computing framework, allowing vast datasets to be stored across the hard drives of a scalable computer cluster rather than on a huge costly hold-it-all device. It persists back to the disk after a map or reduce action. In contrast, does not have its own distributed file system, and it processes data in-memory [@Zaha:etal:2010]. The biggest difference is disk-based computing versus memory-based computing. This is why could work 100 times faster than hadoop for some applications when the data fit in the memory. Some applications such as machine learning or stream processing where data are repeatedly queried makes an ideal framework. For big data that does not fit in memory, ’s operators spill data to disk, allowing it to run well on any sized data. For this purpose, it can be installed on top to take advantage of ’s HDFS. An frontend to is provided in package [[SparkR]{}]{} [@Rpkg:SparkR], which has become part of Apache recently. By using ’s distributed computation engine, the package allows users to run large scale data analysis such as selection, filtering, aggregation from . @Kara:etal:2015 provides a summary of the state-of-the-art on using .

As multicores have become the standard setup for computers today, it is desirable to automatically make use of the cores in implicit parallelism without any explicit requests from the user. The experimental packages [[pnmath]{}]{} and [[pnmath0]{}]{} by Luke Tierney replace a number of internal vector operations in with alternatives that can take advantage of multicores [@Tier:code:2009]. For a serial algorithm such as MCMC, it is desirable to parallelize the computation bottleneck if possible, but this generally involves learning new computing tools and the debugging can be challenging. For instance, @Yan:etal:para:2007 used the parallel linear algebra package (PLAPACK) [@Geji:usin:1997] for the matrix operations (especially the Cholesky decomposition) in a MCMC algorithm for Bayesian spatiotemporal geostatistical models, but the scalability was only moderate.

When random numbers are involved as in the case of simulation, extra care is needed to make sure the parallelized jobs run independent (and preferably reproducible) random-number streams. Package [[rsprng]{}]{} [@Rpkg:rsprng] provides an interface to the Scalable Parallel Random Number Generators (SPRNG) [@Masc:Srin:algo:2000]. Package [[rlecuyer]{}]{} [@Rpkg:rlecuyer] provides an interface to the random number generator with multiple independent streams developed by @L'Ec:etal:2002, the ideas of which are also implemented in the base package [[parallel]{}]{}: make independent streams by separating a single stream with a sufficiently large number of steps apart. Package [[doRNG]{}]{} [@Rpkg:doRNG] provides functions to perform reproducible parallel loops, independent of the parallel environment and associated backend.

From a hardware perspective, many computers have mini clusters of graphics processing units (GPUs) that can help with bottlenecks. GPUs are dedicated numerical processors that were originally designed for rendering three dimensional computer graphics. A GPU has hundreds of processor cores on a single chip and can be programmed to apply the same numerical operations on large data array. @Such:etal:unde:2010 investigated the use of GPUs in massively parallel massive mixture modeling, and showed better performance of GPUs than multicore CPUs, especially for larger samples. To reap the advantage, however, one needs to learn the related tools such as Compute Unified Device Architecture (CUDA), Open Computing Language (OpenCL), and so on, which may be prohibitive. An package [[gputools]{}]{} [@Rpkg:gputools] provides interface to NVidia CUDA toolkit and others.

If one is willing to step out of the comfort zone of and take full control/responsibility of parallel computing, one may program with open source MPI or Open Multi-Processing (OpenMP). MPI is a language-independent communication system designed for programming on parallel computers, targeting high performance, scalability and portability [@Pach:para:1997]. Most MPI implementations are available as libraries from , , and any language that can interface with such libraries, including , or . The interface from can be accessed with package [[Rmpi]{}]{} [@Yu:rmpi:2002] as mentioned earlier. Freely available implementations include OpenMPI (not OpenMP) and MPICH, while others come with license such as Intel MPI. OpenMP is an API that supports multi-platform shared memory multiprocessing programming in and on most processor architectures and operating systems [@chapman2008using]. It is an add on to compilers (e.g., , intel compiler) to take advantage of of shared memory systems such as multicore computers where processors shared the main memory. MPI targets both distributed as well as shared momory systems while OpenMP targets only shared memory systems. MPI provides both process and thread based approach while OpenMP provides only thread based parallilism. OpenMP uses a portable, scalable model that gives programmers a simple and flexible interface for writing multi-threaded programs in and [@Dagu:Enon:open:1998]. Debugging parallel programs can be very challenging.

Commercial Statistical Software {#sec:comm}
===============================

is the core product of Revolution Analytics (formerly Revolution Computing), a company that provides tools, support, and training. focuses on big data, large scale multiprocessor (or high performance) computing, and multicore functionality. Massive datasets are handled via EMA and parallel EMA (PEMA) when multiprocessors or multicores are available. The commercial package [[RevoScaleR]{}]{} [@Rpkg:RevoScaleR] breaks the memory boundary by a special data format that allows efficient storage and retrieval of data. Functions in the package (e.g., for GLM fitting) know to work on a massive dataset one chunk at a time. The computing power boundary is also addressed — functions in the package can exploit multicores or computer clusters. Packages from the aforementioned open source project RHadoop developed by the company provide support for . Other components in allow high speed connection for various types of data sources and threading and inter-process communication for parallel and distributed computing. The same code works on small and big data, and on workstations, servers, clusters, , or in the cloud. The single workstation version of is free for academic use currently, and was used in the case study in Section \[sec:case\].

, one of the most widely used commercial software for statistical analysis, provides big data support through SAS High Performance Analytics. Massive datasets are approached by grid computing, in-database processing, in-memory analytics and connection to . The SAS High Performance Analytics Products include statistics, econometrics, optimization, forecasting, data mining, and text mining, which, respectively, correspond to SAS products STAS, ETS, OR, high-performance forecasting, enterprise miner, and text miner [@Cohe:Rodr:high:2013].

IBM , the Statistical Product and Services Solution, provides big data analytics through SPSS Modeler, SPSS Analytic Server, SPSS Collaboration and Deployment Services, and SPSS Analytic Catalyst [@SPSS]. SPSS Analytic Server is the foundation and it focuses on high performance analytics for data stored in Hadoop-based distributed systems. SPSS modeler is the high-performance data mining workbench, utilizing SPSS Analytic Server to leverage big data in Hadoop environments. Analysts can define analysis in a familiar and accessible workbench to conduct analysis modeling and scoring over high volumes of varied data. SPSS Collaboration and Deployment Services helps manage analytical assets, automate processes and efficiently share results widely and securely. SPSS Analytic Catalyst is the automation of analysis that makes analytics and data more accessible to users.

provides a number of tools to tackle the challenges of big data analytics [@Matlab]. Memory mapped variables map a file or a proportion of a file to a variable in RAM; disk variables direct access to variables from files on disk; datastore allows access to data that do not fit into RAM. Their combination addresses the memory boundary. The computation power boundary is broken by intrinsic multicore math, GPU computing, parallel computing, cloud computing, and support.

A Case Study {#sec:case}
============

The airline on-time performance data from the 2009 ASA Data Expo (<http://stat-computing.org/dataexpo/2009/the-data.html>) is used as a case study to demonstrate a logistic model fitting with a massive dataset that exceeds the RAM of a single computer. The data is publicly available and has been used for demonstration with big data by @Kane:Emer:West:scal:2013 and others. It consists of flight arrival and departure details for all commercial flights within the USA, from October 1987 to April 2008. About 12 million flights were recorded with 29 variables. A compressed version of the pre-processed data set from the bigmemory project (<http://data.jstatsoft.org/v55/i14/Airline.tar.bz2>) is approximately 1.7GB, and it takes 12GB when uncompressed.

The response of the logistic regression is late arrival which was set to 1 if a flight was late by more than 15 minutes and 0 otherwise. Two binary covariates were created from the departure time: night (1 if departure occurred between 8pm and 5am) and weekend (1 if departure occurred on weekends and 0 otherwise). Two continuous covariates were included: departure hour (DepHour, range 0 to 24) and distance from origin to destination (in 1000 miles). In the raw data, the departure time was an integer of the HHmm format. It was converted to minutes first to prepare for DepHour. Three methods are considered in the case study: 1) combination of with package [[bigmemory]{}]{}; 2) combination of with package [[ff]{}]{}; and 3) the academic, single workstation version of . The default settings of [[ff]{}]{} were used. Before fitting the logistic regression, the 12GB raw data needs to be read in from the csv format, and new variables needs to be generated. This leads to a total of $120,748,239$ observations with no missing data. The scripts for the three methods are in the supplementary materials for interested readers.

                        Reading   Transforming   Fitting
  ------------------- --------- -------------- ---------
    [[bigmemory]{}]{}     968.6          105.5    1501.7
           [[ff]{}]{}    1111.3          528.4    1988.0
                          851.7          107.5     189.4

  : Timing results (in seconds) for reading in the whole 12GB data, transforming to create new variables, and fitting the logistic regression with three methods: [[bigmemory]{}]{}, [[ff]{}]{}, and .[]{data-label="tab:timing"}

The scripts were executed in batch mode on a 8-core machine running CenOS (a free Linux distribution functionally compatible with Red Hat Enterprise Linux which is officially supported by ), with Intel Core i7 2.93GHz CPU, and 16GB memory. Table \[tab:timing\] summarizes the timing results of reading in the whole 12GB data, transforming to create new variables, and fitting the logistic regression with the three methods. The chunk sizes were set to be 500,000 observations for all three methods. For , this was set when reading in the data to the format; for the other two methods, this was set at the fitting stage using function . Under the current settings, has a clear advantage in fitting with only 8% of the time used by the other two approaches. This is a result of the joint force of its using all 8 cores implicitly and efficient storage and retrieval of the data; the version of the data is about 1/10 of the size of the external files saved by [[bigmemory]{}]{} or [[ff]{}]{}. Using [[bigmemory]{}]{} and using [[ff]{}]{} in had very similar performance in fitting the logistic regression, but the former took less time in reading, and significantly less time (only about 1/5) in transforming variables of the latter. The [[bigmemory]{}]{} method was quite close to the method in the reading and the transforming tasks. The [[ff]{}]{} method took longer in reading and transforming than the [[bigmemory]{}]{} method, possibly because it used much less memory.

                  Estimate  Std. Error ($\times 10^4$)
  ------------- ---------- ----------------------------
   (Intercept)    $-$2.985            9.470
     DepHour         0.104            0.601
    Distance         0.235            4.032
      Night       $-$0.448            8.173
     Weekend      $-$0.177            5.412

  : Logistic regression results for late arrival.[]{data-label="tab:fit"}

The results of the logistic regression are identical from all methods, and are summarized in Table \[tab:fit\]. Flights with later departure hour or longer distance are more likely to be delayed. Night flights or weekend flights are less likely to be delayed. Given the huge sample size, all coefficients were highly significant. It is possible, however, that p-values can still be useful. A binary covariate with very low rate of event may still have an estimated coefficient with a not-so-low p-value [@Schifano2015], an effect only estimable with big data.

                           1      2     3     4     5     6     7     8
  ------------------- ------ ------ ----- ----- ----- ----- ----- -----
    [[bigmemory]{}]{}   22.1   11.2   7.8   6.9   6.2   6.3   6.4   6.8
           [[ff]{}]{}   21.4   11.0   7.1   6.7   5.8   5.9   6.1   6.8

  : Time results (in seconds) for parallel computing quantiles of departure delay for each day of the week with 1 to 8 cores using [[foreach]{}]{}.[]{data-label="tab:para"}

As an illustration of [[foreach]{}]{} for embarrassingly parallel computing, the example in @Kane:Emer:West:scal:2013 is expanded to include both [[bigmemory]{}]{} and [[ff]{}]{}. The task is to find three quantiles (0.5, 0.9, and 0.99) of departure delays for each day of the week; that is, 7 independent jobs can run on 7 cores separately. To make the task bigger, each job was set to run twice. The resulting 14 jobs were parallelized with on the same Linux machine using 1 to 8 cores for the sake of illustration. The script is included in the supplementary materials. The timing results are summarized in Table \[tab:para\]. There is little difference between the two implementations. When there is no communication overhead, with 14 jobs one would expect the run time to reduce to 1/2, 5/14, 4/14, 3/14, 3/14, 2/14, and 2/14, respectively, with 2, 3, 4, 5, 6, 7 and 8 cores. The impact of communication cost is obvious in Table \[tab:para\]. The time reduction is only closer to the expectation in the ideal case when the number of cores is smaller.

Discussion {#sec:disc}
==========

This article presents a recent snapshot on statistical analysis with big data that exceed the memory and computing capacity of a single computer. Albeit under-appreciated by the general public or even mainstream academic community, computational statisticians have made respectable progress in extending standard statistical analysis to big data, with the most notable achievements in the open source community. Packages [[bigmemory]{}]{} and [[ff]{}]{} make it possible in principle to implement any statistical analysis with their data structure. Nonetheless, for anything that has not been already implemented (e.g., survival analysis, generalized estimating equations, mixed effects model, etc.), one would need to implement an EMA version of the computation task, which may not be straightforward and may involve some steep learning curves. allows easy extension of algorithms that do not require multiple passes of the data, but such analyses are mostly descriptive. An example is visualization, an important tool in exploratory analysis. With big data, the bottleneck is the number of pixels in the screen. The bin-summarize-smooth framework for visualization of large data of @Wick:bin:2014 with package [[bigvis]{}]{} [@Rpkg:bigvis] may be adapted to work with .

Big data present challenges much further beyond the territory of classic statistics, requiring joint workforce with domain knowledge, computing skills, and statistical thinking [@Yu:let:2014]. Statisticians have much to contribute to both the intellectual vitality and the practical utility of big data, but will have to expand their comfort zone to engage high-impact, real world problems which are often less structured or with ambiguity [@Jord:Lin:stat:2014]. Examples are to provide structure for poorly defined problems, or to develop methods/models for new types of data such as image or network. As suggested by @Yu:let:2014, to play a critical role in the arena of big data or own data science, statisticians need to work on real problems and relevant methodology and theory will follow naturally.

Acknowledgement {#acknowledgement .unnumbered}
===============

The authors thank Stephen Archut, Fang Chen, and Joseph Rickert for the big data analytics information on , , and . An earlier version of the manuscript was presented at the “Statistical and Computational Theory and Methodology for Big Data Analysis” workshop in February, 2014, at the Banff International Research Station in Banff, AB, Canada. The discussions and comments from the workshop participants are gratefully acknowledged.

Supplementary Materials {#supplementary-materials .unnumbered}
=======================

Four scripts (and their outputs), along with a descriptive README file are provided for the case study. The first three are the logistic regression with, respectively, combination of [[bigmemory]{}]{} with (), combination of [[ff]{}]{} with (), and (); their output files have extensions. The first two run with , while the third one needs . The fourth script is for the parallel computing with [[foreach]{}]{} combined with [[bigmemory]{}]{} and [[ff]{}]{}, respectively.
