#' Read NONMEM 7.2+ output into an R object.
#'
#' @param index The NONMEM model index, i.e. the numeric part of the filename assuming it follows the convention 'run123.mod'.
#' @param prefix The NONMEM model prefix, assuming it follows the convention 'run123.mod'. The default is \code{"run"}.
#' @param pathNM The path to the NONMEM output. This should not contain a trailing slash.
#' @param ndig Number of significant digits to use. The default is 3.
#' @param ndigB Number of significant digits to use. The default is 3.
#' @param ndigP Number of digits after the decimal point to use for percentages. The default is 1.
#' @param Pci Asymptotic confidence interval to apply when reporting parameter uncertainty. The default is 95.
#' @param ext NONMEM output file extension. The default is \code{".lst"}.
#' @param extmod NONMEM control stream file extension. The default is \code{"mod"}.
#' @param Pvalues Report P-values for parameters? The default is \code{TRUE}.
#' @param RawCI Report confidence intervals without estimate? The default is \code{FALSE}.
#' @param ... Additional arguments.
#'
#' @return A list containing information extracted from the NONMEM output.
#' @details The output list is composed of the following objects:
#' \itemize{
#' \item{"Theta"}{A data frame describing the structural (fixed-effect) parameters, containing parameter name, estimated value, standard error (SE), coefficient of variation (CV), lower and upper confidence limits (CIL and CIU, based on \code{Pci}), and P-value, calculated as \code{2 * (1 - pnorm(abs(theta/theta.se)))}.}
#' \item{"Eta"}{A data frame describing the interindividual random-effects parameters, containing estimated value, standard error (SE), coefficient of variation (CV, calculated as \code{abs(100*(SE/OMEGA))}), coefficient of variation (EtaCV, calculated as \code{100*sqrt(OMEGA)}), and shrinkage.}
#' \item{"Epsilon"}{A data frame describing the residual random-effects parameters, containing estimated value, standard error (SE), coefficient of variation (CV, calculated as \code{abs(100*(SE/OMEGA))}), coefficient of variation (EtaCV, calculated as \code{100*sqrt(SIGMA)}), and shrinkage.}
#' \item{"CorTheta"}{A data frame containing the correlation matrix for fixed effects (\code{"THETA"}).}
#' \item{"CorOmega"}{A data frame containing the correlation matrix for interindividual random effects (\code{"OMEGA"}).}
#' \item{"CorSigma"}{A data frame containing the correlation matrix for residual random effects (\code{"OMEGA"}).}
#' \item{"OmegaMatrix"}{A data frame containing the \code{"OMEGA"} matrix.}
#' \item{"SigmaMatrix"}{A data frame containing the \code{"OMEGA"} matrix.}
#' \item{"CovMatrixTheta"}{A data frame containing the variance-covariance matrix for structural parameters (\code{THETA}).}
#' \item{"CovMatrix"}{A data frame containing the complete variance-covariance matrix.}
#' \item{"OFV"}{The objective function value.}
#' \item{"ThetaString"}{A data frame containing all relevant fixed-effects parameter information, suitable for use in a table of parameter estimates. Contains parameter name, estimate, standard error, coefficient of variation, combined estimate and asymptotic confidence interval, and P-value.}
#' \item{"EtaString"}{A data frame containing all relevant interindivudal random-effects parameter information, suitable for use in a table of parameter estimates. Contains parameter name, estimate (variance), standard error, coefficient of variation, percentage value (calculated as \code{100*sqrt(OMEGA)}), and shrinkage.}
#' \item{"EpsString"}{A data frame containing all relevant residual random-effects parameter information, suitable for use in a table of parameter estimates. Contains parameter name, estimate (variance), standard error, coefficient of variation, percentage value (calculated as \code{100*sqrt(SIGMA)}), and shrinkage.}
#' \item{"RunTime"}{Run time.}
#' \item{"ConditionN"}{Condition number.}
#'
#'
#' }
#'
#' @seealso NONMEM (\url{https://www.iconplc.com/innovation/nonmem/})
#'
#' @author Rik Schoemaker, \email{rik.schoemaker@@occams.com}
#' @author Justin Wilkins, \email{justin.wilkins@@occams.com}
#'
#' @examples
#' \dontrun{
#' nmOutput <- rnm("run315.lst")
#' }
#'
#' @export
rnm <-
function(index,
prefix = "run",
pathNM,
ndig = 3,
ndigB = 3,
ndigP = 1,
Pci = 95,
ext = ".lst",
extmod = ".mod",
Pvalues = TRUE,
RawCI = FALSE,
...) {
FoldChange <- FALSE
fchdig <- 2
string2num <- function(x)
{
oldopts <- options(warn = -1)
on.exit(options(oldopts))
nc <- nchar(x)
tmp <- substring(x, 1:nc, 1:nc)
spc <- tmp == " "
st <- !spc & c(T, spc[-nc])
end <- !spc & c(spc[-1], T)
as.numeric(substring(x, (1:nc)[st], (1:nc)[end]))
}
isReadable <- function(filename)
{
if(!is.null(version$language)) {
return(file.exists(filename)[1])
}
}
FF <- function(x, dig) {
formatC(x, digits = dig, format = "f")
}
FFP <-
function(number, sigdig = 3) {
res <-
formatC(
signif(number, digits = sigdig),
digits = sigdig,
format = "fg",
flag = "#"
)
res <- sub('\\.$', "", res)
res
}
## from old xpose4 read.lst code
readOutput <- function(filename) {
filename <- file.path(pathNM, paste(prefix, index, ext, sep = ""))
listfile <-
scan(
filename,
sep = "\n",
what = character(),
quiet = TRUE,
quote = ""
)
## Find termination messages
minimStart <- grep("#TERM:", listfile)
minimEnd <- grep("#TERE:", listfile)
if (length(minimStart) > 1) {
warning(
"There seems to be more than one estimation step in this run. Only the last has been extracted.\n",
call. = F
)
minimStart <- max(grep("#TERM:", listfile))
minimEnd <- max(grep("#TERE:", listfile))
}
if (length(minimStart) == 0 || length(minimEnd) == 0) {
termes <- NULL
} else {
termes <- listfile[(minimStart + 1):(minimEnd - 1)]
termes <- substring(termes, 2)
}
## Figure out the name of the raw results file
rawfile <- paste(sub("\\.\\w*$", '', filename), ".ext", sep = "")
if (!isReadable(rawfile)) {
stop(
paste0(
"Could not find the raw results file (",
rawfile,
") for output file: ",
filename,
"\n"
)
)
} else {
extfile <-
scan(
rawfile,
sep = "\n",
what = character(),
quiet = TRUE,
quote = ""
)
TableLines <- max(grep("TABLE NO.", extfile))
rawres <- read.table(rawfile, skip = TableLines, header = T)
}
## Extract OFV
ofv <- rawres$OBJ[rawres$ITERATION == -1000000000]
## Extract parameter estimates
## Get lines with relevant info
finalEstLine <- as.numeric(rawres[rawres$ITERATION == -1000000000, ])
## Extract theta estimates
thetas <- finalEstLine[grep("THETA", names(rawres))]
## Extract omega estimates
omindx <- grep("OMEGA", names(rawres))
omega <- list()
seenOM <- NULL
seenOM[1:100] <- 0
for (om in omindx) {
omnam <- names(rawres)[om]
omcol <- as.numeric(sub("OMEGA\\.\\w*\\.", '', omnam, perl = TRUE))
tmp1 <- sub("OMEGA\\.", '', omnam, perl = TRUE)
omrow <- as.numeric(sub("\\.\\w*\\.$", '', tmp1, perl = TRUE))
seenOM[omrow] <- seenOM[omrow] + 1
if (seenOM[omrow] == 1) {
omega[omrow][omcol] <- finalEstLine[om]
} else {
omega[[omrow]] <- c(omega[[omrow]], as.numeric(finalEstLine[om]))
}
}
## Extract sigma estimates
siindx <- grep("SIGMA", names(rawres))
sigma <- list()
seenSI <- NULL
seenSI[1:100] <- 0
if (length(siindx) != 0) {
for (si in siindx) {
sinam <- names(rawres)[si]
sicol <- as.numeric(sub("SIGMA\\.\\w*\\.", '', sinam, perl = TRUE))
tmp1 <- sub("SIGMA\\.", '', sinam, perl = TRUE)
sirow <- as.numeric(sub("\\.\\w*\\.$", '', tmp1, perl = TRUE))
seenSI[sirow] <- seenSI[sirow] + 1
if (seenSI[sirow] == 1) {
sigma[sirow][sicol] <- finalEstLine[si]
} else {
sigma[[sirow]] <- c(sigma[[sirow]], as.numeric(finalEstLine[si]))
}
}
} else {
sigma <- NULL
}
## Extract condition number:
c_line <- as.numeric(rawres[rawres$ITERATION == -1000000003, ])
if (is.na(c_line[1])) {
conditionN <- NULL
} else {
conditionN <- c_line[2]
}
## Extract line with relevant info
seEstLine <- as.numeric(rawres[rawres$ITERATION == -1000000001, ])
if (is.na(seEstLine[1])) {
sethetas <- NULL
seomega <- NULL
sesigma <- NULL
} else {
## Extract theta estimates
sethetas <- seEstLine[grep("THETA", names(rawres))]
for (it in 1:length(sethetas))
if (sethetas[it] == 1.00000E+10)
sethetas[it] <- 0
## Extract omega estimates
omindx <- grep("OMEGA", names(rawres))
seomega <- list()
seeseOM <- NULL
seeseOM[1:100] <- 0
for (om in omindx) {
omnam <- names(rawres)[om]
omcol <- as.numeric(sub("OMEGA\\.\\w*\\.", '', omnam, perl = TRUE))
tmp1 <- sub("OMEGA\\.", '', omnam, perl = TRUE)
omrow <- as.numeric(sub("\\.\\w*\\.$", '', tmp1, perl = TRUE))
#RS modified: in specific cases we can arrive here and still have NaN omega estimates
if (seEstLine[om] == 1.00000E+10 |
is.na(seEstLine[om]))
seEstLine[om] <- 0
seeseOM[omrow] <- seeseOM[omrow] + 1
if (seeseOM[omrow] == 1) {
seomega[omrow][omcol] <- seEstLine[om]
} else {
seomega[[omrow]] <- c(seomega[[omrow]], as.numeric(seEstLine[om]))
}
}
## Extract sigma estimates
siindx <- grep("SIGMA", names(rawres))
sesigma <- list()
seenseSI <- NULL
seenseSI[1:100] <- 0
if (length(siindx) != 0) {
for (si in siindx) {
sinam <- names(rawres)[si]
sicol <-
as.numeric(sub("SIGMA\\.\\w*\\.", '', sinam, perl = TRUE))
tmp1 <- sub("SIGMA\\.", '', sinam, perl = TRUE)
sirow <- as.numeric(sub("\\.\\w*\\.$", '', tmp1, perl = TRUE))
if (seEstLine[si] == 1.00000E+10)
seEstLine[si] <- 0
seenseSI[sirow] <- seenseSI[sirow] + 1
if (seenseSI[sirow] == 1) {
sesigma[sirow][sicol] <- seEstLine[si]
} else {
sesigma[[sirow]] <- c(sesigma[[sirow]], as.numeric(seEstLine[si]))
}
}
} else {
sesigma <- NULL
}
}
## Correlation matrix
rawfile <- paste(sub("\\.\\w*$", '', filename), ".cor", sep = "")
if (!isReadable(rawfile)) {
warning(
paste0(
"Could not find the correlation matrix file (",
rawfile,
") for list file: ",
filename,
"\n"
)
)
cormat <- NULL
npcormat <- NULL
} else {
extfile <-
scan(
rawfile,
sep = "\n",
what = character(),
quiet = TRUE,
quote = ""
)
TableLines <- max(grep("TABLE NO.", extfile))
rawres <- read.table(rawfile, skip = TableLines, header = T)
cormat <-
as.matrix.data.frame(rawres[1:length(thetas), 2:(length(thetas) + 1)])
for (i in 1:length(thetas)) {
cormat[i, i] <- 1
}
npcormat <- NULL
}
## Covariance matrix
rawfile <- paste(sub("\\.\\w*$", '', filename), ".cov", sep = "")
if (!isReadable(rawfile)) {
warning(
paste0(
"Could not find the covariance matrix file (",
rawfile,
") for list file: ",
filename,
"\n"
)
)
covmat <- NULL
fcovmat <- NULL
} else {
extfile <-
scan(
rawfile,
sep = "\n",
what = character(),
quiet = TRUE,
quote = ""
)
TableLines <- max(grep("TABLE NO.", extfile))
rawres <- read.table(rawfile, skip = TableLines, header = T)
covmat <- as.matrix.data.frame(rawres[1:length(thetas), 2:(length(thetas) + 1)])
fcovmat <- as.matrix.data.frame(rawres[1:nrow(rawres), 2:(nrow(rawres)+1)])
}
## Shrinkage
rawfile <- paste(sub("\\.\\w*$", '', filename), ".shk", sep = "")
if (!isReadable(rawfile)) {
warning(paste0(
"Could not find the shrinkage file (",
rawfile,
") for list file: ",
filename,
"\n"
))
shkmat <- NULL
} else {
extfile <-
scan(
rawfile,
sep = "\n",
what = character(),
quiet = TRUE,
quote = ""
)
TableLines <- max(grep("TABLE NO.", extfile))
rawres <- read.table(rawfile, skip = TableLines, header = TRUE)
shkmat <- as.matrix.data.frame(rawres)
shkmat <- rawres
shrink <- shkmat[shkmat$TYPE == 4, -c(1)]
epsshrink <- shkmat[shkmat$TYPE == 5, -c(1, 2)]
epsshrink <- epsshrink[epsshrink != 0]
if (length(epsshrink) == 0) {
epsshrink = NA
}
names(epsshrink) <- paste("EPS", seq(1, length(epsshrink)), sep = ".")
}
ret.list <- list(term = termes, ofv = ofv,
thetas = thetas, omega = omega, sigma = sigma,
sethetas = sethetas, seomegas = seomega, sesigmas = sesigma,
cormat = cormat, npcormat = npcormat, covmat = covmat, fcovmat = fcovmat,
shrink = shrink, epsshrink = epsshrink, conditionN = conditionN)
}
ncrit <- qnorm((1 - (1 - Pci / 100) / 2))
ncrit90 <- qnorm((1 - (1 - 90 / 100) / 2))
fres <-
as.data.frame(scan(
file = file.path(pathNM, paste(prefix, index, ext, sep = "")),
sep = "\t",
what = character(),
quote = "", quiet=T
),
stringsAsFactors = F)[, 1]
fstart <- fres[grep("Run Start Time:", fres)]
fstop <- fres[grep("Run Stop Time:", fres)]
if (length(fstop) > 0) {
Tstart <- chron::chron(times. = substring(fstart, 17, 30),
format = c("h:m:s"))
Tstop <- chron::chron(times. = substring(fstop, 17, 30), format = c("h:m:s"))
Time <- (as.numeric(Tstop - Tstart) * 60 * 60 * 24)
#if the run starts before and ends after midnight, this breaks down :)
if (Time < 0) {
Time <- Time + 60 * 60 * 24
}
} else {
Time <- 0
}
if (exists("THETA"))
rm(THETA, inherits = TRUE)
if (exists("ETAD"))
rm(ETAD, inherits = TRUE)
if (exists("EPSD"))
rm(EPSD, inherits = TRUE)
f <- readOutput(filename = file.path(pathNM, paste(prefix, index, ext, sep = ""))) # from xpose4
name <- paste(prefix, index, sep = "")
fobj <-
paste(
"MINIMUM VALUE OF OBJECTIVE FUNCTION: " ,
formatC(as.numeric(f$ofv), digits = 3, format = "f"),
sep = ""
)
minimisation <- c(f$term, "", fobj)
# Get number of thetas
nth <- length(f$thetas)
# Get number of etas
net <- length(f$omega)
# Get number of epsilons
nep <- length(f$sigma)
#Read names of THETAs
filename <- file.path(pathNM, paste(prefix, index, ext, sep = ""))
rawfile <- paste(sub("\\.\\w*$", '', filename), extmod, sep = "")
if (!isReadable(rawfile)) {
warning(
paste0(
"Could not find the control stream (",
rawfile,
") for list file: ",
filename,
"\n"
)
)
cormat <- NULL
npcormat <- NULL
} else {
NM <-
as.data.frame(scan(
file = file.path(pathNM, paste(name, extmod, sep = "")),
sep = "\n",
what = character(),
quiet = T
), quote = "")
names(NM)[1] <- "V1"
NM$Parameter <-
substr(NM$V1, ifelse(regexpr(";", NM$V1) == -1, 1000, regexpr(";", NM$V1) +
1), 1000)
# lines with $THETA
NM$TS1 <- as.numeric(regexpr("\\$THETA", NM$V1))
f$Parameter <- trimws(NM[NM$TS1 == 1, "Parameter"])
if (f$Parameter[1] == "")
f$Parameter <- paste("THETA", seq(1:nth), sep = "")
if (length(f$Parameter) != nth)
f$Parameter <- paste("THETA", seq(1:nth), sep = "")
}
# omega as matrix:
if (!is.null(f$omega)) {
net <- length(f$omega)
omegaM <- matrix(0, ncol = net, nrow = net)
k <- 0
for (i in 1:net) {
k <- k + 1
for (j in 1:k)
omegaM[i, j] <- f$omega[[i]][j]
}
k <- 1
if (net > 1) {
for (i in 1:(net - 1)) {
k <- k + 1
for (j in k:net)
omegaM[i, j] <- omegaM[j, i]
}
}
omegaM <- as.matrix.data.frame(omegaM)
} else {
omegaM <- NULL
}
# sigma as matrix:
if (!is.null(f$sigma)) {
nep <- length(f$sigma)
sigmaM <- matrix(0, ncol = nep, nrow = nep)
k <- 0
for (i in 1:nep) {
k <- k + 1
for (j in 1:k)
sigmaM[i, j] <- f$sigma[[i]][j]
}
k <- 1
if (nep > 1) {
for (i in 1:(nep - 1)) {
k <- k + 1
for (j in k:nep)
sigmaM[i, j] <- sigmaM[j, i]
}
}
sigmaM <- as.matrix.data.frame(sigmaM)
} else {
sigmaM <- NULL
}
if (!is.null(f$sethetas)) {
f$cvthetas <- abs((f$sethetas / f$thetas) * 100)
f$cithetasl <- f$thetas - ncrit * f$sethetas
f$cithetasu <- f$thetas + ncrit * f$sethetas
f$pvalue <- 2 * (1 - pnorm(abs(f$thetas / f$sethetas)))
THETAs <-
data.frame(
stringsAsFactors = FALSE,
Parameter = f$Parameter,
Estimate = FFP(f$thetas, ndig),
SE = FFP(f$sethetas, ndig),
CV = paste(formatC(
f$cvthetas, format = "f", digits = 1
), "%", sep = ""),
CI = paste(FFP(f$cithetasl, ndig),
"/",
FFP(f$cithetasu, ndig),
sep = ""
),
CIE = paste(
FFP(f$thetas, ndig),
" (",
FFP(f$cithetasl, ndig),
"/",
FFP(f$cithetasu, ndig),
")",
sep = ""
),
CIT = paste(
FFP(10 ** f$thetas, ndigB),
" (",
FFP(10 ** f$cithetasl, ndigB),
"/",
FFP(10 ** f$cithetasu, ndigB),
")",
sep = ""
),
lCIT = paste(
FFP(exp(f$thetas), ndigB),
" (",
FFP(exp(f$cithetasl), ndigB),
"/",
FFP(exp(f$cithetasu), ndigB),
")",
sep = ""
),
lCITp = paste(
FF(100 * (exp(f$thetas) - 1), ndigP),
"% (",
FF(100 * (exp(f$cithetasl) - 1), ndigP),
"%/",
FF(100 * (exp(f$cithetasu) - 1), ndigP),
"%)",
sep = ""
),
fchCss = paste(
FF(1 / (exp(f$thetas)), fchdig),
" (",
FF(1 / (exp(
f$thetas + ncrit90 * f$sethetas
)), fchdig),
"/",
FF(1 / (exp(
f$thetas - ncrit90 * f$sethetas
)), fchdig),
")",
sep = ""
),
logitCIT = paste(FFP((
exp(f$thetas) / (1 + exp(f$thetas))
), ndigB),
" (", FFP((
exp(f$cithetasl) / (1 + exp(f$cithetasl))
), ndigB),
"/", FFP((
exp(f$cithetasu) / (1 + exp(f$cithetasu))
), ndigB), ")", sep = ""),
P = ifelse(
is.na(f$pvalue),
"",
ifelse(f$pvalue < 0.00001, "< 0.00001", FF(f$pvalue, 5))
)
)
THETAs$tCIT <- THETAs$CI
THETAs$tCIT[regexpr("log", THETAs$Parameter) > -1] <- THETAs$lCIT[regexpr("log", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("logit", THETAs$Parameter) > -1] <- THETAs$logitCIT[regexpr("logit", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("10log", THETAs$Parameter) > -1] <- THETAs$CIT[regexpr("10log", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("logP", THETAs$Parameter) > -1] <- THETAs$lCITp[regexpr("logP", THETAs$Parameter) > -1]
THETAs$fchCss[regexpr("logP", THETAs$Parameter) <= -1] <- ""
THETAs$CI <- THETAs$tCIT
THETAs$tCITE <- THETAs$CIE
THETAs$tCITE[regexpr("log", THETAs$Parameter) > -1] <- THETAs$lCIT[regexpr("log", THETAs$Parameter) > -1]
THETAs$tCITE[regexpr("logit", THETAs$Parameter) > -1] <- THETAs$logitCIT[regexpr("logit", THETAs$Parameter) > -1]
THETAs$tCITE[regexpr("10log", THETAs$Parameter) > -1] <- THETAs$CIT[regexpr("10log", THETAs$Parameter) > -1]
THETAs$tCITE[regexpr("logP", THETAs$Parameter) > -1] <- THETAs$lCITp[regexpr("logP", THETAs$Parameter) > -1]
THETAs$fchCss[regexpr("logP", THETAs$Parameter) <= -1] <- ""
THETAs$CIE <- THETAs$tCITE
THETAs$lCIT <- NULL
THETAs$lCITp <- NULL
THETAs$logitCIT <- NULL
THETAs$tCIT <- NULL
THETAs$tCITE <- NULL
THETAs$CIT <- NULL
THETAs$SE[is.na(f$cvthetas)] <- "(Fixed)"
THETAs$CV[is.na(f$cvthetas)] <- " "
THETAs$CI[is.na(f$cvthetas)] <- " "
THETAs$CIE[is.na(f$cvthetas)] <- " "
names(THETAs)[5] <- paste(as.character(Pci), "% CI", sep = "")
NameCI <- paste("Estimate", " (", names(THETAs)[5], ")", sep = "")
names(THETAs)[6] <- NameCI
if (RawCI == FALSE) {
THETAs[, 5] <- NULL
}
THETA <-
data.frame(
Parameter = f$Parameter,
Estimate = f$thetas,
SE = f$sethetas,
CV = f$cvthetas,
CIL = f$cithetasl,
CIU = f$cithetasu,
P = f$pvalue
)
if (Pvalues == FALSE) {
THETAs$P <- NULL
}
THETAs$"Fold change in Css (90% CI)" <- THETAs$fchCss
THETAs$fchCss <- NULL
if (FoldChange == FALSE) {
THETAs$"Fold change in Css (90% CI)" <- NULL
}
}
if (is.null(f$sethetas)) {
empty <- as.vector(mode = "numeric", matrix(data = NA, nrow = nth, ncol = 1))
emptyS <- as.vector(mode = "character", matrix(data = " ", nrow = nth, ncol = 1))
THETA <-
data.frame(
Parameter = f$Parameter,
Estimate = f$thetas,
SE = empty,
CV = empty,
CIL = empty,
CIU = empty
)
THETAs <-
data.frame(
stringsAsFactors = FALSE,
Parameter = f$Parameter,
Estimate = FFP(f$thetas, ndig),
SE = emptyS,
CV = emptyS,
CI = FFP(f$thetas, ndig),
CIT = FFP(10 ** f$thetas, ndigB),
lCIT = FFP(exp(f$thetas), ndigB),
lCITp = paste(FF(100 * (exp(
f$thetas
) - 1), ndigP), "%", sep = ""),
fchCss = FF(1 / (exp(f$thetas)), fchdig),
logitCIT = FFP((exp(f$thetas) / (1 + exp(
f$thetas
))), ndigB),
P = emptyS
)
THETAs$tCIT <- THETAs$CI
THETAs$tCIT[regexpr("log", THETAs$Parameter) > -1] <-
THETAs$lCIT[regexpr("log", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("logit", THETAs$Parameter) > -1] <-
THETAs$logitCIT[regexpr("logit", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("10log", THETAs$Parameter) > -1] <-
THETAs$CIT[regexpr("10log", THETAs$Parameter) > -1]
THETAs$tCIT[regexpr("logP", THETAs$Parameter) > -1] <-
THETAs$lCITp[regexpr("logP", THETAs$Parameter) > -1]
THETAs$fchCss[regexpr("logP", THETAs$Parameter) <= -1] <- ""
THETAs$CIT <- THETAs$tCIT
THETAs$lCIT <- NULL
THETAs$lCITp <- NULL
THETAs$logitCIT <- NULL
THETAs$tCIT <- NULL
THETAs$CI <- " "
names(THETAs)[5] <- paste(as.character(Pci), "% CI", sep = "")
names(THETAs)[6] <- "Back-transformed estimate"
if (RawCI == FALSE) {
THETAs[, 5] <- NULL
}
if (Pvalues == FALSE) {
THETAs$P <- NULL
}
THETAs$"Fold change in Css (90% CI)" <- THETAs$fchCss
THETAs$fchCss <- NULL
if (FoldChange == FALSE) {
THETAs$"Fold change in Css (90% CI)" <- NULL
}
}
f$fetD <-
as.vector(mode = "numeric", matrix(data = NA, nrow = net, ncol = 1))
f$fetSED <- f$fetD
f$fetCVD <- f$fetD
for (i in 1:net) {
f$fetD[i] <- f$omega[[i]][i]
f$fetDsqrt[i] <- 100 * sqrt(f$fetD[i])
}
if (!is.null(f$seomegas)) {
for (i in 1:net) {
f$fetSED[i] <- f$seomegas[[i]][i]
if (f$fetD[i] == 0) {
f$fetCVD[i] <- 0
}
else {
f$fetCVD[i] <- abs(100 * (f$fetSED[i] / f$fetD[i]))
}
}
}
f$fepD <-
as.vector(mode = "numeric", matrix(data = NA, nrow = nep, ncol = 1))
f$fepSED <- f$fepD
f$fepCVD <- f$fepD
for (i in 1:nep) {
f$fepD[i] <- f$sigma[[i]][i]
f$fepDsqrt[i] <- 100 * sqrt(f$fepD[i])
}
if (!is.null(f$sesigmas)) {
for (i in 1:nep) {
f$fepSED[i] <- f$sesigmas[[i]][i]
if (f$fepD[i] == 0) {
f$fepCVD[i] <- 0
}
else {
f$fepCVD[i] <- abs(100 * (f$fepSED[i] / f$fepD[i]))
}
}
}
Shrinkage <- as.data.frame(t(f$shrink))
SubPops <- dim(Shrinkage)[2]
if (SubPops > 1) {
nams <- paste("Shrink Pop", Shrinkage[1,], sep = "")
} else {
nams <- "Shrinkage"
}
Shrinkage <- Shrinkage[-1,]
Shrinkage <- matrix(unlist(Shrinkage), ncol = SubPops)
ShrinkageS <-
paste(formatC(Shrinkage, format = "f", digits = 1), "%", sep = "")
Shrinkage <- as.data.frame(matrix(Shrinkage, ncol = SubPops))
names(Shrinkage) <- nams
ShrinkageS <- as.data.frame(matrix(ShrinkageS, ncol = SubPops))
names(ShrinkageS) <- nams
EShrinkage <- as.data.frame(t(f$epsshrink))
ESubPops <- dim(EShrinkage)[2]
if (ESubPops > 1) {
nams <- paste("Shrink Pop", EShrinkage[1,], sep = "")
} else {
nams <- "Shrinkage"
}
#EShrinkage <- EShrinkage[-1,]
EShrinkage <- matrix(unlist(EShrinkage), ncol = ESubPops)
EShrinkageS <-
paste(formatC(EShrinkage, format = "f", digits = 1), "%", sep = "")
EShrinkage <- as.data.frame(matrix(EShrinkage, ncol = ESubPops))
names(EShrinkage) <- nams
EShrinkageS <- as.data.frame(matrix(EShrinkageS, ncol = ESubPops))
names(EShrinkageS) <- nams
EpsShrinkage <- EShrinkage
if (!is.null(f$seomegas)) {
ETADs <-
data.frame(
stringsAsFactors = FALSE,
Eta = 1:net,
Estimate = FFP(f$fetD, ndig),
SE = FFP(f$fetSED, ndig),
CV = paste(formatC(
f$fetCVD, format = "f", digits = 1
), "%", sep = "")
,
EtaCV = paste(formatC(
f$fetDsqrt, format = "f", digits = 1
), "%", sep = ""),
Shrinkage = ShrinkageS
)
ETADs$SE[f$fetCVD <= 0] <- "(Fixed)"
ETADs$CV[f$fetCVD <= 0] <- " "
ETADs$EtaCV[f$fetCVD <= 0] <- " "
ETAD <-
data.frame(
Eta = 1:net,
Estimate = f$fetD,
SE = f$fetSED,
CV = f$fetCVD,
EtaCV = f$fetDsqrt,
Shrinkage = Shrinkage
)
}
if (is.null(f$seomegas)) {
ETADs <-
data.frame(
Eta = 1:net,
Estimate = FFP(f$fetD, ndig),
SE = f$fetSED,
CV = f$fetCVD,
EtaCV = paste(formatC(
f$fetDsqrt, format = "f", digits = 1
), "%", sep = ""),
Shrinkage = ShrinkageS
)
ETAD <-
data.frame(
Eta = 1:net,
Estimate = f$fetD,
SE = f$fetSED,
CV = f$fetCVD,
EtaCV = f$fetDsqrt,
Shrinkage = Shrinkage
)
}
names(ETADs)[5] <- '100*sqrt(Omega)'
names(ETADs)[seq(6, dim(ETADs)[2])] <- nams
if (!is.null(f$sesigmas)) {
EPSDs <-
data.frame(
stringsAsFactors = FALSE,
Epsilon = 1:nep,
Estimate = FFP(f$fepD, ndig),
SE = FFP(f$fepSED, ndig),
CV = paste(formatC(
f$fepCVD, format = "f", digits = 1
), "%", sep = "")
,
EpsCV = paste(formatC(
f$fepDsqrt, format = "f", digits = 1
), "%", sep = ""),
Shrinkage = EShrinkageS
)
EPSDs$SE[f$fepCVD <= 0] <- "(Fixed)"
EPSDs$CV[f$fepCVD <= 0] <- " "
EPSDs$EpsCV[f$fepCVD <= 0] <- " "
EPSD <-
data.frame(
Epsilon = 1:nep,
Estimate = f$fepD,
SE = f$fepSED,
CV = f$fepCVD,
EpsCV = f$fepDsqrt,
Shrinkage = EShrinkage
)
}
if (is.null(f$sesigmas)) {
EPSDs <-
data.frame(
Epsilon = 1:nep,
Estimate = FFP(f$fepD, ndig),
SE = f$fepSED,
CV = f$fepCVD,
EpsCV = paste(formatC(
f$fepDsqrt, format = "f", digits = 1
), "%", sep = ""),
Shrinkage = EShrinkageS
)
EPSD <-
data.frame(
Epsilon = 1:nep,
Estimate = f$fepD,
SE = f$fepSED,
CV = f$fepCVD,
EpsCV = f$fepDsqrt,
Shrinkage = EShrinkage
)
}
names(EPSDs)[5] <- '100*sqrt(Omega)'
names(EPSDs)[seq(6, dim(EPSDs)[2])] <- nams
offdiag <- FALSE
f$CorOmega <- matrix(data = NA, nrow = net, ncol = net)
for (i in 1:net) {
for (j in 1:i) {
VARi <- sqrt(f$omega[[i]][i])
VARj <- sqrt(f$omega[[j]][j])
if (VARi == 0 | VARj == 0) {
f$CorOmega[i, j] <- 0
}
else {
f$CorOmega[i, j] <- f$omega[[i]][j] / (VARi * VARj)
}
if ((i != j) & (f$CorOmega[i, j] != 0)) {
offdiag <- TRUE
}
}
}
k <- 1
if (net > 1) {
for (i in 1:(net - 1)) {
k <- k + 1
for (j in k:net) {
f$CorOmega[i, j] <- f$CorOmega[j, i]
}
}
}
tmp <- as.data.frame(matrix(data = NA, nrow = net, ncol = net))
for (i in 1:net) {
for (j in 1:net) {
tmp[i, j] <- formatC(f$CorOmega[i, j], format = "f", digits = 3)
}
}
rownames(tmp) <- paste("Eta", 1:net, sep = "")
colnames(tmp) <- paste("Eta", 1:net, sep = "")
f$CorOmega <- tmp
offdiag <- FALSE
f$CorSigma <- matrix(data = NA, nrow = nep, ncol = nep)
for (i in 1:nep) {
for (j in 1:i) {
VARi <- sqrt(f$sigma[[i]][i])
VARj <- sqrt(f$sigma[[j]][j])
if (VARi == 0 | VARj == 0) {
f$CorSigma[i, j] <- 0
}
else {
f$CorSigma[i, j] <- f$sigma[[i]][j] / (VARi * VARj)
}
if ((i != j) & (f$CorSigma[i, j] != 0)) {
offdiag <- TRUE
}
}
}
k <- 1
if (nep > 1) {
for (i in 1:(nep - 1)) {
k <- k + 1
for (j in k:nep) {
f$CorSigma[i, j] <- f$CorSigma[j, i]
}
}
}
tmp <- as.data.frame(matrix(data = NA, nrow = nep, ncol = nep))
for (i in 1:nep) {
for (j in 1:nep) {
tmp[i, j] <- formatC(f$CorSigma[i, j], format = "f", digits = 3)
}
}
rownames(tmp) <- paste("Eps", 1:nep, sep = "")
colnames(tmp) <- paste("Eps", 1:nep, sep = "")
f$CorSigma <- tmp
params <-
list(
Theta = THETA,
Eta = ETAD,
Epsilon = EPSD,
CorTheta = f$cormat,
CorOmega = f$CorOmega,
CorSigma = f$CorSigma,
#NPCorMatrix = f$npcormat,
OmegaMatrix = omegaM,
SigmaMatrix = sigmaM,
CovMatrixTheta = f$covmat,
CovMatrix = f$fcovmat,
OFV = f$ofv,
ThetaString = THETAs,
EtaString = ETADs,
EpsString = EPSDs,
#EpsShrink = f$epsshrink,
RunTime = Time,
conditionN = f$conditionN
)
return(params)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.