1
#' Lack-of-fit test for models fitted to data with replicates
2
#'
3
#' This is a generic function with a method currently only defined for mkinfit
4
#' objects. It fits an anova model to the data contained in the object and
5
#' compares the likelihoods using the likelihood ratio test
6
#' \code{\link[lmtest]{lrtest.default}} from the lmtest package.
7
#'
8
#' The anova model is interpreted as the simplest form of an mkinfit model,
9
#' assuming only a constant variance about the means, but not enforcing any
10
#' structure of the means, so we have one model parameter for every mean
11
#' of replicate samples.
12
#'
13
#' @param object A model object with a defined loftest method
14
#' @param \dots Not used
15
#' @export
16
loftest <- function(object, ...) {
17 2
  UseMethod("loftest")
18
}
19

20
#' @rdname loftest
21
#' @importFrom stats logLik lm dnorm coef
22
#' @seealso lrtest
23
#' @examples
24
#' \dontrun{
25
#' test_data <- subset(synthetic_data_for_UBA_2014[[12]]$data, name == "parent")
26
#' sfo_fit <- mkinfit("SFO", test_data, quiet = TRUE)
27
#' plot_res(sfo_fit) # We see a clear pattern in the residuals
28
#' loftest(sfo_fit)  # We have a clear lack of fit
29
#' #
30
#' # We try a different model (the one that was used to generate the data)
31
#' dfop_fit <- mkinfit("DFOP", test_data, quiet = TRUE)
32
#' plot_res(dfop_fit) # We don't see systematic deviations, but heteroscedastic residuals
33
#' # therefore we should consider adapting the error model, although we have
34
#' loftest(dfop_fit) # no lack of fit
35
#' #
36
#' # This is the anova model used internally for the comparison
37
#' test_data_anova <- test_data
38
#' test_data_anova$time <- as.factor(test_data_anova$time)
39
#' anova_fit <- lm(value ~ time, data = test_data_anova)
40
#' summary(anova_fit)
41
#' logLik(anova_fit) # We get the same likelihood and degrees of freedom
42
#' #
43
#' test_data_2 <- synthetic_data_for_UBA_2014[[12]]$data
44
#' m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"),
45
#'   M1 = list(type = "SFO", to = "M2"),
46
#'   M2 = list(type = "SFO"), use_of_ff = "max")
47
#' sfo_lin_fit <- mkinfit(m_synth_SFO_lin, test_data_2, quiet = TRUE)
48
#' plot_res(sfo_lin_fit) # not a good model, we try parallel formation
49
#' loftest(sfo_lin_fit)
50
#' #
51
#' m_synth_SFO_par <- mkinmod(parent = list(type = "SFO", to = c("M1", "M2")),
52
#'   M1 = list(type = "SFO"),
53
#'   M2 = list(type = "SFO"), use_of_ff = "max")
54
#' sfo_par_fit <- mkinfit(m_synth_SFO_par, test_data_2, quiet = TRUE)
55
#' plot_res(sfo_par_fit) # much better for metabolites
56
#' loftest(sfo_par_fit)
57
#' #
58
#' m_synth_DFOP_par <- mkinmod(parent = list(type = "DFOP", to = c("M1", "M2")),
59
#'   M1 = list(type = "SFO"),
60
#'   M2 = list(type = "SFO"), use_of_ff = "max")
61
#' dfop_par_fit <- mkinfit(m_synth_DFOP_par, test_data_2, quiet = TRUE)
62
#' plot_res(dfop_par_fit) # No visual lack of fit
63
#' loftest(dfop_par_fit)  # no lack of fit found by the test
64
#' #
65
#' # The anova model used for comparison in the case of transformation products
66
#' test_data_anova_2 <- dfop_par_fit$data
67
#' test_data_anova_2$variable <- as.factor(test_data_anova_2$variable)
68
#' test_data_anova_2$time <- as.factor(test_data_anova_2$time)
69
#' anova_fit_2 <- lm(observed ~ time:variable - 1, data = test_data_anova_2)
70
#' summary(anova_fit_2)
71
#' }
72
#' @export
73
loftest.mkinfit <- function(object, ...) {
74

75 2
  name_function <- function(x) {
76 2
    object_name <- paste(x$mkinmod$name, "with error model", x$err_mod)
77 2
    if (length(x$bparms.fixed) > 0) {
78 0
      object_name <- paste(object_name,
79 0
        "and fixed parameter(s)",
80 0
        paste(names(x$bparms.fixed), collapse = ", "))
81
    }
82 2
    return(object_name)
83
  }
84

85
  # Check if we have replicates in the data
86 2
  if (max(aggregate(object$data$observed,
87 2
    by = list(object$data$variable, object$data$time), length)$x) == 1) {
88 2
    stop("Not defined for fits to data without replicates")
89
  }
90

91 2
  data_anova <- object$data
92 2
  data_anova$time <- as.factor(data_anova$time)
93 2
  data_anova$variable <- as.factor(data_anova$variable)
94 2
  if (nlevels(data_anova$variable) == 1) {
95 2
    object_2 <- lm(observed ~ time - 1, data = data_anova)
96
  } else {
97 0
    object_2 <- lm(observed ~ variable:time - 1,
98 0
      data = data_anova)
99
  }
100

101 2
  object_2$mkinmod <- list(name = "ANOVA")
102 2
  object_2$err_mod <- "const"
103 2
  sigma_mle <- sqrt(sum(residuals(object_2)^2)/nobs(object_2))
104 2
  object_2$logLik <- sum(dnorm(x = object_2$residuals,
105 2
      mean = 0, sd = sigma_mle, log = TRUE))
106 2
  object_2$data <- object$data # to make the nobs.mkinfit method work
107 2
  object_2$bparms.optim <- coef(object_2)
108 2
  object_2$errparms <- 1 # We have estimated one error model parameter
109 2
  class(object_2) <- "mkinfit"
110

111 2
  lmtest::lrtest.default(object_2, object, name = name_function)
112
}

Read our documentation on viewing source code .

Loading