Finding threshold for adaptive paradigm via psychometric function

Interfaces brms
#1

I’m on R version 6.0 and using version brms 2.9.0

I have a handful of questions, but I’ll start simple. I’m trying to create a psychometric function to extract threshold values. Reading “Fitting the Psychometric Function” (Bernhard Treutwein): https://link.springer.com/article/10.3758/BF03211951, my understanding of psychometric functions is that fit is enhanced by specifying priors of lapse, guess, threshold, and spread parameters.

I’m trying to fit psychometric functions for an adaptive paradigm of eight tasks, with two conditions (except for one task which has four conditions). The adaptivity component (and independent variable) is response window (increases with incorrect answer, decreases with correct answer). There are roughly 1000 participants (3rd, 5th, and 7th graders). First and foremost, I want to be able to extract a threshold response window for each participant at probability = .8.

I found this link to a posting on this site, which has been helpful: https://discourse.mc-stan.org/c/interfaces/brms.

I’m new to brms but have tried to familiarize myself with the terms, looking through the Cran.r page: https://cran.r-project.org/web/packages/brms/brms.pdf.

I can’t seem to implement threshold or spread parameters in this model. I get the following error:

Error: The parameter 'threshold' is not a valid distributional or non-linear parameter. Did you forget to set 'nl = TRUE'?

So that’s my first issue. My second issue is that when I include more than one participant I get this error for the four chains:

Stan model '547eb5715049d22db6f829f32c7c94cd' does not contain samples.

Here’s my model. I’m not sure what I’m doing wrong. I’ve included a dput of my data below. I’ve kept it to only one condition of one task for simplicity on this site (the dataset is a flanker task and comprises only 11 participants), but I would ideally like to calculate this 80% threshold for all conditions of every task. Is this better done in separate models, or all in the same go?

BF <- bf(
  response~guess + (1-guess-lapse) * Phi(eta),
  eta ~ 1 + pid/rw, #what is eta?,
  #threshold~1,
  #spread ~1,
  guess ~ 1,
  lapse ~ 1, 
  family = bernoulli(link="identity"),
  nl = TRUE
)

priors <- c(
  prior(beta(7, 3),class = "b", nlpar = "eta"),
  #prior(beta(7, 3),class = "b", nlpar = "threshold"),
  #prior(beta(1.4,1.4), nlpar = "spread"),
  prior(beta(.5, 8), nlpar = "lapse", lb = 0, ub = .1),
  prior(beta(1, 5), nlpar = "guess", lb = 0, ub = .1)
)

fit <- brm(
  BF,
  data = flanker.psych,
  inits = 'random', #what is this for?
  control = list(adapt_delta = 0.99), #what is this for?
  prior = priors,
  iter = 1000,
    chains = 4,
    cores = 4
)


Thanks!


Dataset:

