1. The Pair-Clustering Model

Before we introduce MPT models and the ML estimation methods in more detail, we briefly describe the pair-clustering model (e.g., Batchelder & Riefer, Reference Batchelder and Riefer1980, Reference Batchelder and Riefer1986) that we use throughout the article to illustrate the proposed methods. The pair-clustering model can be used to analyze data from a free-recall experiment in which participants are presented with a list of word pairs plus a number of singletons. All words are preselected to be equally difficult in terms of memorizability. Each word pair consists of semantically related words (e.g., rose–tulip), whereas singletons are unrelated to other words in the list. These pairs and the singletons are presented one word at a time, and participants are later asked to recall the items from the list in any order. On the basis of participants’ free recall performance, the studied word pairs and singletons are then assigned to one of six response categories. It is assumed that the probability of each response category can be modeled by two processing trees that include a total of four parameters (see Fig. 1 for a graphical illustration of the model).

Figure 1. The pair-clustering MPT model (adapted from Riefer & Batchelder, Reference Riefer and Batchelder1988, p. 330, Figure 2). Rectangles indicate stimulus classes (left) and observable responses (right). Rectangles with rounded corners represent latent cognitive states. Parameters attached to the branches indicate transition probabilities from left to right, specifically, storing a word pair as a cluster (c), retrieving a stored cluster in free recall (r), storing and retrieving a word from a non-clustered pair in free recall (u), and storing and retrieving a singleton in free recall (a).

Specifically, the model defines four response categories for the word pairs: Category $C_{11}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{11}$$\end{document} includes all cases where both words in a pair are recalled adjacently, $C_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{12}$$\end{document} represents cases in which both words are recalled nonadjacently, $C_{13}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{13}$$\end{document} corresponds to cases where exactly one word is recalled, and finally $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} represents cases in which neither word from the pair is recalled. In addition, the singletons are scored in two categories: singletons recalled successfully ( $C_{21}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{21}$$\end{document} ) versus not successfully ( $C_{22}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{22}$$\end{document} ).

The pair-clustering model proposes four cognitive process parameters that jointly determine the response probabilities: (1) the probability c that a word pair is stored as a cluster in memory, (2) the probability r that a stored cluster is successfully retrieved in free recall, (3) the probability u that one word from a pair that is not stored as part of a cluster is successfully retrieved, and finally, (4) the probability a that a singleton is successfully stored and retrieved. These parameters can be used to derive response probabilities for each of the six categories by calculating the product of all parameters along a branch and then summing up all branch probabilities that refer to the same category. For example, because there is only one branch that terminates in category $C_{11}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{11}$$\end{document} (see Fig. 1), the probability of this category (i.e., both words recalled adjacently) is just the product $c·r$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c \cdot r$$\end{document} (i.e., successful storage as a cluster followed by successful retrieval of this cluster). Applying the same logic to all branches in Fig. 1 results in the following six model equations:

(1) $\begin{array}{cc}P\left(C_{11}\right)& =c·r\\ P\left(C_{12}\right)& =(1-c)·u^2\\ P\left(C_{13}\right)& =(1-c)·u·(1-u)+(1-c)·(1-u)·u\\ P\left(C_{14}\right)& =c·(1-r)+(1-c)·{(1-u)}^2\\ P\left(C_{21}\right)& =a\\ P\left(C_{22}\right)& =1-a\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(C_{11})&= c \cdot r \nonumber \\ P(C_{12})&= (1 - c) \cdot u^2 \nonumber \\ P(C_{13})&= (1 - c) \cdot u \cdot (1 - u) + (1 - c) \cdot (1 - u) \cdot u \nonumber \\ P(C_{14})&= c \cdot (1 - r) + (1 - c) \cdot (1 - u)^2 \nonumber \\ P(C_{21})&= a \nonumber \\ P(C_{22})&= 1 - a \end{aligned}$$\end{document}

The goal of MPT modeling, applied to the pair-clustering model, is to estimate the four cognitive process parameters and to test reasonable hypotheses about these parameters. One example of such a hypothesis would be that unclustered words from a pair behave like singletons in memory, that is, they are stored and retrieved with the same probability (i.e., $u=a$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u = a$$\end{document} ).

In the next section, we begin with a formal introduction of a typical MPT model without random effects. This is followed by an outline of the latent-trait MPT model with some extensions. In the third section, we present three methods of marginal ML estimation for the latent-trait model. Next, we describe the results of a small simulation study in which we compare the different methods with respect to accuracy and speed. In Sect. 3.2, we present an application of our method to an illustrative example using real data. In the final section, we discuss possible extensions and open questions for future research.

2. MPT Models Without Person-Level Random Effects

The pair-clustering model we just described is a specific instance of an MPT model. On a more general level, MPT models are statistical models of response frequencies of mutually exclusive, independent categories. To define the model, we assume that there are $k=1,\dots ,K$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k = 1, \ldots , K$$\end{document} category systems and that each category system consists of $j=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j = 1$$\end{document} , $\dots$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ldots $$\end{document} , $J_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_k$$\end{document} categories $C_{\mathrm{kj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kj}$$\end{document} . In the pair-clustering model, there are $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 2$$\end{document} category systems, one for the word pairs and one for the singletons. The first system consists of four categories: $C_{11}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{11}$$\end{document} , $C_{12}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{12}$$\end{document} , $C_{13}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{13}$$\end{document} , and $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} ; and the singleton system contains two categories: $C_{21}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{21}$$\end{document} and $C_{22}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{22}$$\end{document} . We further assume that data from $t=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = 1$$\end{document} , $\dots$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ldots $$\end{document} , T individuals is available. For each single individual t, the data for category system k are given as a vector of frequencies, $n_{\mathrm{kt}}=(n_{k1t},\dots ,n_{k{J}_{k}t})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{n}_{kt} = (n_{k1t}, \ldots , n_{kJ_{k}t})$$\end{document} . For instance, $n_{1t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{n}_{1t}$$\end{document} contains the four frequencies of person t for the four word pair categories (i.e., k = 1). These frequencies are assumed to stem from a multinomial distribution

(2) $\begin{array}{c}f\left(n_{\mathrm{kt}}\right|\theta )=\left(\begin{array}{c}{N}_{\mathrm{kt}}\\ n_{k1t}...n_{k{J}_{k}t}\end{array}\right)p_{k1t}^{n_{k1t}}·p_{k2t}^{n_{k2t}}\cdots p_{k{J}_{k}t}^{n_{k{J}_{k}t}}=\left(\begin{array}{c}{N}_{\mathrm{kt}}\\ n_{k1t}...n_{k{J}_{k}t}\end{array}\right)\prod _{j=1}^{J_k}{p}_{\mathrm{kjt}}^{n_{\mathrm{kjt}}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}_{kt}|\varvec{\theta }) = \begin{pmatrix} N_{kt} \\ n_{k1t} ... n_{kJ_{k}t} \end{pmatrix} p_{k1t}^{n_{k1t}} \cdot p_{k2t}^{n_{k2t}} \cdots p_{kJ_{k}t}^{n_{kJ_{k}t}} = \begin{pmatrix} N_{kt} \\ n_{k1t} ... n_{kJ_{k}t} \end{pmatrix} \prod _{j = 1}^{J_{k}} p_{kjt}^{n_{kjt}} \end{aligned}$$\end{document}

