Bootstrap Inference for PMM2 Models

Part 1: Bootstrap for Linear Regression

Basic Example: Simple Linear Regression

# Generate data with skewed errors
n <- 200
x <- rnorm(n, mean = 5, sd = 2)

# True parameters
beta_0 <- 10
beta_1 <- 2.5

# Skewed errors: chi-squared distribution
errors <- rchisq(n, df = 4) - 4

# Response
y <- beta_0 + beta_1 * x + errors
dat <- data.frame(y = y, x = x)

# Fit PMM2 model
fit_pmm2 <- lm_pmm2(y ~ x, data = dat, verbose = FALSE)

# Perform bootstrap inference
boot_results <- pmm2_inference(fit_pmm2,
                               formula = y ~ x,
                               data = dat,
                               B = 1000,
                               seed = 123)

# Display summary
summary(boot_results)
#>     Estimate       Std.Error            bias               t.value     
#>  Min.   :2.563   Min.   :0.07524   Min.   :-0.0028871   Min.   :22.84  
#>  1st Qu.:4.352   1st Qu.:0.16282   1st Qu.:-0.0009472   1st Qu.:25.65  
#>  Median :6.142   Median :0.25041   Median : 0.0009928   Median :28.45  
#>  Mean   :6.142   Mean   :0.25041   Mean   : 0.0009928   Mean   :28.45  
#>  3rd Qu.:7.931   3rd Qu.:0.33799   3rd Qu.: 0.0029327   3rd Qu.:31.26  
#>  Max.   :9.721   Max.   :0.42558   Max.   : 0.0048726   Max.   :34.07  
#>     p.value              conf.low       conf.high     
#>  Min.   : 0.000e+00   Min.   :2.416   Min.   : 2.710  
#>  1st Qu.:4.468e-116   1st Qu.:4.033   1st Qu.: 4.672  
#>  Median :8.936e-116   Median :5.651   Median : 6.633  
#>  Mean   :8.936e-116   Mean   :5.651   Mean   : 6.633  
#>  3rd Qu.:1.340e-115   3rd Qu.:7.269   3rd Qu.: 8.594  
#>  Max.   :1.787e-115   Max.   :8.887   Max.   :10.555

Interpretation:

Estimate: PMM2 parameter estimates from original sample
Std.Error: Bootstrap standard error (SD of bootstrap distribution)
CI.Lower / CI.Upper: 95% confidence intervals
t.value: Test statistic for H₀: β = 0
p.value: Two-sided p-value

Understanding Bootstrap Distributions

# Extract bootstrap samples for visualization
# (Note: This requires accessing internals - in practice, use summary output)

# Refit to get bootstrap samples
set.seed(123)
B <- 1000
boot_samples <- matrix(0, nrow = B, ncol = 2)

for(b in 1:B) {
  # Resample residuals
  res_b <- sample(residuals(fit_pmm2), size = n, replace = TRUE)

  # Generate bootstrap data
  y_b <- fitted(fit_pmm2) + res_b
  dat_b <- dat
  dat_b$y <- y_b

  # Refit
  fit_b <- lm_pmm2(y ~ x, data = dat_b, verbose = FALSE)
  boot_samples[b, ] <- coef(fit_b)
}

# Plot bootstrap distributions
par(mfrow = c(2, 2))

# Intercept
hist(boot_samples[, 1], breaks = 30, main = "Bootstrap: Intercept",
     xlab = "Intercept", col = "lightblue", border = "white")
abline(v = beta_0, col = "red", lwd = 2, lty = 2)
abline(v = coef(fit_pmm2)[1], col = "darkblue", lwd = 2)
legend("topright", legend = c("True", "Estimate"),
       col = c("red", "darkblue"), lty = c(2, 1), lwd = 2, cex = 0.8)

# Q-Q plot for intercept
qqnorm(boot_samples[, 1], main = "Q-Q Plot: Intercept")
qqline(boot_samples[, 1], col = "red", lwd = 2)

