R - Standard Deviation and Variance | Application Architect

Key Insights

Standard deviation and variance quantify data dispersion—variance measures average squared deviations from the mean, while standard deviation is its square root in original units
R provides sd() and var() functions with Bessel’s correction (n-1 denominator) by default for sample statistics, requiring manual calculation for population parameters
Understanding when to use sample versus population formulas, handling missing values, and applying these metrics across data frames are critical for real-world statistical analysis

Understanding Variance and Standard Deviation

Variance measures how far data points spread from their mean. It’s calculated by taking the average of squared differences from the mean. Standard deviation is simply the square root of variance, returning the measure to the original units of your data.

For a population:

Variance: σ² = Σ(xi - μ)² / N
Standard deviation: σ = √(σ²)

For a sample:

Variance: s² = Σ(xi - x̄)² / (n-1)
Standard deviation: s = √(s²)

The (n-1) denominator in sample calculations is Bessel’s correction, which provides an unbiased estimate of population variance.

# Create sample data
values <- c(12, 15, 18, 21, 24, 27, 30)

# Built-in functions (sample statistics)
sample_sd <- sd(values)
sample_var <- var(values)

print(paste("Sample SD:", sample_sd))
print(paste("Sample Variance:", sample_var))

# Verify relationship
print(paste("SD squared equals variance:", sample_sd^2 == sample_var))

Output:

[1] "Sample SD: 6.78232998312527"
[1] "Sample Variance: 46"
[1] "SD squared equals variance: TRUE"

Sample vs Population Calculations

R’s built-in functions calculate sample statistics by default. For population statistics, you need custom functions.

# Population standard deviation function
pop_sd <- function(x) {
  sqrt(sum((x - mean(x))^2) / length(x))
}

# Population variance function
pop_var <- function(x) {
  sum((x - mean(x))^2) / length(x)
}

# Compare sample vs population
data <- c(10, 20, 30, 40, 50)

cat("Sample SD:", sd(data), "\n")
cat("Population SD:", pop_sd(data), "\n")
cat("Sample Variance:", var(data), "\n")
cat("Population Variance:", pop_var(data), "\n")

# Show the difference
cat("\nDifference factor:", sd(data) / pop_sd(data), "\n")
cat("This equals sqrt(n/(n-1)):", sqrt(length(data)/(length(data)-1)), "\n")

Output:

Sample SD: 15.81139 
Population SD: 14.14214 
Sample Variance: 250 
Population Variance: 200 

Difference factor: 1.118034 
This equals sqrt(n/(n-1)): 1.118034

Handling Missing Values

Real-world data contains missing values. R’s sd() and var() return NA if any missing values exist unless you specify na.rm = TRUE.

# Data with missing values
incomplete_data <- c(5, 10, NA, 20, 25, NA, 35)

# Default behavior - returns NA
result1 <- sd(incomplete_data)
print(paste("Without na.rm:", result1))

# Remove NA values
result2 <- sd(incomplete_data, na.rm = TRUE)
print(paste("With na.rm:", result2))

# Count missing values
n_missing <- sum(is.na(incomplete_data))
n_complete <- sum(!is.na(incomplete_data))

cat("\nMissing values:", n_missing, "\n")
cat("Complete values:", n_complete, "\n")
cat("SD calculated from", n_complete, "observations\n")

# Custom function with NA handling and reporting
sd_with_report <- function(x) {
  clean_x <- x[!is.na(x)]
  list(
    sd = sd(clean_x),
    n_total = length(x),
    n_used = length(clean_x),
    n_missing = sum(is.na(x))
  )
}

print(sd_with_report(incomplete_data))

Working with Data Frames

Calculating standard deviation across data frame columns is common in exploratory data analysis.

# Create sample dataset
df <- data.frame(
  temperature = c(72, 75, 68, 71, 74, 76, 73),
  humidity = c(45, 52, 48, 50, 55, 58, 51),
  pressure = c(1013, 1015, 1012, 1014, 1016, 1015, 1013)
)

# Calculate SD for each column
sds <- sapply(df, sd)
print(sds)

# Calculate variance for each column
vars <- sapply(df, var)
print(vars)

# Create summary statistics
summary_stats <- data.frame(
  variable = names(df),
  mean = sapply(df, mean),
  sd = sapply(df, sd),
  variance = sapply(df, var),
  cv = sapply(df, function(x) sd(x) / mean(x) * 100)  # coefficient of variation
)

print(summary_stats)

Grouped Calculations

