Data Processing and Analysis with {snow}

The {snow} package in R, which stands for “Simple Network of Workstations,” is designed for parallel computing, allowing users to distribute computational tasks across multiple processors or computers to speed up processing. It provides a high-level interface for parallelizing R code, particularly for tasks that can be split into independent subtasks, such as simulations, bootstrapping, or large-scale data processing. Below is an introduction to {snow}, how it works, and its key components.

The {snow} package, developed by Luke Tierney and others, enables parallel computing in R by leveraging multiple processors or machines in a cluster. It abstracts much of the complexity of parallel programming, making it accessible for users who need to perform computationally intensive tasks. {snow} supports various parallel backends, including MPI (Message Passing Interface), PVM (Parallel Virtual Machine), and socket-based communication, allowing flexibility depending on the system setup.

The package is particularly useful for:

  • Embarrassingly parallel problems: Tasks that can be divided into independent subtasks with no need for communication between them (e.g., Monte Carlo simulations, cross-validation).

  • Cluster computing: Distributing tasks across multiple machines in a network.

  • Multi-core processing: Utilizing multiple cores on a single machine.

{snow} is lightweight, flexible, and integrates well with R’s functional programming paradigm. It is often used alongside other parallel computing packages like {parallel} or {foreach} for specific use cases.

How {snow} Works

{snow} operates by creating a cluster of worker processes (or nodes) that execute tasks in parallel. The main steps in using {snow} are:

  1. Cluster Setup:
    • Initialize a cluster using a chosen communication backend (e.g., MPI, sockets, or PVM).
    • Specify the number of nodes (workers) and their configuration (e.g., local cores or remote machines).
    • Common functions: makeCluster() to create the cluster, specifying the type (e.g., “MPI”, “SOCK”, “PVM”) and number of nodes.
  2. Task Distribution:
    • Divide the computational task into smaller, independent subtasks.
    • Use {snow} functions like clusterApply(), clusterApplyLB(), or parLapply() to distribute these tasks across the cluster’s nodes.
    • Each node processes its assigned task independently, and results are collected back to the master process.
  3. Communication:
    • {snow} handles communication between the master (the main R session) and worker nodes using the chosen backend.
    • Workers execute R functions or scripts sent by the master, and results are returned via the communication protocol (e.g., sockets or MPI).
  4. Cluster Termination:
    • Once tasks are complete, the cluster is shut down using stopCluster() to free resources.

Key Functions in {snow}

Here are the core functions in {snow} and their purposes:

  • makeCluster(spec, type, ...): Creates a cluster of workers. spec specifies the number of nodes or machine names, and type defines the backend (e.g., “SOCK” for sockets, “MPI” for MPI).
  • clusterApply(cl, x, fun, ...): Applies a function fun to each element of x across the cluster cl. Each worker processes a subset of x.
  • clusterApplyLB(cl, x, fun, ...): Similar to clusterApply, but with load balancing, dynamically assigning tasks to workers as they become available.
  • parLapply(cl, x, fun, ...): A parallel version of lapply, distributing tasks across the cluster.
  • clusterExport(cl, list): Exports variables or functions from the master’s environment to the workers’ environments.
  • clusterEvalQ(cl, expr): Evaluates an expression on each worker, useful for initializing worker environments (e.g., loading packages).
  • stopCluster(cl): Shuts down the cluster, freeing resources.

Example Workflow

Here’s a simple example of using {snow} for parallel computation with a socket-based cluster:

Code 
# Load the snow package
library(snow)

# Create a cluster with 4 workers (using sockets)
cl <- makeCluster(4, type = "SOCK")

# Export a variable to all workers
my_data <- 1:10
clusterExport(cl, "my_data")

# Define a function to compute squares
square <- function(x) x^2

# Apply the function in parallel
result <- parLapply(cl, my_data, square)

# Stop the cluster
stopCluster(cl)

# View results
print(unlist(result))
 [1]   1   4   9  16  25  36  49  64  81 100

Output (assuming no errors):

[1]   1   4   9  16  25  36  49  64  81 100

In this example:

  • A cluster of 4 workers is created using sockets.
  • The my_data vector is exported to all workers.
  • The square function is applied to each element of my_data in parallel using parLapply.
  • Results are collected and returned as a list, which is unlisted for display.