# Slope
hist(boot_samples[, 2], breaks = 30, main = "Bootstrap: Slope",
     xlab = "Slope", col = "lightgreen", border = "white")
abline(v = beta_1, col = "red", lwd = 2, lty = 2)
abline(v = coef(fit_pmm2)[2], col = "darkblue", lwd = 2)
legend("topright", legend = c("True", "Estimate"),
       col = c("red", "darkblue"), lty = c(2, 1), lwd = 2, cex = 0.8)

# Q-Q plot for slope
qqnorm(boot_samples[, 2], main = "Q-Q Plot: Slope")
qqline(boot_samples[, 2], col = "red", lwd = 2)


par(mfrow = c(1, 1))

Observations:

Bootstrap distributions are approximately normal (central limit theorem)
Q-Q plots show minor deviations from normality due to error skewness
Bootstrap confidence intervals account for this non-normality

Part 2: Confidence Interval Methods

The pmm2_inference() function provides multiple CI methods:

Percentile Method

The simplest approach: use empirical quantiles of the bootstrap distribution.

\[ CI_{95\%} = \left[ \hat{\beta}^*_{(0.025)}, \hat{\beta}^*_{(0.975)} \right] \]

Advantages: Simple, transformation-respecting Disadvantages: Can be biased if bootstrap distribution is skewed

Bias-Corrected (BC) Method

Adjusts for bias in the bootstrap distribution:

\[ \text{Bias} = \mathbb{E}[\hat{\beta}^*] - \hat{\beta} \]

Advantages: Corrects systematic bias Disadvantages: Requires larger B (≥ 1000)

Comparison with OLS

# Fit OLS for comparison
fit_ols <- lm(y ~ x, data = dat)

# Extract standard errors and CIs
se_ols <- summary(fit_ols)$coefficients[, "Std. Error"]
ci_ols <- confint(fit_ols)

# Bootstrap CIs (from summary)
boot_summary <- boot_results

# Create comparison table
comparison <- data.frame(
  Parameter = c("Intercept", "Slope"),
  True = c(beta_0, beta_1),
  PMM2_Estimate = coef(fit_pmm2),
  OLS_SE = se_ols,
  PMM2_Boot_SE = boot_summary$Std.Error,
  SE_Ratio = boot_summary$Std.Error / se_ols
)

print(comparison, row.names = FALSE, digits = 3)
#>  Parameter True PMM2_Estimate OLS_SE PMM2_Boot_SE SE_Ratio
#>  Intercept 10.0          9.72   0.59       0.4256    0.721
#>      Slope  2.5          2.56   0.11       0.0752    0.686

cat("\n=== Confidence Interval Comparison ===\n\n")
#> 
#> === Confidence Interval Comparison ===
cat("OLS 95% CI (Normal-based):\n")
#> OLS 95% CI (Normal-based):
print(ci_ols)
#>                2.5 %    97.5 %
#> (Intercept) 8.639023 10.967187
#> x           2.330172  2.762934

cat("\nPMM2 95% CI (Bootstrap Percentile):\n")
#> 
#> PMM2 95% CI (Bootstrap Percentile):
print(data.frame(
  Parameter = rownames(ci_ols),
  Lower = boot_summary$conf.low,
  Upper = boot_summary$conf.high
))
#>     Parameter    Lower    Upper
#> 1 (Intercept) 8.886616 10.55485
#> 2           x 2.415535  2.71046

Key Finding: PMM2 bootstrap standard errors are typically 10-30% smaller than OLS standard errors when errors are skewed, reflecting higher statistical efficiency.

Part 3: Hypothesis Testing

Bootstrap enables flexible hypothesis testing without normality assumptions.

Testing Individual Coefficients

# H0: Slope = 0 (no relationship)
# From bootstrap summary
boot_summary <- boot_results

cat("=== Hypothesis Test: Slope = 0 ===\n\n")
#> === Hypothesis Test: Slope = 0 ===
cat("t-statistic:", boot_summary$t.value[2], "\n")
#> t-statistic: 34.06547
cat("p-value:    ", boot_summary$p.value[2], "\n\n")
#> p-value:     2.395759e-254

