R - итеративно применять функцию списка переменных

Question

Моя цель - создать функцию, которая при циклическом переходе по нескольким переменным кадра данных вернет новый фрейм данных, содержащий проценты и доверительные интервалы 95% для каждого уровня каждой переменной.

В качестве примера, если я применил эту функцию «цил» и «утро» из кадра mtcars данных, я хочу это как конечный результат:

variable level    ci.95 
1  cyl  4 34.38 (19.50, 53.11) 
2  cyl  6 21.88 (10.35, 40.45) 
3  cyl  8 43.75 (27.10, 61.94) 
4  am  0 59.38 (40.94, 75.5) 
5  am  1 40.62 (24.50, 59.06) 

Итак, пока у меня есть функция, похоже, работает для одной переменной; тем не менее, у меня есть два вопроса, которые я надеюсь, что сообщество может мне помочь:

1. Генеральный R-исчисления моего кода. Я все еще новичок R. Я читал достаточно сообщений, чтобы знать, что энтузиасты R в целом препятствуют использованию циклов for, но я по-прежнему действительно борюсь с использованием применяемых функций (что в большинстве случаев является альтернативой петлям for).

2. Применение этой функции к списку переменных - в результате создается один кадр данных, содержащий возвращаемые значения из функции для каждого уровня каждой переменной.

Вот где я нахожусь с моим кодом до сих пор:

t1.props <- function(x, data = NULL) { 

    # Grab dataframe and/or variable name 
    if(!missing(data)){ 
    var <- data[,deparse(substitute(x))] 
    } else { 
    var <- x 
    } 

    # Grab variable name for use in ouput 
    var.name <- substitute(x) 

    # Omit observations with missing data 
    var.clean <- na.omit(var) 

    # Number of nonmissing observations 
    n <- length(var.clean) 

    # Grab levels of variable 
    levels <- sort(unique(var.clean)) 

    # Create an empty data frame to store values 
    out <- data.frame(variable = NA, 
        level = NA, 
        ci.95 = NA) 

    # Estimate prop, se, and ci for each level of the variable 
    for(i in seq_along(levels)) { 
    prop <- paste0("prop", i) 
    se <- paste0("se", i) 
    log.prop <- paste0("log.trans", i) 
    log.se <- paste0("log.se", i) 
    log.l <- paste0("log.l", i) 
    log.u <- paste0("log.u", i) 
    lcl <- paste0("lcl", i) 
    ucl <- paste0("ucl", i) 

    # Find the proportion for each level of the variable 
    assign(prop, sum(var.clean == levels[i])/n) 

    # Find the standard error for each level of the variable 
    assign(se, sd(var.clean == levels[i])/
      sqrt(length(var.clean == levels[i]))) 

    # Perform a logit transformation of the original percentage estimate 
    assign(log.prop, log(get(prop)) - log(1 - get(prop))) 

    # Transform the standard error of the percentage to a standard error of its 
    # logit transformation 
    assign(log.se, get(se)/(get(prop) * (1 - get(prop)))) 

    # Calculate the lower and upper confidence bounds of the logit 
    # transformation 
    assign(log.l, 
      get(log.prop) - 
      qt(.975, (length(var.clean == levels[i]) - 1)) * get(log.se)) 
    assign(log.u, 
      get(log.prop) + 
      qt(.975, (length(var.clean == levels[i]) - 1)) * get(log.se)) 

    # Finally, perform inverse logit transformations to get the confidence bounds 
    assign(lcl, exp(get(log.l))/(1 + exp(get(log.l)))) 
    assign(ucl, exp(get(log.u))/(1 + exp(get(log.u)))) 

    # Create a combined 95% CI variable for easy copy/paste into Word tables 
    ci.95 <- paste0(round(get(prop) * 100, 2), " ", 
       "(", sprintf("%.2f", round(get(lcl) * 100, 2)), ",", " ", 
       round(get(ucl) * 100, 2), ")") 

    # Populate the "out" data frame with values 
    out <- rbind(out, c(as.character(var.name), levels[i], ci.95)) 
    } 

    # Remove first (empty) row from out 
    # But only in the first iteration 
    if (is.na(out[1,1])) { 
    out <- out[-1, ] 
    rownames(out) <- 1:nrow(out) 
    } 
    out 
} 

data(mtcars) 
t1.props(cyl, mtcars) 

Я признателен за любую помощь или совет вы должны предложить.

Answer 1

Вы также можете сохранить функцию в основном нетронутыми и использовать lapply над ним:

vars <- c("cyl", "am") 
lapply(vars, t1.props, data=mtcars) 
[[1]] 
    variable level    ci.95 
1  cyl  4 34.38 (19.50, 53.11) 
2  cyl  6 21.88 (10.35, 40.45) 
3  cyl  8 43.75 (27.10, 61.94) 

[[2]] 
    variable level    ci.95 
1  am  0 59.38 (40.94, 75.5) 
2  am  1 40.62 (24.50, 59.06) 

И объединить их все в один кадр данных с:

