0

我有以下代码选择 (4 rows of iris 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 次，为我提供每列的偏差估计分布。这是模仿随机抽样设计。但是，我也想为分层设计执行此操作。所以我使用splitstackshape'sstratified函数。

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 列矩阵可以返回 DT 的元素列表。我认为这可能与 setDT 有关，但我不确定如何修复它。有人知道我哪里出错了吗？

4

1 回答 1

2

我已经为你分成了几个功能

加载 data.table、SimDesign 和 splitstackshape

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

获取n大小分层样本sampsize并返回这些样本的列均值的函数

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，因为您有效地从三个 Species 中的每一个中随机选择一行
   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为每个数据表添加一列之后）。

于 2022-02-24T01:17:24.313 回答