mardi 27 novembre 2018

if-statement in a function using purrr and dplyr (List Column Workflow) in R.

I am trying to estimate differences in means for two hospitals, A and B. Every hospital has different "groups" and I have given them group 1 and 2 in the simulated data set. That is I want to test difference in means between hospital A and B within group 1 and group 2 and in addition I have more than one variable (e.g. value1 and value2). So I have to test value1 between hospital A and B within the groups 1 and 2. My function works fine when I exclude the if-statements in the function body, but when I try to include "method" in the function (to choose how the confidence interval is computed) something goes wrong. I have included the function which works at the end. I am using the infer package for the bootstraping (part of tidyverse or tidymodels).

library(tidyverse)
library(infer)
library(stringr)
library(rlang)

set.seed(1)
A <-data.frame(value1=rnorm(n = 1000, mean = 0.8, sd = 0.2), value2= rnorm(n=10 ,mean=1, sd=0.3)) 
A$hosp <- "A"
A$group <- sample(1:2,nrow(A) , replace=T) 

B= data.frame(value1 = rnorm(n=1200, mean =1 , sd = 0.2), value2= rnorm(n=15, mean=1.1, sd=0.4))
B$hosp <- "B"
B$group <- sample(1:2,nrow(B) , replace=T) 

forskel <- bind_rows(A, B) %>% 
  group_by(group) %>% 
  nest()

rm(A, B)

#function that DOES NOT work
f <- function(dataset, procestid, method=1, reps = 4, alpha = 0.05) { 
  procestid <- enquo(procestid)

  diff_mean <- dataset %>% 
    mutate(diff_means  = map(data, function(.x){.x %>% 
        group_by(hosp) %>% 
        summarise(mean(!!procestid, na.rm=TRUE)) %>% 
        pull() %>% 
        diff() })) %>%
    select(-data) 


  bootstrap <- dataset %>% 
    mutate(distribution =map(data, function(.x){ .x %>% 
        specify(as.formula(paste0(quo_name(procestid), "~ hosp")) ) %>% 
        generate(reps = reps, type = "bootstrap") %>%
        calculate(stat = "diff in means", order = c( "A", "B"))} )) %>% 
    inner_join(diff_mean, by="group")  

  if (method=1) {
    bootstrap2 <- bootstrap %>% mutate(Bias_Corrected_KI=map2(distribution, diff_means, function(.x, .y){ .x %>% 
        summarise( l =quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(alpha/2))),
                   u= quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(1-alpha/2)))    )}))  }
  if (method=2) {
    bootstrap2 <- bootstrap %>% mutate(Percentile_KI = map(distribution, function(.x){.x %>% 
        summarize(l = quantile(stat, alpha/2),
                  u = quantile(stat, 1 - alpha/2))}))  }
  else {
    bootstrap2 <- bootstrap %>% mutate(SD_KI =map2(distribution, diff_means, function(.x,.y){.x %>% 
        get_confidence_interval(level = (1 - alpha), type="se", point_estimate = .y)})) 
  }

}

procestimes <- list("value1", "value2")

The function bellow WORKS, but I would like to choose the bootstraping method with an if statements.

f <- function(dataset, procestid, reps = 4, alpha = 0.05) { 
  procestid <- enquo(procestid)

  diff_mean <- dataset %>% 
    mutate(diff_means  = map(data, function(.x){.x %>% 
        group_by(hosp) %>% 
        summarise(mean(!!procestid, na.rm=TRUE)) %>% 
        pull() %>% 
        diff() })) %>%
    select(-data) 


  bootstrap <- dataset %>% 
    mutate(distribution =map(data, function(.x){ .x %>% 
        specify(as.formula(paste0(quo_name(procestid), "~ hosp")) ) %>% 
        generate(reps = reps, type = "bootstrap") %>%
        calculate(stat = "diff in means", order = c( "A", "B"))} )) %>% 
    inner_join(diff_mean, by="group")  %>% 


   mutate(Bias_Corrected_KI=map2(distribution, diff_means, function(.x, .y){ .x %>% 
        summarise( l =quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(alpha/2))),
                   u= quantile(.x$stat,pnorm(2*qnorm(sum(.x$stat >= .y)/reps) + qnorm(1-alpha/2)))    )}))  %>% 

   mutate(Percentile_KI = map(distribution, function(.x){.x %>% 
        summarize(l = quantile(stat, alpha/2),
                  u = quantile(stat, 1 - alpha/2))}))  %>% 

    mutate(SD_KI =map2(distribution, diff_means, function(.x,.y){.x %>% 
        get_confidence_interval(level = (1 - alpha), type="se", point_estimate = .y)})) 

}

procestimes <- list("value1", "value2")


a <- map(syms(procestimes), f , dataset=forskel, reps=1000)



a[[1]]$Bias_Corrected_KI  

Aucun commentaire:

Enregistrer un commentaire