R/stats_tests.R

In RobWHickman/RobeRt:

# Stats Fun
#
# Functions to help do some stats

#adds in the rolling mean, se and sd for a variable
rowwise_basic_stats <- function(dt, group){
  task_data <- task_data %>%
    .[is.na(task_failure), correct_trials := 1:.N, by = c(group, "block_no")] %>%
    .[is.na(task_failure), rolling_mean := rollapplyr(monkey_final_bid, correct_trials, mean, align = "right"), by = c("offer_value", "date")] %>%
    .[is.na(task_failure), rolling_sd := rollapplyr(monkey_final_bid, correct_trials, sd, align = "right"), by = c("offer_value", "date")] %>%
    .[, rolling_se := rolling_sd / sqrt(correct_trials)]

  return(task_data)
}

#carries out a rowwise welch test (i.e. all trials before the current trial) given a dt
rowwise_welchy <- function(trial, dt, group, group_value, direction){

  #subset by group
  dt_group <- dt[which(dt[[group]] == group_value),] %>%
    .[is.na(task_failure)]

  #how many hypotheses are we testing
  hypotheses <- length(unique(dt$juice_offered))

  #find the most recent trial for each juice offer
  low <- dt_group$monkey_final_bid[which(dt_group$trial <= trial &
                                           dt_group$juice_offered == min(dt_group$juice_offered))]
  mid <- dt_group$monkey_final_bid[which(dt_group$trial <= trial &
                                           dt_group$juice_offered != min(dt_group$juice_offered) &
                                           dt_group$juice_offered != max(dt_group$juice_offered))]
  high <- dt_group$monkey_final_bid[which(dt_group$trial <= trial &
                                            dt_group$juice_offered == max(dt_group$juice_offered))]

  #do a t test with unequal variances and take the p value
  if(length(low) > 1 & length(mid) > 1){
    welch_test_lowmid <- t.test(low, mid, alternative = direction, var.equal = FALSE, conf.level = 1 - (0.05 / hypotheses))$p.value
  } else {welch_test_lowmid <- NA}
  if(length(low) > 1 & length(high) > 1){
    welch_test_lowhigh <- t.test(low, high, alternative = direction, var.equal = FALSE, conf.level = 1 - (0.05 / hypotheses))$p.value
  } else {welch_test_lowhigh <- NA}
  if(length(mid) > 1 & length(high) > 1){
    welch_test_midhigh <- t.test(mid, high, alternative = direction, var.equal = FALSE, conf.level = 1 - (0.05 / hypotheses))$p.value
  } else {welch_test_midhigh <- NA}

  #get the difference between this p value and the bonferroni corrected p level
  lowmid_sig <- welch_test_lowmid - (0.05 / hypotheses)
  lowhigh_sig <- welch_test_lowhigh - (0.05 / hypotheses)
  midhigh_sig <- welch_test_midhigh - (0.05 / hypotheses)

  p_values <- data.frame(block = group_value,
                         trial = trial,
                         lowmid = lowmid_sig,
                         lowhigh = lowhigh_sig,
                         midhigh = midhigh_sig)
  return(p_values)
}