if(boot_summary$p.value[2] < 0.05) {
  cat("Decision: Reject H0 at 5% level\n")
  cat("Conclusion: Significant relationship between x and y\n")
} else {
  cat("Decision: Fail to reject H0\n")
}
#> Decision: Reject H0 at 5% level
#> Conclusion: Significant relationship between x and y

Testing Equality of Coefficients

For multiple predictors, we might test if two slopes are equal.

# Generate multiple regression data
n <- 250
x1 <- rnorm(n)
x2 <- rnorm(n)
errors <- rexp(n, rate = 1) - 1

# True model: similar coefficients
beta <- c(5, 1.5, 1.8, -0.5)
y <- beta[1] + beta[2]*x1 + beta[3]*x2 + beta[4]*x1*x2 + errors

dat_multi <- data.frame(y = y, x1 = x1, x2 = x2)

# Fit model
fit_multi <- lm_pmm2(y ~ x1 + x2 + x1:x2, data = dat_multi, verbose = FALSE)

# Bootstrap
boot_multi <- pmm2_inference(fit_multi,
                             formula = y ~ x1 + x2 + x1:x2,
                             data = dat_multi,
                             B = 1000,
                             seed = 456)

# Test H0: beta_x1 = beta_x2
summary(boot_multi)
#>     Estimate         Std.Error            bias               t.value      
#>  Min.   :-0.5119   Min.   :0.04254   Min.   :-0.0053331   Min.   :-12.03  
#>  1st Qu.: 1.0492   1st Qu.:0.04258   1st Qu.:-0.0012402   1st Qu.: 24.63  
#>  Median : 1.6523   Median :0.04313   Median : 0.0006678   Median : 38.29  
#>  Mean   : 1.9362   Mean   :0.04838   Mean   :-0.0004654   Mean   : 35.26  
#>  3rd Qu.: 2.5392   3rd Qu.:0.04893   3rd Qu.: 0.0014426   3rd Qu.: 48.92  
#>  Max.   : 4.9521   Max.   :0.06474   Max.   : 0.0021360   Max.   : 76.49  
#>     p.value             conf.low         conf.high      
#>  Min.   :0.000e+00   Min.   :-0.5953   Min.   :-0.4286  
#>  1st Qu.:0.000e+00   1st Qu.: 0.9658   1st Qu.: 1.1327  
#>  Median :0.000e+00   Median : 1.5677   Median : 1.7368  
#>  Mean   :5.919e-34   Mean   : 1.8413   Mean   : 2.0310  
#>  3rd Qu.:5.919e-34   3rd Qu.: 2.4433   3rd Qu.: 2.6351  
#>  Max.   :2.368e-33   Max.   : 4.8252   Max.   : 5.0790

Part 4: Bootstrap for Time Series Models

Time series bootstrap requires special care due to temporal dependence.

AR Model Bootstrap

# Generate AR(1) data with skewed innovations
n <- 300
phi <- 0.7
innovations <- rgamma(n, shape = 2, rate = 2) - 1

x <- numeric(n)
x[1] <- innovations[1]
for(i in 2:n) {
  x[i] <- phi * x[i-1] + innovations[i]
}

# Fit AR(1) model
fit_ar <- ar_pmm2(x, order = 1, method = "pmm2", include.mean = FALSE)

# Bootstrap inference for time series
boot_ar <- ts_pmm2_inference(fit_ar,
                             B = 1000,
                             seed = 789)