structure(list(pid = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L), .Label = c("ADMIN-UCSF-bo001", 
"ADMIN-UCSF-bo002", "ADMIN-UCSF-bo004", "ADMIN-UCSF-bo005", "ADMIN-UCSF-bo008", 
"ADMIN-UCSF-bo009", "ADMIN-UCSF-bo010", "ADMIN-UCSF-bo011", 
"ADMIN-UCSF-se104", "ADMIN-UCSF-se105"), class = "factor"), module = structure(c(5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L
), .Label = c("BACKWARDSSPATIALSPAN", "BOXED", "BRT", "FILTER", 
"FLANKER", "ISHIHARA", "SAAT", "SPATIALSPAN", "STROOP", "TASKSWITCH", 
"TNT"), class = "factor"), condition = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("CONGRUENT", 
"Conjunction_12", "Conjunction_4", "Feature_12", "Feature_4", 
"impulsive", "INCONGRUENT", "Left", "R2B0", "R2B2", "R2B4", "R4B0", 
"R4B2", "Right", "Start", "Stay", "sustained", "Switch", "Tap & Trace", 
"Tap Only"), class = "factor"), correct_button = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 
3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), .Label = c("correct", "incorrect", "no_response"), class = "factor"), 
    rw = c(710L, 800L, 820L, 690L, 710L, 760L, 700L, 840L, 680L, 
    780L, 770L, 640L, 720L, 780L, 840L, 830L, 790L, 830L, 760L, 
    700L, 670L, 800L, 740L, 670L, 800L, 600L, 840L, 760L, 920L, 
    910L, 940L, 830L, 900L, 840L, 870L, 320L, 920L, 950L, 820L, 
    440L, 600L, 910L, 790L, 930L, 800L, 880L, 910L, 360L, 800L, 
    280L, 1120L, 1190L, 1030L, 1110L, 1320L, 1250L, 1150L, 1180L, 
    950L, 1030L, 1190L, 1170L, 1180L, 870L, 910L, 1170L, 1240L, 
    960L, 800L, 840L, 880L, 1110L, 1150L, 1210L, 1280L, 950L, 
    890L, 860L, 800L, 870L, 1070L, 820L, 1030L, 1110L, 880L, 
    790L, 810L, 960L, 830L, 1090L, 1100L, 950L, 850L, 780L, 790L, 
    800L, 990L, 880L, 1030L, 840L, 900L, 880L, 1460L, 1480L, 
    910L, 960L, 920L, 940L, 920L, 1440L, 950L, 1470L, 1460L, 
    880L, 880L, 1470L, 1480L, 880L, 840L, 840L, 1000L, 1440L, 
    1160L, 920L, 800L, 1150L, 840L, 2260L, 2250L, 1130L, 1140L, 
    1070L, 1110L, 990L, 2270L, 1150L, 2190L, 2230L, 2110L, 800L, 
    830L, 870L, 910L, 990L, 1030L, 1310L, 1070L, 1630L, 1950L, 
    1990L, 840L, 810L, 850L, 1070L, 830L, 800L, 990L, 980L, 1000L, 
    1060L, 1000L, 960L, 1050L, 990L, 1210L, 770L, 810L, 840L, 
    880L, 920L, 1030L, 1090L, 1200L, 1130L, 1240L, 950L, 790L, 
    790L, 790L, 850L, 940L, 880L, 910L, 780L, 800L, 870L, 820L, 
    670L, 800L, 760L, 630L, 830L, 930L, 710L, 870L, 790L, 840L, 
    810L, 780L, 870L, 830L, 820L, 640L, 840L, 1600L, 360L, 1610L, 
    800L, 880L, 480L, 1620L, 320L, 940L, 960L, 1630L, 920L, 950L, 
    960L, 800L, 800L, 1590L, 400L, 880L, 640L, 840L), rt = c(468.641996383667, 
    489.506006240845, 498.21001291275, 500.387012958527, 513.499021530151, 
    537.150025367737, 567.131042480469, 578.942000865936, 585.861027240753, 
    619.597971439362, 649.703025817871, 651.482999324799, 684.765994548798, 
    699.003994464874, 709.982991218567, 711.242020130157, 721.700012683868, 
    748.438954353333, 901.057004928589, 349.790036678314, 450.733006000519, 
    462.338984012604, 466.818988323212, 950, 950, 487.961947917938, 
    525.457978248596, 586.38596534729, 634.458959102631, 670.267999172211, 
    676.239967346191, 707.340002059937, 711.510002613068, 721.628963947296, 
    724.897027015686, 757.191002368927, 766.102969646454, 795.396983623505, 
    818.078994750977, 848.93000125885, 849.749028682709, 875.111997127533, 
    875.764012336731, 928.068041801453, 951.33101940155, 1017.45498180389, 
    1048.14499616623, 734.741985797882, 800, 1010, 550.595998764038, 
    842.209994792938, 872.821033000946, 900.366961956024, 949.076056480408, 
    975.140988826752, 981.266975402832, 1001.76501274109, 1033.71900320053, 
    1047.33997583389, 1056.77700042725, 1086.42899990082, 1137.69000768661, 
    1149.69998598099, 1214.39599990845, 1260.75398921967, 1376.6970038414, 
    857.559978961945, 840, 870, 990, 1140, 1150, 1400, 1400, 
    342.560946941376, 376.044988632202, 486.966013908386, 488.722026348114, 
    500.590980052948, 585.26599407196, 588.424026966095, 612.955987453461, 
    633.490025997162, 676.449000835419, 689.705014228821, 736.205995082855, 
    786.639034748077, 877.53701210022, 916.054010391235, 943.838000297546, 
    1002.13998556137, 360.725998878479, 398.616015911102, 526.05801820755, 
    585.74903011322, 633.624017238617, 687.473058700562, 702.538013458252, 
    940, 533.69802236557, 584.937989711761, 595.023989677429, 
    595.512986183167, 615.786015987396, 655.584037303925, 697.839021682739, 
    735.269010066986, 768.783986568451, 811.460018157959, 824.137032032013, 
    829.093992710114, 860.583007335663, 915.856957435608, 936.079978942871, 
    1011.21997833252, 1043.23101043701, 73.8150477409363, 483.029007911682, 
    535.443961620331, 699.177026748657, 861.621022224426, 968.819975852966, 
    1001.51002407074, 800, 650.12800693512, 766.891002655029, 
    813.524007797241, 861.293971538544, 946.185052394867, 948.246002197266, 
    1016.33298397064, 1095.98195552826, 1134.80800390244, 1147.3600268364, 
    1160.73602437973, 1278.76400947571, 1411.3729596138, 1446.06500864029, 
    191.150009632111, 1120, 1120, 1140, 1260, 1350, 1350, 1350, 
    1630, 2270, 2270, 590.786993503571, 604.610025882721, 701.962053775787, 
    710.077047348022, 753.491997718811, 807.529985904694, 833.763957023621, 
    849.415004253387, 852.05203294754, 876.634955406189, 962.200999259949, 
    1067.1129822731, 1090.93701839447, 1098.68103265762, 1111.81700229645, 
    840, 880, 950, 950, 950, 1210, 1290, 1290, 1290, 1290, 519.912958145142, 
    529.159009456635, 542.180001735687, 584.949016571045, 597.183048725128, 
    602.801024913788, 612.320005893707, 633.880972862244, 639.142990112305, 
    641.510009765625, 699.303984642029, 728.376030921936, 731.564998626709, 
    756.54000043869, 774.587035179138, 834.501981735229, 887.773990631104, 
    900.014996528625, 948.895037174225, 347.997009754181, 388.738989830017, 
    527.451992034912, 575.237035751343, 632.254958152771, 1120, 
    240.203022956848, 307.247042655945, 416.332960128784, 438.155055046081, 
    529.606997966766, 532.611966133118, 557.883024215698, 616.69796705246, 
    633.761048316956, 646.991014480591, 676.961004734039, 746.246993541718, 
    752.32994556427, 794.799029827118, 814.175963401794, 916.736006736755, 
    934.934973716736, 374.673008918762, 455.820024013519, 539.529025554657, 
    577.70699262619, 581.636965274811, 651.882946491241, 697.731971740723, 
    721.569001674652), response = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 
    0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 
    0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 
    0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
    0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 
    0, 0, 0, 0, 0, 0, 0), age = c(8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L), grade = c("3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
    "3", "3", "3", "3", "3", "3"), gender = c("F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 
    "F", "F", "F", "F", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", "M", 
    "M", "M", "M", "M", "M", "M"), condition_id = c(9L, 1L, 22L, 
    11L, 13L, 4L, 10L, 20L, 7L, 3L, 17L, 6L, 5L, 16L, 24L, 25L, 
    2L, 21L, 18L, 14L, 8L, 23L, 15L, 12L, 19L, 8L, 2L, 7L, 14L, 
    15L, 21L, 3L, 16L, 6L, 24L, 10L, 18L, 20L, 4L, 12L, 13L, 
    23L, 25L, 22L, 5L, 17L, 19L, 11L, 1L, 9L, 5L, 8L, 13L, 12L, 
    25L, 22L, 11L, 9L, 16L, 17L, 18L, 10L, 19L, 14L, 15L, 20L, 
    23L, 4L, 1L, 2L, 3L, 6L, 7L, 21L, 24L, 18L, 24L, 22L, 1L, 
    21L, 16L, 4L, 17L, 13L, 25L, 2L, 5L, 9L, 20L, 15L, 14L, 10L, 
    23L, 3L, 19L, 6L, 11L, 8L, 12L, 7L, 11L, 12L, 20L, 22L, 10L, 
    4L, 7L, 6L, 9L, 25L, 5L, 19L, 24L, 14L, 8L, 23L, 18L, 3L, 
    13L, 2L, 16L, 21L, 17L, 15L, 1L, 7L, 2L, 19L, 20L, 9L, 8L, 
    11L, 10L, 12L, 18L, 15L, 22L, 21L, 23L, 1L, 3L, 4L, 5L, 6L, 
    13L, 16L, 14L, 17L, 24L, 25L, 2L, 4L, 7L, 18L, 3L, 1L, 12L, 
    13L, 15L, 19L, 11L, 14L, 20L, 16L, 23L, 5L, 6L, 8L, 9L, 10L, 
    17L, 21L, 24L, 22L, 25L, 9L, 18L, 2L, 14L, 22L, 10L, 24L, 
    12L, 3L, 5L, 13L, 20L, 16L, 1L, 4L, 15L, 7L, 11L, 17L, 25L, 
    6L, 23L, 21L, 19L, 8L, 3L, 4L, 14L, 2L, 21L, 16L, 25L, 13L, 
    12L, 18L, 24L, 15L, 10L, 8L, 23L, 11L, 9L, 20L, 1L, 5L, 22L, 
    17L, 7L, 19L, 6L)), row.names = c(NA, -225L), class = "data.frame")
