Создание реплицированного фрейма данных без использования rbind, следуйте предыдущей проблеме цикла моделирования - PullRequest
Новая
       5

Создание реплицированного фрейма данных без использования rbind, следуйте предыдущей проблеме цикла моделирования

0 голосов
/ 26 января 2012

Вместо того, чтобы добавлять больше комментариев или увеличивать исходный вопрос, я создал другой вопрос. Я получил отличный совет в предыдущем вопросе ( здесь ), но я недостаточно хорош в 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))
    }
}

Ответ replicate() работал хорошо, когда я думал, что мне нужен только один вывод (mc_amount) из функции: 

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 для тестирования, который составляет всего head() из моего большего набора данных из 7135 наблюдений. Набор 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), но есть две проблемы с результатами. Эти проблемы могут быть связаны, я не понимаю достаточно, чтобы понять, что происходит.

Во-первых, я получаю только 10 столбцов, относящихся к первому RespondentID, поэтому эта функция, похоже, не применяется к 6 из набора данных.

Во-вторых, я получаю 10 столбцов, но результаты RespondentID объединяются в одну ячейку в каждом столбце. Если я добавлю u2 или mc_bca в список возврата, они также будут аналогичным образом объединены в одну ячейку. Я прочитал справку R для return, и она содержит эту строку

значение может быть серией непустых выражений, разделенных запятыми. В этом случае возвращаемое значение представляет собой список оцененных выражений, имена которых установлены в выражения, где это имена объектов R. но я недостаточно разбираюсь в программировании функций R, чтобы понять, насколько это актуально.

Я надеюсь, что есть быстрое и очевидное решение этой проблемы. Мне не удалось найти аналогичную проблему, решение которой я мог бы скопировать, во всех примерах множественных возвратов из функций, которые я обнаружил, использовались переменные, вычисляемые в функции.

Я попробовал альтернативу создания пустого data frame и затем попытался векторизовать результаты в это. Я хуже в векторизации, чем в репликации.

Обновление: пропущено значение min_bca, которое составляет -2,44478269434376

1 Ответ

1 голос
/ 26 января 2012

После еще нескольких правок, надеюсь, вот окончательное решение вашего вопроса. 

    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 для объяснения).

...