# Summary
summary(boot_ar)
#>     Estimate        Std.Error            bias               t.value    
#>  Min.   :0.6787   Min.   :0.03428   Min.   :-2.214e-05   Min.   :19.8  
#>  1st Qu.:0.6787   1st Qu.:0.03428   1st Qu.:-2.214e-05   1st Qu.:19.8  
#>  Median :0.6787   Median :0.03428   Median :-2.214e-05   Median :19.8  
#>  Mean   :0.6787   Mean   :0.03428   Mean   :-2.214e-05   Mean   :19.8  
#>  3rd Qu.:0.6787   3rd Qu.:0.03428   3rd Qu.:-2.214e-05   3rd Qu.:19.8  
#>  Max.   :0.6787   Max.   :0.03428   Max.   :-2.214e-05   Max.   :19.8  
#>     p.value             conf.low        conf.high     
#>  Min.   :3.163e-87   Min.   :0.6115   Min.   :0.7459  
#>  1st Qu.:3.163e-87   1st Qu.:0.6115   1st Qu.:0.7459  
#>  Median :3.163e-87   Median :0.6115   Median :0.7459  
#>  Mean   :3.163e-87   Mean   :0.6115   Mean   :0.7459  
#>  3rd Qu.:3.163e-87   3rd Qu.:0.6115   3rd Qu.:0.7459  
#>  Max.   :3.163e-87   Max.   :0.6115   Max.   :0.7459

# Check if AR coefficient is significantly different from 0.5
cat("\n=== Testing H0: phi = 0.5 ===\n")
#> 
#> === Testing H0: phi = 0.5 ===
boot_summary_ar <- boot_ar
ci_lower <- boot_summary_ar$conf.low[1]
ci_upper <- boot_summary_ar$conf.high[1]

cat("95% CI: [", round(ci_lower, 3), ",", round(ci_upper, 3), "]\n")
#> 95% CI: [ 0.612 , 0.746 ]

if(0.5 >= ci_lower && 0.5 <= ci_upper) {
  cat("Decision: Cannot reject H0 (0.5 is within CI)\n")
} else {
  cat("Decision: Reject H0 (0.5 is outside CI)\n")
}
#> Decision: Reject H0 (0.5 is outside CI)

ARMA Model Bootstrap

# Generate ARMA(1,1) with heavy-tailed innovations
n <- 400
phi <- 0.5
theta <- 0.4
innovations <- rt(n, df = 4)

x <- arima.sim(n = n,
               model = list(ar = phi, ma = theta),
               innov = innovations)

# Fit ARMA model
fit_arma <- arma_pmm2(x, order = c(1, 1), method = "pmm2", include.mean = FALSE)

# Bootstrap
boot_arma <- ts_pmm2_inference(fit_arma,
                               B = 1000,
                               seed = 999,
                               method = "block",
                               block_length = 20)

# Display results
summary(boot_arma)
#>     Estimate        Std.Error            bias             t.value     
#>  Min.   :0.3020   Min.   :0.07011   Min.   :-0.04003   Min.   :3.948  
#>  1st Qu.:0.3682   1st Qu.:0.07170   1st Qu.:-0.03791   1st Qu.:4.982  
#>  Median :0.4344   Median :0.07330   Median :-0.03579   Median :6.016  
#>  Mean   :0.4344   Mean   :0.07330   Mean   :-0.03579   Mean   :6.016  
#>  3rd Qu.:0.5006   3rd Qu.:0.07490   3rd Qu.:-0.03368   3rd Qu.:7.051  
#>  Max.   :0.5668   Max.   :0.07650   Max.   :-0.03156   Max.   :8.085  
#>     p.value             conf.low        conf.high     
#>  Min.   :0.000e+00   Min.   :0.1521   Min.   :0.4519  
#>  1st Qu.:1.970e-05   1st Qu.:0.2214   1st Qu.:0.5150  
#>  Median :3.941e-05   Median :0.2907   Median :0.5781  
#>  Mean   :3.941e-05   Mean   :0.2907   Mean   :0.5781  
#>  3rd Qu.:5.911e-05   3rd Qu.:0.3601   3rd Qu.:0.6411  
#>  Max.   :7.881e-05   Max.   :0.4294   Max.   :0.7042

Key Advantage: Bootstrap inference is robust to innovation distribution and does not require Gaussian assumptions.

Part 5: Choosing Bootstrap Parameters

Number of Bootstrap Samples (B)

The choice of B affects precision and computation time:

