R/analysis_psychometrics.R

In calenwalshe/detectability_maps: What the Package Does (One Line, Title Case)

#' Estimate the psychometric functions
#'
#' @param subject_list
#'
#' @return
#' @export
#' @import ggplot2 tidyr purrr parallel dplyr bbmle
#'
#' @examples
estimate_psychometric <- function(subject_list = c("anqi", "rcw", "can", "arw")) {
  uncertainty_type <- paste0("umap_6cpd_", "Drasdo", "_11_23_20")

  subject_data <- do.call(rbind, map(subject_list, humandetectiondata::format_subject_data))

  all.subjects <- subject_data

  subjects <- unique(all.subjects$SUBJECT)
  estim.params <- lapply(subjects, FUN = function(x) {
    subject_id <- x

    print(x)
    all.subjects.sub <- all.subjects %>% filter(SUBJECT == subject_id)
    split.df         <- split(all.subjects.sub, all.subjects.sub$type)
    df               <- split.df

    #load(paste0("~/Dropbox/Calen/Work/search/detection_experiments/analysis/", uncertainty_type, ".Rdata"))
    #umap6cpdmatrix <<- humandetectiondata::umap6cpdmatrix
    mle.fun        <- {e <- new.env(); e$umap6cpdmatrix <- umap6cpdmatrix; humandetectiondata::f.search.mle(df)}



    testsearch <- mle2(mle.fun,
                       start = list(p_kSp = 1, p_k0 = 1, gamma = 0),
                       method = "L-BFGS-B", fixed = list(p_Beta = 1))

    all.subjects.sub$k_sp  <- testsearch@fullcoef[1]
    all.subjects.sub$k_0   <- testsearch@fullcoef[2]
    all.subjects.sub$beta  <- testsearch@fullcoef[3]
    all.subjects.sub$gamma <- testsearch@fullcoef[4]

    return(all.subjects.sub)
  })

  parameter.df <- do.call(rbind, estim.params)

  #parameter.df <- parameter.df %>% filter(eccentricity == 3) # REMOVE
  save.params <- parameter.df %>%
    group_by(SUBJECT) %>%
    dplyr::select(SUBJECT, k_sp, k_0, beta, gamma) %>%
    unique()

  c_Amp <- parameter.df$tAmp
  c_x   <- parameter.df$X
  c_y   <- parameter.df$Y
  c_deg <- parameter.df$eccentricity
  c_vf  <- parameter.df$condition

  dprime.predictions <- purrr::pmap(parameter.df,
                                    function(tAmp, X, Y, eccentricity, condition,k_sp, k_0, beta,...) {
                                      uncertainty <- humandetectiondata::f.dprime.uni(tAmp, X, Y, beta, k_sp, k_0)
                                      no_uncertainty <- humandetectiondata::f.dprime.no.uni(tAmp, X, Y, beta, k_sp, k_0)

                                      return(data.frame(dprime_type = c("uncertainty", "no_uncertainty"),
                                                        dprime.model = c(uncertainty, no_uncertainty)))
                                    })


  parameter.df$dprime.prediction <- dprime.predictions

  parameter.df <- parameter.df %>%
    unnest(dprime.prediction) %>%
    dplyr::select(-trial)

  # Plotting
  performance.df <- parameter.df %>%
    group_by(tAmp, SUBJECT, eccentricity, experiment_name, dprime_type, condition) %>%
    summarize(hit_r = mean(type == "HIT") / mean(tPresent == 1),
              fa_r = mean(type == "FA") / mean(tPresent == 0),
              n_trials = n(),
              hit_r = ifelse(hit_r == 1, (n_trials - 1) / n_trials, hit_r),
              fa_r = ifelse(fa_r == 0, 1 / n_trials, fa_r),
              dprime.human = qnorm(hit_r) - qnorm(fa_r),
              pc = pnorm(dprime.human/2),
              pc.hat = mean(pnorm(dprime.model/2)),
              dprime.hat = mean(dprime.model))

  performance.df <- performance.df %>%
    group_by(dprime_type) %>%
    gather("measure.type",
           "value",
           c("pc", "pc.hat", "dprime.human", "dprime.hat")) %>%
    unique()

  #fovea
  plot.fov <- ggplot(performance.df %>%
                       filter(measure.type %in% c("pc", "pc.hat"), experiment_name == "fovea"),
                     aes(x = tAmp, y = value, shape = as.factor(SUBJECT), colour = as.factor(measure.type), linetype = dprime_type)) +
    geom_point() +
    geom_line() +
    facet_grid(~condition, scale = "free") +
    expand_limits(y = c(.5, 1)) +
    theme(aspect.ratio = 1)

  #periphery
  plot.periph <- ggplot(performance.df %>%
                          filter(measure.type %in% c("pc", "pc.hat"), experiment_name == "periphery"),
                        aes(x = eccentricity, y = value, shape = as.factor(SUBJECT), colour = as.factor(measure.type), linetype = dprime_type)) +
    geom_point() +
    geom_line() +
    expand_limits(y = c(.5, 1)) +
    facet_grid(~condition, scale = "free") +
    theme(aspect.ratio = 1)

  return(list(periphery_plot = plot.periph, fovea_plot = plot.fov, parameters = save.params, performance = performance.df))
}