Description Usage Arguments Details Value Note Note Author(s) References See Also Examples
fit.mpt.old
function fits binary multinomial processing tree models (MPT models; e.g., Riefer & Batchelder, 1988). However, this function is an old version using the LBFGSB optimization routine. See fit.mpt
for the new version.
1 2 3 4 5 6 7 8 9 10 11 12 13 14  fit.mpt.old(
data,
model.filename,
restrictions.filename = NULL,
n.optim = 5,
fia = NULL,
ci = 95,
starting.values = NULL,
output = c("standard", "fia", "full"),
reparam.ineq = TRUE,
sort.param = TRUE,
model.type = c("easy", "eqn", "eqn2"),
multicore = c("none", "individual", "n.optim"), sfInit = FALSE, nCPU = 2
)

data 
Either a numeric 
model.filename 
A character 
restrictions.filename 

n.optim 
Number of optimization runs. Can be parallelized via 
fia 
Number of random samples to be drawn in the Monte Carlo algorithm to estimate the Fisher Information Approximation (FIA), a minimum description length based measure of model complexity (see Wu, Myung & Batchelder, 2010). The default is 
ci 
A scalar corresponding to the size of the confidence intervals for the parameter estimates. Default is 95 which corresponds to 95% confidence intervals. 
starting.values 
A 
output 
If "fia", 
reparam.ineq 
Logical. Indicates whether or not inequality restrictions (when present in the model file) should be enforced while fitting the model. If 
sort.param 
Logical. If TRUE, parameters are alphabetically sorted in the parameter table. If FALSE, the first parameters in the parameter table are the nonrestricted ones, followed by the restricted parameters. Default is TRUE. 
model.type 
Character vector specifying whether the model file is formatted in the easy way ( 
multicore 
Character vector. If not 
sfInit 
Logical. Relevant if 
nCPU 
Scalar. Only relevant if 
There is a new version of this function using nlminb
and the analytically derived gradient and hessian. See fit.mpt
. We recommend using the new version fit.mpt
, only use this version if you are sure on what to do.
The model file is either of the easy format or the "classical" EQN format (see below).
In the easy format (the default) the model file contains all trees of the model. Trees are separated by at least one empty line. Everything to the right of a hash (#) is ignored (this behavior is new since version 0.9.2). Lines starting with a # are treated as empty. Each line in each tree corresponds to all branches of this tree (concatenated by a +) that correspond to one of the possible response categories. The position of each line must correspond to the position of this response category in the data object (for multiindividual fit to the respective column).
The difference between both types of EQN format ("eqn"
or"eqn2"
) is the way the first line of the model file is treated. If model.file
is set to "eqn"
, MPTinR
will ignore the first line of the model file and will read the rest of the file (as does multiTree; Moshagen, 2010). If model.file
is set to "eqn2"
MPTinR will only read as many lines as indicated in the first line of the EQN model file (as does e.g., HMMTree; Stahl & Klauer, 2007). As default fit.mpt
expects the easy format, but if the filename ends with .eqn or .EQN and model.type
is "easy"
, model.type
is set to "eqn"
For the EQN format consult one of the corresponding papers (see e.g., Moshagen, 2010; Stahl & Klauer, 2007). The positions in the data object (number of column for multiindividual fit) must correspond to the category number in the EQN file.
Note that names of parameters in the model file should not start with hank.
. Variables with these names can lead to unforeseen problems as variables starting with these letters are internally used.
The restrictions file may contain (sequential) equality (i.e., =) and inequality (i.e., <) restrictions and must adhere to the following rules:
1. Inequalities first.
2. If a variable appears in an inequality restriction, it can not be on the LHS of any further restriction.
3. If a variable appears on RHS of an equality restriction, it can not appear on LHS of an equality restriction.
Note that only "<" is supported as inequality operator but not ">"!
Examples of restrictions are (the following could all appear in one restrictions file):
D1 < D2 < D3
D4 = D3
B1 = B3 = 0.3333
X4 = X5 = D3
Restrictions file may contain comments (i.e., everything to the right of a # will be ignored; new behavior since version 0.9.2)
For equality restrictions, the equality restricted parameters are simply exchanged with their restrictions before the fitting.
For inequality restricted parameters, the model is reparameterized so that only the rightmost parameter of an inequality restriction remains the original parameter. Each instance of the other parameters in this restriction is replaced by the product of the rightmost parameter and dummy parameters (see Knapp & Batchelder, 2004). This procedure (which is equivalent to method A described in Knapp & Batchelder, 2004) leads to an equivalent model (although the binary MPT structure is not apparent in the resulting equations).
To prohibit this reparameterization (i.e., if the inequality restrictions hold without reparameterization), you can set reparam.ineq
to FALSE
. This can be useful for obtaining the FIA (see examples in Wu, Myung, & Batchelder, 2010).
The fitting/optimization is achieved via optim
's LBFGSB method by Byrd et al. (1995) with random starting values. To avoid local minima it is useful to set n.optim
> 1. If n.optim
> 1, the summary
of the vector containing the LogLikelihood values returned by each run of optim
is added to the output (to check whether local minima were present). If the model is rather big, n.optim
> 1 can be slow.
To obtain a measure of the model's complexity beyond the number of parameters (and taking inequality restrictions into account), set fia
to a (reasonably high) scalar integer (i.e., a number). Then, fit.mpt
will obtain the Fisher information approximation (FIA), a minimum description based measure of model complexity, using the algorithm provided by Wu, Myung, & Batchelder (2010a, 2010b) ported from Matlab to R. When performing modelselection, this measure is superior to other methods such as the Akaike information criterion (AIC) or Bayesian information criterion (BIC) which basically only take the number of parameters into account.
To get the FIA, fit.mpt.old
performs the following steps:
1. The representation of the model as equations is transformed into the string representation of the model in the contextfree language of MPT models (LBMPT; Purdy & Batchelder, 2009). For this step to be successful it is absolutely necessary that the equations representing the model perfectly map the tree structure of the MPT. That is, the model file is only allowed to contain parameters, their negations (e.g., Dn
and (1  Dn)
) and the operators + and *, but nothing else. Simplifications of the equations will seriously distort this step. This step is achieved by make.mpt.cf
.
2. The context free representation of the model is then fed into the MCMC function computing the FIA (the port of BMPTFIA provided by Wu, Myung & Batchelder (2010a), see bmpt.fia
).
(Actually, both steps are achieved by a call to get.mpt.fia
)
Once again: If one wants to compute the FIA, it is absolutely necessary, that the representation of the model via equations in the model file exactly maps on the structure of the binary MPT (see make.mpt.cf
for more details).
Confidence intervals (CI) are based on the observed Hessian matrix returned by the minimization algorithm (optim
).
For inequality restricted parameters, the CIs are computed using the parameter estimates' variance bounds (see Baldi & Batchelder, 2003; especially equation 19). Note that these bounds represent the "worst case scenario" variances, and can lead to CIs outside parameter boundaries if the set of inequalities is large and/or the variances for the reparameterized model are large (Note that CIs for nonrestricted parameters can be outside the parameter boundaries as well due to large variances).
To set the starting values for the fitting process (e.g., to avoid local minima) one can set starting.values
to a vector of length 2. Then, starting values are randomly drawn from a uniform distribution from starting.values[1]
to starting.values[2]
.
Furthermore, one can specify the starting values individually by supplying a vector with the same length as the number of parameters. Starting values must be ordered according to the alphabetical order of the parameters. Use check.mpt
for a function that returns the alphabetical order of the parameters. If one specifies the starting values like that, n.optim
will be set to 1 as all other values would not make any sense (the optimization routine will produce identical results with identical starting values).
Multicore fitting is achieved via the snowfall
package and needs to be initialized via sfInit
. As initialization needs some time, you can either initialize multicore facilities yourself using sfInit()
and setting the sfInit
argument to FALSE
(the default) or let MPTinR initialize multicore facilities by setting the sfInit
argument to TRUE
. The former is recommended as initializing snowfall
takes some time and only needs to be done once if you run fit.mpt.old
multiple times. If there are any problems with multicore fitting, first try to initialize snowfall
outside MPTinR (e.g., sfInit( parallel=TRUE, cpus=2 )
). If this does not work, the problem is not related to MPTinR but to snowfall (for support and references visit: http://www.imbi.unifreiburg.de/parallel/).
Note that you need to close snowfall via sfStop()
after using MPTinR.
fit.model()
is essentially a copy of fit.mpt.old
that allows the user to specify the upper and lower bounds of the parameters. This function can be used to fit other models than MPT models that can be described in a model file. That is, the model file can contain any type of valid R expressions including R functions (potentially selfwritten) visible in the global environment (i.e., not only +, *, and  as operators). Currently fit.model
should be viewed as experimental.
Note that fit.model()
is usually slower than fit.mpt.old
as there are some more checks in the critical function calculating the likelihood of the model.
The lower.bound
and upper.bound
needs to be of length 1 or equal to the number of free
parameters. If length > 1, parameters are mapped to the bounds in alphabetic order of the parameters. Use check.mpt
to obtain the alphabetical order of parameters for your model.
While it should be possible to specify equality or fixed restrictions it will probably lead to unforeseen consequences to specify inequality restrictions for nonMPT models.
For individual fits (i.e., data
is a vector
) a list
containing one or more of the following components from the best fitting model:
goodness.of.fit 
A 
information.criteria 
A 
model.info 
A 
parameters 
A data.frame containing the parameter estimates and corresponding confidence intervals. If a restriction file was present, the restricted parameters are marked. 
data 
A 
For multiindividual fits (i.e., data
is a matrix
or data.frame
) a list
with similar elements, but the following differences.
The first elements, goodness.of.fit
, information.criteria
, and model.info
, contain the same information as for individual fits, but each are lists
with three elements containing the respective values for: each individual in the list element individual
, the sum of the individual values in the list element sum
, and the values corresponding to the fit for the aggregated data in the list element aggregated
.
parameters
is a list containing:
individual 
A 3dimensional array containing the parameter estimates ([,1,]), confidence intervals [,2:3,], and, if restrictions not 
mean 
A 
aggregated 
A data.frame containing the parameter estimates and corresponding confidence intervals for the aggregated data. If a restriction file was present, the restricted parameters are marked. 
The element data
contains two matrices, one with the observed
, and one with the predicted
data.
If n.optim
> 1, the summary
of the vector (matrix for multiindividual fit) containing the LogLikelihood values returned by each run of optim
is added to the output.
When using R (>= 2.13.0) compiling fit.mpt.old
using compilers
cmpfun
can significantly improve fitting time.
There may be several warnings fit.mpt.old
throws while fitting MPT models. Most of them are not problematic and related to matrix operations needed for confidence intervals. Examples:
In sqrt(var.params) : NaNs produced
In sqrt(min(var.bound.tmp)) : NaNs produced
These warnings are not critical.
Other warnings may relate to the optimization routine (e.g., Optimization routine [...] did not converge successfully
).
In these cases it is recommended to rerun fit.mpt.old
to check if the results are stable.
All (model or restriction) files should end with an empty line, otherwise a warning will be shown.
Henrik Singmann and David Kellen with help from Karl Christoph Klauer and Fabian Hoelzenbein.
Baldi, P. & Batchelder, W. H. (2003). Bounds on variances of estimators for multinomial processing tree models. Journal of Mathematical Psychology, 47, 467470.
Broeder, A., & Schuetz, J. (2009). Recognition ROCs are curvilinearor are they? On premature arguments against the twohighthreshold model of recognition. Journal of Experimental Psychology: Learning, Memory, and Cognition, 35(3), 587. doi:10.1037/a0015279
Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM J. Scientific Computing, 16, 11901208.
Knapp, B. R., & Batchelder, W. H. (2004). Representing parametric order constraints in multitrial applications of multinomial processing tree models. Journal of Mathematical Psychology, 48, 215229.
Moshagen, M. (2010). multiTree: A computer program for the analysis of multinomial processing tree models. Behavior Research Methods, 42, 4254.
Purdy, B. P., & Batchelder, W. H. (2009). A contextfree language for binary multinomial processing tree models. Journal of Mathematical Psychology, 53, 547561.
Riefer, D. M., & Batchelder, W. H. (1988). Multinomial modeling and the measurement of cognitive processes. Psychological Review, 95, 318339.
Stahl, C. & Klauer, K. C. (2007). HMMTree: A computer program for latentclass hierarchical multinomial processing tree models. Behavior Research Methods, 39, 267 273.
Wu, H., Myung, J.I., & Batchelder, W.H. (2010a). Minimum description length model selection of multinomial processing tree models. Psychonomic Bulletin & Review, 17, 275286.
Wu, H., Myung, J.I., & Batchelder, W.H. (2010b). On the minimum description length complexity of multinomial processing trees. Journal of Mathematical Psychology, 54, 291303.
fit.mpt
for the current and recommended function fot fitting MPTs
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156  ## Not run:
# The first example fits the MPT model presented in Riefer and Batchelder (1988, Figure 1)
# to the data presented in Riefer and Batchelder (1988, Table 1)
# Note that Riefer and Batchelder (1988, pp. 328) did some hypotheses tests, that are not done here.
# Rather, we use each condition (i.e., row in Table 1) as a different individual.
# We try to use n.optim = 1 here, but this can lead to local minima
# In general we recommend to set n.optim >= 5
# load the data
data(rb.fig1.data)
#get the character string with the position of the model:
model1 < system.file("extdata", "rb.fig1.model", package = "MPTinR")
model1.eqn < system.file("extdata", "rb.fig1.model.eqn", package = "MPTinR")
# just fit the first "individual":
fit.mpt.old(rb.fig1.data[1,], model1, n.optim = 1)
#fit all "individuals":
fit.mpt.old(rb.fig1.data, model1, n.optim = 1)
#fit all "individuals" using the .EQN model file:
fit.mpt.old(rb.fig1.data, model1.eqn, n.optim = 1)
# The second example fits the MPT model presented in Riefer and Batchelder (1988, Figure 2)
# to the data presented in Riefer and Batchelder (1988, Table 3)
# First, the model without restrictions is fitted: ref.model
# Next, the model with all r set equal is fitted: r.equal
# Then, the model with all c set equal is fitted: c.equal
# Finally, the inferential tests reported by Riefer & Batchelder, (1988, p. 332) are executed.
# Note, that n.optim = 10, because of frequent local minima.
# get the data
data(rb.fig2.data)
# positions of model and restriction files:
model2 < system.file("extdata", "rb.fig2.model", package = "MPTinR")
model2r.r.eq < system.file("extdata", "rb.fig2.r.equal", package = "MPTinR")
model2r.c.eq < system.file("extdata", "rb.fig2.c.equal", package = "MPTinR")
# The full (i.e., unconstrained) model
(ref.model < fit.mpt.old(rb.fig2.data, model2, n.optim = 10))
# All r equal
(r.equal < fit.mpt.old(rb.fig2.data, model2, model2r.r.eq, n.optim = 10))
# All c equal
(c.equal < fit.mpt.old(rb.fig2.data, model2, model2r.c.eq, n.optim = 10))
# is setting all r equal a good idea?
(g.sq.r.equal < r.equal[["goodness.of.fit"]][["G.Squared"]] 
ref.model[["goodness.of.fit"]][["G.Squared"]])
(df.r.equal < r.equal[["goodness.of.fit"]][["df"]] 
ref.model[["goodness.of.fit"]][["df"]])
(p.value.r.equal < pchisq(g.sq.r.equal, df.r.equal , lower.tail = FALSE))
# is setting all c equal a good idea?
(g.sq.c.equal < c.equal[["goodness.of.fit"]][["G.Squared"]] 
ref.model[["goodness.of.fit"]][["G.Squared"]])
(df.c.equal < c.equal[["goodness.of.fit"]][["df"]] 
ref.model[["goodness.of.fit"]][["df"]])
(p.value.c.equal < pchisq(g.sq.c.equal, df.c.equal , lower.tail = FALSE))
# Example from Broeder & Schuetz (2009)
# We fit the data from the 40 individuals from their Experiment 3
# We fit three different models:
# 1. Their 2HTM model: br.2htm
# 2. A restricted 2HTM model with Dn = Do: br.2htm.res
# 3. A 1HTM model (i.e., Dn = 0): br.1htm
# We fit the models with, as well as without, applied inequality restrictions (see Details)
# that is, for some models (.ineq) we impose: G1 < G2 < G3 < G4 < G5
# As will be apparent, the inequality restrictions do not hold for all individuals.
# Finally, we compute the FIA for all models, taking inequalities into account.
# Note: The following examples will take some time (> 1 hour).
data(d.broeder)
m.2htm < system.file("extdata", "5points.2htm.model", package = "MPTinR")
r.2htm < system.file("extdata", "broeder.2htm.restr", package = "MPTinR")
r.1htm < system.file("extdata", "broeder.1htm.restr", package = "MPTinR")
i.2htm < system.file("extdata", "broeder.2htm.ineq", package = "MPTinR")
ir.2htm < system.file("extdata", "broeder.2htm.restr.ineq", package = "MPTinR")
ir.1htm < system.file("extdata", "broeder.1htm.restr.ineq", package = "MPTinR")
# fit the original 2HTM
br.2htm < fit.mpt.old(d.broeder, m.2htm)
br.2htm.ineq < fit.mpt.old(d.broeder, m.2htm, i.2htm)
# do the inequalities hold for all participants?
br.2htm.ineq[["parameters"]][["individual"]][,"estimates",]
br.2htm[["parameters"]][["individual"]][,"estimates",]
# See the difference between forced and nonforced inequality restrictions:
round(br.2htm[["parameters"]][["individual"]][,"estimates",] 
br.2htm.ineq[["parameters"]][["individual"]][,"estimates",],2)
# The same for the other two models
# The restricted 2HTM
br.2htm.res < fit.mpt(d.broeder, m.2htm, r.2htm)
br.2htm.res.ineq < fit.mpt(d.broeder, m.2htm, ir.2htm)
round(br.2htm.res[["parameters"]][["individual"]][,"estimates",] 
br.2htm.res.ineq[["parameters"]][["individual"]][,"estimates",],2)
# The 1HTM
br.1htm < fit.mpt(d.broeder, m.2htm, r.1htm)
br.1htm.ineq < fit.mpt(d.broeder, m.2htm, ir.1htm)
round(br.2htm.res[["parameters"]][["individual"]][,"estimates",] 
br.2htm.res.ineq[["parameters"]][["individual"]][,"estimates",],2)
# These results show that we cannot compute inequality constraints for the non inequality
# imposed models (It would look differently if we excluded critical cases,
# i.e., 2, 6, 7, 10, 18, 21, 25, 29, 32, 34, 35, 37, 38)
# Therefore, we get the FIA for the models as computed above
# WARNING: The following part will take a long time!
br.2htm.fia < fit.mpt.old(d.broeder, m.2htm, fia = 200000)
br.2htm.ineq.fia < fit.mpt.old(d.broeder, m.2htm, i.2htm, fia = 200000)
br.2htm.res.fia < fit.mpt.old(d.broeder, m.2htm, r.2htm, fia = 200000 )
br.2htm.res.ineq.fia < fit.mpt.old(d.broeder, m.2htm, ir.2htm, fia = 200000)
br.1htm.fia < fit.mpt.old(d.broeder, m.2htm, r.1htm, fia = 200000)
br.1htm.ineq.fia < fit.mpt.old(d.broeder, m.2htm, ir.1htm, fia = 200000)
# Model selection using the FIA
(br.select < select.mpt(list(orig.2htm = br.2htm.fia, orig.2htm.ineq = br.2htm.ineq.fia,
res.2htm = br.2htm.res.fia, res.2htm.ineq = br.2htm.res.ineq.fia,
orig.1htm = br.1htm.fia, orig.1htm.ineq = br.1htm.ineq.fia)))
# The same results, ordered by FIA
br.select[order(br.select[,"delta.FIA.sum"]),]
# Compare this with the model selection not using FIA:
select.mpt(list(orig.2htm = br.2htm, orig.2htm.ineq = br.2htm.ineq,
res.2htm = br.2htm.res, res.2htm.ineq = br.2htm.res.ineq,
orig.1htm = br.1htm, orig.1htm.ineq = br.1htm.ineq))
# compare speed of no multicore versus multicore for multiple optimization runs:
require(snowfall)
# change number of CPUs if more are available
nCPU = 2
sfInit( parallel=TRUE, cpus=nCPU, type = "SOCK" )
# NO multicore
system.time(fit.mpt.old(d.broeder, m.2htm))
# multicore:
system.time(fit.mpt.old(d.broeder, m.2htm, multicore = "n.optim"))
sfStop()
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.