where $N_{\mathrm{kt}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{kt}$$\end{document} = $n_{k1t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{k1t}$$\end{document} +... + $n_{k{J}_{k}t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{kJ_{k}t}$$\end{document} and $p_{k1t}+\cdots +p_{k{J}_{k}t}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{k1t} + \cdots + p_{kJ_{k}t} = 1$$\end{document} . This definition requires the vector of all frequencies across all category systems K of person t, that is, $n_t=(n_{1t},...,n_{\mathrm{Kt}})$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{n}_{t} = (\varvec{n}_{1t},..., \varvec{n}_{Kt})$$\end{document} , to follow a product multinomial distribution:

(3) $\begin{array}{c}f\left(n_t\right|\theta )=\prod _{k=1}^K\left(\begin{array}{c}{N}_{\mathrm{kt}}\\ n_{k1t}...n_{k{J}_{k}t}\end{array}\right)\prod _{j=1}^{J_k}{p}_{\mathrm{kjt}}^{n_{\mathrm{kjt}}}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}_{t}|\varvec{\theta }) = \prod _{k = 1}^{K} \begin{pmatrix} N_{kt} \\ n_{k1t} ... n_{kJ_{k}t} \end{pmatrix} \prod _{j = 1}^{J_{k}} p_{kjt}^{n_{kjt}} . \end{aligned}$$\end{document}

The idea behind MPT models is that the probability of the occurrence of an event category can be modeled as a function of the cognitive process parameters. In the standard model (i.e., without person-level random effects), these process parameters $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} are unknown constants that do not vary between individuals. Thus, $p_{\mathrm{kjt}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p_{kjt}$$\end{document} is written as

(4) $\begin{array}{c}{p}_{\mathrm{kjt}}=P\left(C_{\mathrm{kj}}\right|\theta )\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p_{kjt} = P(C_{kj} | \varvec{\theta } ) \end{aligned}$$\end{document}

where $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} contains the S cognitive process parameters and is thus an element of ${[0,1]}^S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0, 1]^S$$\end{document} , $s=1,...,S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s = 1,..., S$$\end{document} . In the example, $\theta =(c,r,u,a)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta } = (c, r, u, a)$$\end{document} , and the probability of $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} (see above) is $p_{14t}=P\left(C_{14}\right|\theta )=c(1-r)+(1-c){(1-u)}^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{14t} = P( C_{14} | \varvec{\theta } ) = c(1-r) + (1-c)(1-u)^2$$\end{document} . To write $P\left(C_{\mathrm{kj}}\right|\theta )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P( C_{kj} | \varvec{\theta } )$$\end{document} more generally, we define $I_{\mathrm{kj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ I_{kj}$$\end{document} to be the number of branches of the tree that terminate in $C_{\mathrm{kj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kj}$$\end{document} , with $i=1,...,I_{\mathrm{kj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i = 1,..., I_{kj}$$\end{document} indexing a specific branch $B_{\mathrm{kji}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{kji}$$\end{document} that terminates in category j of category system k. As described by Hu and Batchelder (Reference Hu and Batchelder1994), we define:

(5) $\begin{array}{c}P\left(B_{\mathrm{kji}}\right|\theta )=\prod _{s=1}^S{\theta }_s^{a_{\mathrm{skji}}}{(1-{\theta }_s)}^{b_{\mathrm{skji}}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P( B_{kji} | \varvec{\theta }) = \prod _{s = 1}^{S} \theta _{s}^{a_{skji}}( 1 -\theta _{s})^{b_{skji}} \end{aligned}$$\end{document}

where $a_{\mathrm{skji}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{skji}$$\end{document} and $b_{\mathrm{skji}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{skji}$$\end{document} indicate how often a process parameter ${\theta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{s}$$\end{document} and its complement $1-{\theta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1-\theta _{s}$$\end{document} , respectively, appear on the branch $B_{\mathrm{kji}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{kji}$$\end{document} . For instance, there are two paths that terminate in $C_{14}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{14}$$\end{document} . To illustrate, we just consider the first of these branches, $B_{141}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{141}$$\end{document} :

(6) $\begin{array}{cc}P\left(B_{141}\right|\theta )& =\prod _{s=1}^4{\theta }_s^{a_{s141}}{(1-{\theta }_s)}^{b_{s141}}\\ & ={\theta }_1^{a_{1141}}{(1-{\theta }_1)}^{b_{1141}}·{\theta }_2^{a_{2141}}{(1-{\theta }_2)}^{b_{2141}}·{\theta }_3^{a_{3141}}{(1-{\theta }_3)}^{b_{3141}}·{\theta }_4^{a_{4141}}{(1-{\theta }_4)}^{b_{1141}}\\ & =c^1{(1-c)}^0·r^0{(1-r)}^1·u^0{(1-u)}^0·a^0{(1-a)}^0=c(1-r)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P( B_{141} | \varvec{\theta })&= \prod _{s = 1}^{4} \theta _{s}^{a_{s141}}( 1 -\theta _{s})^{b_{s141}} \nonumber \\&= \theta _{1}^{a_{1141}}( 1 -\theta _{1})^{b_{1141}} \cdot \theta _{2}^{a_{2141}}( 1 -\theta _{2})^{b_{2141}} \cdot \theta _{3}^{a_{3141}}( 1 -\theta _{3})^{b_{3141}} \cdot \theta _{4}^{a_{4141}}( 1 -\theta _{4})^{b_{1141}}\nonumber \\&= c^{1}(1-c)^{0} \cdot r^{0}(1-r)^{1} \cdot u^{0}(1-u)^{0} \cdot a^{0}(1-a)^{0} = c(1-r) \end{aligned}$$\end{document}

According to these definitions, the probability of $C_{\mathrm{kj}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{kj}$$\end{document} is

(7) $\begin{array}{c}P\left(C_{\mathrm{kj}}\right|\theta )=\sum _{i=1}^{I_k}P\left(B_{\mathrm{kji}}\right|\theta )=\sum _{i=1}^{I_k}\prod _{s=1}^S{\theta }_s^{a_{\mathrm{skji}}}{(1-{\theta }_s)}^{b_{\mathrm{skji}}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(C_{kj} | \varvec{\theta } ) = \sum _{i = 1}^{I_k} P( B_{ kji} | \varvec{\theta }) = \sum _{i = 1}^{I_k} \prod _{s = 1}^{S} \theta _{s}^{a_{skji}}( 1 -\theta _{s})^{b_{skji}} \end{aligned}$$\end{document}

and the multinomial probability of the data for a single individual t is

(8) $\begin{array}{c}f\left(n_t\right|\theta )=\prod _{k=1}^K\left(\begin{array}{c}{N}_{\mathrm{kt}}\\ n_{k1t}...n_{k{J}_{k}t}\end{array}\right)\prod _{j=1}^{J_k}{\left(\sum _{i=1}^{I_{\mathrm{kj}}},\prod _{s=1}^S,{\theta }_s^{a_{\mathrm{skji}}},{(1-{\theta }_s)}^{b_{\mathrm{skji}}}\right)}^{n_{\mathrm{kjt}}}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}_{t}|\varvec{\theta }) = \prod _{k = 1}^{K} \begin{pmatrix} N_{kt} \\ n_{k1t} ... n_{kJ_{k}t} \end{pmatrix} \prod _{j = 1}^{J_{k}} \left[ \sum _{i = 1}^{I_{kj}} \prod _{s = 1}^{S} \theta _{s}^{a_{skji}}( 1 -\theta _{s})^{b_{skji}} \right] ^{n_{kjt}}. \end{aligned}$$\end{document}

Equation 8 can be used to obtain ML estimates of the process parameters $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} . For instance, Hu and Batchelder (Reference Hu and Batchelder1994) showed how the parameters can be estimated with an EM algorithm. Alternatively, methods based on the analytical gradient or finite difference methods can be used for ML parameter estimation. Both approaches are implemented in the R package MPTinR (Singmann and Kellen, Reference Singmann and Kellen2013; Singmann et al., Reference Singmann, Kellen, Gronau, Mueller and Bhel2020) or the package mpt (Wickelmaier and Zeileis, Reference Wickelmaier and Zeileis2011, Reference Wickelmaier and Zeileis2020). Also, a Bayesian approach can be used to estimate the parameters (e.g., Lee & Wagenmakers, Reference Lee and Wagenmakers2014).

3. MPT Models with Person-Level Random Effects

In the previous section, the process parameters contained in $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} were assumed not to differ between individuals. In the latent-trait MPT model, this assumption is relaxed because $R\le S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R {\le } S$$\end{document} elements of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} are written as a function of person-specific random effects that stem from a multivariate normal distribution

(9) $\begin{array}{c}{b}_t\sim MNV(\mu ,\Sigma )\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{b}_t \sim MNV( \varvec{\mu }, \varvec{\Sigma }) \end{aligned}$$\end{document}

where $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} is an $R\times 1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \times 1$$\end{document} vector of expectations and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} is an $R\times R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \times R$$\end{document} covariance matrix. Because each single ${\theta }_{\mathrm{st}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{st}$$\end{document} is a probability parameter, we need to make sure that the resulting coefficient remains in the unit interval [0, 1]. On the basis of the literature on the Generalized Linear Model (GLM), we use the logit-link function (see Coolin et al., Reference Coolin, Erdfelder, Bernstein, Thornton and Thornton2015)

(10) $\begin{array}{c}{\theta }_{\mathrm{st}}=\frac{1}{1+\exp \left(-,(,{\beta }_s,+,b_{\mathrm{st}},)\right)}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta _{st} = \frac{1}{1 + {exp}\left[ -(\beta _s + b_{st})\right] } \end{aligned}$$\end{document}

or the probit-link function (see Klauer, Reference Klauer2010)

(11) $\begin{array}{c}{\theta }_{\mathrm{st}}=\Phi ({\beta }_s+b_{\mathrm{st}})\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta _{st} = \Phi (\beta _s + b_{st}) \end{aligned}$$\end{document}

to map a real-valued person parameter $b_{\mathrm{st}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{st}$$\end{document} on a probability parameter ${\theta }_{\mathrm{st}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{st}$$\end{document} , with ${\beta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _s$$\end{document} serving as a parameter-specific intercept constant (in Eq. 11, $\Phi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi $$\end{document} denotes the cumulative normal distribution). The two approaches typically result in very similar parameter estimates (at least in the GLM framework). To allow for comparisons in the framework of random effects MPT models, we have implemented both link functions in the R package that implements the statistical methods proposed in this article.

We can now define the conditional distribution of the category frequencies of person t given the random effects $b_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{b}_t$$\end{document} :

(12) $\begin{array}{c}f\left(n_t\right|\beta ,b_t)=\prod _{k=1}^K\left(\begin{array}{c}{N}_{\mathrm{kt}}\\ n_{k1t}...n_{k{J}_{k}t}\end{array}\right)\prod _{j=1}^{J_k}{\left(\sum _{i=1}^{I_{\mathrm{kj}}},\prod _{s=1}^S,{\theta }_{\mathrm{st}}^{a_{\mathrm{skji}}},{(1-{\theta }_{\mathrm{st}})}^{b_{\mathrm{skji}}}\right)}^{n_{\mathrm{kjt}}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}_{t}|\varvec{\beta }, \varvec{b}_t) = \prod _{k = 1}^{K} \begin{pmatrix} N_{kt} \\ n_{k1t} ... n_{kJ_{k}t} \end{pmatrix} \prod _{j = 1}^{J_{k}} \left[ \sum _{i = 1}^{I_{kj}} \prod _{s = 1}^{S} \theta _{st}^{a_{skji}}( 1 -\theta _{st})^{b_{skji}} \right] ^{n_{kjt}} \end{aligned}$$\end{document}

where ${\theta }_{\mathrm{st}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{st}$$\end{document} is defined as in Eqs. 10 or 11. Along with the multivariate normal density of the random effects $b_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{b}_t$$\end{document} , Eq. 12 can then be used to obtain ML estimates of the model parameters. We refer to this model as the random effects MPT model.

Before we proceed, we would like to point out that the mean structure in our model is not identified because the expectation of the random effects is $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} and the conditional distribution contains the vector of intercept terms $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} . To identify the mean structure, we define all elements of $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} to be zero when all process parameters are assumed to differ between individuals (i.e., $R=S$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R=S$$\end{document} , where R denotes the number of random effects parameters). In this case, our random effects model corresponds to the model described by Klauer (Reference Klauer2010). However, when only some (but not all) of the process parameters are assumed to be random, we estimate the respective entries in $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} for the $(S-R)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S-R)$$\end{document} fixed parameters while setting the remaining R entries (corresponding to random parameters) to zero and reducing the dimensions of $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} accordingly. Thus, by writing the model as proposed here, we achieve a high degree of flexibility to estimate MPT models with both random and fixed effects. This can be very useful, for example, when MPT models include guessing parameters that need to be estimated and are assumed to be constant across individuals.

Furthermore, we can also extend the model to include person-level covariates contained in the person-specific column vectors $X_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{X}_{t}$$\end{document} . This can be achieved in two ways: First, when the person-level covariates are used to predict R process parameters that are assumed to be random, $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} is person-specific, and we write it as function of the predictor values:

(13) $\begin{array}{c}{\mu }_t=\mu +\Gamma X_t,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\mu }_t = \varvec{\mu } + \varvec{\Gamma }\varvec{X}_{t}, \end{aligned}$$\end{document}

where $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} is an $R\times p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \times p$$\end{document} matrix of weights, and p is the number of person-level predictors used to predict the random process parameters. Note that when there is an assumption that a process parameter is not affected by a predictor in $X_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{X}_{t}$$\end{document} , the respective entry in $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} is simply set to zero. Furthermore, $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} represents the values of the person parameters when the covariates are zero. Second, when a fixed process parameter ${\theta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _s$$\end{document} is assumed to vary as a function of person-level covariates so that the intercept ${\beta }_{\mathrm{st}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{st}$$\end{document} may differ between individuals t, we follow the procedure outlined by Coolin et al. (Reference Coolin, Erdfelder, Bernstein, Thornton and Thornton2015) and write

(14) $\begin{array}{c}{\beta }_{\mathrm{st}}={\beta }_s+{\gamma }_s{X}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta _{st} = \beta _{s} + \varvec{\gamma _{s} } \varvec{X}_{t} \end{aligned}$$\end{document}

where ${\beta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{s}$$\end{document} is the value of the parameter-specific intercept when all person-level covariates are zero and ${\gamma }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma _{s}}$$\end{document} is a row vector containing the weights of the p covariates used to predict ${\theta }_s$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _s$$\end{document} .

3.1. Marginal Maximum Likelihood Estimation

The goal is to estimate the free parameters contained in $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} , $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} , $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} , where $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} and $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} denote the $(S-R)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S-R)$$\end{document} -dimensional vector of intercepts and the $(S-R)\times p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S-R) \times p$$\end{document} matrix of predictor weights, respectively, in the fixed-effects part of our model (corresponding to the model of Coolin et al., Reference Coolin, Erdfelder, Bernstein, Thornton and Thornton2015) while $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} denote to the R-dimensional vector of parameter means, the $R\times p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \times p$$\end{document} matrix of predictor weights, and the $R\times R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \times R$$\end{document} parameter covariance matrix, respectively, in the random-effects part of our model (corresponding to the model of Klauer, Reference Klauer2010). We employ a marginal ML approach that is based on the marginal density of the response frequencies. For a single individual t, this density is

