Fit glm function

library(tidyverse)
library(rlang)
library(broom)
library(quickpsy) # to use qpdat

fit_glm <- function(.data, .response, .predictor) {
  
  predictor <- ensym(.predictor)

  prop <- .data |>
    group_by({{ .predictor }}) |>
    summarise(k = sum({{.response}}), 
              n = n(),
              prop = mean({{ .response }}), .groups = "keep")

  model <- glm(inject(cbind(k, n - k) ~ !!predictor),
                      family = binomial(logit),
                      data = prop)

  param <- tidy(model) |>
    select(term, estimate) |>
    pivot_wider(names_from = term, values_from = estimate) |>
    rename(intercept = `(Intercept)`, sensitivity = {{ .predictor }}) |>
    mutate(bias = -intercept / sensitivity)

  min_max <- .data |>
    distinct({{.predictor}}) |>
    summarise(min = min({{ .predictor }}),
              max = max({{ .predictor }}))

  sequ <- tibble({{.predictor}} := seq(min_max$min,
                                       min_max$max,
                                       length.out = 100))

  psy <- augment(model, newdata = sequ, type.predict = "response") |>
    rename(prop = .fitted)
  
  deviance <- glance(model) |> 
    select(null.deviance, df.null) |> 
    mutate(p_value = pchisq(null.deviance, df = df.null, lower.tail = FALSE)) |>
    rename(null_deviance = null.deviance, df_null = df.null)

  list(prop = prop,
       param = param,
       psy = psy, 
       deviance = deviance)
}

Single group case

fit_all <- fit_glm(qpdat, resp, phase)

ggplot() +
  geom_point(data = fit_all$prop,
             aes(x = phase, y = prop)) +
  geom_line(data = fit_all$psy, 
            aes(x = phase, y = prop)) +
  geom_segment(data = fit_all$param, 
               aes(x = bias, xend = bias, 
                   y = 0, yend = .5))

Multiple groups

fits <- qpdat |> 
  group_by(participant, cond) |> 
  reframe(fit = list(fit_glm(cur_data(), resp, phase))) |> 
  unnest_wider(fit)

proportions <- fits |> 
  select(participant, cond, prop) |> 
  unnest(prop)

psy_functions <- fits |> 
  select(participant, cond, psy) |> 
  unnest(psy)

params <- fits |> 
  select(participant, cond, param) |> 
  unnest(param)

ggplot() +
  facet_wrap(vars(participant)) +
  geom_point(data = proportions,
             aes(x = phase, y = prop, color = cond)) +
  geom_line(data = psy_functions, 
            aes(x = phase, y = prop, color = cond)) +
  geom_segment(data = params, 
               aes(x = bias, xend = bias, 
                   y = 0, yend = .5, color = cond))

All the fits are different from a fit with a flat slope

 fits |> 
  select(participant, cond, deviance) |> 
  unnest(deviance)