/
blrt.R
109 lines (105 loc) · 3.56 KB
/
blrt.R
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
#' @title Conduct Bootstrapped Likelihood Ratio Test
#' @description Conduct Bootstrapped Likelihood Ratio Test to compare two
#' mixture models.
#' @param x An object for which a method exists.
#' @param replications Integer reflecting the number of bootstrapped
#' replications, defaults to `100`.
#' @param ... further arguments to be passed to or from other methods.
#' @return A data.frame.
#' @examples
#' \dontrun{
#' df <- iris[, 1, drop = FALSE]
#' names(df) <- "x"
#' res <- mx_mixture(model = "x ~ m{C}*1
#' x ~~ v{C}*x", classes = 1:2, data = df)
#' BLRT(res, replications = 4)
#' }
#' @export
BLRT <- function(x, replications = 100, ...){
UseMethod("BLRT", x)
}
#' @method BLRT MxModel
#' @export
BLRT.MxModel <- function(x, replications = 100, ...){
dots <- list(...)
object1 <- x
object2 <- dots[[which(sapply(dots, inherits, what = "MxModel"))[1]]]
blrt_internal(object1, object2, replications = replications)
}
# blrt_simple <- function(mod_simple, mod_complex, replications = 100, parallel = TRUE){
#' @method BLRT mixture_list
#' @export
BLRT.mixture_list <- function(x, replications = 100, ...){
df_empty <- data.frame(lr = NA, df = NA, blrt_p = NA, samples = NA)
if(length(x) > 1){
out <- mapply(function(smaller, bigger){
tryCatch({
blrt_internal(smaller, bigger, replications = replications)
},
error = function(e){
df_empty })
}, bigger = x[-1], smaller = x[-length(x)], SIMPLIFY = FALSE)
out <- do.call(rbind, append(out, list(df_empty), 0))
} else {
out <- df_empty
}
out <- data.frame(null = c(NA, sapply(x[-length(x)], function(x){x@name})),
alt = c(NA, sapply(x[-1], function(x){x@name})),
out)[-1, , drop = FALSE]
class(out) <- c("LRT", class(out))
attr(out, "type") <- "Bootstrapped"
return(out)
}
#' @method BLRT list
#' @export
BLRT.list <- BLRT.mixture_list
#' @importFrom future.apply future_lapply
#' @importFrom progressr with_progress progressor
blrt_internal <-
function(mod_simple, mod_complex, replications = 100) {
lrtest <-
mod_simple@output$Minus2LogLikelihood - mod_complex@output$Minus2LogLikelihood
progmsg <- paste0(mod_simple@name, " vs. ", mod_complex@name)
progressr::with_progress({
pgs <- progressr::progressor(steps = replications)
bootres <- future.apply::future_lapply(
X = 1:replications,
future.seed = TRUE,
FUN = function(i) {
pgs(sprintf(progmsg))
tryCatch({
df_sim <- mxGenerateData(mod_simple)
mod_simple@data$observed <-
df_sim
mod_complex@data$observed <-
df_sim
mod_simple <-
mxRun(mod_simple,
silent = TRUE,
suppressWarnings = TRUE)
mod_complex <-
mxRun(mod_complex,
silent = TRUE,
suppressWarnings = TRUE)
c(
mod_simple@output$Minus2LogLikelihood - mod_complex@output$Minus2LogLikelihood,
mod_simple@output$status$code + mod_complex@output$status$code
)
}, error = function(e) {
c(NA, 1)
})
}
)})
bootres <- do.call(rbind, bootres)
isvalid <- bootres[, 2] == 0
lrdist <- bootres[isvalid, 1]
out <- data.frame(
lr = lrtest,
df = length(omxGetParameters(mod_complex)) - length(omxGetParameters(mod_simple)),
blrt_p = sum(lrdist > lrtest) / length(lrdist),
samples = sum(isvalid)
)
class(out) <- c("LRT", class(out))
attr(out, "type") <- "Bootstrapped"
return(out)
}