#2

Welcome to Stan discourse and sorry for the late reply.

Can you write down in mathematical notation what psychometric function exactly you are trying to fit? Then it is easier for me to translate this to brms non-linear syntax.

Perhaps, the following paper about Bayesian IRT modeling with brms may also be of help: https://arxiv.org/abs/1905.09501

#3

Thanks for the response!

Sure–in mathematical notation, the psychometric function I’m trying to fit is

17%20PM

K is the normalization factor given by K = 2^-n.

08%20PM
small k here is the number of different stimulus levels.

I’m still looking through your paper, but it has been helpful in better understanding the brms jargon and in better understanding the brms structure. I’m still figuring out how everything in the formula section comes together.

Below is my improved–but still failing–code. This is a logit function, though ideally it would be a probit, but I can’t quite wrap my head around how the Phi() function fits in with everything else. For instance, if I were to fit as a glmer, the code would be something like: response ~ timepoint*condition * norm + (norm | timepoint/pid/condition), flanker.psych (where norm is the normalized response window per participant, and timepoint is to get threshold values for each of the four timepoints). But I’m confused about how to fit this into brms.

Thanks!

BF <- bf(
  response ~ guess + (1-guess-lapse)/(1+exp(-(norm - threshold)/spread)), #Phi(eta),
  threshold ~ 1 + (norm | pid) + (norm|condition), 
  spread + guess + lapse ~ 1 + (1|pid) + (1|condition), #(1|p|pid) + (1|c|condition),
  family = bernoulli(link = "identity"), #(link = "probit") 
  nl = TRUE)    
  
