HSS: Analysis of Education Expenditure data

library(HSS)
#> Loading required package: lmtest
#> Loading required package: zoo
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
#> Loading required package: BiasedUrn
#> Loading required package: MASS
#> Loading required package: glmnet
#> Loading required package: Matrix
#> Loaded glmnet 4.1-2
#> Loading required package: foreach
#> Loading required package: doParallel
#> Loading required package: iterators
#> Loading required package: parallel
#> Loading required package: itertools

  library(robustbase)

  data(education, package="robustbase")

  k = 3
  nsplit <- 100
  B  <- 99
  p = 3




  EE <- education[-49,] 
  N  <- dim(EE)[1]
  rho = 24 / 49

  X <- EE[,3:5]
  reg <- lm(Y ~ X1+X2+X3, data=EE)
  resid <- reg$residuals

  inclusion_probs <- inclusion_probs_BS <- rep(0,p)
  results_BS <- matrix(NA, nrow=B, ncol=p)

  for (split in 1:nsplit){
    obj <- HSS:::get_pvalue_matrix(X, resid, rho=rho, k=k)
    probs <- obj$probs
    n.ref <- obj$n_ref
    n.ev <- obj$n_ev
    ref <- obj$ref
    ev <- obj$ev

    ## Perform HSS given a p-value matrix
    hrv_subset <- do_hss(probs, n.ref, n.ev,
    ref, ev, k,
    use_alphahat=TRUE,
    alphagrid=FALSE,
    how_compute_expecs = "numeric",
    q=2,
    weight=FALSE)$subset

    inclusion_probs[hrv_subset] <- inclusion_probs[hrv_subset] + 1
  }
  inclusion_probs <- inclusion_probs / nsplit


  for (bb in 1:B){
    inclusion_probs_BS <- rep(0, p)
    residB <- sample(resid, N, replace=TRUE)
    for (split in 1:nsplit){
            obj <- HSS:::get_pvalue_matrix(X, residB, rho=rho, k)
      probs <- obj$probs
      n.ref <- obj$n_ref
      n.ev <- obj$n_ev
      ref <- obj$ref
      ev <- obj$ev

      hrv_subset <- do_hss(probs, n.ref, n.ev,
      ref, ev, k,
      use_alphahat=T,
      alphagrid=F,
      how_compute_expecs="numeric",
      q=2,
      weight=F)$subset
      inclusion_probs_BS[hrv_subset] <- inclusion_probs_BS[hrv_subset] + 1
    }
    inclusion_probs_BS <- inclusion_probs_BS / nsplit
    results_BS[bb, ] <- inclusion_probs_BS
  }

  ecdf_matrix <- apply(results_BS, 2, ecdf)
  (inclusion_p_value <- mapply(function(x, x_ecdf){1-x_ecdf(x)}, inclusion_probs, ecdf_matrix))
#> [1] 0.9797980 0.0000000 0.8686869
  (corrected_inclusion_matrix <- p.adjust(inclusion_p_value, "holm")) 
#> [1] 1 0 1