在不使用 rbind 的情况下创建复制数据帧，继续早期的模拟循环问题

Question

我没有添加更多的 cmets 或延长我的原始问题，而是创建了另一个问题。我在上一个问题 (here) 中收到了很好的建议，但我在 R 方面还不够好，无法在 cmets 中实施这些建议。

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

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 添加到返回列表中，它们也会类似地连接到一个单元格中。我已经阅读了return 的R 帮助，它包含这一行

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

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

我已经尝试过创建一个空的data frame，然后尝试将结果向量化到其中。我在矢量化方面比在复制方面更差。

更新：错过了 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)返回单个值，这使得temp 对于所有RespondentID 都相同。相反，我将其替换为 pmax（有关说明，请参阅 ?max）。