priors <- c(
  prior(beta(7, 3), nlpar = "threshold"),
  prior(beta(1.4, 1.4), nlpar = "spread", lb = .005, ub = .5),
  prior(beta(.5, 8), nlpar = "lapse", lb = 0, ub = .1),
  prior(beta(1, 5), nlpar = "guess", lb = 0, ub = .1)
)  

fit <- brm(
  BF,
  data = flanker.psych, #ace.threshold,
  inits = 0,
  control = list(adapt_delta = 0.99),
  prior = priors
)
#4

Can you just show the math formual of the psychometric function. What you screenshoted seems to be the whole likelihood which we do not need right now.

To model the function via inverse of logit use

response ~ guess + (1-guess-lapse) * inv_logit((norm - threshold) / spread)

for probit simply replace inv_logit by Phi. Another problem is that guess and lapse have to be further restricted so that the whole predictor stays within [0, 1].

I would recommend to first get rid of lapse in order to the the rest working. Once that’s done, you may add lapse again. Basically starting too complicated will likelily lead to failure. Start simple, and gradually work your way up to more complex models.

For the model with just guess and without lapse, you may try

response ~ inv_logit(guess) + (1 - inv_logit(guess)) * Phi((norm - threshold) / spread)

for reasons see Bayesian IRT paper I linked to above.

