Skip to contents

Quantitative Model Comparison

Source: R/QuantitativeAnalysis.R
QuantModelComparison.Rd

subject_modelweights computes the model weights (as probabilities) of individual subjects based on an information criterion (BIC, AIC, or AICc). group_BMS performs a Bayesian model comparison based on marginal likelihoods (alias model evidence), given for different models across different subject on a group level using a fixed effects model and a random effects model on the distribution of model probabilities (see REFERENCE HERE) group_BMS_fits is a wrapper for group_BMS that can be used with the output of fitRTConfModels, i.e. a data frame with information criteria for different models and subjects, using a information criterion to approximate the model evidence.

Usage

subject_modelweights(fits, measure = "BIC")

group_BMS_fits(fits, measure = "BIC", opts = list(), alpha0 = NULL)

group_BMS(mlp, opts = list(), alpha0 = NULL)

Arguments

fits

a data frame as returned by fitRTConfModels. Should contain a column modelindicating the model name, a column subject (alternatively sbj or participant) indicating the grouping structure of the data, and a column with the name given by the measure argument containing the values of the information criterion that should be used to approximate model evidence.

measure

the name of the column indicating the information criterion to approximate model evidence. For outputs of fitRTConfModels, the available measures are 'BIC', 'AIC', and 'AICc'. Any other approximation for the model evidence may be used, the measure is transferred to log model evidence by taking -measure/2.

opts

a list with options for the iteration algorithm to estimate the parameter of the Dirichlet distribution. Following values may be provided:

  • maxiter the maximum number of iterations (Default: 200)

  • tol the tolerance for changes in the free energy approximation to stop the algorithm, if abs(FE(i+1)-FE(i))<tol the algorithm is stopped (Default: 1e-4)

  • eps The number to substitute values of 0 in calls to log (Default: 1e-32)

alpha0

a positive numeric vector representing the parameter of a Dirichlet distribution used as prior over model probabilities. The length should be equal to nrow(mlp) for group_BMS, and equal to the number of unique names in the model column of fits for group_BMS_fits.

mlp

a matrix containing the logarithm of marginal probabilities (i.e. log model evidence) with N columns representing individuals (or any other grouping structure) and K rows representing the models.

Value

subject_modelweights returns a data frame of subject-wise model probabilities with rows for each subject and columns for the models given by name and one column for the subject ID as given in the input. group_BMS and group_BMS_fits return a list with two entries:

  • model_weights: a matrix with rows for each model (row names indicate the model names for group_BMS_fits and for group_BMS if row names are available in mlp), and following columns: alpha (the alpha parameter of the Dirichlet posterior over model probabilities in the population), r (the mean probabilities of each model in the population), ep and pep (exceedance and protected exceedance probabilities for each model), and fx_prop (the posterior model probabilities if a fixed true model is assumed in the population).

  • summary_stats: a vector giving statistics for the Bayesian model comparison that may be used for other analyses: Bayesian omnibus risks: bor (random effects model against the null model), bor_fixed (fixed effects model against the null model), and bor_re_fixed (random effects model against the fixed effects model), and estimations of the Free Energy of the Dirichlet distribution FE (random effects model), FE0 (null model), and FEfixed (fixed effects model)

Details

This set of function can be used for model comparisons on a group level when the models were not fitted hierarchical but by fitting the models independently to different subgroups (e.g. data from different subjects).

The function subject_modelweights computes the model weights for each subject separately to inspect predominant models but also heterogeneity within the population. The functions group_BMS and group_BMS_fits can be used for a Bayesian model selection on the group level following the approach of Rigoux et al. (2014). It therefore aggregates the available subject-level model evidences according to a random effects model, which assumes that model probabilities are distributed in the population according to a Dirichlet distribution. Protected exceedance probability represents the probability that a specific model is more probable than any other model for a random subject out of the population scaled by the Bayesian omnibus risk of assuming a random effects model over a null effect model, in which all models have the same probability.

References

Rigoux, L., Stephan, K. E., Friston, K. J., & Daunizeau, J. (2014). Bayesian model selection for group studies - revisited. NeuroImage, 84, 971–985. doi: 10.1016/j.neuroimage.2013.08.065

Author

Sebastian Hellmann.

Examples