(15) $\begin{array}{c}f\left(n_t\right)={\int }_{b_t}f\left(n_t\right|\beta ,\gamma ,b_t)f\left(b_t\right|\mu ,\Gamma ,\Sigma )d{b}_t.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}_{t}) = \int _{\varvec{b}_t} f(\varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_t) f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) {d}\varvec{b}_t . \end{aligned}$$\end{document}

and for the entire sample, it is

(16) $\begin{array}{c}f\left(n\right)=\prod _{t=1}^{T}f\left(n_t\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(\varvec{n}) = \prod _{t = 1}^{T} f(\varvec{n}_t) \end{aligned}$$\end{document}

Eq. 16 can be used to define the likelihood

(17) $\begin{array}{c}L(\beta ,\gamma ,\mu ,\Gamma ,\Sigma )=\prod _{t=1}^T{\int }_{b_t}f\left(n_t\right|\beta ,\gamma ,b_t)f\left(b_t\right|\mu ,\Gamma ,\Sigma )d{b}_t=\prod _{t=1}^T{\int }_{b_t}{L}_{t}d{b}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L(\varvec{\beta },\varvec{\gamma },\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) = \prod _{t = 1}^{T} \int _{\varvec{b}_t} f(\varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_t) f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) {d}\varvec{b}_t = \prod _{t = 1}^{T} \int _{\varvec{b}_t} L_t {d}\varvec{b}_t \end{aligned}$$\end{document}