I would recommend starting simple that is without any multilevel structure, and once that is working gradually expand the multilevel structure on the non-linear parameters.

#5

So the math formula without the likelihood function would be

19%20PM

Hmmm, I took your advice and started small–without the lapse parameter or random effects, and I’m still getting an error.

Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."

I found a page regarding this error: Rejecting initial value: Log probability evaluates to log(0), i.e. negative infinity. I tried different variations of code and samples from the dataset, but couldn’t get it to run.

Was this the full model you had in mind?

BF <- bf(
  response ~ inv_logit(guess) + (1 - inv_logit(guess)) * inv_logit((norm - threshold) / spread),
  threshold ~ 1, #+ (norm | pid) + (norm|condition), 
  spread~1,
  guess ~ 1, # + (1|pid) + (1|condition), #(1|p|pid) + (1|c|condition),
  family = bernoulli(link = "identity"), #(link = "probit") 
  nl = TRUE)    
  
priors <- c(
  prior(beta(7, 3), nlpar = "threshold"),
  prior(beta(1.4, 1.4), nlpar = "spread", lb = .005, ub = .5),
  #prior(beta(.5, 8), nlpar = "lapse", lb = 0, ub = .1),
  prior(beta(1, 5), nlpar = "guess", lb = 0, ub = .1)
)  

fit <- brm(
  BF,
  data = ace.threshold.t1.samp, #ace.threshold,
  inits = 0,
  control = list(adapt_delta = 0.99),
  prior = priors
)
#6

Now you restrict the range of guess twice, once via lb = 0 and another time via inv_logit. Both work for an intercept only parameter, but only the latter works reliably when you include more predictors or multilevel structure. This model is more or less a 3PL model from item response theory (when adding a lapse it’s a 4PL model). I suggest you carefully read the corresponding sections in https://arxiv.org/abs/1905.09501 and also the relevant parts in https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html

#7

Okay, I think I’m getting this…

I’ve gone through the Item Response Modeling Paper as well as the “Estimating Non-Linear Models With BRMS vignette” with a fine-toothed comb. Then I studied it, and then reviewed it some more. I think I’m getting most of it. I started simple, built up, and was able to construct the 4 parameter model, which I think is correct. (There were still divergent transitions (~900), but I’m running the model now with maxtree depth of 15, so hopefully that resolves the issue; I’ve also only been running one chain since estimating the first 10 participants was taking a while.)

However, I have a couple of questions. Were I to run the model in GLMER, I would model fixed (population-level?) and random (group) effects like this: glmer(response ~ condition * norm + (norm | pid/condition)…). In brms, I thought that I would be able to add the population-level effect of condition (and eventually time, since there are 4 timepoints) as covariates, but I’m getting the error “Error: Factors with more than two levels are not allowed as covariates.” .

It’s possible that I’m not fully comprehending how to move from a mixed-model to a Bayesian one in BRMS, but it still seems like I would obtain the 80% threshold parameter through using fixef() and ranef() to extract estimates (or, if I’m correct, use coef() to extract those together) and then calculate location as (b1/b2) and scale as (1/b2), where b1 is “threshold intercept” (fixed and random), and b2 is norm/response window. Then, I could use qnorm(.8, location, scale) or qlogis (if using logistic) to calculate the 80% threshold. Then I thought that perhaps, since I should be thinking of the “non-linear parameters as placeholders for linear predictor terms , ” I might be able to locate the 80% threshold through extracting individual threshold parameters, but I can’t figure out—or find info on—how to do this. (It looks like Parnames() just extracts the names of the parameters themselves?)

