DstarM (version 0.2.2)

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:

AIC

Akaike Information Criterion.

AICc

AIC with a correction for finite sample sizes.

BIC

Bayesian Information Criterion, also known as Schwarz criterion.

npar

Number of parameters.

logLik

The natural log of the likelihood.

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

# NOT RUN {
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, data = dat)
resObs2 = estObserved(resD2, resND2, data = dat)
resObs3 = estObserved(resD3, resND3, data = dat)
# Compare optimizer fitness
lstObs = list(resObs1, resObs2, resObs3)
sapply(lstObs, function(x) x$fit)          # model 2 performs best!
sapply(lstObs, function(x) x$fit$chisq$sum)# model 2 performs best!
# model 3 has lower Decision fit due to overfitting the decision model
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, data = dat)
resObsC2 = estObserved(resC2, data = dat)
resObsC3 = estObserved(resC3, data = dat)
# Compare optimizer fitness
lstObsC = list(resObsC1, resObsC2, resObsC3)
sapply(lstObsC, function(x) x$fit)
sapply(lstObsC, function(x) x$fit$chisq$sum)
# 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))
# }