不使用 rbind 创建复制数据帧，继续之前的模拟循环问题

Question

我没有添加更多评论或使原来的问题变得更长，而是创建了另一个问题。我在上一个问题中收到了很好的建议（这里），但我在 R 方面还不够好，无法实现评论中的建议。

花了很长时间的原始代码是：

Male.MC <-c()
for (j in 1:100)            {
    for (i in 1:nrow(Male.Distrib))  {
        u2        <- Male.Distrib$stddev_u2[i] * rnorm(1, mean = 0, sd = 1)
        mc_bca    <- Male.Distrib$FixedEff[i] + u2
        temp      <- Lambda.Value*mc_bca+1
        ginv_a    <- temp^(1/Lambda.Value)
        d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
        mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
        z <- data.frame(
        RespondentID = Male.Distrib$RespondentID[i], 
        Subgroup     = Male.Distrib$Subgroup[i], 
        mc_amount    = mc_amount,
        IndvWeight   = Male.Distrib$INDWTS[i]/100
        )
        Male.MC <- as.data.frame(rbind(Male.MC,z))
    }
}

当我认为我只需要函数的一个输出（mc_amount）时，replicate()答案效果很好：

Male.Distrib = read.table('MaleDistrib.txt', check.names=F)

getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778) {
      u2        <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1)
      mc_bca    <- df$FixedEff + u2
      temp      <- Lambda.Value*mc_bca+1
      ginv_a    <- temp^(1/Lambda.Value)
      d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
      mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
      mc_amount
}
replicate(10, getMC(Male.Distrib))

但是，即使进行数据更正后，我得到了意想不到的结果，因此我需要能够查看所有临时计算的值，以确定我的逻辑哪里出了问题。这就是我被困住的地方。我创建了一个名为 tempdata 的较小数据框进行测试，它只是来自 7135 个观测值的较大数据集的 head()。 tempdata 集是：

    RndmEff RespondentID Subgroup RespondentID Replicates IntakeAmt RACE INDWTS    TOTWTS   GRPWTS NUMSUBJECTS TOTSUBJECTS  FixedEff stddev_u2
1  1.343753         9966        6         9966      41067 33.449808    2  41067 120622201 41657878        1466        7135  6.089918  2.645938
2 -5.856516         9967        5         9967       2322  2.533528    3   2322 120622201 22715139        1100        7135  6.755664  2.645938
3 -3.648339         9970        4         9970      17434  9.575439    2  17434 120622201 10520535        1424        7135  7.079757  2.645938
4  2.697533         9972        6         9972      21723 43.340180    2  21723 120622201 41657878        1466        7135  6.089918  2.645938
5  3.531878         9974        3         9974        375 55.660607    3    375 120622201 10791729        1061        7135  6.176319  2.645938
6  6.627767         9976        6         9976      48889 91.480049    2  48889 120622201 41657878        1466        7135  6.089918  2.645938

我使用的更新命令是：

getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778) {
    RespondentID <- df$RespondentID
    u2        <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1)
    mc_bca    <- df$FixedEff + u2
    temp      <- max(Lambda.Value*mc_bca+1,Lambda.Value*Min_bca+1)
    ginv_a    <- temp^(1/Lambda.Value)
    d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
    mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
    return(list(RespondentID, temp, ginv_a, d2ginv_a, mc_amount))
}
Test <- replicate(10, getMC(tempdata))

我为计算变量获得了非常好的布局（temp、ginv_a、 d2ginv_a, mc_amount)，但结果有两个问题。这些问题可能是相关的，我不太了解，无法弄清楚发生了什么。

首先，我只获得与第一个 RespondentID 相关的 10 列，因此该函数似乎不适用于数据集中的 6 列。

其次，我得到 10 列，但 RespondentID 结果连接到每一列的一个单元格中。如果我将 u2 或 mc_bca 添加到返回列表中，它们也会类似地连接到一个单元格中。我已阅读 R 的 return 帮助，它包含这一行

value 可以是一系列用逗号分隔的非空表达式。在这种情况下，返回的值是已计算表达式的列表，名称设置为表达式，其中这些是 R 对象的名称。 但我对 R 函数编程了解不够，不知道这是否相关。

我希望有一个快速且明显的解决方案。我一直无法找到可以复制解决方案的类似问题，我发现的所有函数多次返回的示例都使用了在函数中计算的变量。

我尝试了另一种方法：创建一个空的数据框，然后尝试将结果矢量化为该数据框。我在向量化方面比在复制方面更糟糕。

更新：错过了 min_bca 值，即 -2.44478269434376

Answer 1

经过几次编辑后，希望这是您问题的最终解决方案。

    getMC <- function(df, Lambda.Value=0.4, Male.Resid.Var=12.1029420429778,Min_bca=-2.44478269434376) {
            u2        <- df$stddev_u2 * rnorm(nrow(df), mean = 0, sd = 1)
            mc_bca    <- df$FixedEff + u2
            temp      <- pmax((Lambda.Value*mc_bca+1),(Lambda.Value*Min_bca+1))
            ginv_a    <- temp^(1/Lambda.Value)
            d2ginv_a  <- max(0,(1-Lambda.Value)*temp^(1/Lambda.Value-2))
            mc_amount <- ginv_a + d2ginv_a * Male.Resid.Var / 2
            return(data.frame(RespondentID=df$RespondentID,temp=temp, ginv_a, d2ginv_a, mc_amount))
        }

   data=rep(list(tempdata),10) # change 10 to a higher number of replicates
   result_data=llply(data,getMC, .progress = "text")

一些注意事项：我必须逐行对单个复制中的函数进行故障排除，以找出问题所在（这是您在发布之前应该做的事情，因为上面的问题与此问题无关）。 max(vector1,vector2)返回一个值，该值使得所有RespondentID的temp相同。相反，我用 pmax 替换了它（有关解释，请参阅 ?max）。