Then I just wanted to make sure that I’ve input my group-level factors correctly using the 1|i|condition and pid, since these are likely correlated across parameters. (Also, it seems like overkill to estimate lapse for every condition for every person, so I just did (1|p|pid) for lapse and guess).

Thanks for your help!

Here is my model:

form.thresholds.four <- bf(
  response ~ gamma + (1 - lambda - gamma) * inv_logit((norm - threshold)/spread),
  threshold ~ 1 + (1|i|condition) + (1|p|pid),
  norm ~ 1 + (1|i|condition) + (1|p|pid),
  logitgamma ~ 1 + (1|p|pid),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitlambda ~ 1 + (1|p|pid),
  nlf(lambda ~ inv_logit(logitlambda)),
  spread ~1 + (1|i|condition) + (1|p|pid),
nl = TRUE)

prior <-
  prior(beta(4, 2), class = "b", nlpar = "threshold") + 
  prior(beta(1.4, 1.4), class = "b", nlpar = "spread", lb = .005, ub = .5) +
  prior(beta(.5, 8), nlpar = "logitlambda", lb = 0, ub = .1)+
  prior(beta(1, 5), nlpar = "logitgamma", lb = 0, ub = .1)

fit_thresholds.4pl.norm <- brm(
  formula = form.thresholds.four,
  data = ace.threshold.t1.samp,
  family = bernoulli(link = "identity"),
  prior = prior,
  control = list(adapt_delta = .80, max_treedepth = 10),
  chains = 1,
  cores = 16
)
#8

Estimating both threshold and norm via ~ 1 + (1|i|condition) + (1|p|pid) will very likely yield an unidentified model as you compute norm - threshold. In this case, the model cannot distinquish between what happens in threshold and what happesn in norm.

#9

Okay, so do I wouldn’t need to add condition as a fixed effect?

In that case, how would I extract an 80% threshold from the model? (It seems like there is a quantile function for AsymLaplace, but not for probit or logit functions.

Would something like the following be sufficient; or should I leave norm correlated [(1|i|condition) + (1|p|pid)] while leaving threshold as [1|condition) + (1|pid)]?

form.thresholds.phi.norm.uncor.newparam <- bf(
  response ~ gamma + (1 - lambda - gamma) * Phi((norm - threshold)/spread), #* k
  threshold ~ 1 + (1|condition) + (1|pid),
  norm ~ 1 + (1|condition) + (1|pid),
  logitgamma ~ 1 + (1|p|pid),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitlambda ~ 1 + (1|p|pid),
  nlf(lambda ~ inv_logit(logitlambda)),
  spread ~1 + (1|i|condition) + (1|p|pid),
nl = TRUE)

Thanks!

#10

You should think about what (theoretically) the parameters norm and threshold mean and parameterize them accordingly. Modeling both as depending on both condition and pid is certainly not going to work.

#11

Ah, I see. I don’t actually think I need a parameter on norm.

If I wanted to add condition as a population-level effect, is there a special way I would need to do that? I keep getting the error

“Error: Factors with more than two levels are not allowed as covariates