How It Works Internally

  • Master-Worker Model: The master R session coordinates tasks, while workers (separate R processes) execute them. The master sends tasks and collects results.
  • Communication Backend:
    • SOCK: Uses TCP/IP sockets for communication, suitable for local or networked machines.
    • MPI: Uses the MPI standard (via {Rmpi} package) for high-performance clusters.
    • PVM: Uses the Parallel Virtual Machine system (less common today).
  • Load Balancing: Functions like clusterApplyLB ensure efficient task distribution by assigning new tasks to workers as they complete their current ones.
  • Serialization: R objects (data, functions) are serialized and sent to workers, which deserialize and process them. Results are serialized and sent back to the master.

Advantages of {snow}

  • Simplicity: Provides a high-level interface, hiding much of the complexity of parallel programming.
  • Flexibility: Supports multiple backends, making it adaptable to different hardware setups (local multi-core machines or distributed clusters).
  • Scalability: Works for both small-scale (e.g., multi-core laptops) and large-scale (e.g., HPC clusters) computations.
  • Integration: Works seamlessly with R’s functional programming style and other packages like {foreach}.

Limitations

  • Overhead: Setting up a cluster and communicating data introduces overhead, which may outweigh benefits for small tasks.
  • Dependencies: Some backends (e.g., MPI) require additional software installation and configuration (e.g., {Rmpi}).
  • Memory: Each worker loads a copy of the data and R environment, which can be memory-intensive.
  • Not for Interdependent Tasks: {snow} is best for embarrassingly parallel tasks; tasks requiring frequent communication between workers are better suited for other frameworks (e.g., MPI directly).

Practical Considerations

  • Choosing a Backend: Use “SOCK” for simplicity on local machines or small networks. Use “MPI” for high-performance clusters if {Rmpi} is installed.
  • Environment Setup: Ensure all necessary packages and objects are exported to workers using clusterExport or clusterEvalQ.
  • Task Granularity: Break tasks into appropriately sized chunks to balance computation time and communication overhead.
  • Error Handling: Use try-catch mechanisms in functions to handle errors on workers, as they run independently.

Parallel Data Processing and Basic Statistics with {snow}

This tutorial demonstrates how to use the {snow} package in R for parallel computing, applied to data processing and basic statistical analyses on a large dataset: the NYC Yellow Taxi trip records for January 2023 (yellow_tripdata_2023-01.parquet). The dataset contains approximately 3 million rows of taxi trip data, including variables like pickup/dropoff times, distances, fares, and more. We’ll cover parallel data loading/processing, descriptive statistics, correlation, and simple/multiple linear regression.

{snow} enables distributing tasks across multiple cores or machines, which is ideal for large datasets to reduce computation time. We’ll focus on “embarrassingly parallel” tasks, such as processing data chunks independently and combining results.

Check and Install Required R Packages

Code 
packages<- c(
'tidyverse',
'data.table',
'feather',
'arrow',
'parallel',
'future',
'furrr',
'future.apply',
'nycflights13',
'glue',
'snow',
'matrixStats'
)
#| warning: false
#| error: false

# Install missing packages
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

Verify Instalation

Code 
# Verify 
cat("Installed packages:\n")
Installed packages:
Code 
print(sapply(packages, requireNamespace, quietly = TRUE))
   tidyverse   data.table      feather        arrow     parallel       future 
        TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
       furrr future.apply nycflights13         glue         snow  matrixStats 
        TRUE         TRUE         TRUE         TRUE         TRUE         TRUE

Load R-packages

Code 
# Load packages with suppressed messages
invisible(lapply(packages, function(pkg) {
  suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}))

Check Loaded Packages

Code 
# Check loaded packages
cat("Successfully loaded packages:\n")
Successfully loaded packages:
Code 
print(search()[grepl("package:", search())])
 [1] "package:matrixStats"  "package:glue"         "package:nycflights13"
 [4] "package:future.apply" "package:furrr"        "package:future"      
 [7] "package:parallel"     "package:arrow"        "package:feather"     
[10] "package:data.table"   "package:lubridate"    "package:forcats"     
[13] "package:stringr"      "package:dplyr"        "package:purrr"       
[16] "package:readr"        "package:tidyr"        "package:tibble"      
[19] "package:ggplot2"      "package:tidyverse"    "package:snow"        
[22] "package:stats"        "package:graphics"     "package:grDevices"   
[25] "package:utils"        "package:datasets"     "package:methods"     
[28] "package:base"

