#' Get estimated individual and population parameters
#'
#' Get the individual individual parameters, the population parameters with the population covariates
#' and the population parameters with the individual covariates.
#' @return a list of data frames.
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getEstimatedIndividualParameters2()
#'
#' # r is a list with elements "saem", "conditionalMean", "conditionalSD", "conditionalMode",
#' # "popPopCov" and "popIndCov"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getEstimatedIndividualParameters2 <- function() {
initRsmlx()
ind.param <- mlx.getEstimatedIndividualParameters()
N <- nrow(ind.param$saem)
pop.param <- mlx.getEstimatedPopulationParameters()
pop.param <- pop.param[grep("_pop",names(pop.param))]
df <- as.data.frame(matrix(pop.param,nrow=N,ncol=length(pop.param),byrow=TRUE))
names(df) <- gsub("_pop","",names(pop.param))
ind.param$popPopCov <- data.frame(id=ind.param$saem["id"],df)
rand.eff <- mlx.getEstimatedRandomEffects()
if (!is.null(ind.param$conditionalMean)) {
ip <- ind.param$conditionalMean
re <- rand.eff$conditionalMean
} else if (!is.null(ind.param$conditionalMode)) {
ip <- ind.param$conditionalMode
re <- rand.eff$conditionalMode
} else
stop("The conditional mean or the conditional model should have been computed", call.=FALSE)
ind.dist <- mlx.getIndividualParameterModel()$distribution
var.param <- names(ind.dist)
ind.param$popIndCov <- ind.param$popPopCov
for (nj in var.param) {
dj <- tolower(ind.dist[nj])
yj <- ip[[nj]]
rj <- re[[paste0("eta_",nj)]]
if (dj == "normal") {
yjc <- yj-rj
} else if (dj == "lognormal") {
yjc <- exp(log(yj)-rj)
} else if (dj == "logitnormal") {
yjc <- 1/(1+exp(-log(yj/(1-yj))+rj))
} else if (dj == "probitnormal") {
yjc <- pnorm(qnorm(yj) - rj)
}
ind.param$popIndCov[nj] <- yjc
}
return(ind.param)
}
#--------------------------------------------------------------------
#' Get estimated predictions
#'
#' Get the individual predictions obtained with the estimated individual parameters :
#' @return a list of data frames (one data frame per output).
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getEstimatedPredictions() # r is a list with elements "y1" and "y2"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getEstimatedPredictions <- function() {
initRsmlx()
ip <- getEstimatedIndividualParameters2()
obs.info <- mlx.getObservationInformation()
df <- list()
obsNames <- names(obs.info$mapping)
contNames <- obsNames[obs.info$type[obsNames]=="continuous"]
for (j in seq_along(contNames)) {
name <- contNames[j]
df[j] <- obs.info[name]
df[[j]][name] <- NULL
}
f.pop1 <- mlx.computePredictions(ip$popPopCov)
for (j in seq_along(contNames)) {df[[j]]$popPopCov <- f.pop1[[j]]}
f.pop2 <- mlx.computePredictions(ip$popIndCov)
for (j in seq_along(contNames)) {df[[j]]$popIndCov <- f.pop2[[j]]}
if (!is.null(ip$conditionalMean)) {
f.mean <- mlx.computePredictions(ip$conditionalMean)
for (j in seq_along(contNames)) {df[[j]]$conditionalMean <- f.mean[[j]]}
}
if (!is.null(ip$conditionalMode)) {
f.mode <- mlx.computePredictions(ip$conditionalMode)
for (j in seq_along(contNames)) {df[[j]]$conditionalMode <- f.mode[[j]]}
}
names(df) <- names(f.pop1)
return(df)
}
#--------------------------------------------------------------------
#' Get estimated residuals
#'
#' Get the residuals computed from the individual predictions obtained
#' with the estimated individual parameters:
#' @return a list of data frames (one data frame per output).
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getEstimatedResiduals() # r is a list with elements "y1" and "y2"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getEstimatedResiduals <- function() {
initRsmlx()
df <- getEstimatedPredictions()
obs.info <- mlx.getObservationInformation()
nip <- c("popPopCov", "popIndCov", "conditionalMean", "conditionalMode")
obsNames <- names(obs.info$mapping)
contNames <- obsNames[obs.info$type[obsNames]=="continuous"]
nobs <- length(contNames)
error.model <- mlx.getContinuousObservationModel()$errorModel
error.dist <- mlx.getContinuousObservationModel()$distribution
ep <- error.parameter()
pop.param <- mlx.getEstimatedPopulationParameters()
param.error <- list()
for (j in seq_along(contNames)) {
name <- contNames[j]
dfj <- df[[j]]
ij <- which(names(dfj) %in% nip)
yoj <- replicate(length(ij), obs.info[[name]][[name]])
erj <- tolower(error.dist[j])
if (erj=="normal") {
ypj <- dfj[,ij]
} else if (erj=="lognormal") {
ypj <- log(dfj[,ij])
yoj <- log(yoj)
} else if (erj=="logitnormal") {
limiti <- mlx.getContinuousObservationModel()$limits[[names(error.dist)[j]]]
ypj <- log((dfj[,ij]-limiti[1])/(limiti[2]-dfj[,ij]))
yoj <- log((yoj-limiti[1])/(limiti[2]-yoj))
}
epj <- pop.param[ep[[j]]]
a <- epj[grep("a", names(epj))]
if (length(a)==0)
a <- 0
b <- epj[grep("b", names(epj))]
if (length(b)==0)
b <- 0
c <- epj[grep("c", names(epj))]
if (length(c)==0)
c <- 1
pei <- c(a, b, c)
if (error.model[j]=="combined2")
dfj <- (yoj - ypj)/sqrt(pei[1]^2 + (pei[2]*ypj^pei[3])^2)
else
dfj <- (yoj - ypj)/(pei[1] + pei[2]*ypj^pei[3])
df[[j]][names(dfj)] <- dfj
}
names(df) <- contNames
return(df)
}
#--------------------------------------------------------------------
#' Get simulated predictions
#'
#' Get the individual predictions obtained with the simulated individual parameters :
#' @return a list of data frames (one data frame per output).
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getSimulatedPredictions() # r is a list with elements "Cc" and "E"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getSimulatedPredictions <- function() {
initRsmlx()
sip <- mlx.getSimulatedIndividualParameters()
if (is.null(sip$rep))
sip$rep <- 1
nrep <- max(sip$rep)
obs.info <- mlx.getObservationInformation()
df <- list()
obsNames <- names(obs.info$mapping)
contNames <- obsNames[obs.info$type[obsNames]=="continuous"]
pred <- mlx.getContinuousObservationModel()$prediction
for (j in seq_along(contNames)) {
name <- contNames[j]
df[j] <- obs.info[name]
df[[j]][name] <- NULL
df[[j]][pred[name]] <- 0
df[[j]] <- cbind(rep=1, df[[j]])
}
col.el <- which(!(names(sip) %in% c("rep","id")))
res <- list()
for (irep in seq_len(nrep)) {
parami <- subset(sip, rep==irep)[,col.el]
fi <- mlx.computePredictions(parami)
for (j in seq_along(contNames)) {
name <- contNames[j]
df[[j]][pred[name]] <- fi[[pred[name]]]
df[[j]]["rep"] <- irep
if (irep==1)
res[[j]] <- df[[j]]
else
res[[j]] <- rbind(res[[j]], df[[j]])
}
}
names(res) <- pred[contNames]
return(res)
}
#--------------------------------------------------------------------
#' Get simulated residuals
#'
#' Get the residuals computed from the individual predictions obtained
#' with the simulated individual parameters:
#' @return a list of data frames (one data frame per output).
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getSimulatedResiduals() # r is a list with elements "y1" and "y2"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getSimulatedResiduals <- function() {
initRsmlx()
df <- getSimulatedPredictions()
obs.info <- mlx.getObservationInformation()
obsNames <- names(obs.info$mapping)
contNames <- obsNames[obs.info$type[obsNames]=="continuous"]
error.model <- mlx.getContinuousObservationModel()$errorModel
error.dist <- mlx.getContinuousObservationModel()$distribution
ep <- error.parameter()
pop.param <- mlx.getEstimatedPopulationParameters()
param.error <- list()
nrep <- max(df[[1]]["rep"])
for (j in seq_along(contNames)) {
name <- contNames[j]
dfj <- df[[j]]
ij <- which(names(dfj) == names(df)[j])
yoj <- rep(obs.info[[name]][[name]],nrep)
erj <- tolower(error.dist[j])
if (erj=="normal") {
ypj <- dfj[,ij]
} else if (erj=="lognormal") {
ypj <- log(dfj[,ij])
yoj <- log(yoj)
} else if (erj=="logitnormal") {
limiti <- mlx.getContinuousObservationModel()$limits[[names(error.dist)[j]]]
ypj <- log((dfj[,ij]-limiti[1])/(limiti[2]-dfj[,ij]))
yoj <- log((yoj-limiti[1])/(limiti[2]-yoj))
}
epj <- pop.param[ep[[j]]]
a <- epj[grep("a", names(epj))]
if (length(a)==0)
a <- 0
b <- epj[grep("b", names(epj))]
if (length(b)==0)
b <- 0
c <- epj[grep("c", names(epj))]
if (length(c)==0)
c <- 1
pei <- c(a, b, c)
if (error.model[j]=="combined2")
dfj <- (yoj - ypj)/sqrt(pei[1]^2 + (pei[2]*ypj^pei[3])^2)
else
dfj <- (yoj - ypj)/(pei[1] + pei[2]*ypj^pei[3])
df[[j]][,ij] <- dfj
names(df[[j]])[ij] <- "residual"
}
names(df) <- contNames
return(df)
}
#--------------------------------------------------------------------
#' Get estimated covariance and correlation matrices
#'
#' @return a list of two matrices.
#' @examples
#' \dontrun{
#' # Assume that the Monolix project "warfarinPKPD.mlxtran" has been loaded
#' r = getEstimatedCovarianceMatrix() # r is a list with elements "cor.matrix" and "cov.matrix"
#'
#' # See http://monolix.lixoft.com/rsmlx/newconnectors/ for more detailed examples
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#' }
#' @export
getEstimatedCovarianceMatrix <- function() {
initRsmlx()
param <- mlx.getEstimatedPopulationParameters()
pname <- names(param)
i.omega <- grep("^omega_",pname)
if (length(i.omega)>0) {
oest <- 1
oname <- gsub("^omega_","",pname[i.omega])
omega <- param[i.omega]
} else {
i.omega <- grep("^omega2_",pname)
oname <- gsub("^omega2_","",pname[i.omega])
omega <- sqrt(param[i.omega])
oest <- 2
}
i.corr <- grep("^corr_",pname)
d <- length(i.omega)
c <- param[i.corr]
R <- diag(rep(1,d))
rownames(R) <- colnames(R) <- oname
if (length(c)>0) {
for (j in 1:length(c)) {
cj <- names(c)[j]
sj <- strsplit(cj,"_")[[1]]
R[sj[3],sj[2]] <- R[sj[2],sj[3]] <- c[j]
}
}
C <- diag(omega)%*%R%*%diag(omega)
return(list(cor.matrix=R, cov.matrix=C))
}
#--------------------------------------------------------------------
error.parameter <- function(project=NULL) {
if (is.null(project)) {
dp <- mlx.getProjectSettings()$directory
if (!is.null(dp))
project <- paste0(dp,".mlxtran")
# project <- paste0(basename(dp),".mlxtran")
}
if (!file.exists(project))
stop("Enter a valid project", call.=FALSE)
con = file(project, open = "r")
lines = readLines(con, warn=FALSE)
close(con)
lines <- lines[grep("errorModel", lines)]
r <- list()
if (length(lines)>0) {
lines <- gsub(" ","",lines)
for (k in (1: length(lines))) {
lk <- lines[k]
i1 <- regexpr("\\(",lk)
i2 <- regexpr("\\)",lk)
r[[k]] <- strsplit(substr(lk,(i1+1),(i2-1)),",")[[1]]
}
}
return(r)
}
#-------------------------------------------------
mlx.getFit <- function() {
project <- mlx.getProjectSettings()$project
con <- file(project, open = "r")
lines <- readLines(con, warn=FALSE)
close(con)
lines <- gsub("\\;.*","",lines)
l.fit <- grep("<FIT>", lines)
lines <- lines[l.fit:length(lines)]
l.data <- lines[grep("data", lines)[1]]
l.model <- lines[grep("model", lines)[1]]
ll <- gsub(" ","",l.data)
ll <- gsub("data=","",ll)
if (grepl("\\{",ll)) {
i1 <- regexpr("\\{",ll)
i2 <- regexpr("\\}",ll)
ll <- substr(ll,i1+1,i2-1)
data.names <- strsplit(ll,",")[[1]]
} else {
data.names <- ll
}
ll <- gsub(" ","",l.model)
ll <- gsub("model=","",ll)
if (grepl("\\{",ll)) {
i1 <- regexpr("\\{",ll)
i2 <- regexpr("\\}",ll)
ll <- substr(ll,i1+1,i2-1)
model.names <- strsplit(ll,",")[[1]]
} else {
model.names <- ll
}
return(list(data=data.names, model=model.names))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.