and the log-likelihood of the data

(18) $\begin{array}{c}ll(\beta ,\gamma ,\mu ,\Gamma ,\Sigma )=\sum _{t=1}^T\text{log}\left({\int }_{b_t},L_t,d,b_t\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} ll(\varvec{\beta },\varvec{\gamma }, \varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) = \sum _{t = 1}^{T} \text {log}\left( \int _{\varvec{b}_t} L_t {d}\varvec{b}_t\right) . \end{aligned}$$\end{document}

One problem with the marginal ML approach is that there is no analytical solution for the values of $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} , $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} , $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} that maximize Eq. 18. Rather, numerical approximations have to be used to estimate the model parameters. For these approximations, we use the analytical gradient of the log-likelihood function given by

(19) $\begin{array}{c}\frac{\partial ll}{\partial {\tau }_k}=\sum _{t=1}^T\frac{1}{f\left(n_t\right)}{\int }_{b_t}{L}_t\frac{\partial \text{log}\left(L_t\right)}{\partial {\tau }_k}d{b}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial ll}{\partial \tau _k} = \sum _{t = 1}^{T} \frac{1}{ f (\varvec{n}_{t}) } \int _{\varvec{b}_t} L_t \frac{\partial \text {log}(L_t)}{\partial \tau _k} {d}\varvec{b}_t \end{aligned}$$\end{document}

where $\tau$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} denotes an element of $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} , $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} , $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , or $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} . To implement the gradient, one needs the first derivative of the log-likelihood function for a single person t with regard to a parameter. Note that $\text{log}\left(L_t\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {log}(L_t)$$\end{document} is

(20) $\begin{array}{c}\text{log}f\left(n_t,|\beta ,\gamma ,,b_t\right)+\text{log}f\left(b_t,|\mu ,\Gamma ,\Sigma \right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {log}f\left( \varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_t\right) + \text {log}f \left( \varvec{b}_t|\varvec{\mu }, \varvec{\Gamma }, \varvec{\Sigma }\right) \end{aligned}$$\end{document}

Thus, for a single element ${\sigma }_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{l}$$\end{document} contained in $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} , the derivative is

(21) $\begin{array}{c}\frac{\partial \text{log}{L}_t}{\partial {\sigma }_l}=-\frac{1}{2}tr\left({\Sigma }^{-1},\frac{\partial \Sigma }{\partial {\sigma }_l}\right)+\frac{1}{2}{\left(b_t,-,{\mu }_t\right)}^{{}^{\prime }}{\Sigma }^{-1}\frac{\partial \Sigma }{\partial {\sigma }_l}{\Sigma }^{-1}(b_t-{\mu }_t)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \text {log}L_t}{\partial \sigma _l} = -\frac{1}{2} tr\left( \varvec{\Sigma }^{-1} \frac{\partial \varvec{\Sigma }}{ \partial \sigma _l }\right) + \frac{1}{2} \left( \varvec{b}_{t} - \varvec{\mu }_{t}\right) ^{'} \varvec{\Sigma }^{-1} \frac{\partial \varvec{\Sigma }}{ \partial \sigma _l } \varvec{\Sigma }^{-1} (\varvec{b}_{t} - \varvec{\mu }_{t}) \end{aligned}$$\end{document}

with

(22) $\begin{array}{c}\frac{\partial \Sigma }{\partial {\sigma }_l}=1_l+1_l^{{}^{\prime }}-1_l·1_l^{{}^{\prime }}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \varvec{\Sigma }}{ \partial \sigma _l } = \varvec{1}_{l} + \varvec{1}_{l}^{'} - \varvec{1}_{l} \cdot \varvec{1}_{l}^{'} \end{aligned}$$\end{document}

and $1_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{1}_{l}$$\end{document} is an R x R matrix that contains a 1 in the position of ${\sigma }_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _l$$\end{document} and zeroes in all other positions. For elements in $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} or $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , the derivative is

(23) $\begin{array}{c}\frac{\partial \text{log}{L}_t}{\partial {\mu }_{\mathrm{tl}}}={(b_t-{\mu }_t)}^{{}^{\prime }}{\Sigma }^{-1}\frac{\partial {\mu }_t}{\partial {\mu }_{\mathrm{tl}}}\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \text {log}L_t}{\partial \mu _{tl}} = (\varvec{b}_t - \varvec{\mu }_t)^{'} \varvec{\Sigma }^{-1} \frac{\partial \varvec{\mu }_{t}}{ \partial \mu _{tl} } \end{aligned}$$\end{document}

where

(24) $\begin{array}{c}\frac{{\mu }_t}{\partial {\mu }_l}=1_l\text{and}\frac{{\mu }_t}{\partial {\gamma }_l}=1_l{X}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\varvec{\mu }_{t}}{ \partial \mu _l} = \varvec{1}_l \text { and } \frac{\varvec{\mu }_{t}}{ \partial \gamma _l} = \varvec{1}_l \varvec{X}_t \end{aligned}$$\end{document}

and $1_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{1}_{l}$$\end{document} is an R x 1 vector or an R x p matrix that contains a 1 in the position of the parameter that is being estimated and zeroes in all other positions. Finally, the derivative of $\text{log}\left(L_t\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {log}(L_t)$$\end{document} with regard to the elements in $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} is

(25) $\begin{array}{c}\frac{\partial \text{log}{L}_t}{\partial {\beta }_l}\sum _{t=1}^T\sum _{k=1}^K\sum _{j=1}^{J_k}\frac{n_{\mathrm{kjt}}}{P\left(C_{\mathrm{kj}}\right|{\theta }_t)}\sum _{i=1}^{I_{\mathrm{kj}}}\left(a_{\mathrm{skji}},·,{\left(\frac{\partial {\theta }_{\mathrm{st}}}{\partial {\beta }_l}\right)}^{-1},-,b_{\mathrm{skji}},·,{\left(1,-,\frac{\partial {\theta }_{\mathrm{st}}}{\partial {\beta }_l}\right)}^{-1}\right)·P\left(B_{\mathrm{kji}}\right|{\theta }_t)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \text {log}L_t}{\partial \beta _{l}} \sum _{t=1}^{T} \sum _{k=1}^{K} \sum _{j=1}^{J_{k}} \frac{ n_{kjt} }{ P(C_{kj} | \varvec{\theta }_t ) } \sum _{i = 1}^{I_{kj}} \left[ a_{skji} \cdot \left( \frac{\partial \theta _{st}}{\partial \beta _{l}}\right) ^{-1} - b_{skji} \cdot {\left( 1 - \frac{\partial \theta _{st}}{\partial \beta _{l}}\right) }^{-1} \right] \cdot P(B_{kji}|\varvec{\theta }_t) \end{aligned}$$\end{document}

where $\frac{\partial {\theta }_{\mathrm{st}}}{\partial {\beta }_l}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial \theta _{st}}{\partial \beta _{l}}$$\end{document} is the derivative of the chosen link function with respect to the intercept parameter ${\beta }_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{l}$$\end{document} . The same formula can be used for the parameters in $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} , but the derivative of the link function has to be computed for the element ${\gamma }_{\mathrm{sl}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _{sl}$$\end{document} rather than ${\beta }_l$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{l}$$\end{document} (i.e., $\frac{\partial {\theta }_{\mathrm{st}}}{\partial {\gamma }_{\mathrm{sl}}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\partial \theta _{st}}{\partial \gamma _{sl}}$$\end{document} ).

A second problem with using the marginal ML approach is that the likelihood or log-likelihood function, respectively, and the gradient involve integrals that are not tractable. However, one can approximate the integrals using a number of numerical techniques (see Tuerlinckx et al., Reference Tuerlinckx, Rijmen, Verbeke and De Boeck2006; Nestler, Reference Nestler2020, Reference Nestler2021). In the R package that implements marginal ML estimation methods, we implemented three approaches: the Laplace approximation, the Adaptive Gauss–Hermite Quadrature (AGHQ), and Quasi Monte Carlo (QMC) Sampling.

Laplace Approximation The basic idea behind the Laplace approximation is that it can be used to replace the function within the integral with another function that has a closed-form expression. Imagine that the modes of the random effects $b_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{b}_t$$\end{document} of

(26) $\begin{array}{c}l\left(b_t\right)=f\left(n_t\right|\beta ,\gamma ,b_t)f\left(b_t\right|\mu ,\Gamma ,\Sigma )\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l( \varvec{b}_t) = f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_t) f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) \end{aligned}$$\end{document}