Loading and Exploring the Data Sequentially

Load the data using {arrow} for memory efficiency. Convert to a data.frame for easier manipulation (note: this collects the data into memory).

Code 
# Read the Parquet file
dataFolder<-"/home/zia207/WebSites/R_Website/High_Performance_Computing_R/R_Markdown/Data/"
# Load January 2023 data
df <- read_parquet(paste0(dataFolder, "yellow_tripdata_2023-01.parquet"))
# Basic exploration
dim(df)  # Dimensions: ~3,066,766 rows, 19 columns
[1] 3066766      19

The dataset is large, so sequential processing can be slow for intensive tasks. We’ll use {snow} to parallelize.

Setting Up the {snow} Cluster

Create a cluster using available cores (e.g., 4 workers). Use the “SOCK” type for local multi-core processing.

Always stop the cluster when done to free resources:

Code 
# At the end of your script/session
stopCluster(cl)

Parallel Data Processing: Cleaning and Feature Engineering

  • For large data, split into chunks, process each in parallel (e.g., add features, filter invalid rows), and recombine.

  • First, split the data into a list of chunks (e.g., 4 chunks for 4 workers).

Code 
# Split data into chunks
n_chunks <- n_cores
chunk_indices <- rep(1:n_chunks, each = ceiling(nrow(df) / n_chunks), length.out = nrow(df))
df_chunks <- split(df, chunk_indices)

Define a processing function (e.g., compute trip duration in minutes, filter invalid trips: positive distance, fare, duration > 0, reasonable passenger count).

Code 
# Processing function for each chunk
process_chunk <- function(chunk) {
  # Compute trip duration (in minutes)
  chunk$trip_duration <- as.numeric(difftime(chunk$tpep_dropoff_datetime, 
                                             chunk$tpep_pickup_datetime, units = "mins"))
  
  # Filter invalid rows
  chunk <- chunk[chunk$trip_distance > 0 & 
                 chunk$fare_amount > 0 & 
                 chunk$trip_duration > 0 & 
                 chunk$passenger_count >= 1 & 
                 chunk$passenger_count <= 6, ]  # Reasonable passenger range
  
  return(chunk)
}

Apply in parallel and recombine:

This parallelizes the cleaning, speeding up for large datasets.

Code 
# At the end of your script/session
# Safely stop the cluster
if (!is.null(cl)) {
  tryCatch({
    stopCluster(cl)
  }, error = function(e) {
    cat("Error stopping cluster:", conditionMessage(e), "\n")
  })
  cl <- NULL
}

Parallel Descriptive Statistics

  • For descriptive stats (mean, median, sd, min, max), compute on chunks and aggregate.

  • Focus on numeric variables: trip_distance, fare_amount, tip_amount, trip_duration.

  • Define functions to compute local stats per chunk, then combine.

Code 
# Create cluster with error handling
n_cores <- parallel::detectCores() - 1
cl <- NULL
tryCatch({
  cl <- makeCluster(n_cores, type = "SOCK")
}, error = function(e) {
  cat("Error creating cluster:", conditionMessage(e), "\n")
  stop("Cluster creation failed.")
})

# Load required libraries on workers

