为分层抽样准备抽样分布的最佳 R 包函数

Question

我正在尝试在 R 中准备一个演示，说明小群体的重复分层随机抽样如何导致均值接近正态抽样分布。例如，考虑下面的 R 代码（它可以工作，但由于循环而非常慢）。

#Dummy population made up of dice throws - 18 per row
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P1 <- as.data.frame(c(5,6,5,1,6,4,2,2,4,4,6,6,5,2,3,5,1,6))
P1$Zn <- 1
names(P1) <- c('Die','Zn')
Dt <- P1

P2 <- as.data.frame(c(2,5,4,5,5,5,3,3,2,5,6,1,2,5,4,3,6,1))
P2$Zn <- 2
names(P2) <- c('Die','Zn')
Dt <- rbind(Dt,P2)

# Empty dataframe to hold random draws
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Smps <- data.frame(Die = numeric(), Zn= numeric(),Drw = numeric())

# Draw stratifed samples one from each row
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
print(paste('Start','at',Sys.time()))
n <- 10000          # number of draws
r <- 2              # number of rows (the strata)
for (j in 1:n){
  # for a 2 strata
  for (i in 1:r){
    #sub set strata
    x <- subset(Dt, Dt$Zn == i)
    # random sample
    y <- x[sample(1:18,1),]
    y$Drw <- j
    #append sample
    Smps <- rbind(Smps,y)
  }
  # report progress
  if(right(j,3) == '000'){
    print(paste(j,'at',Sys.time()))
    flush.console()
  }
}

# Compute the sample means
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Mns <-aggregate(Smps[, 1], list(Smps$Drw), mean)

# Density plot of means
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
d <- density(Mns$x)
plot(d,xlab = 'Means', las=1, main = '')
polygon(d, col="blue", border="blue")

我期待有一个 R 包具有执行这种类型的分层抽样的功能，但我正在努力寻找以我能理解的方式工作的包。输入带有分组字段的数据帧和要从每个组中抽取的样本数的东西是我所期望的，它已经被编写为允许一组重复采样。任何指向有效示例的指针都将不胜感激。理想情况下，我想准备说来自具有更多层的已知人群的 100,000 个分层样本，然后绘制均值的分布（但很快）

Answer 1

在解决这个问题一段时间后，我发现了一个名为“fifer”（https://www.rdocumentation.org/packages/fifer/versions/1.1）的包，它似乎在包中包含一个分层函数，但不幸的是，这个包不适用于最新版本的 R。我做到了但是，请从 Ananda Mahto (https://gist.github.com/mrdwab/6424112) 中找到一个巧妙的分层函数，该函数运行良好，但代价是脚本中的函数相当长，而不是加载包的一行。我使用此功能解决上述问题的方法如下。

#Dummy population made up of dice throws - 18 per row
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P1 <- as.data.frame(c(5,6,5,1,6,4,2,2,4,4,6,6,5,2,3,5,1,6))
P1$Zn <- 1
names(P1) <- c('Die','Zn')
Dt <- P1

P2 <- as.data.frame(c(2,5,4,5,5,5,3,3,2,5,6,1,2,5,4,3,6,1))
P2$Zn <- 2
names(P2) <- c('Die','Zn')
Dt <- rbind(Dt,P2)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Stratfed function from web
# https://gist.github.com/mrdwab/6424112
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

stratified <- function(df, group, size, select = NULL, 
                       replace = FALSE, bothSets = FALSE) {
  if (is.null(select)) {
    df <- df
  } else {
    if (is.null(names(select))) stop("'select' must be a named list")
    if (!all(names(select) %in% names(df)))
      stop("Please verify your 'select' argument")
    temp <- sapply(names(select),
                   function(x) df[[x]] %in% select[[x]])
    df <- df[rowSums(temp) == length(select), ]
  }
  df.interaction <- interaction(df[group], drop = TRUE)
  df.table <- table(df.interaction)
  df.split <- split(df, df.interaction)
  if (length(size) > 1) {
    if (length(size) != length(df.split))
      stop("Number of groups is ", length(df.split),
           " but number of sizes supplied is ", length(size))
    if (is.null(names(size))) {
      n <- setNames(size, names(df.split))
      message(sQuote("size"), " vector entered as:\n\nsize = structure(c(",
              paste(n, collapse = ", "), "),\n.Names = c(",
              paste(shQuote(names(n)), collapse = ", "), ")) \n\n")
    } else {
      ifelse(all(names(size) %in% names(df.split)),
             n <- size[names(df.split)],
             stop("Named vector supplied with names ",
                  paste(names(size), collapse = ", "),
                  "\n but the names for the group levels are ",
                  paste(names(df.split), collapse = ", ")))
    }
  } else if (size < 1) {
    n <- round(df.table * size, digits = 0)
  } else if (size >= 1) {
    if (all(df.table >= size) || isTRUE(replace)) {
      n <- setNames(rep(size, length.out = length(df.split)),
                    names(df.split))
    } else {
      message(
        "Some groups\n---",
        paste(names(df.table[df.table < size]), collapse = ", "),
        "---\ncontain fewer observations",
        " than desired number of samples.\n",
        "All observations have been returned from those groups.")
      n <- c(sapply(df.table[df.table >= size], function(x) x = size),
             df.table[df.table < size])
    }
  }
  temp <- lapply(
    names(df.split),
    function(x) df.split[[x]][sample(df.table[x],
                                     n[x], replace = replace), ])
  set1 <- do.call("rbind", temp)

  if (isTRUE(bothSets)) {
    set2 <- df[!rownames(df) %in% rownames(set1), ]
    list(SET1 = set1, SET2 = set2)
  } else {
    set1
  }
}


# Empty dataframe to hold random draws
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Smps <- data.frame(Die = numeric(), Zn = numeric())

# Right function for reporting progress
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
right = function(text, num_char) {
  substr(text, nchar(text) - (num_char-1), nchar(text))
}

# Draw stratifed samples one from each row
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n <- 10000          # number of draws
for (j in 1:n){
    y <- stratified(Dt,"Zn",1)
    y <- cbind(y,j)
    Smps <- rbind(Smps,y)
  # report progress
  if(right(j,3) == '000'){
    print(paste(j,'at',Sys.time()))
    flush.console()
  }
}

# Compute the sample means
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Mns <-aggregate(Smps[, 1], list(Smps$j), mean)

# Density plot of means
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
d <- density(Mns$x)
plot(d,xlab = 'Means', las=1, main = '')
polygon(d, col="blue", border="blue")