# Define a data frame with information criteria from model fits
# (this is a sub-data.frame from an output of fitRTConfModels with
# 8 subjects, three models and rounded information criteria)
fits <- data.frame(
  participant = rep(1:8, each=3),
  model = rep(c("dynaViTE", "2DSD", "PCRMt"), 8),
  BIC = c(5318, 5665, 1659, 3856, 5508, 3982, 3950, 3998,
          4114, 4216, 4314, 4419, 3170, 3489, 3256, 1950,
          1934, 2051, 3194, 3317, 3359, 9656, 10161, 4024),
  AIC = c(5211, 5577, 1577, 3750, 5420, 3899, 3843, 3911,
          4031, 4109, 4226, 4337, 3063, 3401, 3173, 1844,
          1847, 1969, 3087, 3229, 3277, 9549, 10074, 3942),
  AICc = c(5212, 5578, 1577, 3751, 5421, 3900, 3844, 3911,
           4032, 4110, 4227, 4337, 3064, 3402, 3174, 1845,
           1848, 1970, 3088, 3230, 3277, 9550, 10074, 3942))
# Compute subject-wise model probabitilities based on different ICs
subject_modelweights(fits, measure = "BIC")
#>           2DSD        PCRMt     dynaViTE participant
#> 1 0.000000e+00 1.000000e+00 0.0000000000           1
#> 2 0.000000e+00 4.359610e-28 1.0000000000           2
#> 3 3.775135e-11 2.442601e-36 1.0000000000           3
#> 4 5.242886e-22 8.300611e-45 1.0000000000           4
#> 5 5.370691e-70 2.115131e-19 1.0000000000           5
#> 6 9.996646e-01 3.923080e-26 0.0003353501           6
#> 7 1.953842e-27 1.481512e-36 1.0000000000           7
#> 8 0.000000e+00 1.000000e+00 0.0000000000           8
subject_modelweights(fits, measure = "AIC")
#>           2DSD        PCRMt  dynaViTE participant
#> 1 0.000000e+00 1.000000e+00 0.0000000           1
#> 2 0.000000e+00 4.416326e-33 1.0000000           2
#> 3 1.713908e-15 1.500786e-41 1.0000000           3
#> 4 3.924396e-26 3.093350e-50 1.0000000           4
#> 5 4.020060e-74 1.299581e-24 1.0000000           5
#> 6 1.824255e-01 5.876547e-28 0.8175745           6
#> 7 1.462486e-31 5.521082e-42 1.0000000           7
#> 8 0.000000e+00 1.000000e+00 0.0000000           8
subject_modelweights(fits, measure = "AICc")
#>           2DSD        PCRMt  dynaViTE participant
#> 1 0.000000e+00 1.000000e+00 0.0000000           1
#> 2 0.000000e+00 4.416326e-33 1.0000000           2
#> 3 2.825757e-15 1.500786e-41 1.0000000           3
#> 4 3.924396e-26 5.100072e-50 1.0000000           4
#> 5 4.020060e-74 1.299581e-24 1.0000000           5
#> 6 1.824255e-01 5.876547e-28 0.8175745           6
#> 7 1.462486e-31 9.102726e-42 1.0000000           7
#> 8 0.000000e+00 1.000000e+00 0.0000000           8
# Conduct group-level Bayesian model selection based on BIC
group_BMS_fits(fits, measure="BIC")
#> $model_weights
#>             alpha         r     ep    pep fx_prop
#> 2DSD     1.998461 0.1816782 0.0267 0.0267       0
#> PCRMt    1.500000 0.1363636 0.0135 0.0135       1
#> dynaViTE 7.501539 0.6819581 0.9598 0.9598       0
#> 
#> $summary_stats
#>          bor           FE          FE0      FEfixed    bor_fixed bor_re_fixed 
#>        0.000     1715.178   -13010.289   -13433.099        1.000        0.000 
#> 


## General group-level Bayesian model selection based on any marginal log-probabilities
# Compute marginal log-likelihood based on BIC from fits
mlp <- matrix(NA, ncol=8, nrow=3)
for (i in 1:8) mlp[,i] <- fits[(i-1)*3 + 1:3, "BIC"]
mlp <- - mlp/(2)
rownames(mlp) <- c("dynaViTE", "2DSD", "PCRMt")
# conduct group BMS:
group_BMS(mlp)
#> $model_weights
#>             alpha         r     ep    pep fx_prop
#> dynaViTE 7.501539 0.6819581 0.9634 0.9634       0
#> 2DSD     1.998461 0.1816782 0.0252 0.0252       0
#> PCRMt    1.500000 0.1363636 0.0114 0.0114       1
#> 
#> $summary_stats
#>          bor           FE          FE0      FEfixed    bor_fixed bor_re_fixed 
#>        0.000     1715.178   -13010.289   -13433.099        1.000        0.000 
#>