clusterEvalQ(cl, {
  library(arrow)
  library(matrixStats)
})
[[1]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[2]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[3]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[4]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[5]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[6]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[7]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[8]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[9]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[10]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[11]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[12]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[13]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[14]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[15]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[16]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[17]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[18]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[19]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[20]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[21]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[22]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[23]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[24]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[25]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[26]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[27]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[28]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[29]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[30]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[31]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[32]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[33]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[34]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[35]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[36]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[37]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[38]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[39]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[40]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[41]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[42]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[43]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[44]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[45]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[46]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[47]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[48]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[49]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[50]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[51]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[52]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[53]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[54]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[55]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[56]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[57]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[58]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[59]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[60]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[61]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[62]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[63]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[64]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[65]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[66]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[67]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[68]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[69]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[70]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[71]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[72]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[73]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[74]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[75]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[76]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[77]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[78]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[79]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"
Code 
# Variables of interest
vars <- c("trip_distance", "fare_amount", "tip_amount", "trip_duration")

# Function to compute local sums, counts, etc., for mean/sd
local_stats <- function(chunk, vars) {
  list(
    n = nrow(chunk),
    sums = colSums(chunk[, vars], na.rm = TRUE),
    sum_squares = colSums(chunk[, vars]^2, na.rm = TRUE),
    mins = apply(chunk[, vars], 2, min, na.rm = TRUE),
    maxs = apply(chunk[, vars], 2, max, na.rm = TRUE),
    medians = apply(chunk[, vars], 2, median, na.rm = TRUE)
  )
}

# Split processed data into chunks
n_chunks <- n_cores
df_chunks <- split(df_processed, rep(1:n_chunks, each = ceiling(nrow(df_processed) / n_chunks), length.out = nrow(df_processed)))

# Export and compute local stats in parallel
clusterExport(cl, c("local_stats", "vars"))
local_results <- parLapply(cl, df_chunks, local_stats, vars = vars)

# Combine results
total_n <- sum(sapply(local_results, function(x) x$n))
global_sums <- rowSums(sapply(local_results, function(x) x$sums))
global_sum_squares <- rowSums(sapply(local_results, function(x) x$sum_squares))
global_means <- global_sums / total_n
global_vars <- (global_sum_squares / total_n) - (global_means)^2
global_sds <- sqrt(global_vars)

# Compute global min/max using rowMins/rowMaxs
global_mins <- matrixStats::rowMins(sapply(local_results, function(x) x$mins), na.rm = TRUE)
global_maxs <- matrixStats::rowMaxs(sapply(local_results, function(x) x$maxs), na.rm = TRUE)

# For medians, compute sequentially on full data (or approximate)
global_medians <- apply(df_processed[, vars], 2, median, na.rm = TRUE)

# Display in a table
descriptive_stats <- data.frame(
  Variable = vars,
  Mean = global_means,
  SD = global_sds,
  Min = global_mins,
  Median = global_medians,
  Max = global_maxs
)
print(descriptive_stats)
                   Variable      Mean        SD        Min Median      Max
trip_distance trip_distance  3.397158 21.922633 0.01000000   1.80 14098.55
fare_amount     fare_amount 18.142215 17.112374 0.01000000  12.80  1160.10
tip_amount       tip_amount  3.326246  3.780658 0.00000000   2.79   380.80
trip_duration trip_duration 15.439935 42.659988 0.01666667  11.55  6179.40
Code 
# Safely stop the cluster
if (!is.null(cl)) {
  tryCatch({
    stopCluster(cl)
  }, error = function(e) {
    cat("Error stopping cluster:", conditionMessage(e), "\n")
  })
  cl <- NULL
}

This approach avoids loading the full dataset on each worker and computes aggregates efficiently.

Parallel Correlation Analysis

  • For correlation (e.g., Pearson), compute on the full numeric subset sequentially if memory allows, or approximate via chunks.

  • For parallel: Compute local covariance matrices, then combine.

  • Select numeric vars.

Code 
# Create cluster with error handling
n_cores <- parallel::detectCores() - 1
cl <- NULL
tryCatch({
  cl <- makeCluster(n_cores, type = "SOCK")
}, error = function(e) {
  cat("Error creating cluster:", conditionMessage(e), "\n")
  stop("Cluster creation failed.")
})


# Load required libraries on workers
clusterEvalQ(cl, {
  library(arrow)
  library(matrixStats)
})
[[1]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[2]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[3]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[4]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[5]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[6]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[7]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[8]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[9]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[10]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[11]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[12]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[13]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[14]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[15]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[16]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[17]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[18]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[19]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[20]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[21]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[22]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[23]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[24]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[25]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[26]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[27]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[28]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[29]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[30]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[31]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[32]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[33]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[34]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[35]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[36]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[37]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[38]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[39]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[40]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[41]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[42]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[43]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[44]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[45]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[46]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[47]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[48]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[49]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[50]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[51]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[52]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[53]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[54]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[55]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[56]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[57]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[58]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[59]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[60]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[61]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[62]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[63]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[64]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[65]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[66]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[67]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[68]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[69]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[70]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[71]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[72]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[73]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[74]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[75]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[76]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[77]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[78]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"       

[[79]]
 [1] "matrixStats" "arrow"       "snow"        "stats"       "graphics"   
 [6] "grDevices"   "utils"       "datasets"    "methods"     "base"
Code 
# Numeric data: select variables and remove NAs/infinite values
vars <- c("trip_distance", "fare_amount", "tip_amount", "trip_duration")
num_df <- df_processed[, vars, drop = FALSE]

# Check for non-numeric columns
if (!all(sapply(num_df, is.numeric))) {
  stop("Non-numeric columns detected in num_df. Please ensure all variables are numeric.")
}

# Remove rows with NA or infinite values
num_df <- num_df[complete.cases(num_df) & apply(num_df, 1, function(x) all(is.finite(x))), ]
if (nrow(num_df) == 0) {
  stop("No valid rows remain after removing NAs and infinite values.")
}

# Function for local covariance components
local_cov <- function(chunk) {
  # Ensure chunk has enough rows
  if (nrow(chunk) < 2) {
    return(list(n = 0, means = rep(0, length(vars)), sum_prods = matrix(0, nrow = length(vars), ncol = length(vars))))
  }
  
  # Remove NAs and infinite values
  chunk <- chunk[complete.cases(chunk) & apply(chunk, 1, function(x) all(is.finite(x))), , drop = FALSE]
  n <- nrow(chunk)
  
  if (n < 2) {
    return(list(n = 0, means = rep(0, length(vars)), sum_prods = matrix(0, nrow = length(vars), ncol = length(vars))))
  }
  
  # Compute means and covariance components
  means <- colMeans(chunk, na.rm = TRUE)
  sum_prods <- crossprod(scale(chunk, center = TRUE, scale = FALSE)) / (n - 1) * (n - 1)
  
  list(
    n = n,
    means = means,
    sum_prods = sum_prods
  )
}

# Split numeric data into chunks
n_chunks <- n_cores
num_chunks <- split(num_df, rep(1:n_chunks, each = ceiling(nrow(num_df) / n_chunks), length.out = nrow(num_df)))

# Export and compute local covariance in parallel
clusterExport(cl, c("local_cov", "vars"))
local_covs <- parLapply(cl, num_chunks, local_cov)

# Check local results for issues
cat("Local chunk sizes:\n")
Local chunk sizes:
Code 
print(sapply(local_covs, function(x) x$n))
    1     2     3     4     5     6     7     8     9    10    11    12    13 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   14    15    16    17    18    19    20    21    22    23    24    25    26 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   27    28    29    30    31    32    33    34    35    36    37    38    39 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   40    41    42    43    44    45    46    47    48    49    50    51    52 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   53    54    55    56    57    58    59    60    61    62    63    64    65 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   66    67    68    69    70    71    72    73    74    75    76    77    78 
36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 36509 
   79 
36451
Code 
# Combine results
total_n <- sum(sapply(local_covs, function(x) x$n))
if (total_n < 2) {
  stop("Not enough valid observations across chunks to compute correlations.")
}

global_means <- rowSums(sapply(local_covs, function(x) x$means * x$n)) / total_n
global_cov <- Reduce(`+`, lapply(local_covs, function(x) {
  if (x$n == 0) {
    matrix(0, nrow = length(vars), ncol = length(vars))
  } else {
    x$sum_prods + outer(x$means - global_means, x$means - global_means) * x$n
  }
})) / total_n

# Compute correlation from covariance
global_sd <- sqrt(diag(global_cov))
if (any(global_sd == 0)) {
  warning("Zero standard deviation detected for some variables. Correlation will be NA for those.")
}
global_cor <- global_cov / outer(global_sd, global_sd)

# Ensure row and column names
dimnames(global_cor) <- list(vars, vars)

# Display correlation matrix
print(round(global_cor, 3))
              trip_distance fare_amount tip_amount trip_duration
trip_distance         1.000       0.194      0.120         0.043
fare_amount           0.194       1.000      0.608         0.222
tip_amount            0.120       0.608      1.000         0.126
trip_duration         0.043       0.222      0.126         1.000
Code 
# Safely stop the cluster
if (!is.null(cl)) {
  tryCatch({
    stopCluster(cl)
  }, error = function(e) {
    cat("Error stopping cluster:", conditionMessage(e), "\n")
  })
  cl <- NULL
}

Parallel Linear Regression

Code 
# Create cluster with error handling
n_cores <- parallel::detectCores() - 1
cl <- NULL
tryCatch({
  cl <- makeCluster(n_cores, type = "SOCK")
}, error = function(e) {
  cat("Error creating cluster:", conditionMessage(e), "\n")
  stop("Cluster creation failed.")
})

# Ensure data is clean
vars <- c("trip_distance", "fare_amount", "passenger_count", "trip_duration")
reg_df <- df_processed[, vars, drop = FALSE]
reg_df <- reg_df[complete.cases(reg_df) & apply(reg_df, 1, function(x) all(is.finite(x))), ]

# Function to fit model on a chunk and return coefficients and sample size
fit_model <- function(data, formula) {
  if (nrow(data) < 2) {
    return(list(coefficients = NULL, n = 0))
  }
  model <- lm(formula, data = data)
  list(coefficients = coef(model), n = nrow(data))
}

# Split data into chunks
n_chunks <- n_cores
chunks <- split(reg_df, rep(1:n_chunks, each = ceiling(nrow(reg_df) / n_chunks), length.out = nrow(reg_df)))

Simple Linear Regression (SLR)

Code 
# --- OLS: fare_amount ~ trip_distance ---
# Define formula explicitly
ols_formula <- fare_amount ~ trip_distance

# Export and fit OLS in parallel
clusterExport(cl, c("fit_model", "reg_df", "ols_formula"))
ols_results <- parLapply(cl, chunks, function(chunk) fit_model(chunk, ols_formula))

# Combine OLS results (weighted by sample size)
ols_n <- sapply(ols_results, function(x) x$n)
total_n <- sum(ols_n)
if (total_n < 2) stop("Not enough valid observations to fit model.")
ols_coefs <- Reduce("+", lapply(ols_results, function(x) if (is.null(x$coefficients)) rep(0, 2) else x$coefficients * x$n)) / total_n

# Display OLS results
cat("OLS Coefficients (Parallel):\n")
OLS Coefficients (Parallel):
Code 
names(ols_coefs) <- c("(Intercept)", "trip_distance")
print(round(ols_coefs, 3))
  (Intercept) trip_distance 
        8.176         3.014
Code 
# Compare with sequential fit
ols_seq <- lm(fare_amount ~ trip_distance, data = reg_df)
cat("\nOLS Coefficients (Sequential):\n")

OLS Coefficients (Sequential):
Code 
print(round(coef(ols_seq), 3))
  (Intercept) trip_distance 
       18.025         0.150

Multiple Linear Regression (MLR)

Code 
# --- MLR: fare_amount ~ trip_distance + passenger_count + trip_duration ---
# Define formula explicitly
mlr_formula <- fare_amount ~ trip_distance + passenger_count + trip_duration

# Export and fit MLR in parallel
clusterExport(cl, c("fit_model", "reg_df", "mlr_formula"))
mlr_results <- parLapply(cl, chunks, function(chunk) fit_model(chunk, mlr_formula))

# Combine MLR results
mlr_n <- sapply(mlr_results, function(x) x$n)
total_n <- sum(mlr_n)
if (total_n < 2) stop("Not enough valid observations to fit model.")
mlr_coefs <- Reduce("+", lapply(mlr_results, function(x) if (is.null(x$coefficients)) rep(0, 4) else x$coefficients * x$n)) / total_n

# Display MLR results
cat("\nMLR Coefficients (Parallel):\n")

MLR Coefficients (Parallel):
Code 
names(mlr_coefs) <- c("(Intercept)", "trip_distance", "passenger_count", "trip_duration")
print(round(mlr_coefs, 3))
    (Intercept)   trip_distance passenger_count   trip_duration 
          7.552           3.000           0.186           0.025
Code 
# Compare with sequential fit
mlr_seq <- lm(fare_amount ~ trip_distance + passenger_count + trip_duration, data = reg_df)
cat("\nMLR Coefficients (Sequential):\n")

MLR Coefficients (Sequential):
Code 
print(round(coef(mlr_seq), 3))
    (Intercept)   trip_distance passenger_count   trip_duration 
         16.053           0.143           0.475           0.085
Code 
# Safely stop the cluster
if (!is.null(cl)) {
  tryCatch({
    stopCluster(cl)
  }, error = function(e) {
    cat("Error stopping cluster:", conditionMessage(e), "\n")
  })
  cl <- NULL
}

Simple and Multiple Linear Regression with Parallel Bootstrapping

  • Bootstrapping estimates variability of coefficients by resampling data with replacement.
  • This is embarrassingly parallel, as each bootstrap sample can be processed independently.
  • We’ll use {snow} to parallelize bootstrap resampling and model fitting.
Code 
# Create cluster with error handling
n_cores <- parallel::detectCores() - 1
cl <- NULL
tryCatch({
  cl <- makeCluster(n_cores, type = "SOCK")
}, error = function(e) {
  cat("Error creating cluster:", conditionMessage(e), "\n")
  stop("Cluster creation failed.")
})

# Ensure data is clean
vars <- c("trip_distance", "fare_amount", "passenger_count", "trip_duration")
reg_df <- df_processed[, vars, drop = FALSE]
reg_df <- reg_df[complete.cases(reg_df) & apply(reg_df, 1, function(x) all(is.finite(x))), ]

# Bootstrap function for regression
bootstrap_lm <- function(data, idx, formula) {
  model <- lm(formula, data = data[idx, ])
  coef(model)
}

# Generate bootstrap indices
n_boot <- 1000
boot_indices <- lapply(1:n_boot, function(i) sample(1:nrow(reg_df), replace = TRUE))

Simple Linear Regression (SLR)

Code 
# --- OLS: fare_amount ~ trip_distance ---
# Sequential fit for reference
ols_model <- lm(fare_amount ~ trip_distance, data = reg_df)
cat("OLS Coefficients:\n")
OLS Coefficients:
Code 
print(summary(ols_model)$coefficients)
                Estimate   Std. Error   t value Pr(>|t|)
(Intercept)   18.0253620 0.0099879609 1804.7089        0
trip_distance  0.1499316 0.0004452982  336.6992        0
Code 
# Define OLS formula explicitly
ols_formula <- fare_amount ~ trip_distance

# Parallel bootstrap for OLS
clusterExport(cl, c("bootstrap_lm", "reg_df", "ols_formula"))
boot_results_ols <- parLapplyLB(cl, boot_indices, function(idx) bootstrap_lm(reg_df, idx, ols_formula))
boot_coefs_ols <- do.call(rbind, boot_results_ols)
ols_ci <- apply(boot_coefs_ols, 2, quantile, probs = c(0.025, 0.975))
cat("\nOLS Bootstrap 95% CIs:\n")

OLS Bootstrap 95% CIs:
Code 
print(round(ols_ci, 3))
      (Intercept) trip_distance
2.5%       17.601         0.094
97.5%      18.215         0.274
Code 
# Convert boot_coefs_ols to a data frame for easier plotting
ols_boot_df <- as.data.frame(boot_coefs_ols)
colnames(ols_boot_df) <- c("(Intercept)", "trip_distance")

# Sequential OLS fit for reference
ols_seq_coefs <- coef(ols_model)

# Base R: Histograms for OLS coefficients
par(mfrow = c(1, 2))  # Set up a 1x2 plot grid
for (coef_name in colnames(ols_boot_df)) {
  hist(ols_boot_df[[coef_name]], main = paste("Bootstrap Distribution:", coef_name),
       xlab = "Coefficient Estimate", col = "lightblue", breaks = 30)
  abline(v = ols_seq_coefs[coef_name], col = "red", lwd = 2, lty = 2)  # Sequential estimate
}

Code 
par(mfrow = c(1, 1))  # Reset plot grid

Multiple Linear Regression (MLR)

Code 
# --- MLR: fare_amount ~ trip_distance + passenger_count + trip_duration ---
# Sequential fit for reference
mlr_model <- lm(fare_amount ~ trip_distance + passenger_count + trip_duration, data = reg_df)
cat("\nMLR Coefficients:\n")

MLR Coefficients:
Code 
print(summary(mlr_model)$coefficients)
                   Estimate   Std. Error   t value Pr(>|t|)
(Intercept)     16.05262060 0.0181965390 882.17988        0
trip_distance    0.14270084 0.0004348296 328.17643        0
passenger_count  0.47490616 0.0108395463  43.81237        0
trip_duration    0.08481967 0.0002237588 379.06741        0
Code 
# Define MLR formula explicitly
mlr_formula <- fare_amount ~ trip_distance + passenger_count + trip_duration

# Parallel bootstrap for MLR
clusterExport(cl, c("bootstrap_lm", "reg_df", "mlr_formula"))
boot_results_mlr <- parLapplyLB(cl, boot_indices, function(idx) bootstrap_lm(reg_df, idx, mlr_formula))
boot_coefs_mlr <- do.call(rbind, boot_results_mlr)
mlr_ci <- apply(boot_coefs_mlr, 2, quantile, probs = c(0.025, 0.975))
cat("\nMLR Bootstrap 95% CIs:\n")

MLR Bootstrap 95% CIs:
Code 
print(round(mlr_ci, 3))
      (Intercept) trip_distance passenger_count trip_duration
2.5%       15.704         0.090           0.449         0.081
97.5%      16.215         0.261           0.498         0.089
Code 
# Convert boot_coefs_mlr to a data frame
mlr_boot_df <- as.data.frame(boot_coefs_mlr)
colnames(mlr_boot_df) <- c("(Intercept)", "trip_distance", "passenger_count", "trip_duration")

# Sequential MLR fit for reference
mlr_seq_coefs <- coef(mlr_model)

# Base R: Histograms for MLR coefficients
par(mfrow = c(2, 2))  # Set up a 2x2 plot grid
for (coef_name in colnames(mlr_boot_df)) {
  hist(mlr_boot_df[[coef_name]], main = paste("Bootstrap Distribution:", coef_name),
       xlab = "Coefficient Estimate", col = "lightgreen", breaks = 30)
  abline(v = mlr_seq_coefs[coef_name], col = "red", lwd = 2, lty = 2)  # Sequential estimate
}

Code 
par(mfrow = c(1, 1))  # Reset plot grid

# ggplot2: Density plots for MLR coefficients (optional)
mlr_boot_long <- reshape2::melt(mlr_boot_df, variable.name = "Coefficient", value.name = "Estimate")
ggplot(mlr_boot_long, aes(x = Estimate, fill = Coefficient)) +
  geom_density(alpha = 0.5) +
  geom_vline(data = data.frame(Coefficient = names(mlr_seq_coefs), Estimate = mlr_seq_coefs),
             aes(xintercept = Estimate), color = "red", linetype = "dashed") +
  facet_wrap(~ Coefficient, scales = "free") +
  theme_minimal() +
  labs(title = "MLR Bootstrap Coefficient Distributions", x = "Coefficient Estimate", y = "Density")

Code 
# Safely stop the cluster
if (!is.null(cl)) {
  tryCatch({
    stopCluster(cl)
  }, error = function(e) {
    cat("Error stopping cluster:", conditionMessage(e), "\n")
  })
  cl <- NULL
}

Cleanup and Best Practices

  • Always stop the cluster: stopCluster(cl)
  • Monitor memory: For very large data, use smaller chunks or sample first.
  • Compare timings: Use system.time() for sequential vs. parallel to see speedups.
  • Error handling: Wrap functions in tryCatch for robustness.

Summary and Conclusion

The {snow} package is a powerful tool for parallel computing in R, making it easier to speed up computations by distributing tasks across multiple cores or machines. Its simplicity and flexibility make it ideal for embarrassingly parallel problems, though users must consider overhead and setup requirements for optimal use. By leveraging functions like parLapply and clusterApply, users can significantly reduce computation time for large-scale analyses. This tutorial demonstrated how to set up a {snow} cluster, process large datasets in parallel, and perform statistical analyses efficiently. Always remember to stop the cluster after use to free up resources.

References

  1. {snow} Package Documentation (CRAN)
    • Relevance: 95%
    • Link: CRAN: snow
    • Citation: Tierney, L., et al. (2023). snow: Simple Network of Workstations. R package version 0.4-4.
  2. “Parallel R” by McCallum and Weston
    • Relevance: 85%
    • Link: O’Reilly
    • Citation: McCallum, Q. E., & Weston, S. (2011). Parallel R. O’Reilly Media.
  3. “Applied Linear Statistical Models” by Kutner et al.
    • Relevance: 75%
    • Citation: Kutner, M. H., et al. (2005). Applied Linear Statistical Models (5th ed.). McGraw-Hill.