Thank you both for your help!
I have also tried computing the PSE “manually” from the model’s estimates (see some messy code below), but problems come when considering the random/varying effects… here I am a bit uncertain about my implementation, and as @zacho said, it would be good to know whether there is a more programmatic way of achieving this.
Best wishes,
Ladislas
# retrieving the fixed/constant effects
fixed_effects <- fixef(trial_model)
# retrieving the random/varying effects for dyad and ppt
# remember that (1 | dyad / ppt) is equivalent to (1 | dyad) + (1 | dyad:ppt)
random_effects_dyad <- ranef(trial_model)$dyad
random_effects_ppt <- ranef(trial_model)$`dyad:ppt`
# creating a data frame with test, dyad, participant (ppt), trial
results <- df_trial %>%
select(-ppt) %>%
rename(ppt = `dyad:ppt`) %>%
distinct(test, dyad, ppt, trial) %>%
mutate(dyad = as.character(dyad) ) %>%
mutate(trial = as.numeric(trial), test = ifelse(test == "individual", -0.5, 0.5) ) %>%
rowwise() %>%
mutate(
# fixed/constant effects
intercept = fixed_effects["Intercept", "Estimate"],
test_effect = fixed_effects["test1", "Estimate"] * test,
trial_effect = fixed_effects["trial", "Estimate"] * trial,
test_trial_effect = fixed_effects["test1:trial", "Estimate"] * test * trial,
stim_effect = fixed_effects["stim", "Estimate"],
test_stim_effect = fixed_effects["test1:stim", "Estimate"] * test,
stim_trial_effect = fixed_effects["stim:trial", "Estimate"] * trial,
test_stim_trial_effect = fixed_effects["test1:stim:trial", "Estimate"] * test * trial,
# random/varying effects (dyad-level)
random_intercept_dyad = random_effects_dyad[, , "Intercept"][dyad, "Estimate"],
random_test_dyad = random_effects_dyad[, , "test1"][dyad, "Estimate"] * test,
random_trial_dyad = random_effects_dyad[, , "trial"][dyad, "Estimate"] * trial,
random_stim_dyad = random_effects_dyad[, , "stim"][dyad, "Estimate"],
random_test_stim_dyad = random_effects_dyad[, , "test1:stim"][dyad, "Estimate"] * test,
random_stim_trial_dyad = random_effects_dyad[, , "stim:trial"][dyad, "Estimate"] * trial,
random_test_stim_trial_dyad = random_effects_dyad[, , "test1:stim:trial"][dyad, "Estimate"] * test * trial,
# random-varying effects (ppt-level, nested in dyad)
random_intercept_ppt = random_effects_ppt[, , "Intercept"][ppt, "Estimate"],
random_test_ppt = random_effects_ppt[, , "test1"][ppt, "Estimate"] * test,
random_trial_ppt = random_effects_ppt[, , "trial"][ppt, "Estimate"] * trial,
random_stim_ppt = random_effects_ppt[, , "stim"][ppt, "Estimate"],
random_test_stim_ppt = random_effects_ppt[, , "test1:stim"][ppt, "Estimate"] * test,
random_stim_trial_ppt = random_effects_ppt[, , "stim:trial"][ppt, "Estimate"] * trial,
random_test_stim_trial_ppt = random_effects_ppt[, , "test1:stim:trial"][ppt, "Estimate"]* test * trial,
# combining fixed and random effects for the numerator
numerator = intercept +
test_effect + trial_effect + test_trial_effect +
random_intercept_dyad + random_test_dyad + random_trial_dyad +
random_stim_dyad + random_test_stim_dyad + random_stim_trial_dyad + random_test_stim_trial_dyad +
random_intercept_ppt + random_test_ppt + random_trial_ppt +
random_stim_ppt + random_test_stim_ppt + random_stim_trial_ppt + random_test_stim_trial_ppt,
# combining fixed and random effects for the denominator
denominator = stim_effect +
test_stim_effect + stim_trial_effect + test_stim_trial_effect +
random_stim_dyad + random_test_stim_dyad + random_stim_trial_dyad + random_test_stim_trial_dyad +
random_stim_ppt + random_test_stim_ppt + random_stim_trial_ppt + random_test_stim_trial_ppt,
# computing stim at which p(resp) = 0.5
stim_at_p_0.5 = -numerator / denominator
) %>%
ungroup()