将 splitstackshape 合并到循环中

Incorporating splitstackshape into loop

我有以下代码选择(4 行虹膜 x 1000)*100 并计算每列的偏差。

library(SimDesign)
library(data.table)

do.call(rbind,lapply(1:100, function(x) {
  bias(
    setDT(copy(iris))[as.vector(sapply(1:1000, function(X) sample(1:nrow(iris),4)))][
      , lapply(.SD, mean), by=rep(c(1:1000),4), .SDcols=c(1:4)][,c(2:5)],
    parameter=c(5,3,2,1), #parameter is the true population value used to calculate bias
    type='relative' #denotes the type of bias being calculated 
  )
}))

这需要 4 行的 1000 个样本,通过样本 # 计算平均值,给我 1000 个平均值。为每一列找到 1000 均值的偏差,然后再进行 99 次,为我提供每一列的偏差估计分布。这是在模仿随机抽样设计。但是,我也想为分层设计这样做。所以我使用splitstackshapestratified函数。

do.call(rbind,lapply(1:100, function(x) {
  bias(
    setDT(copy(iris))[as.vector(sapply(1:1000, function(X) stratified(iris,group="Species", size=1)))][
      , lapply(.SD, mean), by=rep(c(1:1000),4), .SDcols=c(1:4)][,c(2:5)],
    parameter=c(5,3,2,1), 
    type='relative'
  )
}))

我原以为这只是换出函数的问题,但我一直收到错误 (i is invalid type (matrix))。也许将来 2 列矩阵可以 return DT 的元素列表。我认为这可能与 setDT 有关,但我不确定如何修复它。有人知道我哪里错了吗?

我已经为你分成了几个函数

加载 data.table、SimDesign 和 splitstackshape

library(SimDesign)
library(data.table)
library(splitstackshape)

获取 n 大小为 sampsize 的分层样本和这些样本的 return 列均值的函数

get_samples <- function(n, sampsize=4) {
  rbindlist(lapply(1:n, function(x) {
    splitstackshape::stratified(iris, group="Species",sampsize)[, id:=x]
  }))[, lapply(.SD, mean), by=.(Species, id)]
}

现在,让我们得到这些样本的 y 此类迭代中的偏差分布

get_bias_distribution <- function(y=100, samples_per_iter=50, size_per_iter=4) {
  rbindlist(lapply(1:y, function(y) {
    samples = get_samples(samples_per_iter, sampsize=size_per_iter)[, id:=NULL]
    samples[, as.list(bias(
      estimate=.SD,parameter=c(5,3,2,1),type="relative")*100),
      by=.(Species)][, iter:=y]  
  }))
}

用法(使用默认值)

get_bias_distribution()    

输出:

        Species Sepal.Length Sepal.Width Petal.Length Petal.Width iter
  1:     setosa    -1.236667    22.61833    -26.70000   -39.69667    1
  2: versicolor    46.476667   -11.99500    115.12833    16.82167    1
  3:  virginica    80.596667    -0.20000    180.21833    53.89000    1
  4:     setosa    -1.513333    20.87000    -27.46167   -38.83667    2
  5: versicolor    45.333333   -11.34833    112.84833    17.84500    2
 ---                                                                  
296: versicolor    48.250000   -12.26833    113.37000    17.71167   99
297:  virginica    77.366667    -2.87000    175.60000    53.07167   99
298:     setosa    -1.005000    22.67500    -27.02833   -39.69500  100
299: versicolor    47.921667   -10.28333    110.97833    16.86833  100
300:  virginica    76.153333    -2.44000    174.46167    52.62167  100

关于上面问题的一些评论

  1. 当您调用 stratified(iris,group="Species", size=1) 时,您将获得第 3 行 data.table,因为您实际上是从三个物种中的每一个中随机选择一行
   Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1:          4.9         3.6          1.4         0.1     setosa
2:          6.3         2.5          4.9         1.5 versicolor
3:          7.7         2.8          6.7         2.0  virginica
  1. 当你将它包装在 sapply(1:1000, function(x)...) 中时,你会得到 5 x 1000 列矩阵,其中每列包含 5 个长度为 3 的列表。下面,我将向你展示如果你这样做的话会是什么样子sapply(1:6, function(x)...)
             [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     
Sepal.Length numeric,3 numeric,3 numeric,3 numeric,3 numeric,3 numeric,3
Sepal.Width  numeric,3 numeric,3 numeric,3 numeric,3 numeric,3 numeric,3
Petal.Length numeric,3 numeric,3 numeric,3 numeric,3 numeric,3 numeric,3
Petal.Width  numeric,3 numeric,3 numeric,3 numeric,3 numeric,3 numeric,3
Species      factor,3  factor,3  factor,3  factor,3  factor,3  factor,3

这并不是您真正想要的,因为您无法 lapply 按照您当时的预期方式处理这些内容。您要做的是使用 lapply(1:1000, function(x) ...) 创建此类 3 行数据表的列表,然后将它们绑定在一起(在为每个数据表添加 id 列之后)。