bootstrapping of 2 metrics in a data.table by group

Question

I'm trying to boostrap 2 metrics in a data.table, by group.

library(data.table)
library(boot)
library(moments)

mystat2 <- function(data, i) {
  m <- mean(data[i])
  s <- skewness(data[i])
  list(m = m, s = s)
}


AAval_len <- 10
ABval_len <- 15
BDval_len <- 8
AAval <- rnorm(AAval_len)
ABval <- rnorm(ABval_len, mean = 2)
BDval <- rnorm(BDval_len, mean = 3)
dt <- data.table(val = c(AAval, ABval, BDval ),
                 group = c(rep("A", AAval_len),  rep("A", ABval_len), rep("B", BDval_len)),
                 subgroup = c(rep("A", AAval_len),  rep("B", ABval_len), rep("D", BDval_len)))


dt[, boot.ci(boot(val, mystat2, R = 1e3), cconf = 0.95, type = "perc"), by = list(group, subgroup)]

Expected results is a data.table similar to bellow:

expected <- data.table(group = c("A", "A", "B" ),
                       subgroup = c("A", "B", "D"),
                       level = c(0.95, 0.95, 0.95),
                       mean_percentile_low = c(-0.1, 1.9, 2.9),
                       mean_percentile_high = c(0.1, 2.1, 3.1),
                       skew_percentile_low = c(-0.6, -0.6, -0.6),
                       skew_percentile_high = c(0.6, 0.6, 0.6))

So I meet two issues.
First, the boot function has a problem with my statistic function returning 2 metrics instead of 1:

> dt[, boot.ci(boot(val, mystat2, R = 1e3), cconf = 0.95, type = "perc"), by = list(group, subgroup)]
Error in t.star[r, ] <- res[[r]] : 
  incorrect number of subscripts on matrix

Second, my data.table has a problem aggregating the results:

> dt[, boot.ci(boot(val, mystat2, R = 1e3), cconf = 0.95, type = "perc"), by = list(group, subgroup)]
Error in `[.data.table`(dt, , boot.ci(boot(val, mystat2, R = 1000), cconf = 0.95,  : 
  All items in j=list(...) should be atomic vectors or lists. If you are trying something like j=list(.SD,newcol=mean(colA)) then use := by group instead (much quicker), or cbind or merge afterwards.

Pretty sure I can handle it with a For loop but would love to have a more elegant solution...

EDIT:

I'm close but not there, with the code bellow i can look at only one statistic and don't get the columns named properly...

mystat3 <- function(data, i) {
  m <- mean(data[i])
  s <- skewness(data[i])
  list(m = m, s = s)
  m
}
dt[, { 
  tmp <-  .SD[, boot.ci(boot(val, mystat3, R = 1e3), conf = 0.95, type = "perc")]
  l <- list(level = tmp$percent[1], low <- tmp$percent[4], high <- tmp$percent[5])
  }, by = list(group, subgroup)]

user_4447540 · Answer

Here are some adjustments you can make:

Your function mystat2 should return a vector of statistics, not a list. Change as below

mystat2 <- function(data, i) {
  c(mean(data[i]), skewness(data[i]))
}

You can use the boot.ci on each statistic separately, using the index parameter, and then collect the conf level and the low and high values, and return all five in a named list. You can do all of this, including the initial bootstrap, in a function like this one:

get_boot_ci <- function(val) {
  b_out = boot(val,mystat2,R=1e3)
  ci_mean = boot.ci(b_out, type="perc", index=1)$percent[c(1,4:5)]
  ci_skew = boot.ci(b_out, type="perc", index=2)$percent[c(1,4:5)]
  list("level" = ci_mean[1],
       "mean_percentile_low" = ci_mean[2],
       "mean_percentile_high" = ci_mean[3],
       "skew_percentile_low" = ci_skew[2],
       "skew_percentile_high" = ci_skew[3]
  )
}

Now, just apply the function above to the column val, by your groups of interest

dt[, get_boot_ci(val), by=.(group, subgroup)]

Output:

   group subgroup level mean_percentile_low mean_percentile_high skew_percentile_low skew_percentile_high
1:     A        A  0.95          -0.7117178            0.2915394          -1.5501161            1.4951826
2:     A        B  0.95           1.4666828            2.4285060          -2.2816065            0.5033563
3:     B        D  0.95           2.4914297            3.3688430          -0.5834545            1.6190134

bootstrapping of 2 metrics in a data.table by group

1 Answers