Purpose	Recommended B	Computation Time
Quick check	200-500	Seconds
Standard inference	1000-2000	Minutes
Publication-quality	5000-10000	Hours

# Assess stability of bootstrap SEs with different B
B_values <- c(100, 500, 1000, 2000)
se_results <- data.frame(B = B_values, SE_Intercept = NA, SE_Slope = NA)

for(i in seq_along(B_values)) {
  boot_temp <- pmm2_inference(fit_pmm2, formula = y ~ x, data = dat,
                              B = B_values[i], seed = 123)
  summary_temp <- boot_temp
  se_results$SE_Intercept[i] <- summary_temp$Std.Error[1]
  se_results$SE_Slope[i] <- summary_temp$Std.Error[2]
}

print(se_results)

Rule of Thumb: Use B ≥ 1000 for confidence intervals, B ≥ 2000 for hypothesis tests.

Parallel Bootstrap

For large datasets or complex models, use parallel computing:

# Enable parallel bootstrap (requires 'parallel' package)
boot_parallel <- pmm2_inference(fit_pmm2,
                                formula = y ~ x,
                                data = dat,
                                B = 2000,
                                seed = 123,
                                parallel = TRUE,
                                cores = 4)  # Use 4 CPU cores

# Check speedup
# Typical speedup: 3-4x on 4 cores

Performance Tip: Parallel bootstrap is most beneficial when:

B ≥ 1000
Model fitting is computationally intensive (ARMA, ARIMA)
Multiple cores available

Part 6: Diagnostic Checks

Convergence of Bootstrap Estimates

Check if bootstrap replicates converged successfully:

# Rerun bootstrap with verbose output for diagnostics
boot_verbose <- pmm2_inference(fit_pmm2,
                               formula = y ~ x,
                               data = dat,
                               B = 500,
                               seed = 321)

# Check for failed bootstrap samples
# (In production code, inspect convergence attribute)
summary(boot_verbose)
#>     Estimate       Std.Error            bias               t.value     
#>  Min.   :2.563   Min.   :0.07911   Min.   :-0.0114838   Min.   :21.12  
#>  1st Qu.:4.352   1st Qu.:0.17442   1st Qu.:-0.0083642   1st Qu.:23.94  
#>  Median :6.142   Median :0.26974   Median :-0.0052445   Median :26.76  
#>  Mean   :6.142   Mean   :0.26974   Mean   :-0.0052445   Mean   :26.76  
#>  3rd Qu.:7.931   3rd Qu.:0.36505   3rd Qu.:-0.0021249   3rd Qu.:29.58  
#>  Max.   :9.721   Max.   :0.46037   Max.   : 0.0009947   Max.   :32.40  
#>     p.value             conf.low       conf.high     
#>  Min.   :0.000e+00   Min.   :2.408   Min.   : 2.718  
#>  1st Qu.:1.443e-99   1st Qu.:4.011   1st Qu.: 4.694  
#>  Median :2.885e-99   Median :5.613   Median : 6.671  
#>  Mean   :2.885e-99   Mean   :5.613   Mean   : 6.671  
#>  3rd Qu.:4.328e-99   3rd Qu.:7.216   3rd Qu.: 8.647  
#>  Max.   :5.771e-99   Max.   :8.818   Max.   :10.623

cat("\nAll bootstrap samples should converge successfully.\n")
#> 
#> All bootstrap samples should converge successfully.
cat("If many fail, consider:\n")
#> If many fail, consider:
cat("- Increasing max_iter in lm_pmm2()\n")
#> - Increasing max_iter in lm_pmm2()
cat("- Checking for outliers or influential points\n")
#> - Checking for outliers or influential points
cat("- Simplifying the model\n")
#> - Simplifying the model

Bootstrap Distribution Visualization

# For demonstration, regenerate bootstrap samples
set.seed(123)
B <- 1000
sample_size <- nrow(dat)
boot_coefs <- matrix(0, nrow = B, ncol = 2)