are available for each of the T individuals. One can then show (see, e.g., Pinheiro & Bates, Reference Pinheiro and Bates1995) that the likelihood function given in Eq. 17 can be approximated by

(27) $\begin{array}{c}L(\beta ,\gamma ,\mu ,\Gamma ,\Sigma )\approx \prod _{t=1}^T{\left(2\pi \right)}^{S/2}{\left|\hat{\Omega }\right|}^{1/2}f\left(n_t\right|\beta ,\gamma ,{\hat{b}}_t)f\left({\hat{b}}_t\right|\Sigma ,\mu ,\Gamma )\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L(\varvec{\beta },\varvec{\gamma },\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) \approx \prod _{t = 1}^{T} (2\pi )^{S/2}|\hat{\varvec{\Omega }}|^{1/2} f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \hat{\varvec{b}}_t) f( \hat{\varvec{b}}_t|\varvec{\Sigma },\varvec{\mu },\varvec{\Gamma }) \end{aligned}$$\end{document}

where $\widehat{b_t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}_t}$$\end{document} denotes the modes, and $\hat{\Omega }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{\Omega }}$$\end{document} is the asymptotic covariance matrix of these modes.

A problem with using Eq. 27 is that one has to estimate the modes $\widehat{b_t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}_t}$$\end{document} for all T individuals given the (unknown) parameters contained in $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} , $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} , $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} , $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} , and $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} . In practical implementations, this problem is circumvented by first estimating the modes given the current parameter estimates. Then, the model parameters are estimated by maximizing Eq. 27 given the current modes ${\hat{b}}_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}}_t$$\end{document} . This two-step procedure is repeated until the algorithm converges. The Laplace approximation is very fast (compared with the other methods), and, in the well-known R package lme4 (Bates et al., Reference Bates, Mächler, Bolker and Walker2015), it is the default method for estimating the parameters of a Generalized Linear Mixed Model (GLMM).

Gauss–Hermite Quadrature The basic idea behind quadrature approaches is that they can be used to approximate the numerical value of the integral. When one assumes that the random effects are normally distributed (as we did), one first generates M vectors of size R x 1 of Gauss–Hermite (GH) nodes $x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}$$\end{document} and weights $w$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{w}$$\end{document} . For each of the node vectors (e.g., $x_m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_m$$\end{document} ), one then computes $f\left(n_t\right|\beta ,\gamma ,b_t=x_m)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_t = \varvec{x}_m )$$\end{document} . The integral for person t can then be approximated by a weighted sum

(28) $\begin{array}{c}{L}_t(\beta ,\gamma ,\mu ,\Gamma ,\Sigma )\approx \sum _{n_1=1}^M\cdots \sum _{n_S=1}^{M}f\left(n_t\right|\beta ,\gamma ,b_{n_1\cdots n_R})w_{n_1}\cdots w_{n_R}.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L_{t}(\varvec{\beta },\varvec{\gamma },\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) \approx \sum _{n_{1} = 1}^{M} \cdots \sum _{n_{S} = 1}^{M} f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_{n_{1} \cdots n_{R}} ) \varvec{w}_{n_{1}} \cdots \varvec{w}_{n_{R}} . \end{aligned}$$\end{document}

The main problem with this GH quadrature is that the number of points M increases exponentially with the size of R. For example, when M = 10 and R = 4, there are ${10}^4=1000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4 = 1000$$\end{document} single node vectors. This is problematic because $f\left(n_t\right|\beta ,\gamma ,b_t=x_m)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_t = \varvec{x}_m )$$\end{document} needs to be computed for each vector $x_m$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_m$$\end{document} and each person t. Thus, when M is large, the computational burden associated with approximating the integral is too large to be feasible. For instance, when M = 10 and R = 10, one would need to compute ${10}^{10}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{10}$$\end{document} vector values for each person t. For this reason, we use Adaptive GH quadrature in our implementation (AGHQ, Rabe-Hesketh et al., Reference Rabe-Hesketh, Skrondal and Pickles2002; Tuerlinckx et al., Reference Tuerlinckx, Rijmen, Verbeke and De Boeck2006). In AGHQ, the GH quadrature points $x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}$$\end{document} for each individual are first centered and scaled with the individual’s modes ${\hat{b}}_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}}_t$$\end{document} . These scaled nodes and weights are then used to compute a weighted sum that is similar to the one provided in Eq. 28. By scaling the node vectors with the modes, fewer nodes are required to achieve a precise approximation of the integrals.

Quasi Monte Carlo Integration Even when using AGHQ, the computational burden is too high for high-dimensional random effect distributions. An alternative one could use is Monte Carlo (MC) integration. MC integration rests on the observation that the integral in Eq. 15 can be seen as an expectation of the function $f\left(n_t\right|\beta ,\gamma ,b_t)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_t)$$\end{document} with respect to the random effect distribution $f\left(b_t\right|\mu ,\Gamma ,\Sigma )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma })$$\end{document} :

(29) $\begin{array}{c}{\int }_{b_t}f\left(n_t\right|\beta ,\gamma ,b_t)f\left(b_t\right|\mu ,\Gamma ,\Sigma )d{b}_t=\text{E}\left(f,\left(,n_t,\right|,\beta ,,,\gamma ,,,b_t,)\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{\varvec{b}_t} f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_t) f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) {d}\varvec{b}_t = \text {E} \left[ f(\varvec{n}_{t}|\varvec{\beta }, \varvec{\gamma }, \varvec{b}_t) \right] . \end{aligned}$$\end{document}

One can thus draw M random samples from $f\left(b_t\right|\mu ,\Gamma ,\Sigma )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma })$$\end{document} to approximate the integral (Robert and Casella, Reference Robert and Casella2010) with

(30) $\begin{array}{c}{L}_t(\beta ,\gamma ,\mu ,\Gamma ,\Sigma )\approx \frac{1}{M}\sum _{m=1}^{M}f\left(n_t\right|\beta ,\gamma ,b_{\mathrm{tm}}).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L_{t}(\varvec{\beta },\varvec{\gamma }, \varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) \approx \frac{1}{M} \sum _{m = 1}^{M} f(\varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_{tm}) . \end{aligned}$$\end{document}

where one inserts the m-th random draw for $b_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{b}_{t}$$\end{document} . The advantage of MC integration is that the number of draws M does not have to increase when another random effect is added. However, one disadvantage of the technique is that M must be sufficiently large to be precise. Furthermore, because samples are randomly drawn from the random effect distribution, a (Monte Carlo) sampling error is introduced into the estimation. For these two reasons, here, we use Quasi Monte Carlo (QMC) integration, which builds on deterministic sequences of points instead of randomly drawn points in MC integration (hence the name “Quasi”). There are different ways to generate such sequences (see González et al., Reference González, Tuerlinckx, De Boeck and Cools2006, for an introduction). In accordance with the GLMM literature (e.g., Crowther, Reference Crowther2017; González et al., Reference González, Tuerlinckx, De Boeck and Cools2006), we use Halton sequences in our implementation because smaller numbers of draws M are required to achieve a precise approximation. To further decrease the computational burden, we additionally scale the Halton numbers with $\widehat{b_t}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}_t}$$\end{document} and $\hat{\Omega }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\Omega }$$\end{document} .