form.thresholds.phi.norm.uncor.newparam <- bf(
  response ~ gamma + (1 - lambda - gamma) * condition * Phi((norm - threshold)/spread), #* k
  threshold ~ 1 + (1|i|condition) + (1|p|pid),
  logitgamma ~ 1 + (1|p|pid),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitlambda ~ 1 + (1|p|pid),
  nlf(lambda ~ inv_logit(logitlambda)),
  spread ~1 + (1|i|condition) + (1|p|pid),
#12

Inside the linear formulas for the non-linear parameters.

#13

Interesting…I tried that before and got the following warning:

Error: Duplicated group-level effects are not allowed.
Occured for effect 'Intercept' of group 'pid'.

So then I returned to a simpler model, with only an intercept on the spread parameter. That worked. Then I added random slopes of response window (norm). That worked. Then I added group effect of participant to spread, then condition. It all worked. Surprising. I have to deal with quite a few divergent transitions, but this is exciting!

Would I have to interpolate the 80% threshold, or is there a way to pull out an 80% threshold parameter? If I use coef or fixef/ranef I’ll only get a threshold value for the prior I set, right? If I were using GLMER, I could just use fixef/ranef to get the location and scale and then use qnorm(), but that’s without having estimated parameters.

form.thresholds.phi.nonormparam.cond <- bf(
  response ~ gamma + (1 - lambda - gamma) * Phi((norm - threshold)/spread), #*k=
  threshold ~ 1 + condition + (norm|p|pid) + (norm|condition),
  logitgamma ~ 1 + (1|p|pid),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitlambda ~ 1 + (1|p|pid),
  nlf(lambda ~ inv_logit(logitlambda)),
  spread ~ 1 + (1|p|pid) + (1|condition),
nl = TRUE)
#14

In order to find the 80% threshold, you probably have to do a lot of manual work of transforming the parameters yourself in the right way.

I don’t see what you mean by “If I use coef or fixef/ranef I’ll only get a threshold value for the prior I set, right?” though? Neither is coef (or the others) returning a threshold per se (it may be of course depending on the model, so it may be true in your case) nore are those prior values. Instead, they are the posterior estimates of the group-level coeffiecients.

#15

I thought I had found a way to do this is Matlab, but alas the model wasn’t hierarchical, which is what I’m pretty sure I need given that the tasks are adaptive and there are few trials per participant.

I’ve read through the item response examples again, which was helpful review. To clarify the last question I asked many months ago–if I wanted to locate the mean of the psychometric function for each participant for every task (each participant would have 14, since there are 14 conditions), I should be able to get these means by adding the intercept, with the mean for each of the conditions, with the mean for each participant, right? But when I call coef() there are values for the four lapse, guess, spread, and threshold parameters. What would give me the intercept and mean difference values for the actual function vs the parameters? Does that make sense?

If I were to create a psychometric function equivalent in glmer without the guess/lapse, etc. parameters, it would be something like

glmer(response ~  condition * norm + (norm | pid/condition)

(where norm is the normalized response window). To be doing the equivalent with brms should I be including that interaction of condition*norm in the psychometric function itself (before the gamma variable) such that it would be response ~ condition * norm + gamma + (1-lambda…), or should I be putting this interaction into the threshold parameter below the function? In your two-item response example, you add the covariates directly into the model, but I wasn’t sure if that is what I should be doing here?

Thanks a lot for your help! Hope you’ve been well.

fit.form <- bf(
  response ~ gamma + (1 - lambda - gamma) * Phi((norm - threshold)/spread),
  threshold ~ 1 + (norm|p|pid) + (norm|condition),
  logitgamma ~ 1 + (1|p|pid),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitlambda ~ 1 + (1|p|pid),
  nlf(lambda ~ inv_logit(logitlambda)),
  spread ~ 1 + (1|p|pid) + (1|condition),
nl = TRUE)

prior <-
  prior(beta(9, 3), class = "b", nlpar = "threshold", lb = 0, ub = 1) +
  prior(beta(1.4, 1.4), class = "b", nlpar = "spread", lb = .005, ub = .5) +
  prior(beta(.5, 8), nlpar = "logitlambda", lb = 0, ub = .1)+
  prior(beta(1, 5), nlpar = "logitgamma", lb = 0, ub = .1)

fit_form <- brm(
  formula = form.thresholds.phi.nonormparam.cond,
  data = ace.threshold.t1.samp,
  family = bernoulli(link = "identity"),
  prior = prior,
  control = list(adapt_delta = .85, max_treedepth = 15),
  inits = 0,
  chains = 1,
  cores = 16
)
#16

Hi, sorry for the late response. I am currently unsure what exactly you mean. In general your assumption is right, that you can just add coefficients together, on a per-draw basis, that is, when coef(..., summary = FALSE). On which of the non-linear parameters you want to do this transformation I cannot tell you because each of them has a different meaning and it depends on in what parameter you are interested in.