[R-sig-ME] Hierarchical Psychometric Function in BRMS

Mon Mar 16 10:10:44 CET 2020

Hi James,

since I am working with brms and glmer, I feel I should be able to give a
response (although addressing Paul in the Stan-Forum might be
a better option), there seem to be two questions, and some missing details,
that might lead to even more questions.... let's begin....

My questions:
1. "14 executive functions". Does this mean every participant completed
each of 14 tasks supposed to measure different facets of the general
construct "executive functions in working memory"? (If not, please
clarify). What term is this in the model "condition" or "norm"? (Given that
you have random slopes for "norm" it seems to be "norm" ?) Then what is
condition?

2. "adaptive tasks with 25 to 40 trials" Does this mean "tailored testing"?
(I.e., the trial that comes next within the task depends on the decisions
(their error) from all previous trials?)

3. "Goal: disentangle the response window at which participants reach a
70%", - if you have tailored testing (I am not sure), which already is
designed to sort trials to meander around 75% accuracy for maximum
information/variance , this threshold seems a bit unmotivated, can you give
more background?

4. "four different time points" , I suppose these are four sessions, in
each the participants have completed subsets of the 14 tasks

Your (secondary) questions (I ignore points 1 to 3 now, but they need
clarification):
"I'm not sure whether the four timepoints can be fit at once because
probability distributions for random factor of participant are already used
to account for repeated measures of participant completing 14 conditions)."
My answer:
- Regardless of the technical details:  First, "time points"  has only four
levels, thus, it would not make sense to separate their "random" intercepts
from other variance sources in the design, no matter which. Computing
standard deviations of a distribution for which you only have 4
observations/levels is problematic. Second, nonetheless assuming that "time
points" (e.g., increasing ability over time) has an effect, then
controlling for it is pretty legit, so, it makes sense to include "time
points" into the fixed effects. Also legit.

5. "The other problem I'm having is using coef() or fixef()/ranef() to
withdraw (or locate) the overall intercept and slope such that I can use
the qlogis() function to determine the psychometric threshold at 70% (since
I don't think it would be accurate to directly pull the 70% threshold
estimate from the parameter itself?)."
My answer:
- Do you mean, by 70% threshold, the "location" on the predictor(s) (the
logit) at which the predicted probably of the response is 70%? (Please keep
in mind, that you have two interacting predictors in your model, which
means getting these estimates for one predictor requires to either ignore
variance of the other predictor, which needs theoretical clarification if
you want to interpret this; or taking it into account - see below.) Anyway,
the "manual" way to do this, is to make predictions, based on the
coefficients, and then search the point of crossing 70%. For this you want
to use the "emmeans" package which works for both glmer and brms (but I am
not sure whether it works also for the non-linear models; if not, you need
to ask Paul Buerkner in the Stan forum how to do it ;)); it sure works with
standard hierarchical regression output from brms.) . In the emmeans
package you find the function "emmip", which is what you desire.