3.2. Standard Errors, Goodness-of-Fit Tests, and Random Effects

Once the ML estimates have been determined, the estimates and their corresponding standard errors can be used to compute z statistics and confidence intervals. From standard ML theory, it follows that the parameters are asymptotically normally distributed with a covariance matrix that is obtained by calculating the inverse of the information matrix. This matrix is the negative of the matrix of second derivatives that is given by

(31) $\begin{array}{c}\frac{\partial ll}{\partial {\tau }_k{\tau }_j}=\sum _{t=1}^T\frac{1}{f{\left(n_t\right)}^2}\left(f,\left(n_t\right),\frac{\partial d_t}{\partial {\tau }_j},-,d_t,d_t^T\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial ll}{\partial \tau _k\tau _j} = \sum _{t = 1}^{T} \frac{1}{ f(\varvec{n}_{t})^2 } \left[ f(\varvec{n}_{t}) \frac{\partial d_t}{\partial \tau _j} - d_{t} d^{T}_{t} \right] \end{aligned}$$\end{document}

where

(32) $\begin{array}{c}{d}_t={\int }_{b_t}{L}_t\frac{\partial \text{log}\left(L_t\right)}{\partial {\tau }_k}d{b}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} d_{t} = \int _{\varvec{b}_t} L_t \frac{\partial \text {log}(L_t)}{\partial \tau _k} {d}\varvec{b}_t \end{aligned}$$\end{document}

and

(33) $\begin{array}{c}\frac{\partial d_t}{\partial {\tau }_j}={\int }_{b_t}\left(\frac{{\partial }^2\text{log}\left(L_t\right)}{\partial {\tau }_k{\tau }_j},+,\frac{\partial \text{log}\left(L_t\right)}{\partial {\tau }_k},\frac{\partial \text{log}\left(L_t\right)}{\partial {\tau }_j}\right)L_{t}d{b}_t\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial d_t}{\partial \tau _j} = \int _{\varvec{b}_t} \left[ \frac{\partial ^2 \text {log}(L_t)}{\partial \tau _k \tau _j} + \frac{\partial \text {log}(L_t)}{\partial \tau _k} \frac{\partial \text {log}(L_t)}{\partial \tau _j} \right] L_t {d}\varvec{b}_t \end{aligned}$$\end{document}

Again, due to the integrals involved in Eqs. 32 and 33, the matrix of second derivatives is hard to compute. In our implementation of the ML approach users can therefore choose whether standard errors are based on a numerical approximation of the Hessian with finite-difference methods using the analytical gradient (see Eq. 19) or on the exact hessian computed with Eq. 31.

Furthermore, before researchers interpret MPT parameter estimates, they need to show that their model actually fits the observed data. In addition, researchers very often want to compare the fit of a current model with the fit of a more restricted model that imposes psychologically motivated constraints on the parameters (e.g., equality constraints or parameter fixations). Goodness-of-fit tests and model comparisons can both be conducted by means of a likelihood ratio test

(34) $\begin{array}{c}LR=-2·(l{l}_r-l{l}_u)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} LR = -2\cdot ( ll_r - ll_u ) \end{aligned}$$\end{document}

where $l{l}_r$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ll_r$$\end{document} ( $l{l}_u$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ll_u$$\end{document} ) is the log-likelihood value of the restricted (unrestricted) model. Under certain regularity conditions (see Reed & Cressie, Reference Reed and Cressie1988), the LR test statistic asymptotically follows a central ${\chi }^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2$$\end{document} distribution with degrees of freedom equal to the difference between the parameters of the unrestricted and restricted models if the more restricted model actually holds. The likelihood ratio test can be used to compare the fit of two models that differ in the mean structure (e.g., u = a in the pair-clustering model) or in the covariance structure (e.g., ${\sigma }_{\mathrm{ua}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{ua}$$\end{document} = 0).

The likelihood ratio test is based on the assumption that the models being compared are nested. For non-nested model comparisons, we suggest using the Akaike or the Bayesian Information Criterion (AIC or BIC, respectively):

(35) $\begin{array}{cc}\mathrm{AIC}& =-2·ll+2·df\\ BIC& =-2·ll+\text{log}\left(n\right)·df,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} AIC&= -2 \cdot ll + 2 \cdot df \nonumber \\ BIC&= -2 \cdot ll + \text {log}(n)\cdot df, \end{aligned}$$\end{document}

where $n={\sum }_{t=1}^T{n}_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = \sum _{t=1}^{T} n_t$$\end{document} and $df=R+(S-R)+p_{\Gamma }+p_{\gamma }+\text{dim}\left(\Sigma \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$df = R + (S - R) + p_{\varvec{\Gamma }} + p_{\varvec{\gamma }} + \text {dim}(\varvec{\Sigma })$$\end{document} . Here, R is the number of estimated cognitive process parameters set to be random (i.e., parameters in $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\mu }$$\end{document} ), $S-R$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S - R$$\end{document} is the number of fixed cognitive process parameters (i.e., parameters in $\beta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\beta }$$\end{document} ), $p_{\Gamma }+p_{\gamma }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\varvec{\Gamma }} + p_{\varvec{\gamma }}$$\end{document} gives the total number of to-be-estimated weights of the person-level covariates in the random and the fixed effects part of the model (i.e., parameters in $\Gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Gamma }$$\end{document} and $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\gamma }$$\end{document} ), respectively, and $\text{dim}\left(\Sigma \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {dim}(\varvec{\Sigma })$$\end{document} represents the number of estimated covariance parameters (see Bates et al., Reference Bates, Mächler, Bolker and Walker2015).

Finally, it can also be of interest to obtain estimates of the individual random effects for each participant t. As a random effects estimator, ${\hat{b}}_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\varvec{b}}_t$$\end{document} , it is possible to use a participant’s mode, which can be obtained by maximizing Eq. 26, while treating the final model parameter estimates as fixed. Fortunately, these values are already required when estimating the model so that they do not have to be estimated again. Another choice would be the empirical Bayes estimator (Bock and Aitken, Reference Bock and Aitken1981):

(36) $\begin{array}{c}{\hat{b}}_t=\prod _{t=1}^T\frac{1}{f\left(n_t\right)}{\int }_{b_t}{b}_t\text{log}\left(f\left(n_t\right|\beta ,\gamma ,b_t)\right)f\left(b_t\right|\mu ,\Gamma ,\Sigma )d{b}_t.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{b}_{t} = \prod _{t = 1}^{T} \frac{1}{ f(\varvec{n}_{t}) } \int _{\varvec{b}_t} \varvec{b}_t \text {log}(f(\varvec{n}_{t}|\varvec{\beta },\varvec{\gamma }, \varvec{b}_t)) f(\varvec{b}_t|\varvec{\mu },\varvec{\Gamma },\varvec{\Sigma }) {d}\varvec{b}_t . \end{aligned}$$\end{document}

In this approach, we approximate the integrals using one of the integral approximation methods described above.

4. Simulation Study

We performed a simulation study to assess the frequentist properties of the suggested ML estimation approaches and also compare it with the performance of a Bayesian approach. Specifically, we examined the effect of the number of participants and the number of responses per participant on the bias of the parameter estimates and the coverage rate of the corresponding confidence or credibility intervals, respectively. In addition, we examined for the ML estimator how the number of quadrature points in the AGHQ approach and the size of the Halton sequence in the QMC method affect the adequacy of the model parameter estimates.

Population Model and Simulation Conditions. We used the pair-clustering model for this simulation. This MPT model comprises two category systems, including four and two categories, respectively, with response probabilities described by four process parameters c, r, u, and a. In our simulation study, the population covariance matrix of these four parameters (describing how they vary across participants) was set to

