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
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/loftest.R
\name{loftest}
\alias{loftest}
\alias{loftest.mkinfit}
\title{Lack-of-fit test for models fitted to data with replicates}
\usage{
loftest(object, ...)
\method{loftest}{mkinfit}(object, ...)
}
\arguments{
\item{object}{A model object with a defined loftest method}
\item{\dots}{Not used}
}
\description{
This is a generic function with a method currently only defined for mkinfit
objects. It fits an anova model to the data contained in the object and
compares the likelihoods using the likelihood ratio test
\code{\link[lmtest]{lrtest.default}} from the lmtest package.
}
\details{
The anova model is interpreted as the simplest form of an mkinfit model,
assuming only a constant variance about the means, but not enforcing any
structure of the means, so we have one model parameter for every mean
of replicate samples.
}
\examples{
\dontrun{
test_data <- subset(synthetic_data_for_UBA_2014[[12]]$data, name == "parent")
sfo_fit <- mkinfit("SFO", test_data, quiet = TRUE)
plot_res(sfo_fit) # We see a clear pattern in the residuals
loftest(sfo_fit) # We have a clear lack of fit
#
# We try a different model (the one that was used to generate the data)
dfop_fit <- mkinfit("DFOP", test_data, quiet = TRUE)
plot_res(dfop_fit) # We don't see systematic deviations, but heteroscedastic residuals
# therefore we should consider adapting the error model, although we have
loftest(dfop_fit) # no lack of fit
#
# This is the anova model used internally for the comparison
test_data_anova <- test_data
test_data_anova$time <- as.factor(test_data_anova$time)
anova_fit <- lm(value ~ time, data = test_data_anova)
summary(anova_fit)
logLik(anova_fit) # We get the same likelihood and degrees of freedom
#
test_data_2 <- synthetic_data_for_UBA_2014[[12]]$data
m_synth_SFO_lin <- mkinmod(parent = list(type = "SFO", to = "M1"),
M1 = list(type = "SFO", to = "M2"),
M2 = list(type = "SFO"), use_of_ff = "max")
sfo_lin_fit <- mkinfit(m_synth_SFO_lin, test_data_2, quiet = TRUE)
plot_res(sfo_lin_fit) # not a good model, we try parallel formation
loftest(sfo_lin_fit)
#
m_synth_SFO_par <- mkinmod(parent = list(type = "SFO", to = c("M1", "M2")),
M1 = list(type = "SFO"),
M2 = list(type = "SFO"), use_of_ff = "max")
sfo_par_fit <- mkinfit(m_synth_SFO_par, test_data_2, quiet = TRUE)
plot_res(sfo_par_fit) # much better for metabolites
loftest(sfo_par_fit)
#
m_synth_DFOP_par <- mkinmod(parent = list(type = "DFOP", to = c("M1", "M2")),
M1 = list(type = "SFO"),
M2 = list(type = "SFO"), use_of_ff = "max")
dfop_par_fit <- mkinfit(m_synth_DFOP_par, test_data_2, quiet = TRUE)
plot_res(dfop_par_fit) # No visual lack of fit
loftest(dfop_par_fit) # no lack of fit found by the test
#
# The anova model used for comparison in the case of transformation products
test_data_anova_2 <- dfop_par_fit$data
test_data_anova_2$variable <- as.factor(test_data_anova_2$variable)
test_data_anova_2$time <- as.factor(test_data_anova_2$time)
anova_fit_2 <- lm(observed ~ time:variable - 1, data = test_data_anova_2)
summary(anova_fit_2)
}
}
\seealso{
lrtest
}
|
