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))
}
2
  • 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, 2019 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, 2019 at 22:06

2 Answers 2

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
5
  • 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, 2019 at 22:17
  • @Biotechgeek You can convert any tibble to dataframe using as.data.frame().
    – Shree
    Aug 20, 2019 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, 2019 at 22:35
  • Should there be a sd(growth, na.rm = T)^2
    – Rspacer
    Aug 20, 2019 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, 2019 at 22:43
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.)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.