for(b in 1:B) {
  res_b <- sample(residuals(fit_pmm2), sample_size, replace = TRUE)
  y_b <- fitted(fit_pmm2) + res_b
  dat_b <- dat
  dat_b$y <- y_b

  fit_b <- lm_pmm2(y ~ x, data = dat_b, verbose = FALSE)
  boot_coefs[b, ] <- coef(fit_b)
}

# Scatter plot of bootstrap estimates
par(mfrow = c(1, 2))

# Bivariate distribution
plot(boot_coefs[, 1], boot_coefs[, 2],
     pch = 16, col = rgb(0, 0, 1, 0.1),
     xlab = "Intercept", ylab = "Slope",
     main = "Joint Bootstrap Distribution")
points(coef(fit_pmm2)[1], coef(fit_pmm2)[2],
       pch = 19, col = "red", cex = 2)

# Correlation
cor_boot <- cor(boot_coefs[, 1], boot_coefs[, 2])
cat("Bootstrap correlation between intercept and slope:", round(cor_boot, 3), "\n")
#> Bootstrap correlation between intercept and slope: -0.87

# Density plot
plot(density(boot_coefs[, 2]), main = "Slope Bootstrap Density",
     xlab = "Slope", lwd = 2, col = "blue")
abline(v = beta_1, col = "red", lwd = 2, lty = 2)
abline(v = coef(fit_pmm2)[2], col = "darkblue", lwd = 2)
legend("topright", legend = c("True", "Estimate", "Bootstrap"),
       col = c("red", "darkblue", "blue"),
       lty = c(2, 1, 1), lwd = 2)


par(mfrow = c(1, 1))

Part 7: Practical Guidelines

When to Use Bootstrap for PMM2

Recommended:

Always for PMM2 inference—it’s the primary method
When error distribution is clearly non-Gaussian (check Q-Q plots)
For moderate to large samples (n ≥ 100)
When you need robust confidence intervals

Alternative methods:

Analytical standard errors exist but assume asymptotic normality
May underestimate uncertainty with skewed errors

Common Pitfalls and Solutions

Problem	Symptom	Solution
Too few bootstrap samples	Wide variation in CIs across runs	Increase B to ≥ 1000
Non-convergence	Many failed bootstrap fits	Increase max_iter, check data quality
Extreme outliers	Very wide bootstrap CIs	Investigate outliers, consider robust methods
Small sample size	Unstable bootstrap distribution	Consider analytical methods, increase n if possible

Part 8: Advanced Topics

Studentized Bootstrap

For even better coverage, use studentized (pivotal) bootstrap:

# Not yet implemented in EstemPMM
# Future feature: studentized CIs for improved small-sample performance

Block Bootstrap for Time Series

For time series with strong autocorrelation, block bootstrap may be needed:

# Not yet implemented in EstemPMM
# Current approach: residual bootstrap (adequate for most cases)
# Future: moving block bootstrap for highly dependent series

Summary

Bootstrap inference for PMM2 provides:

Distribution-free confidence intervals without normality assumptions
Robust standard errors accounting for error distribution skewness
Flexible hypothesis testing for linear and time series models
Practical diagnostics to assess estimation quality

Best Practices Checklist

Use B ≥ 1000 for standard inference
Check convergence of bootstrap samples
Compare PMM2 vs OLS standard errors to quantify efficiency
Visualize bootstrap distributions to assess normality
Use parallel computing for large B or complex models
Report both percentile and bias-corrected CIs

References

Efron, B., & Tibshirani, R. J. (1994). An Introduction to the Bootstrap. CRC Press.
Davison, A. C., & Hinkley, D. V. (1997). Bootstrap Methods and Their Application. Cambridge University Press.
Zabolotnii, S., Warsza, Z.L., Tkachenko, O. (2018). “Polynomial Estimation of Linear Regression Parameters for the Asymmetric PDF of Errors”. Springer, AUTOMATION 2018.
Package functions: ?pmm2_inference, ?ts_pmm2_inference
GitHub: https://github.com/SZabolotnii/EstemPMM