Calculate statistics by groups using base R or dplyr.

# Create grouped data
set.seed(42)
grouped_data <- data.frame(
  category = rep(c("A", "B", "C"), each = 10),
  value = c(rnorm(10, 100, 15), rnorm(10, 120, 20), rnorm(10, 90, 10))
)

# Base R approach
aggregate(value ~ category, data = grouped_data, FUN = sd)
aggregate(value ~ category, data = grouped_data, FUN = var)

# Using tapply
tapply(grouped_data$value, grouped_data$category, sd)

# dplyr approach (if available)
if (require(dplyr, quietly = TRUE)) {
  grouped_data %>%
    group_by(category) %>%
    summarise(
      n = n(),
      mean = mean(value),
      sd = sd(value),
      variance = var(value),
      se = sd(value) / sqrt(n())  # standard error
    )
}

Weighted Standard Deviation

When observations have different weights or frequencies, use weighted calculations.

# Weighted variance function
weighted_var <- function(x, w) {
  weighted_mean <- sum(x * w) / sum(w)
  sum(w * (x - weighted_mean)^2) / sum(w)
}

# Weighted SD function
weighted_sd <- function(x, w) {
  sqrt(weighted_var(x, w))
}

# Example: grade distribution
grades <- c(85, 90, 78, 92, 88)
frequencies <- c(5, 8, 3, 6, 4)  # number of students with each grade

# Unweighted (treats each unique grade equally)
unweighted_sd <- sd(grades)

# Weighted (accounts for frequency)
weighted_sd_result <- weighted_sd(grades, frequencies)

cat("Unweighted SD:", unweighted_sd, "\n")
cat("Weighted SD:", weighted_sd_result, "\n")

# Expand data to verify
expanded_grades <- rep(grades, frequencies)
cat("Expanded data SD:", sd(expanded_grades), "\n")

Rolling Standard Deviation

Calculate moving standard deviation for time series analysis.

# Rolling SD function
rolling_sd <- function(x, window) {
  n <- length(x)
  result <- rep(NA, n)
  
  for (i in window:n) {
    result[i] <- sd(x[(i - window + 1):i])
  }
  
  return(result)
}

# Time series example
set.seed(123)
time_series <- cumsum(rnorm(50, 0, 1))
window_size <- 10

rolling_sds <- rolling_sd(time_series, window_size)

# Display results
results_df <- data.frame(
  index = 1:length(time_series),
  value = time_series,
  rolling_sd = rolling_sds
)

print(tail(results_df, 15))

# Using zoo package for comparison (if available)
if (require(zoo, quietly = TRUE)) {
  zoo_rolling_sd <- rollapply(time_series, width = window_size, 
                               FUN = sd, align = "right", fill = NA)
  cat("\nZoo package produces same results:", 
      all.equal(rolling_sds, as.numeric(zoo_rolling_sd), 
                check.attributes = FALSE))
}

Coefficient of Variation

The coefficient of variation (CV) expresses standard deviation as a percentage of the mean, enabling comparison across different scales.

# CV function
cv <- function(x, na.rm = FALSE) {
  (sd(x, na.rm = na.rm) / mean(x, na.rm = na.rm)) * 100
}

# Compare variability across different scales
dataset1 <- c(100, 110, 120, 130, 140)  # mean = 120
dataset2 <- c(10, 11, 12, 13, 14)       # mean = 12

cat("Dataset 1 - SD:", sd(dataset1), "CV:", cv(dataset1), "%\n")
cat("Dataset 2 - SD:", sd(dataset2), "CV:", cv(dataset2), "%\n")

# Both have same relative variability despite different SDs

Performance Considerations

For large datasets, consider computational efficiency.

# Compare methods for large datasets
n <- 1000000
large_data <- rnorm(n)

# Built-in function
system.time({
  result1 <- sd(large_data)
})

# Manual calculation
system.time({
  result2 <- sqrt(sum((large_data - mean(large_data))^2) / (n - 1))
})

# Verify results match
print(all.equal(result1, result2))

# For matrices, use colSds from matrixStats package
if (require(matrixStats, quietly = TRUE)) {
  mat <- matrix(rnorm(1000000), ncol = 100)
  
  system.time({
    base_result <- apply(mat, 2, sd)
  })
  
  system.time({
    fast_result <- colSds(mat)
  })
}

Standard deviation and variance are fundamental for understanding data distribution, detecting outliers, and performing statistical tests. Master these calculations in R to build a solid foundation for advanced analytics and modeling work.