1

I have the following data-set (dat) with 8 unique treatment groups. I want to sample 3 points from each unique group and store their mean and variance. I want to do this 1000 times over (sample with replacement) using a loop to store all the values in output. I tried to do this loop and I keep running into unexpected '=' in:"output[i] <- summarise(group_by(new_df[i], fertilizer,crop, level),mean[i]="

Any suggestions on how to fix it, or make it more

fertilizer <- c("N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P","N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P")

crop <- c("alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group","alone","group")

level <- c("low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","low")

growth <- c(0,0,1,2,90,5,2,5,8,55,1,90,2,4,66,80,1,90,2,33,56,70,99,100,66,80,1,90,2,33,0,0,1,2,90,5,2,2,5,8,55,1,90,2,4,66,0,0)

dat <- data.frame(fertilizer, crop, level, growth)

library(dplyr)

for(i in 1:1000){
  new_df[i] <- dat %>% 
                  group_by(fertilizer, crop, level) %>% 
                  sample_n(3)
  output[i] <- summarise(
                  group_by(new_df[i], fertilizer, crop, level),
                  mean[i] = mean(growth), 
                  var[i] = sd(growth) * sd(growth))
}
Rspacer
  • 2,369
  • 1
  • 14
  • 40
  • `summarize(..., mean[i]=...)` is bad for a few reasons: (1) `summarize` doesn't take indexed assignment like that (while it works with simple vectors on the REPL); (2) my opinion, naming a variable the same as a (common) function is perhaps bad form, but that just my two cents. Mostly the first. – r2evans Aug 20 '19 at 22:05
  • 1
    I fixed the typo in your code, please check the code you give us before asking the question. – r2evans Aug 20 '19 at 22:06

2 Answers2

3

Try this:

replicate(2, {
  dat %>%
    group_by(fertlizer, crop, level) %>%
    sample_n(3) %>%
    summarize(mu = mean(growth), sigma2 = sd(growth)^2) %>%
    ungroup()
}, simplify = FALSE)
# [[1]]
# # A tibble: 8 x 5
#   fertlizer crop  level    mu  sigma2
#   <fct>     <fct> <fct> <dbl>   <dbl>
# 1 N         alone high   1       1   
# 2 N         alone low   30.7  2641.  
# 3 N         group high  33.3  2408.  
# 4 N         group low   56     553   
# 5 P         alone high  22.7  1409.  
# 6 P         alone low    2.33    2.33
# 7 P         group high  40.3  1336.  
# 8 P         group low   23    1387   
# [[2]]
# # A tibble: 8 x 5
#   fertlizer crop  level    mu sigma2
#   <fct>     <fct> <fct> <dbl>  <dbl>
# 1 N         alone high   30.3  2670.
# 2 N         alone low    52.7  2069.
# 3 N         group high   61.7  2408.
# 4 N         group low    20     925 
# 5 P         alone high   35.3  3042.
# 6 P         alone low    19.7   990.
# 7 P         group high   14.3   270.
# 8 P         group low    32    2524.  

(Replace 2 with your 1000.)

r2evans
  • 141,215
  • 6
  • 77
  • 149
3

I don't think you need a loop. You can do this faster by sampling 3*1000 values per group at once, assign sample_id and add it to grouping variables, and finaly summarize to get desired values. This way you are calling all functions only once. -

dat %>% 
  group_by(fertilizer, crop, level) %>% 
  sample_n(3*1000, replace = T) %>% 
  mutate(sample_id = rep(1:1000, each = 3)) %>% 
  group_by(sample_id, add = TRUE) %>% 
  summarise(
    mean = mean(growth, na.rm = T),
    var = sd(growth)^2
  ) %>% 
  ungroup()

# A tibble: 8,000 x 6
   fertilizer crop  level sample_id  mean      var
   <chr>      <chr> <chr>     <int> <dbl>    <dbl>
 1 N          alone high          1 30.7  2640.   
 2 N          alone high          2  1       0    
 3 N          alone high          3 60.3  2640.   
 4 N          alone high          4  1.33    0.333
 5 N          alone high          5  1.33    0.333
 6 N          alone high          6 60.3  2640.   
 7 N          alone high          7  1.33    0.333
 8 N          alone high          8 30.3  2670.   
 9 N          alone high          9  1.33    0.333
10 N          alone high         10 60.7  2581.   
# ... with 7,990 more rows
Shree
  • 10,835
  • 1
  • 14
  • 36
  • Is there a way to get the final table as a data frame than a tibble? Tibbles are harder to maniuplate later – Rspacer Aug 20 '19 at 22:17
  • @Biotechgeek You can convert any tibble to dataframe using `as.data.frame()`. – Shree Aug 20 '19 at 22:19
  • When I try applying the same concept to my original data, I get an error with `Error in summarise_impl(.data, dots) : Evaluation error: dims [product 8000] do not match the length of object [1].` Why do you think? – Rspacer Aug 20 '19 at 22:35
  • Should there be a `sd(growth, na.rm = T)^2` – Rspacer Aug 20 '19 at 22:41
  • Hard to say but seems like whatever summary functions you are using are not producing one value per group. Yes, `sd(growth, na.rm = T)^2` is better. – Shree Aug 20 '19 at 22:43