lst <- lapply(vars, t1.props, data=mtcars) 
do.call(rbind,lst) 

данных

You должен упростить задания var и var.name:

t1.props <- function(x, data = NULL) { 

    # Grab dataframe and/or variable name 
    if(!missing(data)){ 
    var <- data[,x] 
    } else { 
    var <- x 
    } 

    # Grab variable name for use in ouput 
    var.name <- x 

    # Omit observations with missing data 
    var.clean <- na.omit(var) 

    # Number of nonmissing observations 
    n <- length(var.clean) 

    # Grab levels of variable 
    levels <- sort(unique(var.clean)) 

    # Create an empty data frame to store values 
    out <- data.frame(variable = NA, 
        level = NA, 
        ci.95 = NA) 

    # Estimate prop, se, and ci for each level of the variable 
    for(i in seq_along(levels)) { 
    prop <- paste0("prop", i) 
    se <- paste0("se", i) 
    log.prop <- paste0("log.trans", i) 
    log.se <- paste0("log.se", i) 
    log.l <- paste0("log.l", i) 
    log.u <- paste0("log.u", i) 
    lcl <- paste0("lcl", i) 
    ucl <- paste0("ucl", i) 

    # Find the proportion for each level of the variable 
    assign(prop, sum(var.clean == levels[i])/n) 

    # Find the standard error for each level of the variable 
    assign(se, sd(var.clean == levels[i])/
      sqrt(length(var.clean == levels[i]))) 

    # Perform a logit transformation of the original percentage estimate 
    assign(log.prop, log(get(prop)) - log(1 - get(prop))) 

    # Transform the standard error of the percentage to a standard error of its 
    # logit transformation 
    assign(log.se, get(se)/(get(prop) * (1 - get(prop)))) 

    # Calculate the lower and upper confidence bounds of the logit 
    # transformation 
    assign(log.l, 
      get(log.prop) - 
      qt(.975, (length(var.clean == levels[i]) - 1)) * get(log.se)) 
    assign(log.u, 
      get(log.prop) + 
      qt(.975, (length(var.clean == levels[i]) - 1)) * get(log.se)) 

    # Finally, perform inverse logit transformations to get the confidence bounds 
    assign(lcl, exp(get(log.l))/(1 + exp(get(log.l)))) 
    assign(ucl, exp(get(log.u))/(1 + exp(get(log.u)))) 

    # Create a combined 95% CI variable for easy copy/paste into Word tables 
    ci.95 <- paste0(round(get(prop) * 100, 2), " ", 
        "(", sprintf("%.2f", round(get(lcl) * 100, 2)), ",", " ", 
        round(get(ucl) * 100, 2), ")") 

    # Populate the "out" data frame with values 
    out <- rbind(out, c(as.character(var.name), levels[i], ci.95)) 
    } 

    # Remove first (empty) row from out 
    # But only in the first iteration 
    if (is.na(out[1,1])) { 
    out <- out[-1, ] 
    rownames(out) <- 1:nrow(out) 
    } 
    out 
} 

Answer 2

Хорошая вещь обо всех функциях, которые вы используете, состоит в том, что они уже векторизованы (кроме sd и qt, но вы можете легко их векторизовать для конкретных аргументов с помощью Vectorize). Это означает, что вы можете передавать им векторы без необходимости писать один цикл. Я оставил части вашей функции, которые занимаются подготовкой ввода и получением результата.

t1.props <- function(var, data=mtcars) { 
    N <- nrow(data) 
    levels <- names(table(data[,var])) 
    count <- unclass(table(data[,var]))  # counts 
    prop <- count/N       # proportions 
    se <- sqrt(prop * (1-prop)/(N-1))   # standard errors of props. 
    lprop <- log(prop) - log(1-prop)   # logged prop 
    lse <- se/(prop*(1-prop))    # logged se 
    stat <- Vectorize(qt, "df")(0.975, N-1) # tstats 
    llower <- lprop - stat*lse     # log lower 
    lupper <- lprop + stat*lse     # log upper 
    lower <- exp(llower)/(1 + exp(llower)) # lower ci 
    upper <- exp(lupper)/(1 + exp(lupper)) # upper ci 

    data.frame(variable=var, 
       level=levels, 
       perc=100*prop, 
       lower=100*lower, 
       upper=100*upper) 
} 

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

## Apply your function to two variables 
do.call(rbind, lapply(c("cyl", "am"), t1.props)) 
# variable level perc lower upper 
# 4  cyl  4 34.375 19.49961 53.11130 
# 6  cyl  6 21.875 10.34883 40.44691 
# 8  cyl  8 43.750 27.09672 61.94211 
# 0  am  0 59.375 40.94225 75.49765 
# 1  am  1 40.625 24.50235 59.05775 

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

Я думаю, что самое важное, что нужно изменить в вашем коде, это использование assign и get. Вместо этого вы можете хранить переменные в списках или другую структуру данных и использовать setNames, names<- или names(...) <-, чтобы назвать компоненты, когда это необходимо.