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
|
#' @importFrom nlme intervals
#' @export
nlme::intervals
#' Confidence intervals for parameters in saem.mmkin objects
#'
#' @param object The fitted saem.mmkin object
#' @param level The confidence level. Must be the default of 0.95 as this is what
#' is available in the saemix object
#' @param backtransform In case the model was fitted with mkin transformations,
#' should we backtransform the parameters where a one to one correlation
#' between transformed and backtransformed parameters exists?
#' @param \dots For compatibility with the generic method
#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the
#' class attribute
#' @export
intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...)
{
if (!identical(level, 0.95)) {
stop("Confidence intervals are only available for a level of 95%")
}
mod_vars <- names(object$mkinmod$diffs)
pnames <- names(object$mean_dp_start)
# Confidence intervals are available in the SaemixObject, so
# we just need to extract them and put them into a list modelled
# after the result of nlme::intervals.lme
conf.int <- object$so@results@conf.int
rownames(conf.int) <- conf.int$name
colnames(conf.int)[2] <- "est."
confint_trans <- as.matrix(conf.int[pnames, c("lower", "est.", "upper")])
# Fixed effects
# In case objects were produced by earlier versions of saem
if (is.null(object$transformations)) object$transformations <- "mkin"
if (object$transformations == "mkin" & backtransform) {
bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
object$transform_rates, object$transform_fractions)
bpnames <- names(bp)
# Transform boundaries of CI for one parameter at a time,
# with the exception of sets of formation fractions (single fractions are OK).
f_names_skip <- character(0)
for (box in mod_vars) { # Figure out sets of fractions to skip
f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
n_paths <- length(f_names)
if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
}
confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
dimnames = list(bpnames, colnames(confint_trans)))
confint_back[, "est."] <- bp
for (pname in pnames) {
if (!pname %in% f_names_skip) {
par.lower <- confint_trans[pname, "lower"]
par.upper <- confint_trans[pname, "upper"]
names(par.lower) <- names(par.upper) <- pname
bpl <- backtransform_odeparms(par.lower, object$mkinmod,
object$transform_rates,
object$transform_fractions)
bpu <- backtransform_odeparms(par.upper, object$mkinmod,
object$transform_rates,
object$transform_fractions)
confint_back[names(bpl), "lower"] <- bpl
confint_back[names(bpu), "upper"] <- bpu
}
}
confint_ret <- confint_back
} else {
confint_ret <- confint_trans
}
attr(confint_ret, "label") <- "Fixed effects:"
# Random effects
sdnames <- intersect(rownames(conf.int), paste("SD", pnames, sep = "."))
corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
ranef_ret <- as.matrix(conf.int[c(sdnames, corrnames), c("lower", "est.", "upper")])
sdnames_ret <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
corrnames_ret <- gsub("Corr\\.(.*)\\.(.*)", "corr(\\1,\\2)", corrnames)
rownames(ranef_ret) <- c(sdnames_ret, corrnames_ret)
attr(ranef_ret, "label") <- "Random effects:"
# Error model
enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1")
err_ret <- as.matrix(conf.int[enames, c("lower", "est.", "upper")])
res <- list(
fixed = confint_ret,
random = ranef_ret,
errmod = err_ret
)
class(res) <- c("intervals.saemix.mmkin", "intervals.lme")
attr(res, "level") <- level
return(res)
}
|
