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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
#' Fit nonlinear mixed-effects models built from one or more kinetic
#' degradation models and one or more error models
#'
#' The name of the methods expresses that (**m**ultiple) **h**ierarchichal
#' (also known as multilevel) **m**ulticompartment **kin**etic models are
#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear
#' mixed-effects model fitting functions.
#'
#' @param objects A list of [mmkin] objects containing fits of the same
#' degradation models to the same data, but using different error models.
#' @param backend The backend to be used for fitting. Currently, only saemix is
#' supported
#' @param algorithm The algorithm to be used for fitting (currently not used)
#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
#' model fitting function.
#' @param cores The number of cores to be used for multicore processing. This
#' is only used when the \code{cluster} argument is \code{NULL}. On Windows
#' machines, cores > 1 is not supported, you need to use the \code{cluster}
#' argument to use multiple logical processors. Per default, all cores detected
#' by [parallel::detectCores()] are used, except on Windows where the default
#' is 1.
#' @param cluster A cluster as returned by [makeCluster] to be used for
#' parallel execution.
#' @importFrom parallel mclapply parLapply detectCores
#' @return A two-dimensional [array] of fit objects and/or try-errors that can
#' be indexed using the degradation model names for the first index (row index)
#' and the error model names for the second index (column index), with class
#' attribute 'mhmkin'.
#' @author Johannes Ranke
#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
#' @export
mhmkin <- function(objects, backend = "saemix", algorithm = "saem", ...) {
UseMethod("mhmkin")
}
#' @export
#' @rdname mhmkin
mhmkin.list <- function(objects, backend = "saemix",
...,
cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
{
call <- match.call()
dot_args <- list(...)
backend_function <- switch(backend,
saemix = "saem"
)
deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod)
names(deg_models) <- dimnames(objects[[1]])$model
n.deg <- length(deg_models)
ds <- lapply(objects[[1]][1, ], function(x) x$data)
for (other in objects[-1]) {
# Check if the degradation models in all objects are the same
for (deg_model_name in names(deg_models)) {
if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec,
deg_models[[deg_model_name]]$spec))
{
stop("The mmkin objects have to be based on the same degradation models")
}
}
# Check if they have been fitted to the same dataset
other_object_ds <- lapply(other[1, ], function(x) x$data)
for (i in 1:length(ds)) {
if (!all.equal(ds[[i]][c("time", "variable", "observed")],
other_object_ds[[i]][c("time", "variable", "observed")]))
{
stop("The mmkin objects have to be fitted to the same datasets")
}
}
}
n.o <- length(objects)
error_models = sapply(objects, function(x) x[[1]]$err_mod)
n.e <- length(error_models)
n.fits <- n.deg * n.e
fit_indices <- matrix(1:n.fits, ncol = n.e)
dimnames(fit_indices) <- list(degradation = names(deg_models),
error = error_models)
fit_function <- function(fit_index) {
w <- which(fit_indices == fit_index, arr.ind = TRUE)
deg_index <- w[1]
error_index <- w[2]
mmkin_row <- objects[[error_index]][deg_index, ]
res <- try(do.call(backend_function, args = c(list(mmkin_row), dot_args)))
return(res)
}
fit_time <- system.time({
if (is.null(cluster)) {
results <- parallel::mclapply(as.list(1:n.fits), fit_function,
mc.cores = cores, mc.preschedule = FALSE)
} else {
results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
}
})
attributes(results) <- attributes(fit_indices)
attr(results, "call") <- call
attr(results, "time") <- fit_time
class(results) <- "mhmkin"
return(results)
}
#' Subsetting method for mhmkin objects
#'
#' @param x An [mhmkin] object.
#' @param i Row index selecting the fits for specific models
#' @param j Column index selecting the fits to specific datasets
#' @param drop If FALSE, the method always returns an mhmkin object, otherwise
#' either a list of fit objects or a single fit object.
#' @return An object of class \code{\link{mhmkin}}.
#' @rdname mhmkin
#' @export
`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
class(x) <- NULL
x_sub <- x[i, j, drop = drop]
if (!drop) {
class(x_sub) <- "mhmkin"
}
return(x_sub)
}
#' Print method for mhmkin objects
#'
#' @rdname mhmkin
#' @export
print.mhmkin <- function(x, ...) {
cat("<mhmkin> object\n")
cat("Status of individual fits:\n\n")
print(convergence(x))
}
#' @export
convergence.mhmkin <- function(object, ...) {
all_summary_warnings <- character()
if (inherits(object[[1]], "saem.mmkin")) {
test_func <- function(fit) {
if (inherits(fit$so, "try-error")) return("E")
else return("OK")
}
} else {
stop("Only mhmkin objects containing saem.mmkin objects currently supported")
}
result <- lapply(object, test_func)
result <- unlist(result)
dim(result) <- dim(object)
dimnames(result) <- dimnames(object)
class(result) <- "convergence.mhmkin"
return(result)
}
#' @export
print.convergence.mhmkin <- function(x, ...) {
class(x) <- NULL
print(x, quote = FALSE)
cat("\n")
if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
if (any(x == "E")) cat("E: Error\n")
}
#' @export
AIC.mhmkin <- function(object, ..., k = 2) {
res <- sapply(object, function(x) {
if (inherits(x, "try-error")) return(NA)
else return(AIC(x$so, k = k))
})
dim(res) <- dim(object)
dimnames(res) <- dimnames(object)
return(res)
}
#' @export
BIC.mhmkin <- function(object, ...) {
res <- sapply(object, function(x) {
if (inherits(x, "try-error")) return(NA)
else return(BIC(x$so))
})
dim(res) <- dim(object)
dimnames(res) <- dimnames(object)
return(res)
}
|
