DstarM (version 0.1.0)

calcIC: Calculate Information Criteria

Description

Calculate information criteria for D*M models. This function can be called with the resObserved argument or the resDecision and, in case of D*M analyses, also the nondecision analysis. If resObserved is not supplied a call to estObserved will be made.

Usage

calcIC(resObserved, resDecision, resND, data, npar)

Arguments

resObserved

an object of class D*M, specifically output of the function estObserved.

resDecision

an object of class D*M, specifically output of the function estDstarM.

resND

an object of class D*M, specifically output of the function estND.

data

The data used to estimate resDecision.

npar

The number of parameters estimated.

Value

a list (S3 object of class 'D*M') that contains:

Details

Calculates several information criteria.

Calculate information criteria for D*M models. The results should be interpreted carefully! D*M models are not estimated by maximizing a likelihood, and therefore more complex models do not always necessarily have a higher likelihood than less complex models. Also, adding parameters to the decision model might improve the fit of the decision model, but can worsen the fit of the nondecision model. As a result the total fit decreases when adding unneeded parameters. This is shown in the example. Caution is adviced when interpreting likelihood based information criteria.

Examples

rm(list = ls())
set.seed(42)
# Simulate data, fit 3 models, compare models.
# simulate data with three stimuli of different difficulty.
# this implies different drift rates across conditions.
# define a time grid. A more reasonable stepsize is .01; this is just for speed.
tt = seq(0, 5, .1)
pars = c(.8, 2, .5, .5, .5, # condition 1
        .8, 3, .5, .5, .5, # condition 2
        .8, 4, .5, .5, .5) # condition 3
pdfND = dbeta(tt, 10, 30)
# simulate data
dat = simData(n = 3e3, pars = pars, tt = tt, pdfND = pdfND)
# define restriction matrices
restr1 = restr2 = restr3 = matrix(1:5, 5, 3)
## restr1 allows nothing to differ over conditions
restr2[2, 2:3] = 6:7 # allow drift rates to differ
restr3[1:3, 2:3] = 6:11 # allow all decision model parameters to differ
# fix variance parameters
fixed = matrix(c('sz1', .5, 'sv1', .5), 2, 2)
## Not run: 
# # Run D*M analysis
# resD1 = estDstarM(dat = dat, tt = tt, restr = restr1, fixed = fixed, Optim = list(parallelType = 1))
# resD2 = estDstarM(dat = dat, tt = tt, restr = restr2, fixed = fixed, Optim = list(parallelType = 1))
# resD3 = estDstarM(dat = dat, tt = tt, restr = restr3, fixed = fixed, Optim = list(parallelType = 1))
# # Estimate nondecision density - The warnings can be ignored.
# # Quantiles for nondecision densities cannot be calculated due to the big time grid (by .1).
# resND1 = estND(resD1, Optim = list(parallelType = 1))
# resND2 = estND(resD2, Optim = list(parallelType = 1))
# resND3 = estND(resD3, Optim = list(parallelType = 1))
# # Estimate observed densities
# resObs1 = estObserved(resD1, resND1)
# resObs2 = estObserved(resD2, resND2)
# resObs3 = estObserved(resD3, resND3)
# # Compare optimizer fitness
# lstObs = list(resObs1, resObs2, resObs3)
# # fit$total is the sum of fit$Decision and
# # fit$ND. Model 3 has the lowest fit
# sapply(lstObs, function(x) x$fit$total)
# # model 3 has lowest Decision fit due to overfitting
# sapply(lstObs, function(x) x$fit$Decision)
# # However, model 3 has worse nondecision fit compared to model 2 due to overfitting!
# sapply(lstObs, function(x) x$fit$ND)
# # Calculate information criteria
# IC1 = calcIC(resObs1, resD1, data = dat)
# IC2 = calcIC(resObs2, resD2, data = dat)
# IC3 = calcIC(resObs3, resD3, data = dat)
# sapply(list(IC1, IC2, IC3), function(x) x$logLik)
# # Do likelihood ratio tests
# anova(IC1, IC2, IC3) # unorthodox output, model 2 is best?
# which.min(c(IC1$AIC, IC2$AIC, IC3$AIC)) # model 2 is best
# which.min(c(IC1$BIC, IC2$BIC, IC3$BIC)) # model 2 is best
# # We can do the same for traditional analyses
# # likelihood based methods do work for traditional analyses
# # define restriction matrices
# restrC1 = restrC2 = restrC3 = matrix(1:7, 7, 3)
# # restr1 allows nothing to differ over conditions
# restrC2[2, 2:3] = 8:9 # allow drift rates to differ
# restrC3[c(1:2, 4), 2:3] = 8:13 # allow all decision model parameters to differ
# # Do traditional analyses
# resC1 = estDstarM(dat = dat, tt = tt, restr = restrC1, fixed = fixed,
# DstarM = FALSE, Optim = list(parallelType = 1))
# resC2 = estDstarM(dat = dat, tt = tt, restr = restrC2, fixed = fixed,
# DstarM = FALSE, Optim = list(parallelType = 1))
# resC3 = estDstarM(dat = dat, tt = tt, restr = restrC3, fixed = fixed,
# DstarM = FALSE, Optim = list(parallelType = 1))
# # resObs does not need to be calculated now since resC* contains the full model
# # Estimate observed densities
# resObsC1 = estObserved(resC1)
# resObsC2 = estObserved(resC2)
# resObsC3 = estObserved(resC3)
# # Compare optimizer fitness
# lstObsC = list(resObsC1, resObsC2, resObsC3)
# sapply(lstObsC, function(x) x$fit$total)
# # Calculate information criteria
# ICC1 = calcIC(resObserved = resObsC1, data = dat)
# ICC2 = calcIC(resDecision = resC2, data = dat) # both input methods are possible
# ICC3 = calcIC(resObserved = resObsC3, data = dat)
# sapply(list(ICC1, ICC2, ICC3), function(x) x$logLik)
# # Do likelihood ratio tests
# # correct model is retrieved by all methods
# anova(ICC1, ICC2, ICC3)
# which.min(c(ICC1$AIC, ICC2$AIC, ICC3$AIC))
# which.min(c(ICC1$BIC, ICC2$BIC, ICC3$BIC))
# ## End(Not run)