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
|
##' This function computes Bayes factors corresponding to restrictions on a full model.
##'
##' See the help for \code{\link{anovaBF}} and \code{\link{anovaBF}} or details.
##'
##' Models, priors, and methods of computation are provided in Rouder et al.
##' (2012) and Liang et al (2008).
##'
##' @title Function to compute Bayes factors for general designs
##' @param formula a formula containing the full model for the analysis
##' (see Examples)
##' @param data a data frame containing data for all factors in the formula
##' @param whichRandom a character vector specifying which factors are random
##' @param whichModels which set of models to compare; see Details
##' @param neverExclude a character vector containing a regular expression (see
##' help for \link{regex} for details) that indicates which terms to always keep
##' in the analysis
##' @param iterations How many Monte Carlo simulations to generate, if relevant
##' @param progress if \code{TRUE}, show progress with a text progress bar
##' @param rscaleFixed prior scale for standardized, reduced fixed effects. A
##' number of preset values can be given as strings; see Details.
##' @param rscaleRandom prior scale for standardized random effects
##' @param rscaleCont prior scale for standardized slopes
##' @param rscaleEffects A named vector of prior settings for individual factors,
##' overriding rscaleFixed and rscaleRandom. Values are scales, names are factor names.
##' @param multicore if \code{TRUE} use multiple cores through the \code{doMC}
##' package. Unavailable on Windows.
##' @param method approximation method, if needed. See \code{\link{nWayAOV}} for
##' details.
##' @param noSample if \code{TRUE}, do not sample, instead returning NA.
##' @return An object of class \code{BFBayesFactor}, containing the computed
##' model comparisons
##' @param callback callback function for third-party interfaces
##' @author Richard D. Morey (\email{richarddmorey@@gmail.com})
##' @export
##' @references
##' Rouder, J. N., Morey, R. D., Speckman, P. L., Province, J. M., (2012)
##' Default Bayes Factors for ANOVA Designs. Journal of Mathematical
##' Psychology. 56. p. 356-374.
##'
##' Liang, F. and Paulo, R. and Molina, G. and Clyde, M. A. and
##' Berger, J. O. (2008). Mixtures of g-priors for Bayesian Variable
##' Selection. Journal of the American Statistical Association, 103, pp.
##' 410-423
##'
##' @note The function \code{generalTestBF} can compute Bayes factors for all
##' restrictions of a full model against the null
##' hypothesis that all effects are 0. The total number of tests
##' computed -- if all tests are requested -- will be \eqn{2^K-1}{2^K - 1}
##' for \eqn{K} factors or covariates.
##' This number increases very quickly with the number of tested predictors. An option is included to
##' prevent testing too many models: \code{options('BFMaxModels')}, which defaults to 50,000, is
##' the maximum number of models that will be analyzed at once. This can
##' be increased by increased using \code{\link{options}}.
##'
##' It is possible to reduce the number of models tested by only testing the
##' most complex model and every restriction that can be formed by removing
##' one factor or interaction using the \code{whichModels} argument. See the
##' help for \code{\link{anovaBF}} for details.
##'
##' @examples
##' ## Puzzles example: see ?puzzles and ?anovaBF
##' data(puzzles)
##' ## neverExclude argument makes sure that participant factor ID
##' ## is in all models
##' result = generalTestBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
##' neverExclude="ID", progress=FALSE)
##' result
##'
##' @keywords htest
##' @seealso \code{\link{lmBF}}, for testing specific models, and
##' \code{\link{regressionBF}} and \code{anovaBF} for other functions for
##' testing multiple models simultaneously.
generalTestBF <-
function(formula, data, whichRandom = NULL,
whichModels = "withmain", neverExclude=NULL, iterations = 10000, progress = getOption('BFprogress', interactive()),
rscaleFixed = "medium", rscaleRandom = "nuisance", rscaleCont="medium", rscaleEffects = NULL, multicore = FALSE, method="auto",
noSample=FALSE, callback=function(...) as.integer(0))
{
data <- marshallTibble(data)
checkFormula(formula, data, analysis = "lm")
# pare whichRandom down to terms that appear in the formula
whichRandom <- whichRandom[whichRandom %in% fmlaFactors(formula, data)[-1]]
dataTypes <- createDataTypes(formula, whichRandom, data, analysis = "lm")
models = enumerateGeneralModels(formula, whichModels, neverExclude,
includeBottom = whichModels!="top",
data = data)
if(length(models)>getOption('BFMaxModels', 50000)) stop("Maximum number of models exceeded (",
length(models), " > ",getOption('BFMaxModels', 50000) ,"). ",
"The maximum can be increased by changing ",
"options('BFMaxModels').")
if(multicore){
message("Note: Progress bars and callbacks are suppressed when running multicore.")
if(!requireNamespace("doMC", quietly = TRUE)){
stop("Required package (doMC) missing for multicore functionality.")
}
doMC::registerDoMC()
if(foreach::getDoParWorkers()==1){
warning("Multicore specified, but only using 1 core. Set options(cores) to something >1.")
}
bfs <- foreach::"%dopar%"(
foreach::foreach(gIndex=models, .options.multicore=mcoptions),
lmBF(gIndex,data = data, whichRandom = whichRandom,
rscaleFixed = rscaleFixed, rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont, rscaleEffects = rscaleEffects, iterations = iterations, method=method,
progress=FALSE,noSample=noSample)
)
}else{ # Single core
checkCallback(callback,as.integer(0))
bfs = NULL
myCallback <- function(prgs){
frac <- (i - 1 + prgs/1000)/length(models)
ret <- callback(frac*1000)
return(as.integer(ret))
}
if(progress){
pb = txtProgressBar(min = 0, max = length(models), style = 3)
}else{
pb = NULL
}
for(i in 1:length(models)){
oneModel <- lmBF(models[[i]],data = data, whichRandom = whichRandom,
rscaleFixed = rscaleFixed, rscaleRandom = rscaleRandom,
rscaleCont = rscaleCont, rscaleEffects = rscaleEffects, iterations = iterations,
progress = FALSE, method = method,noSample=noSample,callback=myCallback)
if(inherits(pb,"txtProgressBar")) setTxtProgressBar(pb, i)
bfs = c(bfs,oneModel)
}
if(inherits(pb,"txtProgressBar")) close(pb)
checkCallback(callback,as.integer(1000))
}
# combine all the Bayes factors into one BFBayesFactor object
bfObj = do.call("c", bfs)
if(whichModels=="top") bfObj = BFBayesFactorTop(bfObj)
return(bfObj)
}
|