(37) $\begin{array}{c}\Sigma =\left(\begin{array}{cccc}0.50& 0.08& 0.04& 0.00\\ 0.08& 0.35& 0.03& 0.00\\ 0.04& 0.03& 0.20& 0.07\\ 0.00& 0.00& 0.07& 0.20\\ \end{array}\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\Sigma } = \begin{pmatrix} 0.50 &{}\quad 0.08 &{}\quad 0.04 &{}\quad 0.00 \\ 0.08 &{}\quad 0.35 &{}\quad 0.03 &{}\quad 0.00 \\ 0.04 &{}\quad 0.03 &{}\quad 0.20 &{}\quad 0.07 \\ 0.00 &{}\quad 0.00 &{}\quad 0.07 &{}\quad 0.20 \\ \end{pmatrix} \end{aligned}$$\end{document}

where the non-zero covariance terms reflect correlations of 0.30, 0.20, and 0.10. The mean values of the parameters were set to ${\mu }_c=0.50$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{c} = 0.50$$\end{document} , ${\mu }_r=0.40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{r} = 0.40$$\end{document} , ${\mu }_u=0.25$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{u} = 0.25$$\end{document} , and ${\mu }_a=0.15$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{a} = 0.15$$\end{document} (see Batchelder & Riefer, Reference Batchelder and Riefer1986).

We manipulated the number of simulated participants (75 vs. 125) and the number of simulated responses per participant (25 vs. 75 vs. 125). In accordance with Batchelder and Riefer (Reference Batchelder and Riefer1986), about 80% of the responses were assigned to the first category system (i.e., 20 vs. 60 vs. 100), leaving 20% for the second system (i.e., 5 vs. 15 vs. 25). The R package mvtnorm and the function rmultinorm were used to generate the samples. We drew 500 samples from the population for each of the six simulation conditions.

Estimators All the functions that were required to estimate the parameters with ML were implemented in R (R Core Team, 2020; a working version of the package can be downloaded from https://osf.io/w97m5/). To examine the properties of the ML estimation procedures, the number of quadrature points in the AGHQ method was set to 3, 4, or 5, resulting in $4^3=64$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^3 = 64$$\end{document} , $4^4=256$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^4 = 256$$\end{document} or $4^5=1,024$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^5 = 1,024$$\end{document} node vectors per participant. For the QMC method, the size of the Halton sequence was set to 500, 1000, or 2000. To obtain Bayesian estimates, we employed the TreeBUGS package (Heck et al., Reference Heck, Arnold and Arnold2018a). TreeBUGS uses the JAGS-MCMC sampler (Plummer, Reference Plummer2003) to approximate the posterior distribution of the model’s parameters. For each replication, we fitted the model with the traitMPT function using the default settings. That is, weakly informative priors were specified for all parameters (i.e., normal distributions for the means and a scaled Wishart distributions for the covariance matrix). Furthermore, three chains of 20,000 samples were generated from the posterior distributions, whereby the first 2000 samples were discarded for parameter estimation (i.e., burn-in period). Since TreeBUGS by default provides the means of the posterior distributions as parameter estimates, we decided to used them as the Bayes estimates.

Dependent Measures We used the relative bias (RB) of the parameter estimates and the coverage rate (CR) to investigate the statistical performances of the ML approaches and the Bayes estimator. For the relative bias, we first computed the average parameter estimate in a simulation condition. We then computed the difference between this average and the true parameter and thereafter divided the difference by the true parameter. We consider relative biases below 10% to be acceptable, biases of 10–20% to be substantial, and biases above 20% to be unacceptable (e.g., Forero et al., Reference Forero, Maydeu-Olivares and Gallardo-Pujol2009; Morris et al., Reference Morris, White and Crowther2019). For ML, the confidence interval with the standard errorFootnote 1 of an estimate was computed in each replication to determine the observed coverage of the 95% confidence intervals. The coverage was then coded 1 if the true parameter value was included in the interval and 0 if the true parameter was not. We used the same approach to determine the observed coverage of the 95% credibility intervals, but employed the 95% credibility intervals as provided by TreeBUGS as the basis for the coding.

Table 1 Relative frequencies of converged replications (CR) and average computation time (in s), depending on the estimator, the type of ML approximation method, the number of individuals T, and the number of responses N per individual.

AGHQ = Adaptive Gauss–Hermit quadrature with 3, 4, or 5 nodes; QMC = Quasi Monte Carlo integration with 500, 1000, or 2000 points. In case of Bayes, convergence rates were calculated using $\hat{R}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}$$\end{document} values. Computation times were determined on an Intel Core i7-6700 with four cores and 16GB RAM.

Results We first examined the percentage of samples in which the estimation algorithm converged for the two estimators, the different approximation methods, the different numbers of simulated participants, and the different numbers of simulated responses per participant. In case of ML, a sample was counted as converged when there were neither inadmissible estimates nor undefined standard errors in the final solution. For Bayes, judging the convergence is a more difficult issue (Hoff, Reference Hoff2009). Here, we decided to use $\hat{R}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{R}$$\end{document} > 1.05 as the criterion, because it can be used well in simulation studies and it is suggested in the literature (e.g., Lynch, Reference Lynch2007). However, we acknowledge that a more or less stringent cutoff may lead to different results and that this should be taken into account in the interpretation of our findings.

As can be seen in Table 1, convergence rates for both estimators increased with the number of participants and the number of responses. In case of ML, the larger the number of points used by a method to approximate the integrals, the better the convergence of the ML estimator. When the number of responses was 75 or 125, convergence rates were acceptable for the AGHQ and QMC methods. In conditions with 25 responses, convergence rates were highest for conditions with the largest size of points. In this case, convergence rates were also similar to the convergence rate of the Bayes estimator. With regard to the Laplace approximation, we found that convergence rates were always below 70% (and even near 0% when number of responses was 25). Further analyses showed that all convergence failures occurred because some (or all) of the standard errors were undefined for the final estimates. These estimates were also very biased. We think that this bias can be explained by noting that the precision of the Laplace approximation depends on whether the shape of the function that is being integrated resembles a multivariate normal distribution. Hence, the method does not perform well for highly non-normal cases (e.g., when the responses are Bernoulli distributed; see Engel, Reference Engel1998), and we suspect that similar performance problems occur for the MPT model. Finally, Table 1 also contains the average run times of the different methods. We note that comparing the computation times across ML and Bayes is difficult, because they depend on how the approaches are implemented in R (i.e., using C++ in the background or parallelization), how many chains are generated etc. In case of Bayes, all methods were slower for larger numbers of participants. Similarly, ML methods became slower the larger numbers of points used to approximate the integrals and the larger the number of participants.

Table 2 Relative bias of parameter estimates in percent, depending on the approximation method, the number of individuals T, and the number of responses N per individual.

AGHQ = Adaptive Gauss–Hermit quadrature with 3, 4, or 5 nodes; QMC = Quasi Monte Carlo integration with 500, 1000, or 2000 points.

In the following, we drop the Laplace approximation from further consideration when we discuss the precision of the ML estimates (i.e., RB) and the confidence intervals (i.e., CR). Furthermore, to facilitate the interpretation of the results, we decided to average the results for the indices per parameter group (i.e., for the cognitive process parameter mean values c, r, u, and a, called $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ; the variance parameters in $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} , termed ${\sigma }_b^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{b}^{2}$$\end{document} ; and the covariance parameters in $\Sigma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }$$\end{document} , termed ${\sigma }_{\mathrm{bb}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{bb}$$\end{document} ). The resulting values for the relative bias are displayed in Table 2. When we consider relative biases below 10% to be acceptable and biases of 10–20% to be substantial (Forero et al., Reference Forero, Maydeu-Olivares and Gallardo-Pujol2009; Morris et al., Reference Morris, White and Crowther2019), we found for both ML methods that relative biases decreased as the numbers of node vectors, sample sizes, and numbers of responses per participant increased. Importantly, relative biases were generally low and acceptable in all simulation conditions for the cognitive process parameters and the variance parameters. For the latter, however, biases were somewhat larger but still acceptable in case of AGHQ when the number of responses was 25 and the number of persons was 75. For the covariance parameters, the relative biases were larger and more substantial the smaller the number of points used and the smaller the number of responses. For the Bayes estimator, the relative bias was negligible for the cognitive process parameters and the variance parameters in all conditions. However, replicating the results of Klauer (Reference Klauer2010), the covariance parameters were unacceptably and substantially biased when the number of responses was 25 or 75. When the number of responses was 125, the relative biases were still large but acceptable.

