dev/se_variance.R

In friendly/heplots: Visualizing Hypothesis Tests in Multivariate Linear Models

# Function to calculate the standard error of variance
se_variance <- function(X) {
  # Remove any NA values
  X <- X[!is.na(X)]
  
  # Calculate sample variance
  sample_var <- var(X)
  
  # Calculate sample size
  n <- length(X)
  
  # Calculate standard error of variance
  # Using the formula: SE(s²) = sqrt((2σ⁴) / (n-1))
  se <- sqrt((2 * (sample_var^2)) / (n - 1))
  
  return(se)
}

# Function to calculate the fourth sample moment
fourth_sample_moment <- function(X) {
  # Remove any NA values
  X <- X[!is.na(X)]
  
  # Calculate the mean of X
  mean_X <- mean(X)
  
  # Calculate the fourth moment
  # Fourth moment is the average of (x - mean)^4
  fourth_moment <- mean((X - mean_X)^4)
  
  return(fourth_moment)
}
# Example usage:
# x <- c(1, 2, 3, 4, 5)
# se_variance(x)

# Function to calculate standard error of variance using bootstrap
bootstrap_variance_se <- function(X, num_bootstrap = 1000, seed = NULL) {
  # Set seed for reproducibility if provided
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  # Remove NA values
  X <- X[!is.na(X)]
  n <- length(X)
  
  # Bootstrap resampling
  bootstrap_variances <- sapply(1:num_bootstrap, function(i) {
    # Resample with replacement
    bootstrap_sample <- sample(X, size = n, replace = TRUE)
    
    # Calculate variance of the bootstrap sample
    var(bootstrap_sample)
  })
  
  # Calculate standard error as standard deviation of bootstrap variances
  bootstrap_se <- sd(bootstrap_variances)
  
  # Additional outputs for detailed analysis
  results <- list(
    standard_error = bootstrap_se,
    original_variance = var(X),
    bootstrap_variances = bootstrap_variances
  )
  
  return(results)
}

# Example usage function
example_bootstrap_variance_se <- function() {
  # Generate some sample data
  set.seed(123)
  sample_data <- rnorm(100, mean = 10, sd = 2)
  
  # Calculate bootstrap standard error of variance
  bootstrap_result <- bootstrap_variance_se(sample_data, num_bootstrap = 2000)
  
  # Print results
  cat("Original Sample Variance:", bootstrap_result$original_variance, "\n")
  cat("Bootstrap Standard Error of Variance:", bootstrap_result$standard_error, "\n")
}