用于 plogis 的基本 R 代码的 Tidyverse 版本？

Question

问题

我正在将一些使用基本 R 的 R 代码转换为使用 tidyverse 的代码，但我有点卡住了。这是我正在尝试转换的代码的相关 sn-p：

J <- c(2, 3, 7, 3, 50) # male or not, eth, age, income level, state
  poststrat <- as.data.frame(array(NA, c(prod(J), length(J)+1))) # Columns of post-strat matrix, plus one for size
  colnames(poststrat) <- c("male", "eth", "age","income", "state",'N')
  count <- 0
  for (i1 in 1:J[1]){
    for (i2 in 1:J[2]){
      for (i3 in 1:J[3]){
        for (i4 in 1:J[4]){
          for (i5 in 1:J[5]){
            count <- count + 1
            # Fill them in so we know what category we are referring to
            poststrat[count, 1:5] <- c(i1-1, i2, i3,i4,i5) 
          }
        }
      }
    }
  }
  
  coef_male <- c(0,-0.3)
  coef_eth <- c(0, 0.6, 0.9)
  coef_age <- c(0,-0.2,-0.3, 0.4, 0.5, 0.7, 0.8)
  coef_income <- c(0,-0.2, 0.6)
  coef_state <- c(0, round(rnorm(49, 0, 1), 1))
  coef_age_male <- t(cbind(c(0, .1, .23, .3, .43, .5, .6),
                           c(0, -.1, -.23, -.5, -.43, -.5, -.6)))
  true_popn <- data.frame(poststrat[, 1:5], cat_pref = rep(NA, prod(J)))
  for (j in 1:prod(J)) {
    true_popn$cat_pref[j] <- plogis(
      coef_male[poststrat[j, 1] + 1] +
        coef_eth[poststrat[j, 2]] + coef_age[poststrat[j, 3]] +
        coef_income[poststrat[j, 4]] + coef_state[poststrat[j, 5]] +
        coef_age_male[poststrat[j, 1] + 1, poststrat[j, 3]]
    )
  }

这正是我正在努力解决的代码部分：

  true_popn <- data.frame(poststrat[, 1:5], cat_pref = rep(NA, prod(J)))
  for (j in 1:prod(J)) {
    true_popn$cat_pref[j] <- plogis(
      coef_male[poststrat[j, 1] + 1] +
        coef_eth[poststrat[j, 2]] + coef_age[poststrat[j, 3]] +
        coef_income[poststrat[j, 4]] + coef_state[poststrat[j, 5]] +
        coef_age_male[poststrat[j, 1] + 1, poststrat[j, 3]]
    )
  }

我知道这段代码的作用——它遍历所有分配的系数，并通过plogis 函数将对数赔率转换为概率。我不清楚的是如何使用 tidyverse 做到这一点。

这是我尝试过的方法

我的直觉是使用 expand_grid 函数，然后将 plogis 函数应用于行总和：

  coef_male <- c(0,-0.3)
  coef_eth <- c(0, 0.6, 0.9)
  coef_age <- c(0,-0.2,-0.3, 0.4, 0.5, 0.7, 0.8)
  coef_income <- c(0,-0.2, 0.6)
  coef_state <- c(0, round(rnorm(49, 0, 1), 1))
  coef_age_male <- t(cbind(c(0, .1, .23, .3, .43, .5, .6),
                           c(0, -.1, -.23, -.5, -.43, -.5, -.6)))
  expand_grid(coef_male,
              coef_eth,
              coef_age,
              coef_income,
              coef_state)

此代码有效，因为我可以获得我想要的 6,300 个单元格，但它无法合并我想要的 coef_age_male 变量。使用带有此变量的扩展网格会生成一个大于我需要的 6,300 的网格。

我的问题

有人可以帮我用 plogis 函数将循环替换为 dplyr 吗？我认为它会更容易被同事理解和阅读。

Answer 1

下面是不包含原始代码的系数data.frame（例如male：0、1）。顺序不同，但如有必要，我们可以使用不同的顺序，arrange。

library(tidyverse)

coef_male <- c(0,-0.3)
coef_eth <- c(0, 0.6, 0.9)
coef_age <- c(0,-0.2,-0.3, 0.4, 0.5, 0.7, 0.8)
coef_income <- c(0,-0.2, 0.6)
coef_state <- c(0, round(rnorm(49, 0, 1), 1))
coef_age_male <- t(cbind(c(0, .1, .23, .3, .43, .5, .6),
                         c(0, -.1, -.23, -.5, -.43, -.5, -.6)))
# lets turn the matrix above into a vector
coef_age_male2 = c(0, .1, .23, .3, .43, .5, .6, 
                   0, -.1, -.23, -.5, -.43, -.5, -.6) 
# now we can create a look-up data.frame
look_up <- crossing(male = coef_male,
                    age = coef_age) %>% 
  arrange(desc(male)) %>% 
  mutate(age_male = coef_age_male2)

# create a matrix of coefficients and join together
cross_df(list(
         male = coef_male, 
         eth = coef_eth,
         age = coef_age,
         income = coef_income,
         state = coef_state)) %>% 
  left_join(look_up, by = c("male", "age")) %>% 
  mutate(cat_pref = select(.,male:age_male) %>% 
                      rowSums() %>%
                      plogis()
         )
#> # A tibble: 6,300 x 7
#>     male   eth   age income state age_male cat_pref
#>    <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
#>  1   0     0     0        0     0     0.23    0.557
#>  2  -0.3   0     0        0     0    -0.23    0.371
#>  3   0     0.6   0        0     0     0.23    0.696
#>  4  -0.3   0.6   0        0     0    -0.23    0.517
#>  5   0     0.9   0        0     0     0.23    0.756
#>  6  -0.3   0.9   0        0     0    -0.23    0.591
#>  7   0     0    -0.2      0     0     0.1     0.475
#>  8  -0.3   0    -0.2      0     0    -0.1     0.354
#>  9   0     0.6  -0.2      0     0     0.1     0.622
#> 10  -0.3   0.6  -0.2      0     0    -0.1     0.5  
#> # … with 6,290 more rows

^{由reprex package (v0.3.0) 于 2021-01-16 创建}