Table 3 Coverage rate of parameter estimates, depending on the approximation method, the number of individuals T, and the number of responses N per individual.

AGHQ = Adaptive Gauss–Hermit quadrature with 3, 4, or 5 nodes, QMC = Quasi Monte Carlo integration with 500, 1000, or 2000 points.

With regard to the CR (see Table 3), we found that for all three types of parameters, the CR moved closer to the nominal value as the sample size and the number of responses increased. For AGHQ, the CR was near its nominal value when the number of points was $4^4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^4$$\end{document} and the number of responses at least 75. When the number of responses is 25, the CR was near the nominal value when $4^5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4^5$$\end{document} points were used. A similar pattern of results was observed for the QMC approximation, although the amount of undercoverage was greater than for the AGHQ approximation. Furthermore, when the number of responses was 25, undercoverage occurred irrespective of how many points were investigated in our simulation. In case of Bayes, we found that the coverage was nominal for the cognitive process parameters and the variance parameters. However, overcoverage occurred for the covariance parameters although this tendency disappeared the larger the number of responses.

To summarize, the simulation study reveals that the AGHQ method works better than the QMC method when the number of nodes is at least four. Relative biases were low even for small sample sizes and few responses per participant and confidence interval estimates were close to nominal. When the number of responses was 75, the QMC also provided acceptable results, at least when the number of points was 1000. Finally, when the number of responses was small (i.e., R = 25), AGHQ with five nodes yield estimates that are at least as good as the estimates of a Bayesian approach.

5. Illustrative Example

To illustrate the proposed ML approach, we analyzed neuropsychological data from Schilken (Reference Schilken1998), who employed the pair-clustering model to analyze and compare the memory performance of 22 epileptic patients with a right-temporal focus of epileptic seizures, 21 epileptic patients with a left-temporal focus, and 20 healthy controls matched with respect to age and intelligence. Each participant learned word lists consisting of 10 semantically strongly related word pairs (e.g., armchair-sofa) and 5 singletons that were not related to other words in the list (e.g., sunflower). All words were comparable in terms of difficulty and memorizability when considered in isolation. In addition, three words were added to the beginning and another three words to the end of the word list to absorb primacy and recency effects in free recall. Because these primacy and recency buffer words were excluded from further analysis, $N=15$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 15$$\end{document} responses per participant and study-test cycle remained for analysis-10 responses for the pairs and 5 for the singletons. The studying of the word list and the subsequent free recall test were repeated 6 times to examine learning effects on the c, r, and u parameters of the pair-clustering model, resulting in a total of six study-test trials per participant.

The same data set was also analyzed by Klauer (Reference Klauer2006, Reference Klauer2010), thus providing us with the opportunity to compare our results with Klauer’s, which were based on the Bayesian Latent-Trait model (cf. Klauer, 2010, pp. 86/87). To maximize heterogeneity between participants, we followed the procedure outlined by Klauer (Reference Klauer2006, Reference Klauer2010) and analyzed all $T=22+21+20=63$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 22 + 21 + 20 = 63$$\end{document} participants conjointly in our first model. Also following Klauer’s guidelines, we restricted our attention to the first two study-test cycles for each participant. This left us with $T=63$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 63$$\end{document} participants, $K=4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 4$$\end{document} category systems (i.e., those for word pairs and singletons in Trials 1 and 2) with $J_1=J_3=4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_1 = J_3 = 4$$\end{document} and $J_2=J_4=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2 = J_4 = 2$$\end{document} categories (for word pairs and singletons, respectively), and $N_1=N_3=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1 = N_3 = 10$$\end{document} as well as $N_2=N_4=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_2 = N_4 = 5$$\end{document} responses per participant within the four category systems. Finally, again in line with Klauer (Reference Klauer2010)’s suggestions, we imposed the restriction that unclustered words in a pair must match the singletons in terms of trial-specific storage and retrieval probabilities (i.e., $u=a$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u = a$$\end{document} ). Hence, there were three parameters to estimate per trial, $c^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(1)}$$\end{document} , $r^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(1)}$$\end{document} , $u^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{(1)}$$\end{document} for Trial 1 and $c^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(2)}$$\end{document} , $r^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(2)}$$\end{document} , $u^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{(2)}$$\end{document} for Trial 2.

Table 4 Parameter estimates of the mean structure for the illustrative data example.

Parameters were obtained with AGHQ with 4 nodes. Est = Parameter estimates; CI = confidence interval; TE = probability-transformed estimates; Bayes = Bayesian Estimates as reported in Table 2 in Klauer (Reference Klauer2010). D1 = First dummy variable, that is, 1 for right-temporal epileptic patients (all other participants are coded zero). D2 = Second dummy variable, that is, 1 for left-temporal epileptic patients (all other participants are coded zero). ${\Delta }_{D.}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _{D.}$$\end{document} is the estimate of the parameter for the group coded in the respective dummy variable (i.e., ${\delta }_{.,D_{.}}+{\mu }_{.}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{.,D_{.}} + \mu _{.}$$\end{document} ). For model 2, probability transformed estimates (column TE) indicate parameter means in the three groups (control group, right-temporal epileptics, left-temporal epileptics, respectively).

To test the goodness of fit, we estimated the model that we just specified and a more general model with four parameters $c^{\left(d\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(d)}$$\end{document} , $r^{\left(d\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(d)}$$\end{document} , $u^{\left(d\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{(d)}$$\end{document} , and $a^{\left(d\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^{(d)}$$\end{document} per trial. For both models, we used the AGHQ approximation method with 4 nodes to fit the two models. The LR test statistic comparing the original and the more general model was LR = 9.02, which is not significantly different from zero, ${\chi }_{\mathrm{crit}}^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^{{2}}_{{crit}}$$\end{document} = 30.1, p =.98, df = 19. Table 4 shows the parameters of the mean structure (i.e., $\mu$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} ) in the restricted model. As can be seen, parameter values increased from Trial 1 to Trial 2, and the probability-transformed parameters were very similar to the probability-transformed parameters reported in Klauer (Reference Klauer2010). Variance and correlation parameters were also quite similar across the two approaches (see Table 5). A notable exception was the variance of $r^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(1)}$$\end{document} , where the ML estimate was considerably larger than the Bayesian estimates. Also, the correlations involving this parameter were smaller for ML compared with Bayes.

Table 5 Parameter estimates of the covariance structure for the illustrative example with real data.

The first two columns present variance estimates based on the AGHQ method with 4 nodes (Est) and the corresponding Bayesian estimates (Bayes) as reported by Klauer (Reference Klauer2010, Table 2), respectively. The correlation matrix displays the Bayesian estimates of Klauer (Reference Klauer2010) above the diagonal and the corresponding AGHQ estimates using 4 nodes below the diagonal.

We decided to go one step further than Klauer (Reference Klauer2010) and also analyze effects of the clinical group on the parameter estimates. For this purpose, we added two dummy variables as covariates to our model, the first one coded “1” for the right-temporal epileptic patients (patients in the two other groups were coded 0) and the second dummy variable coded “1” for the left-temporal epileptic patients (again, all other patients were coded “0”). This model provided us with the opportunity to estimate the average (negative) effect of each clinical group relative to the control group on each of the model parameters, that is, $c^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(1)}$$\end{document} , $r^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(1)}$$\end{document} , $u^{\left(1\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{(1)}$$\end{document} for Trial 1 and $c^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(2)}$$\end{document} , $r^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^{(2)}$$\end{document} , $u^{\left(2\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{(2)}$$\end{document} for Trial 2. The results are also shown in Table 4. The results indicate that estimates for the $c^{(.)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^{(.)}$$\end{document} parameters were lower for epileptic patients compared with control patients whereas this pattern is less clear for the remaining parameters.