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 L-BFGS-B 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 )
Either a numeric
Number of optimization runs. Can be parallelized via
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
A scalar corresponding to the size of the confidence intervals for the parameter estimates. Default is 95 which corresponds to 95% confidence intervals.
Logical. Indicates whether or not inequality restrictions (when present in the model file) should be enforced while fitting the model. If
Logical. If TRUE, parameters are alphabetically sorted in the parameter table. If FALSE, the first parameters in the parameter table are the non-restricted ones, followed by the restricted parameters. Default is TRUE.
Character vector specifying whether the model file is formatted in the easy way (
Character vector. If not
Logical. Relevant if
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 (see http://www.psychologie.uni-freiburg.de/Members/singmann/R/mptinr) 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 multi-individual fit to the respective column).
The difference between both types of EQN format (
"eqn2") is the way the first line of the model file is treated. If
model.file is set to
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 set to
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 multi-individual 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
FALSE. This can be useful for obtaining the FIA (see examples in Wu, Myung, & Batchelder, 2010).
The fitting/optimization is achieved via
optim's L-BFGS-B 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 Log-Likelihood 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 model-selection, 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 context-free language of MPT models (L-BMPT; 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.,
(1 - Dn)) and the operators + and *, but nothing else. Simplifications of the equations will seriously distort this step. This step is achieved by
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
(Actually, both steps are achieved by a call to
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 (
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 non-restricted 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
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.uni-freiburg.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 self-written) visible in the global environment (i.e., not only +, *, and - as operators). Currently
fit.model should be viewed as experimental.
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.
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 non-MPT models.
For individual fits (i.e.,
data is a
list containing one or more of the following components from the best fitting model:
A data.frame containing the parameter estimates and corresponding confidence intervals. If a restriction file was present, the restricted parameters are marked.
For multi-individual fits (i.e.,
data is a
list with similar elements, but the following differences.
The first elements,
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
parameters is a list containing:
A 3-dimensional array containing the parameter estimates ([,1,]), confidence intervals [,2:3,], and, if restrictions not
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.
data contains two matrices, one with the
observed, and one with the
n.optim > 1, the
summary of the vector (matrix for multi-individual fit) containing the Log-Likelihood values returned by each run of
optim is added to the output.
When using R (>= 2.13.0) compiling
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, 467-470.
Broeder, A., & Schuetz, J. (2009). Recognition ROCs are curvilinear-or are they? On premature arguments against the two-high-threshold 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, 1190-1208.
Knapp, B. R., & Batchelder, W. H. (2004). Representing parametric order constraints in multi-trial applications of multinomial processing tree models. Journal of Mathematical Psychology, 48, 215-229.
Moshagen, M. (2010). multiTree: A computer program for the analysis of multinomial processing tree models. Behavior Research Methods, 42, 42-54.
Purdy, B. P., & Batchelder, W. H. (2009). A context-free language for binary multinomial processing tree models. Journal of Mathematical Psychology, 53, 547-561.
Riefer, D. M., & Batchelder, W. H. (1988). Multinomial modeling and the measurement of cognitive processes. Psychological Review, 95, 318-339.
Stahl, C. & Klauer, K. C. (2007). HMMTree: A computer program for latent-class 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, 275-286.
Wu, H., Myung, J.I., & Batchelder, W.H. (2010b). On the minimum description length complexity of multinomial processing trees. Journal of Mathematical Psychology, 54, 291-303.
fit.mpt for the current and recommended function fot fitting MPTs
http://www.psychologie.uni-freiburg.de/Members/singmann/R/mptinr for additional information on model files and restriction files
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 non-forced 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.