#assuming this is your model with a continuous predictor ("continuous") and
a factorial predictor ("factor"):
model<-glmer(response ~  continuous * factor + (continuous | pid))
emmip(model,~continuous,at = list(continuous = c(1,2,3,4,5,6),
type="response",CIs=TRUE, engine="ggplot" )
# this gives you the probability predictions for "continuous" from 1 to 6
(you can make these as "fine" as you want), while ignoring "factor"
# if you want it "by factor" (taking the interaction into account) you can
write:
emmip(model,~continuous|factor ,at = list(continuous = c(1,2,3,4,5,6),
type="response",CIs=TRUE, engine="ggplot" )
#All you have to do is search for the point crossing 70% then :) .

However, as noted, non-linear brms models might not directly translate to
the emmeans architecture (I don't know), and there is a more elegant
solution anyway:

1. A standard logistic function predicts 50% when the logit becomes 0
(before applying the exponential ratio rule; I ignore the fact that your
gamma and lambda model terms absolutely destroy this property... :))
2. The "intercept" shifts the whole logit statically (or by factorial
conditions), such that it indicates "where" 50% is predicted (in a given
condition). For example, in standard models
1/(1+exp(intercept+varyingeffects)) the intercept says for which value of
varyingeffects  the term becomes 0).
3. You can "make the intercept" to indicate a 70% prediction instead of a
50% prediction, if you add a constant on the logit level; that is:
1/(1+exp(-.8477)) = (about) 70%; and
 1/(1+exp(-.8477+intercept+varyingeffects)) shifts the intercept by this
constant, such that it now indicates the value of varyingeffects which
predicts 70%. I guess. .. :)) There could be more detail to that (which I
don't see right now), but it sure is a starting point.

Hope this helps, with your actual questions.
The rest seems to be a different matter.... (e.g., taking dependencies of
tailored testing into account etc).

But one final note: I have once tried to fit simpler models with
constructing the logit myself, like you do, and then setting,  family =
bernoulli(link = "identity"), which never worked (it never converged). ...
Just saying: I think Paul makes some points about the identifiability of
those models in his vignettes, which you should check, if your model fails
converging.
(But dropping lambda and gamma, might be worth considering in any case. If
you simulate logistic functions hierarchically, then they do not
approximate 100% on average (which would be the reason you use gamma and
lambda), but the limited growth approximates e.g., 80 % depending on the
individual variations in the slope parameters of the logistic function.
This means, you don't need "maximum performance" parameters, but can
approximate this behavior by the assumption of hierarchically clustered
variance. Which also makes the model simpler... , and identifiable, and you
could use the "elegant" way of determining 70%).

Best, Ree

Am Mo., 16. März 2020 um 04:28 Uhr schrieb Ades, James <
jades using health.ucsd.edu>:

> Hi all,
>
> Given that this is a mixed-model listserv, I'm hoping that a BRMS question
> might fit within that purview.
>
> A quick synopsis of the dataset: there are 14 different conditions of
> executive function tasks ( ~1000 3rd, 5th, 7th graders). Given that these
> tasks use an adaptive paradigm (tasks might have anywhere from 25 to 40
> trials), I'm trying to disentangle the response window at which
> participants reach a 70% performance threshold. There are four separate
> timepoints. (I'm not sure whether the four timepoints can be fit at once
> because probability distributions for random factor of participant are
> already used to account for repeated measures of participant completing 14
> conditions, but that question is secondary to ensuring that I'm fitting one
> time point correctly and adequately extracting those the intercept/slope
> parameters).
>
> If I were to only input this into glmer without the priors, I'd write the
> model as:
> ```
> glmer(response ~  condition * norm + (norm | pid/condition)
> ```
> (In a glmer model, I can extract intercept/slope parameters fine).
>
> My current model is below. My question isn't so much with the psychometric
> function or the priors, which, besides the threshold, I've borrowed from
> Treutwein and Strasburger:
> https://link.springer.com/article/10.3758/BF03211951--though if there are
> contentions with any of the those, feel free to raise them--as it is
> whether I've correctly structured the non-linear parameters. The reason for
> modeling all four parameters is to minimize bias, but threshold is the only
> estimate that I'm concerned with. So regarding the multi-level structure,
> I've created parameters for lapse, guess, spread, and threshold. It seems
> reasonable to expect that threshold and spread will vary for every
> participant for every condition, while lapse and guessing (forced yes/no)
> will likely not differ much from condition to condition within participant
> (though if there are arguments that it would make for an improved model,
> I'm fine including lapse and guess parameters for every condition as well).
>
> The other problem I'm having is using coef() or fixef()/ranef() to
> withdraw (or locate) the overall intercept and slope such that I can use
> the qlogis() function to determine the psychometric threshold at 70% (since
> I don't think it would be accurate to directly pull the 70% threshold
> estimate from the parameter itself?).
>
> Does all of that make sense? This is all a little bit over my head and
> though I've culled Buerkner's item-response vignettes (Here:
> https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html
> and here: https://arxiv.org/pdf/1905.09501.pdf, they're similar but
> fundamentally different, so they only get me so far).
>
> I've included a small sample of ~five participants here:
> https://drive.google.com/file/d/1YFnQRSjnp5hVziQx5wQzaIhn75KigaGx/view?usp=sharing
>
> Thanks in advance for any and all help! Hope everyone is staying healthy!
>
> James
>
>
> ```
> thresholds <- bf(
>   response ~ (gamma + (1 - lambda - gamma) * Phi((norm -
> threshold)/spread)),
>   threshold ~ 1 + (1|p|pid) + (1|c|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|c|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_thresholds <- brm(
>   formula = thresholds,
>   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
> )
> ```
>
>
>
> [
> https://media.springernature.com/w110/springer-static/cover/journal/13414.jpg
> ]<https://link.springer.com/article/10.3758/BF03211951>
> Fitting the psychometric function | SpringerLink<
> https://link.springer.com/article/10.3758/BF03211951>
> A constrained generalized maximum likelihood routine for fitting
> psychometric functions is proposed, which determines optimum values for the
> complete parameter set—that is, threshold and slopeas well as for guessing
> and lapsing probability. The constraints are realized by Bayesian prior
> distributions for each of these parameters. The fit itself results from
> maximizing the posterior ...
> link.springer.com
>
> Abstract R arXiv:1905.09501v2 [stat.CO] 20 Jul 2019<
> https://arxiv.org/pdf/1905.09501.pdf>
> Paul-Christian B urkner 3 dictions via a nested non-linear formula syntax,
> the implementation of several distributions designed for response times
> data, and extentions of distributions for ordinal data, for example
> arxiv.org
>
> Estimating Non-Linear Models with brms<
> https://cran.r-project.org/web/packages/brms/vignettes/brms_nonlinear.html
> >
> Introduction. This vignette provides an introduction on how to fit
> non-linear multilevel models with brms.Non-linear models are incredibly
> flexible and powerful, but require much more care with respect to model
> specification and priors than typical generalized linear models.
> cran.r-project.org
>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]