## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(tempodisco)

## -----------------------------------------------------------------------------
data("td_bc_study")
# For speed, only look at the first 10 participants
ids <- unique(td_bc_study$id)[1:10]
study <- subset(td_bc_study, id %in% ids)

## -----------------------------------------------------------------------------
# Check for nonsystematic discounting
is_nonsys <- by(study, study$id, function(ptpt) { # Stratify by participant ID
  mod <- td_bcnm(ptpt, discount_function = 'model-free') # Fit a "model-free" discount function to the participant-level data
  any(nonsys(mod)) # Check whether either non-systematic criterion is met
}, simplify = T)
# Keep data from systematic discounters
keep_ids <- rownames(is_nonsys)[!is_nonsys]
study <- subset(study, id %in% keep_ids)

## -----------------------------------------------------------------------------
rows <- by(study, study$id, function(ptpt) { # Stratify by participant ID
  # Test several discount functions and use the best-fitting one
  mod <- td_bcnm(ptpt, discount_function = c('hyperbolic', 'exponential'))
  # Return a one-row dataframe containing model-agnostic measures of discounting
  data.frame(
    id = ptpt$id[1],
    AUC = AUC(mod),
    ED50 = ED50(mod)
  )
}, simplify = F)
discount_measures <- do.call(rbind, rows)
cor.test(discount_measures$AUC, log(discount_measures$ED50))

## -----------------------------------------------------------------------------
# Compare Kirby scoring to nonlinear modeling
k_vals <- by(study, study$id, function(ptpt) { # Stratify by participant ID
  # Get k from Kirby MCQ scoring
  mod <- kirby_score(ptpt)
  k_kirby <- coef(mod)['k']
  # Get k from nonlinear model fitting
  mod <- td_bcnm(ptpt, discount_function = 'hyperbolic')
  k_nm <- coef(mod)['k']
  # Return a one-row dataframe containing both
  data.frame(
    k_nm = k_nm,
    k_kirby = k_kirby
  )
}, simplify = F)
# Stack all the rows to produce a single dataframe
df <- do.call(rbind, k_vals)
# What's the correlation?
cor.test(log(df$k_nm), log(df$k_kirby))

## -----------------------------------------------------------------------------
best_mods <- by(study, study$id, function(ptpt) {
  # Test both the exponential and hyperbolic discount functions
  mod <- td_bcnm(ptpt, discount_function = c('hyperbolic', 'exponential'))
  # Return the name of the best-fitting function
  mod$config$discount_function$name
}, simplify = T)
table(best_mods)