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
|
#' @importFrom lmtest lrtest
#' @export
lmtest::lrtest
#' Likelihood ratio test for mkinfit models
#'
#' Compare two mkinfit models based on their likelihood. If two fitted
#' mkinfit objects are given as arguments, it is checked if they have been
#' fitted to the same data. It is the responsibility of the user to make sure
#' that the models are nested, i.e. one of them has less degrees of freedom
#' and can be expressed by fixing the parameters of the other.
#'
#' Alternatively, an argument to mkinfit can be given which is then passed
#' to \code{\link{update.mkinfit}} to obtain the alternative model.
#'
#' The comparison is then made by the \code{\link[lmtest]{lrtest.default}}
#' method from the lmtest package. The model with the higher number of fitted
#' parameters (alternative hypothesis) is listed first, then the model with the
#' lower number of fitted parameters (null hypothesis).
#'
#' The method for mmkin objects only works for column objects with two members.
#'
#' @importFrom stats logLik update
#' @param object An \code{\link{mkinfit}} object
#' @param object_2 Optionally, another mkinfit object fitted to the same data.
#' @param \dots Argument to \code{\link{mkinfit}}, passed to
#' \code{\link{update.mkinfit}} for creating the alternative fitted object.
#' @examples
#' \dontrun{
#' test_data <- subset(synthetic_data_for_UBA_2014[[12]]$data, name == "parent")
#' sfo_fit <- mkinfit("SFO", test_data, quiet = TRUE)
#' dfop_fit <- mkinfit("DFOP", test_data, quiet = TRUE)
#' lrtest(dfop_fit, sfo_fit)
#' lrtest(sfo_fit, dfop_fit)
#'
#' # The following two examples are commented out as they fail during
#' # generation of the static help pages by pkgdown
#' #lrtest(dfop_fit, error_model = "tc")
#' #lrtest(dfop_fit, fixed_parms = c(k2 = 0))
#'
#' # However, this equivalent syntax works for static help pages
#' lrtest(dfop_fit, update(dfop_fit, error_model = "tc"))
#' lrtest(dfop_fit, update(dfop_fit, fixed_parms = c(k2 = 0)))
#' }
#' @export
lrtest.mkinfit <- function(object, object_2 = NULL, ...) {
name_function <- function(x) {
object_name <- paste(x$mkinmod$name, "with error model", x$err_mod)
if (length(x$bparms.fixed) > 0) {
object_name <- paste(object_name,
"and fixed parameter(s)",
paste(names(x$bparms.fixed), collapse = ", "))
}
return(object_name)
}
if (is.null(object_2)) {
object_2 <- update(object, ...)
} else {
data_object <- object$data[c("time", "variable", "observed")]
data_object_2 <- object_2$data[c("time", "variable", "observed")]
if (!identical(data_object, data_object_2)) {
stop("It seems that the mkinfit objects have not been fitted to the same data")
}
}
if (attr(logLik(object), "df") > attr(logLik(object_2), "df")) {
lmtest::lrtest.default(object, object_2, name = name_function)
} else {
lmtest::lrtest.default(object_2, object, name = name_function)
}
}
#' @rdname lrtest.mkinfit
#' @export
lrtest.mmkin <- function(object, ...) {
if (nrow(object) != 2 | ncol(object) > 1) stop("Only works for a column containing two mkinfit objects")
lrtest(object[[1, 1]], object[[2, 1]])
}
|
