Nothing
R/StatisticalModel.R
In PFIM: Population Fisher Information Matrix
#################################################################################
#' Class "StatisticalModel"
#'
#' @description Class \code{StatisticalModel} represents a statistical model.
#'
#' @name StatisticalModel-class
#' @aliases StatisticalModel
#' @docType class
#' @include ModelEquations.R
#' @export StatisticalModel
#' @exportClass StatisticalModel
#' @section Mathematical description of the statistical model:
#' \describe{
#' \item{-}{\deqn{y = f(x, \theta) + g*\epsilon}, this part is considered in class \linkS4class{Response}}
#' \item{-}{\deqn{f(x, \theta)}, this part is considered in class \linkS4class{Response}}
#' \item{-}{\deqn{\theta={(\mu, \omega)}}, this part is considered in class \linkS4class{ModelParameter}}
#' \item{-}{\deqn{\epsilon}, this is a slot of the object \code{NormalDistribution} with mean = 0 and covariate_matrix = I}
#' }
#'
#' @section Objects from the class: \code{StatisticalModel} objects are typically created by calls to \code{StatisticalModel} and contain the following slots:
#'
#' @section Slots for \code{StatisticalModel} objects:
#' \describe{
#' \item{\code{modelEquations}:}{An object from the class \code{ModelEquations}}
#' \item{\code{responses}:}{A list of objects of type Responses -> f(x, theta)}
#' \item{\code{correlations}:}{A list giving all the covariables.}
#' \item{\code{model_parameters}:}{A list giving all the parameters of the models.}
#' }
#################################################################################
StatisticalModel<-setClass(
Class = "StatisticalModel",
representation = representation(
modelEquations = "ModelEquations",
responses = "list", # list of object of type Response
correlations = "list", # list of all covariables
model_parameters = "list", # list of all parameters of the models
computeFIM = "logical",
parametersForSolverOde = "list",
variables="list"),
prototype = prototype(
parametersForSolverOde=list( .relStep = .Machine$double.eps^(1/3),
atol = 1e-6,
rtol = 1e-6),
computeFIM = TRUE),
validity = function(object)
{
return(TRUE)
}
)
# Initialize method
setMethod(
f = "initialize",
signature = "StatisticalModel",
definition = function(.Object, modelEquations, responses, correlations,
model_parameters, computeFIM, parametersForSolverOde, variables)
{
if(!missing(modelEquations))
.Object@modelEquations <- modelEquations
if(!missing(responses))
.Object@responses <- responses
if(!missing(correlations))
.Object@correlations <- correlations
if(!missing(model_parameters))
.Object@model_parameters <- model_parameters
if(!missing(computeFIM))
.Object@computeFIM <- computeFIM
if(!missing(parametersForSolverOde))
.Object@parametersForSolverOde <- parametersForSolverOde
if(!missing(variables))
.Object@variables <- variables
validObject(.Object)
return(.Object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Set parameters for the ode solver
#'
#' @name setParametersOdeSolver
#' @param object A \code{StatisticalModel} object.
#' @param value A list giving the values of the parameters.
#' @return The \code{StatisticalModel} object with the new parameters for the ode solver.
setGeneric("setParametersOdeSolver",
function(object, value)
{
standardGeneric("setParametersOdeSolver")
}
)
setMethod( f="setParametersOdeSolver",
signature="StatisticalModel",
definition = function(object, value)
{
if ( length ( value ) !=0 )
{
if( length( value$.relStep ) !=0 ){
object@parametersForSolverOde$.relStep = value$.relStep
}
if( length( value$atol ) !=0 ){
object@parametersForSolverOde$atol = value$atol
}
if( length( value$rtol ) !=0 ){
object@parametersForSolverOde$rtol = value$rtol
}
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get parameters for the ode solver
#'
#' @name getParametersOdeSolver
#' @param object A \code{StatisticalModel} object.
#' @return The parameters for the ode solver.
setGeneric("getParametersOdeSolver",
function(object)
{
standardGeneric("getParametersOdeSolver")
}
)
setMethod( f="getParametersOdeSolver",
signature="StatisticalModel",
definition = function(object)
{
parametersOdeSolver = list( .relStep = object@parametersForSolverOde$.relStep,
atol = object@parametersForSolverOde$ato,
rtol = object@parametersForSolverOde$rtol )
return( parametersOdeSolver )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Add a response to a statistical model.
#'
#' @name addResponse
#' @param object A \code{StatisticalModel} object.
#' @param value A character string giving the name of the response to add.
#' @return The \code{StatisticalModel} object with the added response.
setGeneric("addResponse",
function(object, value)
{
standardGeneric("addResponse")
}
)
setMethod( f="addResponse",
signature="StatisticalModel",
definition = function(object, value)
{
object@responses[[ value@name ]] <- value
validObject(object)
return(object)
}
)
# ---------------------------------------------------------------------------------------------------------------
#' Add responses to a statistical model.
#'
#' @name addResponses
#' @param object A \code{StatisticalModel} object.
#' @param listOfResponses A list of character string giving the names of the responses to add.
#' @return The \code{StatisticalModel} object with the added responses.
setGeneric("addResponses",
function(object, listOfResponses)
{
standardGeneric("addResponses")
}
)
setMethod( f="addResponses",
signature="StatisticalModel",
definition = function(object, listOfResponses)
{
for( i in 1:length( listOfResponses ) ){
response = listOfResponses[[i]]
object@responses[[ response@name ]] <- response
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Set the correlation.
#'
#' @name defineCorrelation
#' @param object A \code{StatisticalModel} object.
#' @param correlationlist ...
#' @return Return correlationlist
#
setGeneric("defineCorrelation",
function(object, correlationlist)
{
standardGeneric("defineCorrelation")
}
)
setMethod(f="defineCorrelation",
signature= "StatisticalModel",
definition=function(object, correlationlist )
{
object@correlations = correlationlist
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Compute the residual variance thanks to the function g of the model error.
#'
#' @name CalculatedResidualVariance
#' @param objectStatisticalModel A \code{StatisticalModel} object.
#' @param objectModelError A \code{ModelError} object.
#' @param x_i variable x_i of the model error.
#' #' @return CalculatedResidualVariance
setGeneric("CalculatedResidualVariance",
function(objectStatisticalModel, objectModelError, x_i ) #
{
standardGeneric("CalculatedResidualVariance")
}
)
setMethod(f="CalculatedResidualVariance",
signature= "StatisticalModel",
definition=function(objectStatisticalModel, objectModelError, x_i)
{
calculated_residual_variance = g(objectModelError, x_i)
return(calculated_residual_variance)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Check the parameters in the model equations.
#'
#' @name checkParameterInEquations
#' @param object A \code{StatisticalModel} object.
#' @return The \code{ModelParameter} objects of the model parameters in the model equation.
setGeneric("checkParameterInEquations",
function(object)
{
standardGeneric("checkParameterInEquations")
}
)
setMethod(f="checkParameterInEquations",
signature= "StatisticalModel",
definition=function(object)
{
modelParameters = getModelParameters( object )
modelEquations = getEquationsStatisticalModel( object )
parameterInEquation = list()
dC = getDerivatives( modelEquations )
allParametersinEquation = unlist(lapply( dC, function(x) all.vars( x )))
nameParameters = names( modelParameters )
for ( nameParameter in nameParameters )
{
if ( nameParameter %in% allParametersinEquation)
{
parameterInEquation[[nameParameter]] = "TRUE"
}
else
{
modelParameters[[nameParameter]] = NULL
}
}
return( modelParameters )
})
# -------------------------------------------------------------------------------------------------------------------
#' Set the parameters for the evaluation of the model.
#'
#' @rdname setParametersForEvaluateModel
#' @param object A \code{StatisticalModel} object.
#' @param administrations An \code{Administration} object.
#' @param sampling_times A \code{SamplingTimes} object.
#' @param cond_init A list for the initial conditions of the \code{StatisticalModel} object.
#' @return A list containing the parameters used for the model evaluation.
setGeneric("setParametersForEvaluateModel",
function(object, administrations, sampling_times, cond_init)
{
standardGeneric("setParametersForEvaluateModel")
}
)
setMethod(f = "setParametersForEvaluateModel",
signature="StatisticalModel",
definition= function(object, administrations, sampling_times, cond_init)
{
# model parameters
samplingTimesModel = list()
sampleTimeResponse = list()
inputsModel = list()
### get model equations
modelEquations = getEquationsStatisticalModel( object )
classModel = class( modelEquations )
responses = getResponsesStatisticalModel( object )
### responses name and number
responseNames = names( responses )
numberOfResponse = length( responseNames )
numberOfResponsePK = length( administrations)
numberOfResponsePD = numberOfResponse - numberOfResponsePK
responsePKNames = responseNames[1:numberOfResponsePK]
responsePDNames = setdiff(responseNames,responsePKNames)
inputsModel$responsePKNames = responsePKNames
inputsModel$responsePDNames = responsePDNames
inputsModel$responseNames = c( responsePKNames, responsePDNames )
### get sampling times for each responses
for ( responseName in responseNames )
{
sampleTimeResponse[[responseName ]] = getSampleTime( sampling_times[[ responseName ]] )
samplingTimesModel = append( samplingTimesModel, sampleTimeResponse[[ responseName ]] )
}
samplingTimesModel = sort( unique( unlist( samplingTimesModel ) ) )
inputsModel$sampleTimeResponse = sampleTimeResponse
### parameters
Tinf = list()
dose = list()
time_dose = list()
tau = list()
timeDoseForOde = c()
for ( response in responsePKNames)
{
Tinf[[response]] = getTinf( administrations[[response]])
dose[[response]] = getAmountDose( administrations[[response]])
time_dose[[response]] = getTimeDose( administrations[[response]])
tau[[response]] = getTau( administrations[[response]])
if ( length( tau[[response]] )==0 )
{
tau[[response]] = 0
}
timeDoseForOde = c(timeDoseForOde, time_dose[[response]])
inputsModel$Tinf[[response]] = Tinf[[response]]
inputsModel$dose[[response]] = dose[[response]]
inputsModel$time_dose[[response]] = time_dose[[response]]
inputsModel$tau[[response]] = tau[[response]]
inputsModel$timeDoseForOde[[response]] = sort(unique(c( timeDoseForOde, max(samplingTimesModel))))
}
###### model parameters
## check case model parameter in cond_init but no inequations
modelParameters = getModelParameters(object)
if (classModel == "ModelODEquations")
{
modelParametersInEquation = checkParameterInEquations( object )
}else
{
modelParametersInEquation = modelParameters
}
parameterInEquation = list()
for( parameter in modelParameters )
{
name = getNameModelParameter( parameter )
inputsModel[[name]] = getMu( parameter )
assign(name,inputsModel[[name]] )
}
inputsModel$model_parameters = modelParameters
# set initial conditions
# cond_init & case for doses bolus model
# dose in cond init and not in ode equations
cond_init_ode = list()
if (classModel == "ModelODEquations")
{
dC = getDerivatives(modelEquations)
doseRespPKInEquation = list()
doseRespPKInCondinit = list()
for ( nameResponsePK in responsePKNames )
{
nameDose = paste0("dose_",nameResponsePK)
doseRespPKInEquation[[nameDose]] = sapply( dC, function(x) all( grepl( nameDose, x) ) )
doseRespPKInCondinit[[nameDose]] = vapply( cond_init, function(x) all( grepl( nameDose, x) ),logical(1) )
for ( i in 1:length( doseRespPKInEquation[[nameDose]] ) )
{
for ( j in 1:length( doseRespPKInCondinit[[nameDose]] ) )
{
if ( ( doseRespPKInEquation[[nameDose]][i] == FALSE ) & ( doseRespPKInCondinit[[nameDose]][j] == TRUE ) )
{
assign( nameDose, dose[[nameResponsePK]] )
}
}
}
}
}
for ( nameCondiInit in names( cond_init ) )
{
cond_init_ode[[nameCondiInit]] = eval( cond_init[[nameCondiInit]] )
}
cond_init_ode = c( unlist( cond_init_ode ) )
# model equations
inputsModel$modelEquations = modelEquations
# variable of the model
inputsModel$variablesNames = names( cond_init )
# get parameters for ode solver
parametersOdeSolver = getParametersOdeSolver( object )
inputsModel$parametersOdeSolver = parametersOdeSolver
return(list(samplingTimesModel = samplingTimesModel,
cond_init_ode = cond_init_ode,
inputsModel = inputsModel,
modelParametersInEquation = modelParametersInEquation ) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Parameters used for computing the model gradient by finite-differences.
#'
#' @rdname parametersForComputingGradient
#' @param object A \code{StatisticalModel} object.
#' @return A list containing the parameters used for computing the gradient of the model.
setGeneric("parametersForComputingGradient",
function(object)
{
standardGeneric("parametersForComputingGradient")
}
)
setMethod(f = "parametersForComputingGradient",
signature="StatisticalModel",
definition= function(object)
{
minAbsPar = 0
parametersOdeSolver = getParametersOdeSolver( object )
.relStep = parametersOdeSolver$.relStep
listPars = getModelParameters( object )
valuePars = unlist(lapply(listPars, slot, name = "mu"))
.relStep = pmax(abs(valuePars), minAbsPar) * .relStep
pars <- as.numeric(valuePars)
npar <- length(pars)
incr <- pmax(abs(pars), minAbsPar) * .relStep
baseInd <- diag(npar)
frac <- c(1, incr, incr^2)
cols <- list(0, baseInd, -baseInd)
for ( i in seq_along(pars)[ -npar ] ) {
cols <- c( cols, list( baseInd[ , i ] + baseInd[ , -(1:i) ] ) )
frac <- c( frac, incr[ i ] * incr[ -(1:i) ] )
}
indMat <- do.call( "cbind", cols)
shifted <- pars + incr * indMat
indMat <- t(indMat)
Xcols <- list(1, indMat, indMat^2)
for ( i in seq_along(pars)[ - npar ] ) {
Xcols <- c( Xcols, list( indMat[ , i ] * indMat[ , -(1:i) ] ) )
}
return( list( Xcols = Xcols, shifted = shifted, frac = frac ) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Evaluation for the model, analytic, ode, infusion
#'
#' @rdname EvaluationModel
#' @param object A \code{StatisticalModel} object.
#' @param samplingTimesModel A vector containing the sampling times of the model.
#' @param cond_init_ode A vector containing the initial conditions for an ODE model.
#' @param inputsModel A list containing the parameters used for the model evaluation.
#' @param parametersGradient A list containing the parameters used for the evaluation of the gradient of the model
#' @param computeFIM A boolean giving TRUE id the FIM is compute, FALSE otherwise.
#' @return A list contaning a dataframe giving the results of the evaluation and a list giving the gradient of the model.
setGeneric("EvaluationModel",
function(object, samplingTimesModel,cond_init_ode, inputsModel, parametersGradient, computeFIM )
{
standardGeneric("EvaluationModel")
})
setMethod(f="EvaluationModel",
signature= "StatisticalModel",
definition=function(object, samplingTimesModel,cond_init_ode, inputsModel, parametersGradient, computeFIM )
{
modelEquations = getEquationsStatisticalModel( object )
classModel = class( modelEquations )
if ( classModel == "ModelInfusionODEquations" )
{
evaluationModel = EvaluateModelODEInfusion( modelEquations, samplingTimesModel,cond_init_ode,
inputsModel, parametersGradient, computeFIM )
} else if ( classModel == "ModelODEquations" ) {
evaluationModel = EvaluateModelODE( modelEquations, samplingTimesModel,cond_init_ode,
inputsModel, parametersGradient, computeFIM )
} else if ( classModel == "ModelEquations" ) {
evaluationModel = EvaluateModel( modelEquations, samplingTimesModel,inputsModel, computeFIM )
}else if ( classModel == "ModelInfusionEquations" ) {
evaluationModel = EvaluateModelInfusion( modelEquations, samplingTimesModel,inputsModel, computeFIM )
}
return( evaluationModel )
})
# -------------------------------------------------------------------------------------------------------------------
#' Evaluate an \code{StatisticalModel} object.
#'
#' @rdname Evaluate
#' @param object A \code{StatisticalModel} object.
#' @param administrations An \code{Administration} object.
#' @param sampling_times A \code{SamplingTimes} object.
#' @param cond_init A list for the initial conditions of the \code{StatisticalModel} object.
#' @param fim \code{FIM} object.
#' @return A \code{fim} object giving the Fisher Information Matrix of the \code{StatisticalModel} object.
setGeneric("Evaluate",
function(object, administrations, sampling_times, cond_init, fim )
{
standardGeneric("Evaluate")
})
setMethod(f="Evaluate",
signature= "StatisticalModel",
definition=function(object, administrations, sampling_times, cond_init, fim )
{
# initial parameters
MF_var = NA
V_total = NA
position = 1
errorVariances = list()
sigmaDerivatives = list()
modelParameters = getModelParameters( object )
modelEquations = getEquationsStatisticalModel( object )
modelResponses = getResponsesStatisticalModel( object )
# set parameters for evaluating the model
parametersGradient = parametersForComputingGradient( object )
parametersForEvaluateModel = setParametersForEvaluateModel( object,
administrations,
sampling_times, cond_init )
samplingTimesModel = parametersForEvaluateModel$samplingTimesModel
cond_init_ode = parametersForEvaluateModel$cond_init_ode
inputsModel = parametersForEvaluateModel$inputsModel
# for parameters not in equation but in the initial conditions
modelParametersInEquation = parametersForEvaluateModel$modelParametersInEquation
indexModelParametersInEquation = which( modelParameters %in% modelParametersInEquation )
indexModelParametersNotInEquation = which( modelParameters %in% modelParametersInEquation == FALSE )
# set model parameters for the evaluation
object@model_parameters <- inputsModel$model_parameters[indexModelParametersInEquation]
# -----------------------------------
# evaluate model
# -----------------------------------
# sampling time for plot
if ( object@computeFIM == FALSE ){
minSamplingTimesModel = 0
maxSamplingTimesModel = max( samplingTimesModel )
samplingTimesModel = sort( unique( c( samplingTimesModel,
linspace( minSamplingTimesModel, maxSamplingTimesModel, 500 ) ) ) )
evaluationModel = EvaluationModel( object,
samplingTimesModel,
cond_init_ode,
inputsModel,
parametersGradient, object@computeFIM )
predictedResponses = evaluationModel$evaluationResponse
responseGradient = evaluationModel$responseGradient
# gradient of the parameters in model equations
for ( response in modelResponses )
{
nameResponse = response@name
if( length( modelParametersInEquation ) == 1 )
{
responseGradient[[nameResponse]] = responseGradient[[nameResponse]][indexModelParametersInEquation]
}else{
responseGradient[[nameResponse]] = responseGradient[[nameResponse]][,indexModelParametersInEquation]
}
}
return( list( fim = fim,
predictedResponses = predictedResponses, sensitivityIndicesModel = responseGradient ) )
}
else
{
# sampling time of the designs
evaluationModel = EvaluationModel( object,
samplingTimesModel,
cond_init_ode,
inputsModel,
parametersGradient, object@computeFIM )
# predicted responses and gradients
predictedResponses = evaluationModel$evaluationResponse
responseGradient = evaluationModel$responseGradient
responseAllGradient = evaluationModel$responseAllGradient
nTotalSamplingTimes = length( unlist( inputsModel$sampleTimeResponse ) )
# gradient of the parameters in model equations
for ( response in modelResponses )
{
nameResponse = response@name
responseGradient[[nameResponse]] = responseAllGradient[,indexModelParametersInEquation]
}
responseAllGradient = responseAllGradient[,indexModelParametersInEquation]
modelParameters = getModelParameters( object )
modelParameters = modelParameters[indexModelParametersInEquation]
responseAllGradient = as.matrix( responseAllGradient )
## gradient for parameter not in eqs but in initial conditions
# parameter not in equations
modelParametersNotInEquation = inputsModel$model_parameters[indexModelParametersNotInEquation]
variablesNames = inputsModel$variablesNames
responsePKNames = inputsModel$responsePKNames
valueGradient = list()
# assign values parameters not in equation & PK doses
for( modelParameterNotInEquation in modelParametersNotInEquation )
{
name = getNameModelParameter( modelParameterNotInEquation )
assign( name, getMu( modelParameterNotInEquation ) )
}
for ( responsePKName in responsePKNames )
{
assign( paste0("dose_",responsePKName ), inputsModel$dose[[responsePKName]] )
}
# compute and eval the gradients for variables in initial conditions
for ( variableName in variablesNames )
{
equationCondInit = cond_init[[variableName]]
gradientExpression = sapply( names( modelParametersNotInEquation ),
function(x){ D( equationCondInit, x ) } )
valueGradient = unlist( sapply( gradientExpression, function(x){eval(x)} ) )
}
if ( length( valueGradient ) !=0 )
{
gradientParameterNotInEquations = rep(0, dim( responseAllGradient )[1])
gradientParameterNotInEquations[1] = valueGradient
responseAllGradient = cbind( responseAllGradient, gradientParameterNotInEquations )
}
# evaluate FIM: individual, population, Bayesian
object@model_parameters <- inputsModel$model_parameters
modelParameters = getModelParameters( object )
for ( response in modelResponses )
{
nameResponse = response@name
samplingTimesResponse = getSampleTime( sampling_times[[nameResponse]] )
predResp = predictedResponses[predictedResponses$time %in% samplingTimesResponse,nameResponse]
errorModelDerivatives = EvaluateErrorModelDerivatives( response, predResp )
errorVariances = append( errorVariances, bdiag( errorModelDerivatives$errorVariance ) )
errorVariances = bdiag( errorVariances )
for( errorModelParameter in errorModelDerivatives$sigmaDerivatives )
{
emptysigmaDerivatives = matrix( 0, ncol = nTotalSamplingTimes, nrow = nTotalSamplingTimes )
nTime = length( inputsModel$sampleTimeResponse[[ nameResponse ]] )
range <- position:( position + nTime - 1 )
emptysigmaDerivatives[ range, range ] <- errorModelParameter
# Add the MF_beta matrix for the parameters for each response
# add the MF_var in different blocks for each response
sigmaDerivatives = c( sigmaDerivatives, list( emptysigmaDerivatives ) )
}
position <- position + nTime
}
# get fixed parameters
dataFixedParameters = getFixedParameters( object )
indexMuNoFixedParameters = dataFixedParameters$indexMuNoFixedParameters
indexMuFixedParameters = dataFixedParameters$indexMuFixedParameters
indexOmegaNoFixedParameters = dataFixedParameters$indexOmegaNoFixedParameters
indexOmegaFixedParameters = dataFixedParameters$indexOmegaFixedParameters
indexZeroOmegaParameters = dataFixedParameters$indexZeroOmegaParameters
# names for the parameters
namesMuNoFixedParameters = names( modelParameters[indexMuNoFixedParameters] )
namesOmegaNoFixedParameters = names( modelParameters[indexOmegaNoFixedParameters] )
# names for the columns/rows for the FIM
muNames = paste0( "\u03bc_", namesMuNoFixedParameters )
omegaNames = paste0( "\u03c9\u00B2_", namesOmegaNoFixedParameters )
indexMuOmegaNoFixedParameters = sort( unique ( c( indexMuNoFixedParameters ) ) )
modelParameters = getModelParameters( object )
modelParametersNoFixedMu = modelParameters
# vector for sigma of error model
sigmaNames = c()
parameterNames = c()
# PopulationFim and IndividualFim
for (response in modelResponses )
{
if ( is( fim, "PopulationFim" ) )
{
MF = PopulationFIMEvaluateVariance(response,
modelEquations,
modelParametersNoFixedMu,
administrations,
sampling_times,
responseAllGradient,
errorVariances, sigmaDerivatives )
parameterNames <- c( muNames, omegaNames )
}
else
if ( is( fim, "IndividualFim" ) )
{
MF = IndividualFIMEvaluateVariance(response,
modelEquations,
modelParametersNoFixedMu,
administrations,
sampling_times,
responseAllGradient,
errorVariances, sigmaDerivatives )
parameterNames <- c( muNames )
}
### Add the MF_beta matrix for the parameters for each response
### add the MF_var in different blocks for each response
if ( is.matrix( MF_var ))
{
MF_var = bdiag( MF_var, MF$MF_var )
V_total = bdiag( V_total, MF$V )
}
else
{
### Create the first matrix and the first block
MF_var = bdiag( MF$MF_var )
V_total = bdiag( MF$V )
}
sigmaNames <- c( sigmaNames, getSigmaNames( response ) )
}
# compute MF_beta
if ( length( indexMuNoFixedParameters )!=0 )
{
responseAllGradient = responseAllGradient[, indexMuNoFixedParameters ]
}
MF_beta = t( responseAllGradient ) %*% solve( V_total ) %*% responseAllGradient
# Bayesian FIM
if ( is( fim, "BayesianFim" ) )
{
mu = c()
omega2 = c()
nParameter = 0
for (parameter in modelParameters )
{
omegaValue = getOmega( parameter )
omega2 = c( omega2, omegaValue^2 )
distribution = getDistribution( parameter )
if ( class( distribution )[1] == "LogNormalDistribution" )
{
muValue = getMu( parameter )
mu = c( mu, muValue )
}
else
if ( class( distribution )[1] == "NormalDistribution" )
{
mu = c( mu, 1 )
}
else
print("That distribution type is not supported for bayesian FIM")
}
# case: parameters fixed=FALSE & omega != 0
indexParameters = intersect ( indexOmegaNoFixedParameters, indexMuNoFixedParameters )
if (length(omega2)==1)
{
Omega = as.matrix(omega2[indexParameters])
meanParametersMatrix = as.matrix( mu[indexParameters] )
}else
{
Omega = diag( omega2[indexParameters] )
meanParametersMatrix = diag( mu[indexParameters] )
}
# subsetting beta matrix: no fixed mu and no fixed omega
colnames( MF_beta ) <- c( muNames )
rownames( MF_beta ) <- c( muNames )
muNames = paste0( "\u03bc_",names(modelParameters[indexParameters]))
MF_beta = as.matrix(MF_beta[muNames,muNames ])
# compute MF_beta
if( dim( MF_beta )[1] == 0 )
{
identity = c(1.0)
}else{
identity = diag( dim( MF_beta )[1] )
}
MF_beta = t( identity ) %*% MF_beta %*% identity + solve( meanParametersMatrix %*% Omega %*% meanParametersMatrix )
# set col and row names and Fisher matrix
colnames( MF_beta ) <- c( muNames )
rownames( MF_beta ) <- c( muNames )
fim <- `setMfisher<-`( fim, MF_beta )
# parameters indices used for plot
fim@parametersIndicesMuFIM = as.vector( indexParameters )
# set mu and omega : no fixed no omega=0
fim<-setOmega( fim, Omega )
fim<-setMu( fim, meanParametersMatrix )
return( list( fim = fim,
predictedResponses = predictedResponses,
sensitivityIndicesModel = responseGradient ) )
}
# case : PopulationFim & parameters with omega^2=0
# case parameter omega == 0 fixed and mu fixed
MF_i_total = bdiag( MF_beta, MF_var )
if ( is( fim, "PopulationFim" ) )
{
if ( length( indexOmegaFixedParameters ) !=0 ){
MF_var = MF_var[-c(indexOmegaFixedParameters),-c(indexOmegaFixedParameters)]
}
}
MF_i_total = bdiag( MF_beta, MF_var )
parameterNames = c( parameterNames, sigmaNames )
colnames( MF_i_total ) <- c( parameterNames )
rownames( MF_i_total ) <- c( parameterNames )
fim <- `setMfisher<-`( fim, MF_i_total )
# parameters indices used for plot
fim@parametersIndicesMuFIM = indexMuNoFixedParameters
fim@parametersIndicesOmegaFIM = indexOmegaNoFixedParameters
return( list( fim = fim,
predictedResponses = predictedResponses,
sensitivityIndicesModel = responseGradient ) )
}
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the fixed and non fixed model parameters.
#'
#' @name getFixedParameters
#' @param object A \code{StatisticalModel} object.
#' @return A list that contains the name and indices of the fixed parameters.
setGeneric("getFixedParameters",
function(object)
{
standardGeneric("getFixedParameters")
}
)
setMethod("getFixedParameters",
"StatisticalModel",
function(object)
{
modelParameters = getModelParameters( object )
# indices parameters mu and omega fixed or not
indexMuFixedParameters = which( lapply( modelParameters, slot, "fixedMu" ) == TRUE )
indexOmegaFixedParameters =which( lapply( modelParameters, slot, "fixedOmega" ) == TRUE )
indexMuNoFixedParameters = which( lapply( modelParameters, slot, "fixedMu" ) == FALSE )
indexOmegaNoFixedParameters = which( lapply( modelParameters, slot, "fixedOmega" ) == FALSE )
indexZeroOmegaParameters = which( lapply( modelParameters, slot, "omega" ) == 0 )
indexZeroMuParameters = which( lapply( modelParameters, slot, "mu" ) == 0 )
indexNoZeroOmegaParameters = which( lapply( modelParameters, slot, "omega" ) != 0 )
indexNoZeroMuParameters = which( lapply( modelParameters, slot, "mu" ) != 0 )
indexMuFixedParameters = sort( unique( c( indexMuFixedParameters, indexZeroMuParameters ) ) )
indexOmegaFixedParameters = sort( unique( c( indexOmegaFixedParameters, indexZeroOmegaParameters ) ) )
return ( output = list( indexMuFixedParameters = indexMuFixedParameters,
indexOmegaFixedParameters = indexOmegaFixedParameters,
indexMuNoFixedParameters = indexMuNoFixedParameters,
indexOmegaNoFixedParameters = indexOmegaNoFixedParameters,
indexZeroOmegaParameters = indexZeroOmegaParameters,
indexZeroMuParameters = indexZeroMuParameters,
indexNoZeroOmegaParameters = indexNoZeroOmegaParameters,
indexNoZeroMuParameters = indexNoZeroMuParameters) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Define model equations
#'
#' @name defineModelEquations
#' @param object A \code{StatisticalModel} object.
#' @param equations An expression giving the equations of the model.
#' @return The \code{StatisticalModel} object with the equations.
setGeneric("defineModelEquations",
function(object, equations)
{
standardGeneric("defineModelEquations")
}
)
setMethod("defineModelEquations",
"StatisticalModel",
function(object, equations )
{
if ( class( equations ) %in% c( "PKModel","PKPDModel" ) ) {
object@modelEquations = getEquations( equations )
}else{
object@modelEquations = equations
}
return( object )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define a parameter of a statistical model.
#'
#' @name defineParameter
#' @param object A \code{StatisticalModel} object.
#' @param parameter An expression giving a parameter of the \code{StatisticalModel} object.
#' @return Return \code{StatisticalModel} object with new parameters.
setGeneric("defineParameter",
function(object, parameter)
{
standardGeneric("defineParameter")
}
)
setMethod("defineParameter",
"StatisticalModel",
function(object, parameter )
{
parameterName = getNameModelParameter( parameter )
object@model_parameters[[ parameterName ]] = parameter
object@model_parameters = object@model_parameters[ order( names( object@model_parameters ) ) ]
return( object )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define the parameters of a statistical model.
#'
#' @name defineParameters
#' @param object A \code{StatisticalModel} object.
#' @param listOfParameters A list of string giving the parameters of the \code{StatisticalModel} object.
#' @return Return \code{StatisticalModel} object with new parameters.
setGeneric("defineParameters",
function(object, listOfParameters)
{
standardGeneric("defineParameters")
}
)
setMethod("defineParameters",
"StatisticalModel",
function(object, listOfParameters )
{
listOfParameters = listOfParameters[order(names(listOfParameters))]
for (i in 1:length( listOfParameters)) {
parameter = listOfParameters[[i]]
parameterName = getNameModelParameter( parameter )
object@model_parameters[[ parameterName ]] = parameter
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Show the content of a \code{StatisticalModel} object.
#'
#' @rdname show
#' @param object \code{StatisticalModel} object.
#' @return Display the responses name of the model equations,
#' the ordinary derivatives of the model equations (for an ODE model), the parameters of the model.
setMethod("show",
signature = "StatisticalModel",
definition = function(object)
{
modelEquations = getEquationsStatisticalModel(object)
cat("*********************** Statistical model ******************************\n")
# ------------------------------------------------------------------------
# ModelEquations
# ------------------------------------------------------------------------
if( class(object@modelEquations ) %in% c("ModelEquations")){
cat( "\n **** Model equation:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
equation = getEquation( object@modelEquations, nameResponse )
cat( paste0(nameResponse ," = ", equation,"\n\n" ) )
}
cat( " **** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat( " Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat( "\n" )
}
}
# ------------------------------------------------------------------------
# ModelInfusionEquations
# ------------------------------------------------------------------------
else if( class(object@modelEquations ) %in% c("ModelInfusionEquations")){
cat( "\n **** Model equation:","\n\n" )
responses = getResponsesStatisticalModel(object)
namesResponses = names( responses )
for (nameResponse in namesResponses )
{
nameEquations = paste0( c( "DuringInfusion_","AfterInfusion_" ), nameResponse )
for (nameEquation in nameEquations )
{
equation = getEquation( object@modelEquations, nameEquation )
cat( paste0( paste0( nameEquation ) ," = ", equation,"\n\n" ) )
}
}
cat( " **** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat( " Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat( "\n" )
}
}
# ------------------------------------------------------------------------
# ModelODEquations
# ------------------------------------------------------------------------
else if(class(object@modelEquations) %in% c("ModelODEquations"))
{
derivatives <- getDerivatives( object@modelEquations )
cat("\n")
cat("**** Model responses:","\n\n" )
for( response in object@responses)
{
nameResponse = getNameResponse( response )
equation = getEquation( object@modelEquations, nameResponse )
cat( paste0( nameResponse ," = ", equation,"\n" ) )
}
cat("\n")
cat("**** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat("Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat("\n")
}
cat("**** Ordinary derivatives of model equation:\n\n")
for(i in 1:length(derivatives))
{
cat( names(derivatives)[i]," = ")
print(derivatives[[i]][[1]])
}
cat( "\n" )
}
# ------------------------------------------------------------------------
# ModelInfusionODEquations
# ------------------------------------------------------------------------
else if(class(modelEquations) %in% c("ModelInfusionODEquations"))
{
responses = getResponsesStatisticalModel(object)
namesResponses = names( responses )
cat("\n")
cat("**** Model responses:","\n\n" )
equationsDuringInfusion = modelEquations@duringInfusionResponsesEquations
equationsAfterInfusion = modelEquations@afterInfusionResponsesEquations
for (nameResponse in namesResponses )
{
nameEquationsDuringInfusion = paste0( c( "DuringInfusion_" ), nameResponse )
nameEquationsAfterInfusion = paste0( c( "AfterInfusion_" ), nameResponse )
cat( paste0( nameEquationsDuringInfusion ," = ", equationsDuringInfusion[[nameEquationsDuringInfusion]],"\n" ) )
cat( paste0( nameEquationsAfterInfusion ," = ", equationsAfterInfusion[[nameEquationsAfterInfusion]],"\n" ) )
}
cat("\n")
cat("**** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat("Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat("\n")
}
cat("\n")
cat("**** Ordinary derivatives of model equation:\n\n")
derivatives = getDerivatives( modelEquations )
namesDerivatives = names( derivatives )
for (namesDerivative in namesDerivatives )
{
cat( paste0( namesDerivative ," = ", derivatives[[namesDerivative]],"\n" ) )
}
}
# ------------------------------------------------------------------------
# parameters
# ------------------------------------------------------------------------
nameCol <- c()
muCol <- c()
omegaCol <- c()
distributionCol <- c()
for(parameter in object@model_parameters)
{
nameCol <- c( nameCol, getNameModelParameter(parameter))
muCol <- c(muCol, getMu(parameter))
omegaCol <- c(omegaCol, getOmega(parameter))
distributionCol <- c(distributionCol, class(getDistribution(parameter))[1])
}
parameterFrame <- matrix(c( as.character(muCol),
as.character(omegaCol),
distributionCol),ncol=3)
colnames(parameterFrame) <- c("\u03bc", "\u03c9", "Distribution" )
rownames(parameterFrame) <- nameCol
cat("\n")
cat("**** Model parameters: \n")
print(parameterFrame)
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the model parameters of a statistical model.
#'
#' @name getModelParameters
#' @param object \code{getModelParameters} object.
#' @return A list giving the model parameters.
setGeneric("getModelParameters",
function(object)
{
standardGeneric("getModelParameters")
}
)
setMethod("getModelParameters",
"StatisticalModel",
function(object)
{
return(object@model_parameters)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get the responses of a statistical model.
#'
#' @name getResponsesStatisticalModel
#' @param object A \code{getResponsesStatisticalModel} object.
#' @return A list giving the responses of a statistical model.
setGeneric("getResponsesStatisticalModel",
function(object)
{
standardGeneric("getResponsesStatisticalModel")
}
)
setMethod("getResponsesStatisticalModel",
"StatisticalModel",
function(object)
{
return(object@responses)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get the equations of a statistical model.
#'
#' @name getEquationsStatisticalModel
#' @param object A \code{StatisticalModel} object.
#' @return A \code{ModelEquationsg} object giving the equations of the \code{StatisticalModel} object.
setGeneric(
"getEquationsStatisticalModel",
function(object) {
standardGeneric("getEquationsStatisticalModel")
})
setMethod("getEquationsStatisticalModel",
signature("StatisticalModel"),
function(object)
{
return( object@modelEquations)
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the SE and RSE of the parameters.
#'
#' @name getErrorModelStandardErrors
#' @param object A \code{StatisticalModel} object.
#' @param fim A \code{Fim} object giving the Fisher Information Matrix.
#' @return A dataframe giving he SE and RSE of the parameters.
setGeneric("getErrorModelStandardErrors",
function(object, fim)
{
standardGeneric("getErrorModelStandardErrors")
}
)
setMethod("getErrorModelStandardErrors",
"StatisticalModel",
function(object, fim)
{
se <- getSE(fim)
sigmaCol <- c()
rseCol <- c()
responses = getResponsesStatisticalModel( object )
for(response in responses)
{
modelError <- getModelError(response)
sigmaCol <- c( sigmaCol, getSigmaValues(modelError))
}
if (length(se) >1 ){
if (length(se) != length(sigmaCol) ){
se <- se[-c(1:(length(se)-length(sigmaCol)))]
}else{
se <- (sort(se))
}}
rseCol <- c(rseCol, se/sigmaCol*100 )
parameterFrame <- data.frame( value = sigmaCol, se = se, "rse"= rseCol, row.names=NULL, check.names=FALSE)
return(parameterFrame)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define a variable in a statistical model.
#'
#' @name defineVariable
#' @param .Object A \code{StatisticalModel} object.
#' @param variable A character string giving the variable to be defined.
#' @return The \code{StatisticalModel} object with the new variable.
setGeneric("defineVariable",
function(.Object, variable )
{
standardGeneric("defineVariable")
}
)
setMethod(f="defineVariable",
signature= "StatisticalModel",
definition=function(.Object, variable)
{
variableName = variable@name
.Object@variables[[ variableName ]] = variable
return(.Object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define variables in a statistical model.
#'
#' @name defineVariables
#' @param .Object A \code{StatisticalModel} object.
#' @param listOfVariables A list of character string giving the variables to be defined.
#' @return The \code{StatisticalModel} object with the new variables.
setGeneric( "defineVariables",
function(.Object, listOfVariables ) #
{
standardGeneric( "defineVariables" )
}
)
setMethod(f = "defineVariables",
signature = "StatisticalModel",
definition = function(.Object, listOfVariables )
{
for (i in 1:length( listOfVariables)) {
modelVariable = ModelVariable( listOfVariables[[i]] )
variableName = modelVariable@name
.Object@variables[[ variableName ]] = modelVariable
}
return(.Object)
}
)
##########################################################################################################
# END Class "StatisticalModel"
##########################################################################################################
Try the PFIM package in your browser
Any scripts or data that you put into this service are public.
PFIM documentation built on June 24, 2022, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#################################################################################
#' Class "StatisticalModel"
#'
#' @description Class \code{StatisticalModel} represents a statistical model.
#'
#' @name StatisticalModel-class
#' @aliases StatisticalModel
#' @docType class
#' @include ModelEquations.R
#' @export StatisticalModel
#' @exportClass StatisticalModel
#' @section Mathematical description of the statistical model:
#' \describe{
#' \item{-}{\deqn{y = f(x, \theta) + g*\epsilon}, this part is considered in class \linkS4class{Response}}
#' \item{-}{\deqn{f(x, \theta)}, this part is considered in class \linkS4class{Response}}
#' \item{-}{\deqn{\theta={(\mu, \omega)}}, this part is considered in class \linkS4class{ModelParameter}}
#' \item{-}{\deqn{\epsilon}, this is a slot of the object \code{NormalDistribution} with mean = 0 and covariate_matrix = I}
#' }
#'
#' @section Objects from the class: \code{StatisticalModel} objects are typically created by calls to \code{StatisticalModel} and contain the following slots:
#'
#' @section Slots for \code{StatisticalModel} objects:
#' \describe{
#' \item{\code{modelEquations}:}{An object from the class \code{ModelEquations}}
#' \item{\code{responses}:}{A list of objects of type Responses -> f(x, theta)}
#' \item{\code{correlations}:}{A list giving all the covariables.}
#' \item{\code{model_parameters}:}{A list giving all the parameters of the models.}
#' }
#################################################################################
StatisticalModel<-setClass(
Class = "StatisticalModel",
representation = representation(
modelEquations = "ModelEquations",
responses = "list", # list of object of type Response
correlations = "list", # list of all covariables
model_parameters = "list", # list of all parameters of the models
computeFIM = "logical",
parametersForSolverOde = "list",
variables="list"),
prototype = prototype(
parametersForSolverOde=list( .relStep = .Machine$double.eps^(1/3),
atol = 1e-6,
rtol = 1e-6),
computeFIM = TRUE),
validity = function(object)
{
return(TRUE)
}
)
# Initialize method
setMethod(
f = "initialize",
signature = "StatisticalModel",
definition = function(.Object, modelEquations, responses, correlations,
model_parameters, computeFIM, parametersForSolverOde, variables)
{
if(!missing(modelEquations))
.Object@modelEquations <- modelEquations
if(!missing(responses))
.Object@responses <- responses
if(!missing(correlations))
.Object@correlations <- correlations
if(!missing(model_parameters))
.Object@model_parameters <- model_parameters
if(!missing(computeFIM))
.Object@computeFIM <- computeFIM
if(!missing(parametersForSolverOde))
.Object@parametersForSolverOde <- parametersForSolverOde
if(!missing(variables))
.Object@variables <- variables
validObject(.Object)
return(.Object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Set parameters for the ode solver
#'
#' @name setParametersOdeSolver
#' @param object A \code{StatisticalModel} object.
#' @param value A list giving the values of the parameters.
#' @return The \code{StatisticalModel} object with the new parameters for the ode solver.
setGeneric("setParametersOdeSolver",
function(object, value)
{
standardGeneric("setParametersOdeSolver")
}
)
setMethod( f="setParametersOdeSolver",
signature="StatisticalModel",
definition = function(object, value)
{
if ( length ( value ) !=0 )
{
if( length( value$.relStep ) !=0 ){
object@parametersForSolverOde$.relStep = value$.relStep
}
if( length( value$atol ) !=0 ){
object@parametersForSolverOde$atol = value$atol
}
if( length( value$rtol ) !=0 ){
object@parametersForSolverOde$rtol = value$rtol
}
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get parameters for the ode solver
#'
#' @name getParametersOdeSolver
#' @param object A \code{StatisticalModel} object.
#' @return The parameters for the ode solver.
setGeneric("getParametersOdeSolver",
function(object)
{
standardGeneric("getParametersOdeSolver")
}
)
setMethod( f="getParametersOdeSolver",
signature="StatisticalModel",
definition = function(object)
{
parametersOdeSolver = list( .relStep = object@parametersForSolverOde$.relStep,
atol = object@parametersForSolverOde$ato,
rtol = object@parametersForSolverOde$rtol )
return( parametersOdeSolver )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Add a response to a statistical model.
#'
#' @name addResponse
#' @param object A \code{StatisticalModel} object.
#' @param value A character string giving the name of the response to add.
#' @return The \code{StatisticalModel} object with the added response.
setGeneric("addResponse",
function(object, value)
{
standardGeneric("addResponse")
}
)
setMethod( f="addResponse",
signature="StatisticalModel",
definition = function(object, value)
{
object@responses[[ value@name ]] <- value
validObject(object)
return(object)
}
)
# ---------------------------------------------------------------------------------------------------------------
#' Add responses to a statistical model.
#'
#' @name addResponses
#' @param object A \code{StatisticalModel} object.
#' @param listOfResponses A list of character string giving the names of the responses to add.
#' @return The \code{StatisticalModel} object with the added responses.
setGeneric("addResponses",
function(object, listOfResponses)
{
standardGeneric("addResponses")
}
)
setMethod( f="addResponses",
signature="StatisticalModel",
definition = function(object, listOfResponses)
{
for( i in 1:length( listOfResponses ) ){
response = listOfResponses[[i]]
object@responses[[ response@name ]] <- response
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Set the correlation.
#'
#' @name defineCorrelation
#' @param object A \code{StatisticalModel} object.
#' @param correlationlist ...
#' @return Return correlationlist
#
setGeneric("defineCorrelation",
function(object, correlationlist)
{
standardGeneric("defineCorrelation")
}
)
setMethod(f="defineCorrelation",
signature= "StatisticalModel",
definition=function(object, correlationlist )
{
object@correlations = correlationlist
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Compute the residual variance thanks to the function g of the model error.
#'
#' @name CalculatedResidualVariance
#' @param objectStatisticalModel A \code{StatisticalModel} object.
#' @param objectModelError A \code{ModelError} object.
#' @param x_i variable x_i of the model error.
#' #' @return CalculatedResidualVariance
setGeneric("CalculatedResidualVariance",
function(objectStatisticalModel, objectModelError, x_i ) #
{
standardGeneric("CalculatedResidualVariance")
}
)
setMethod(f="CalculatedResidualVariance",
signature= "StatisticalModel",
definition=function(objectStatisticalModel, objectModelError, x_i)
{
calculated_residual_variance = g(objectModelError, x_i)
return(calculated_residual_variance)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Check the parameters in the model equations.
#'
#' @name checkParameterInEquations
#' @param object A \code{StatisticalModel} object.
#' @return The \code{ModelParameter} objects of the model parameters in the model equation.
setGeneric("checkParameterInEquations",
function(object)
{
standardGeneric("checkParameterInEquations")
}
)
setMethod(f="checkParameterInEquations",
signature= "StatisticalModel",
definition=function(object)
{
modelParameters = getModelParameters( object )
modelEquations = getEquationsStatisticalModel( object )
parameterInEquation = list()
dC = getDerivatives( modelEquations )
allParametersinEquation = unlist(lapply( dC, function(x) all.vars( x )))
nameParameters = names( modelParameters )
for ( nameParameter in nameParameters )
{
if ( nameParameter %in% allParametersinEquation)
{
parameterInEquation[[nameParameter]] = "TRUE"
}
else
{
modelParameters[[nameParameter]] = NULL
}
}
return( modelParameters )
})
# -------------------------------------------------------------------------------------------------------------------
#' Set the parameters for the evaluation of the model.
#'
#' @rdname setParametersForEvaluateModel
#' @param object A \code{StatisticalModel} object.
#' @param administrations An \code{Administration} object.
#' @param sampling_times A \code{SamplingTimes} object.
#' @param cond_init A list for the initial conditions of the \code{StatisticalModel} object.
#' @return A list containing the parameters used for the model evaluation.
setGeneric("setParametersForEvaluateModel",
function(object, administrations, sampling_times, cond_init)
{
standardGeneric("setParametersForEvaluateModel")
}
)
setMethod(f = "setParametersForEvaluateModel",
signature="StatisticalModel",
definition= function(object, administrations, sampling_times, cond_init)
{
# model parameters
samplingTimesModel = list()
sampleTimeResponse = list()
inputsModel = list()
### get model equations
modelEquations = getEquationsStatisticalModel( object )
classModel = class( modelEquations )
responses = getResponsesStatisticalModel( object )
### responses name and number
responseNames = names( responses )
numberOfResponse = length( responseNames )
numberOfResponsePK = length( administrations)
numberOfResponsePD = numberOfResponse - numberOfResponsePK
responsePKNames = responseNames[1:numberOfResponsePK]
responsePDNames = setdiff(responseNames,responsePKNames)
inputsModel$responsePKNames = responsePKNames
inputsModel$responsePDNames = responsePDNames
inputsModel$responseNames = c( responsePKNames, responsePDNames )
### get sampling times for each responses
for ( responseName in responseNames )
{
sampleTimeResponse[[responseName ]] = getSampleTime( sampling_times[[ responseName ]] )
samplingTimesModel = append( samplingTimesModel, sampleTimeResponse[[ responseName ]] )
}
samplingTimesModel = sort( unique( unlist( samplingTimesModel ) ) )
inputsModel$sampleTimeResponse = sampleTimeResponse
### parameters
Tinf = list()
dose = list()
time_dose = list()
tau = list()
timeDoseForOde = c()
for ( response in responsePKNames)
{
Tinf[[response]] = getTinf( administrations[[response]])
dose[[response]] = getAmountDose( administrations[[response]])
time_dose[[response]] = getTimeDose( administrations[[response]])
tau[[response]] = getTau( administrations[[response]])
if ( length( tau[[response]] )==0 )
{
tau[[response]] = 0
}
timeDoseForOde = c(timeDoseForOde, time_dose[[response]])
inputsModel$Tinf[[response]] = Tinf[[response]]
inputsModel$dose[[response]] = dose[[response]]
inputsModel$time_dose[[response]] = time_dose[[response]]
inputsModel$tau[[response]] = tau[[response]]
inputsModel$timeDoseForOde[[response]] = sort(unique(c( timeDoseForOde, max(samplingTimesModel))))
}
###### model parameters
## check case model parameter in cond_init but no inequations
modelParameters = getModelParameters(object)
if (classModel == "ModelODEquations")
{
modelParametersInEquation = checkParameterInEquations( object )
}else
{
modelParametersInEquation = modelParameters
}
parameterInEquation = list()
for( parameter in modelParameters )
{
name = getNameModelParameter( parameter )
inputsModel[[name]] = getMu( parameter )
assign(name,inputsModel[[name]] )
}
inputsModel$model_parameters = modelParameters
# set initial conditions
# cond_init & case for doses bolus model
# dose in cond init and not in ode equations
cond_init_ode = list()
if (classModel == "ModelODEquations")
{
dC = getDerivatives(modelEquations)
doseRespPKInEquation = list()
doseRespPKInCondinit = list()
for ( nameResponsePK in responsePKNames )
{
nameDose = paste0("dose_",nameResponsePK)
doseRespPKInEquation[[nameDose]] = sapply( dC, function(x) all( grepl( nameDose, x) ) )
doseRespPKInCondinit[[nameDose]] = vapply( cond_init, function(x) all( grepl( nameDose, x) ),logical(1) )
for ( i in 1:length( doseRespPKInEquation[[nameDose]] ) )
{
for ( j in 1:length( doseRespPKInCondinit[[nameDose]] ) )
{
if ( ( doseRespPKInEquation[[nameDose]][i] == FALSE ) & ( doseRespPKInCondinit[[nameDose]][j] == TRUE ) )
{
assign( nameDose, dose[[nameResponsePK]] )
}
}
}
}
}
for ( nameCondiInit in names( cond_init ) )
{
cond_init_ode[[nameCondiInit]] = eval( cond_init[[nameCondiInit]] )
}
cond_init_ode = c( unlist( cond_init_ode ) )
# model equations
inputsModel$modelEquations = modelEquations
# variable of the model
inputsModel$variablesNames = names( cond_init )
# get parameters for ode solver
parametersOdeSolver = getParametersOdeSolver( object )
inputsModel$parametersOdeSolver = parametersOdeSolver
return(list(samplingTimesModel = samplingTimesModel,
cond_init_ode = cond_init_ode,
inputsModel = inputsModel,
modelParametersInEquation = modelParametersInEquation ) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Parameters used for computing the model gradient by finite-differences.
#'
#' @rdname parametersForComputingGradient
#' @param object A \code{StatisticalModel} object.
#' @return A list containing the parameters used for computing the gradient of the model.
setGeneric("parametersForComputingGradient",
function(object)
{
standardGeneric("parametersForComputingGradient")
}
)
setMethod(f = "parametersForComputingGradient",
signature="StatisticalModel",
definition= function(object)
{
minAbsPar = 0
parametersOdeSolver = getParametersOdeSolver( object )
.relStep = parametersOdeSolver$.relStep
listPars = getModelParameters( object )
valuePars = unlist(lapply(listPars, slot, name = "mu"))
.relStep = pmax(abs(valuePars), minAbsPar) * .relStep
pars <- as.numeric(valuePars)
npar <- length(pars)
incr <- pmax(abs(pars), minAbsPar) * .relStep
baseInd <- diag(npar)
frac <- c(1, incr, incr^2)
cols <- list(0, baseInd, -baseInd)
for ( i in seq_along(pars)[ -npar ] ) {
cols <- c( cols, list( baseInd[ , i ] + baseInd[ , -(1:i) ] ) )
frac <- c( frac, incr[ i ] * incr[ -(1:i) ] )
}
indMat <- do.call( "cbind", cols)
shifted <- pars + incr * indMat
indMat <- t(indMat)
Xcols <- list(1, indMat, indMat^2)
for ( i in seq_along(pars)[ - npar ] ) {
Xcols <- c( Xcols, list( indMat[ , i ] * indMat[ , -(1:i) ] ) )
}
return( list( Xcols = Xcols, shifted = shifted, frac = frac ) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Evaluation for the model, analytic, ode, infusion
#'
#' @rdname EvaluationModel
#' @param object A \code{StatisticalModel} object.
#' @param samplingTimesModel A vector containing the sampling times of the model.
#' @param cond_init_ode A vector containing the initial conditions for an ODE model.
#' @param inputsModel A list containing the parameters used for the model evaluation.
#' @param parametersGradient A list containing the parameters used for the evaluation of the gradient of the model
#' @param computeFIM A boolean giving TRUE id the FIM is compute, FALSE otherwise.
#' @return A list contaning a dataframe giving the results of the evaluation and a list giving the gradient of the model.
setGeneric("EvaluationModel",
function(object, samplingTimesModel,cond_init_ode, inputsModel, parametersGradient, computeFIM )
{
standardGeneric("EvaluationModel")
})
setMethod(f="EvaluationModel",
signature= "StatisticalModel",
definition=function(object, samplingTimesModel,cond_init_ode, inputsModel, parametersGradient, computeFIM )
{
modelEquations = getEquationsStatisticalModel( object )
classModel = class( modelEquations )
if ( classModel == "ModelInfusionODEquations" )
{
evaluationModel = EvaluateModelODEInfusion( modelEquations, samplingTimesModel,cond_init_ode,
inputsModel, parametersGradient, computeFIM )
} else if ( classModel == "ModelODEquations" ) {
evaluationModel = EvaluateModelODE( modelEquations, samplingTimesModel,cond_init_ode,
inputsModel, parametersGradient, computeFIM )
} else if ( classModel == "ModelEquations" ) {
evaluationModel = EvaluateModel( modelEquations, samplingTimesModel,inputsModel, computeFIM )
}else if ( classModel == "ModelInfusionEquations" ) {
evaluationModel = EvaluateModelInfusion( modelEquations, samplingTimesModel,inputsModel, computeFIM )
}
return( evaluationModel )
})
# -------------------------------------------------------------------------------------------------------------------
#' Evaluate an \code{StatisticalModel} object.
#'
#' @rdname Evaluate
#' @param object A \code{StatisticalModel} object.
#' @param administrations An \code{Administration} object.
#' @param sampling_times A \code{SamplingTimes} object.
#' @param cond_init A list for the initial conditions of the \code{StatisticalModel} object.
#' @param fim \code{FIM} object.
#' @return A \code{fim} object giving the Fisher Information Matrix of the \code{StatisticalModel} object.
setGeneric("Evaluate",
function(object, administrations, sampling_times, cond_init, fim )
{
standardGeneric("Evaluate")
})
setMethod(f="Evaluate",
signature= "StatisticalModel",
definition=function(object, administrations, sampling_times, cond_init, fim )
{
# initial parameters
MF_var = NA
V_total = NA
position = 1
errorVariances = list()
sigmaDerivatives = list()
modelParameters = getModelParameters( object )
modelEquations = getEquationsStatisticalModel( object )
modelResponses = getResponsesStatisticalModel( object )
# set parameters for evaluating the model
parametersGradient = parametersForComputingGradient( object )
parametersForEvaluateModel = setParametersForEvaluateModel( object,
administrations,
sampling_times, cond_init )
samplingTimesModel = parametersForEvaluateModel$samplingTimesModel
cond_init_ode = parametersForEvaluateModel$cond_init_ode
inputsModel = parametersForEvaluateModel$inputsModel
# for parameters not in equation but in the initial conditions
modelParametersInEquation = parametersForEvaluateModel$modelParametersInEquation
indexModelParametersInEquation = which( modelParameters %in% modelParametersInEquation )
indexModelParametersNotInEquation = which( modelParameters %in% modelParametersInEquation == FALSE )
# set model parameters for the evaluation
object@model_parameters <- inputsModel$model_parameters[indexModelParametersInEquation]
# -----------------------------------
# evaluate model
# -----------------------------------
# sampling time for plot
if ( object@computeFIM == FALSE ){
minSamplingTimesModel = 0
maxSamplingTimesModel = max( samplingTimesModel )
samplingTimesModel = sort( unique( c( samplingTimesModel,
linspace( minSamplingTimesModel, maxSamplingTimesModel, 500 ) ) ) )
evaluationModel = EvaluationModel( object,
samplingTimesModel,
cond_init_ode,
inputsModel,
parametersGradient, object@computeFIM )
predictedResponses = evaluationModel$evaluationResponse
responseGradient = evaluationModel$responseGradient
# gradient of the parameters in model equations
for ( response in modelResponses )
{
nameResponse = response@name
if( length( modelParametersInEquation ) == 1 )
{
responseGradient[[nameResponse]] = responseGradient[[nameResponse]][indexModelParametersInEquation]
}else{
responseGradient[[nameResponse]] = responseGradient[[nameResponse]][,indexModelParametersInEquation]
}
}
return( list( fim = fim,
predictedResponses = predictedResponses, sensitivityIndicesModel = responseGradient ) )
}
else
{
# sampling time of the designs
evaluationModel = EvaluationModel( object,
samplingTimesModel,
cond_init_ode,
inputsModel,
parametersGradient, object@computeFIM )
# predicted responses and gradients
predictedResponses = evaluationModel$evaluationResponse
responseGradient = evaluationModel$responseGradient
responseAllGradient = evaluationModel$responseAllGradient
nTotalSamplingTimes = length( unlist( inputsModel$sampleTimeResponse ) )
# gradient of the parameters in model equations
for ( response in modelResponses )
{
nameResponse = response@name
responseGradient[[nameResponse]] = responseAllGradient[,indexModelParametersInEquation]
}
responseAllGradient = responseAllGradient[,indexModelParametersInEquation]
modelParameters = getModelParameters( object )
modelParameters = modelParameters[indexModelParametersInEquation]
responseAllGradient = as.matrix( responseAllGradient )
## gradient for parameter not in eqs but in initial conditions
# parameter not in equations
modelParametersNotInEquation = inputsModel$model_parameters[indexModelParametersNotInEquation]
variablesNames = inputsModel$variablesNames
responsePKNames = inputsModel$responsePKNames
valueGradient = list()
# assign values parameters not in equation & PK doses
for( modelParameterNotInEquation in modelParametersNotInEquation )
{
name = getNameModelParameter( modelParameterNotInEquation )
assign( name, getMu( modelParameterNotInEquation ) )
}
for ( responsePKName in responsePKNames )
{
assign( paste0("dose_",responsePKName ), inputsModel$dose[[responsePKName]] )
}
# compute and eval the gradients for variables in initial conditions
for ( variableName in variablesNames )
{
equationCondInit = cond_init[[variableName]]
gradientExpression = sapply( names( modelParametersNotInEquation ),
function(x){ D( equationCondInit, x ) } )
valueGradient = unlist( sapply( gradientExpression, function(x){eval(x)} ) )
}
if ( length( valueGradient ) !=0 )
{
gradientParameterNotInEquations = rep(0, dim( responseAllGradient )[1])
gradientParameterNotInEquations[1] = valueGradient
responseAllGradient = cbind( responseAllGradient, gradientParameterNotInEquations )
}
# evaluate FIM: individual, population, Bayesian
object@model_parameters <- inputsModel$model_parameters
modelParameters = getModelParameters( object )
for ( response in modelResponses )
{
nameResponse = response@name
samplingTimesResponse = getSampleTime( sampling_times[[nameResponse]] )
predResp = predictedResponses[predictedResponses$time %in% samplingTimesResponse,nameResponse]
errorModelDerivatives = EvaluateErrorModelDerivatives( response, predResp )
errorVariances = append( errorVariances, bdiag( errorModelDerivatives$errorVariance ) )
errorVariances = bdiag( errorVariances )
for( errorModelParameter in errorModelDerivatives$sigmaDerivatives )
{
emptysigmaDerivatives = matrix( 0, ncol = nTotalSamplingTimes, nrow = nTotalSamplingTimes )
nTime = length( inputsModel$sampleTimeResponse[[ nameResponse ]] )
range <- position:( position + nTime - 1 )
emptysigmaDerivatives[ range, range ] <- errorModelParameter
# Add the MF_beta matrix for the parameters for each response
# add the MF_var in different blocks for each response
sigmaDerivatives = c( sigmaDerivatives, list( emptysigmaDerivatives ) )
}
position <- position + nTime
}
# get fixed parameters
dataFixedParameters = getFixedParameters( object )
indexMuNoFixedParameters = dataFixedParameters$indexMuNoFixedParameters
indexMuFixedParameters = dataFixedParameters$indexMuFixedParameters
indexOmegaNoFixedParameters = dataFixedParameters$indexOmegaNoFixedParameters
indexOmegaFixedParameters = dataFixedParameters$indexOmegaFixedParameters
indexZeroOmegaParameters = dataFixedParameters$indexZeroOmegaParameters
# names for the parameters
namesMuNoFixedParameters = names( modelParameters[indexMuNoFixedParameters] )
namesOmegaNoFixedParameters = names( modelParameters[indexOmegaNoFixedParameters] )
# names for the columns/rows for the FIM
muNames = paste0( "\u03bc_", namesMuNoFixedParameters )
omegaNames = paste0( "\u03c9\u00B2_", namesOmegaNoFixedParameters )
indexMuOmegaNoFixedParameters = sort( unique ( c( indexMuNoFixedParameters ) ) )
modelParameters = getModelParameters( object )
modelParametersNoFixedMu = modelParameters
# vector for sigma of error model
sigmaNames = c()
parameterNames = c()
# PopulationFim and IndividualFim
for (response in modelResponses )
{
if ( is( fim, "PopulationFim" ) )
{
MF = PopulationFIMEvaluateVariance(response,
modelEquations,
modelParametersNoFixedMu,
administrations,
sampling_times,
responseAllGradient,
errorVariances, sigmaDerivatives )
parameterNames <- c( muNames, omegaNames )
}
else
if ( is( fim, "IndividualFim" ) )
{
MF = IndividualFIMEvaluateVariance(response,
modelEquations,
modelParametersNoFixedMu,
administrations,
sampling_times,
responseAllGradient,
errorVariances, sigmaDerivatives )
parameterNames <- c( muNames )
}
### Add the MF_beta matrix for the parameters for each response
### add the MF_var in different blocks for each response
if ( is.matrix( MF_var ))
{
MF_var = bdiag( MF_var, MF$MF_var )
V_total = bdiag( V_total, MF$V )
}
else
{
### Create the first matrix and the first block
MF_var = bdiag( MF$MF_var )
V_total = bdiag( MF$V )
}
sigmaNames <- c( sigmaNames, getSigmaNames( response ) )
}
# compute MF_beta
if ( length( indexMuNoFixedParameters )!=0 )
{
responseAllGradient = responseAllGradient[, indexMuNoFixedParameters ]
}
MF_beta = t( responseAllGradient ) %*% solve( V_total ) %*% responseAllGradient
# Bayesian FIM
if ( is( fim, "BayesianFim" ) )
{
mu = c()
omega2 = c()
nParameter = 0
for (parameter in modelParameters )
{
omegaValue = getOmega( parameter )
omega2 = c( omega2, omegaValue^2 )
distribution = getDistribution( parameter )
if ( class( distribution )[1] == "LogNormalDistribution" )
{
muValue = getMu( parameter )
mu = c( mu, muValue )
}
else
if ( class( distribution )[1] == "NormalDistribution" )
{
mu = c( mu, 1 )
}
else
print("That distribution type is not supported for bayesian FIM")
}
# case: parameters fixed=FALSE & omega != 0
indexParameters = intersect ( indexOmegaNoFixedParameters, indexMuNoFixedParameters )
if (length(omega2)==1)
{
Omega = as.matrix(omega2[indexParameters])
meanParametersMatrix = as.matrix( mu[indexParameters] )
}else
{
Omega = diag( omega2[indexParameters] )
meanParametersMatrix = diag( mu[indexParameters] )
}
# subsetting beta matrix: no fixed mu and no fixed omega
colnames( MF_beta ) <- c( muNames )
rownames( MF_beta ) <- c( muNames )
muNames = paste0( "\u03bc_",names(modelParameters[indexParameters]))
MF_beta = as.matrix(MF_beta[muNames,muNames ])
# compute MF_beta
if( dim( MF_beta )[1] == 0 )
{
identity = c(1.0)
}else{
identity = diag( dim( MF_beta )[1] )
}
MF_beta = t( identity ) %*% MF_beta %*% identity + solve( meanParametersMatrix %*% Omega %*% meanParametersMatrix )
# set col and row names and Fisher matrix
colnames( MF_beta ) <- c( muNames )
rownames( MF_beta ) <- c( muNames )
fim <- `setMfisher<-`( fim, MF_beta )
# parameters indices used for plot
fim@parametersIndicesMuFIM = as.vector( indexParameters )
# set mu and omega : no fixed no omega=0
fim<-setOmega( fim, Omega )
fim<-setMu( fim, meanParametersMatrix )
return( list( fim = fim,
predictedResponses = predictedResponses,
sensitivityIndicesModel = responseGradient ) )
}
# case : PopulationFim & parameters with omega^2=0
# case parameter omega == 0 fixed and mu fixed
MF_i_total = bdiag( MF_beta, MF_var )
if ( is( fim, "PopulationFim" ) )
{
if ( length( indexOmegaFixedParameters ) !=0 ){
MF_var = MF_var[-c(indexOmegaFixedParameters),-c(indexOmegaFixedParameters)]
}
}
MF_i_total = bdiag( MF_beta, MF_var )
parameterNames = c( parameterNames, sigmaNames )
colnames( MF_i_total ) <- c( parameterNames )
rownames( MF_i_total ) <- c( parameterNames )
fim <- `setMfisher<-`( fim, MF_i_total )
# parameters indices used for plot
fim@parametersIndicesMuFIM = indexMuNoFixedParameters
fim@parametersIndicesOmegaFIM = indexOmegaNoFixedParameters
return( list( fim = fim,
predictedResponses = predictedResponses,
sensitivityIndicesModel = responseGradient ) )
}
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the fixed and non fixed model parameters.
#'
#' @name getFixedParameters
#' @param object A \code{StatisticalModel} object.
#' @return A list that contains the name and indices of the fixed parameters.
setGeneric("getFixedParameters",
function(object)
{
standardGeneric("getFixedParameters")
}
)
setMethod("getFixedParameters",
"StatisticalModel",
function(object)
{
modelParameters = getModelParameters( object )
# indices parameters mu and omega fixed or not
indexMuFixedParameters = which( lapply( modelParameters, slot, "fixedMu" ) == TRUE )
indexOmegaFixedParameters =which( lapply( modelParameters, slot, "fixedOmega" ) == TRUE )
indexMuNoFixedParameters = which( lapply( modelParameters, slot, "fixedMu" ) == FALSE )
indexOmegaNoFixedParameters = which( lapply( modelParameters, slot, "fixedOmega" ) == FALSE )
indexZeroOmegaParameters = which( lapply( modelParameters, slot, "omega" ) == 0 )
indexZeroMuParameters = which( lapply( modelParameters, slot, "mu" ) == 0 )
indexNoZeroOmegaParameters = which( lapply( modelParameters, slot, "omega" ) != 0 )
indexNoZeroMuParameters = which( lapply( modelParameters, slot, "mu" ) != 0 )
indexMuFixedParameters = sort( unique( c( indexMuFixedParameters, indexZeroMuParameters ) ) )
indexOmegaFixedParameters = sort( unique( c( indexOmegaFixedParameters, indexZeroOmegaParameters ) ) )
return ( output = list( indexMuFixedParameters = indexMuFixedParameters,
indexOmegaFixedParameters = indexOmegaFixedParameters,
indexMuNoFixedParameters = indexMuNoFixedParameters,
indexOmegaNoFixedParameters = indexOmegaNoFixedParameters,
indexZeroOmegaParameters = indexZeroOmegaParameters,
indexZeroMuParameters = indexZeroMuParameters,
indexNoZeroOmegaParameters = indexNoZeroOmegaParameters,
indexNoZeroMuParameters = indexNoZeroMuParameters) )
})
# -------------------------------------------------------------------------------------------------------------------
#' Define model equations
#'
#' @name defineModelEquations
#' @param object A \code{StatisticalModel} object.
#' @param equations An expression giving the equations of the model.
#' @return The \code{StatisticalModel} object with the equations.
setGeneric("defineModelEquations",
function(object, equations)
{
standardGeneric("defineModelEquations")
}
)
setMethod("defineModelEquations",
"StatisticalModel",
function(object, equations )
{
if ( class( equations ) %in% c( "PKModel","PKPDModel" ) ) {
object@modelEquations = getEquations( equations )
}else{
object@modelEquations = equations
}
return( object )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define a parameter of a statistical model.
#'
#' @name defineParameter
#' @param object A \code{StatisticalModel} object.
#' @param parameter An expression giving a parameter of the \code{StatisticalModel} object.
#' @return Return \code{StatisticalModel} object with new parameters.
setGeneric("defineParameter",
function(object, parameter)
{
standardGeneric("defineParameter")
}
)
setMethod("defineParameter",
"StatisticalModel",
function(object, parameter )
{
parameterName = getNameModelParameter( parameter )
object@model_parameters[[ parameterName ]] = parameter
object@model_parameters = object@model_parameters[ order( names( object@model_parameters ) ) ]
return( object )
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define the parameters of a statistical model.
#'
#' @name defineParameters
#' @param object A \code{StatisticalModel} object.
#' @param listOfParameters A list of string giving the parameters of the \code{StatisticalModel} object.
#' @return Return \code{StatisticalModel} object with new parameters.
setGeneric("defineParameters",
function(object, listOfParameters)
{
standardGeneric("defineParameters")
}
)
setMethod("defineParameters",
"StatisticalModel",
function(object, listOfParameters )
{
listOfParameters = listOfParameters[order(names(listOfParameters))]
for (i in 1:length( listOfParameters)) {
parameter = listOfParameters[[i]]
parameterName = getNameModelParameter( parameter )
object@model_parameters[[ parameterName ]] = parameter
}
return(object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Show the content of a \code{StatisticalModel} object.
#'
#' @rdname show
#' @param object \code{StatisticalModel} object.
#' @return Display the responses name of the model equations,
#' the ordinary derivatives of the model equations (for an ODE model), the parameters of the model.
setMethod("show",
signature = "StatisticalModel",
definition = function(object)
{
modelEquations = getEquationsStatisticalModel(object)
cat("*********************** Statistical model ******************************\n")
# ------------------------------------------------------------------------
# ModelEquations
# ------------------------------------------------------------------------
if( class(object@modelEquations ) %in% c("ModelEquations")){
cat( "\n **** Model equation:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
equation = getEquation( object@modelEquations, nameResponse )
cat( paste0(nameResponse ," = ", equation,"\n\n" ) )
}
cat( " **** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat( " Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat( "\n" )
}
}
# ------------------------------------------------------------------------
# ModelInfusionEquations
# ------------------------------------------------------------------------
else if( class(object@modelEquations ) %in% c("ModelInfusionEquations")){
cat( "\n **** Model equation:","\n\n" )
responses = getResponsesStatisticalModel(object)
namesResponses = names( responses )
for (nameResponse in namesResponses )
{
nameEquations = paste0( c( "DuringInfusion_","AfterInfusion_" ), nameResponse )
for (nameEquation in nameEquations )
{
equation = getEquation( object@modelEquations, nameEquation )
cat( paste0( paste0( nameEquation ) ," = ", equation,"\n\n" ) )
}
}
cat( " **** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat( " Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat( "\n" )
}
}
# ------------------------------------------------------------------------
# ModelODEquations
# ------------------------------------------------------------------------
else if(class(object@modelEquations) %in% c("ModelODEquations"))
{
derivatives <- getDerivatives( object@modelEquations )
cat("\n")
cat("**** Model responses:","\n\n" )
for( response in object@responses)
{
nameResponse = getNameResponse( response )
equation = getEquation( object@modelEquations, nameResponse )
cat( paste0( nameResponse ," = ", equation,"\n" ) )
}
cat("\n")
cat("**** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat("Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat("\n")
}
cat("**** Ordinary derivatives of model equation:\n\n")
for(i in 1:length(derivatives))
{
cat( names(derivatives)[i]," = ")
print(derivatives[[i]][[1]])
}
cat( "\n" )
}
# ------------------------------------------------------------------------
# ModelInfusionODEquations
# ------------------------------------------------------------------------
else if(class(modelEquations) %in% c("ModelInfusionODEquations"))
{
responses = getResponsesStatisticalModel(object)
namesResponses = names( responses )
cat("\n")
cat("**** Model responses:","\n\n" )
equationsDuringInfusion = modelEquations@duringInfusionResponsesEquations
equationsAfterInfusion = modelEquations@afterInfusionResponsesEquations
for (nameResponse in namesResponses )
{
nameEquationsDuringInfusion = paste0( c( "DuringInfusion_" ), nameResponse )
nameEquationsAfterInfusion = paste0( c( "AfterInfusion_" ), nameResponse )
cat( paste0( nameEquationsDuringInfusion ," = ", equationsDuringInfusion[[nameEquationsDuringInfusion]],"\n" ) )
cat( paste0( nameEquationsAfterInfusion ," = ", equationsAfterInfusion[[nameEquationsAfterInfusion]],"\n" ) )
}
cat("\n")
cat("**** Model error:","\n\n" )
for( response in object@responses )
{
nameResponse = getNameResponse( response )
cat("Model error",nameResponse, "\n" )
modelError = getModelError( response )
show( modelError )
cat("\n")
}
cat("\n")
cat("**** Ordinary derivatives of model equation:\n\n")
derivatives = getDerivatives( modelEquations )
namesDerivatives = names( derivatives )
for (namesDerivative in namesDerivatives )
{
cat( paste0( namesDerivative ," = ", derivatives[[namesDerivative]],"\n" ) )
}
}
# ------------------------------------------------------------------------
# parameters
# ------------------------------------------------------------------------
nameCol <- c()
muCol <- c()
omegaCol <- c()
distributionCol <- c()
for(parameter in object@model_parameters)
{
nameCol <- c( nameCol, getNameModelParameter(parameter))
muCol <- c(muCol, getMu(parameter))
omegaCol <- c(omegaCol, getOmega(parameter))
distributionCol <- c(distributionCol, class(getDistribution(parameter))[1])
}
parameterFrame <- matrix(c( as.character(muCol),
as.character(omegaCol),
distributionCol),ncol=3)
colnames(parameterFrame) <- c("\u03bc", "\u03c9", "Distribution" )
rownames(parameterFrame) <- nameCol
cat("\n")
cat("**** Model parameters: \n")
print(parameterFrame)
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the model parameters of a statistical model.
#'
#' @name getModelParameters
#' @param object \code{getModelParameters} object.
#' @return A list giving the model parameters.
setGeneric("getModelParameters",
function(object)
{
standardGeneric("getModelParameters")
}
)
setMethod("getModelParameters",
"StatisticalModel",
function(object)
{
return(object@model_parameters)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get the responses of a statistical model.
#'
#' @name getResponsesStatisticalModel
#' @param object A \code{getResponsesStatisticalModel} object.
#' @return A list giving the responses of a statistical model.
setGeneric("getResponsesStatisticalModel",
function(object)
{
standardGeneric("getResponsesStatisticalModel")
}
)
setMethod("getResponsesStatisticalModel",
"StatisticalModel",
function(object)
{
return(object@responses)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Get the equations of a statistical model.
#'
#' @name getEquationsStatisticalModel
#' @param object A \code{StatisticalModel} object.
#' @return A \code{ModelEquationsg} object giving the equations of the \code{StatisticalModel} object.
setGeneric(
"getEquationsStatisticalModel",
function(object) {
standardGeneric("getEquationsStatisticalModel")
})
setMethod("getEquationsStatisticalModel",
signature("StatisticalModel"),
function(object)
{
return( object@modelEquations)
})
# -------------------------------------------------------------------------------------------------------------------
#' Get the SE and RSE of the parameters.
#'
#' @name getErrorModelStandardErrors
#' @param object A \code{StatisticalModel} object.
#' @param fim A \code{Fim} object giving the Fisher Information Matrix.
#' @return A dataframe giving he SE and RSE of the parameters.
setGeneric("getErrorModelStandardErrors",
function(object, fim)
{
standardGeneric("getErrorModelStandardErrors")
}
)
setMethod("getErrorModelStandardErrors",
"StatisticalModel",
function(object, fim)
{
se <- getSE(fim)
sigmaCol <- c()
rseCol <- c()
responses = getResponsesStatisticalModel( object )
for(response in responses)
{
modelError <- getModelError(response)
sigmaCol <- c( sigmaCol, getSigmaValues(modelError))
}
if (length(se) >1 ){
if (length(se) != length(sigmaCol) ){
se <- se[-c(1:(length(se)-length(sigmaCol)))]
}else{
se <- (sort(se))
}}
rseCol <- c(rseCol, se/sigmaCol*100 )
parameterFrame <- data.frame( value = sigmaCol, se = se, "rse"= rseCol, row.names=NULL, check.names=FALSE)
return(parameterFrame)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define a variable in a statistical model.
#'
#' @name defineVariable
#' @param .Object A \code{StatisticalModel} object.
#' @param variable A character string giving the variable to be defined.
#' @return The \code{StatisticalModel} object with the new variable.
setGeneric("defineVariable",
function(.Object, variable )
{
standardGeneric("defineVariable")
}
)
setMethod(f="defineVariable",
signature= "StatisticalModel",
definition=function(.Object, variable)
{
variableName = variable@name
.Object@variables[[ variableName ]] = variable
return(.Object)
}
)
# -------------------------------------------------------------------------------------------------------------------
#' Define variables in a statistical model.
#'
#' @name defineVariables
#' @param .Object A \code{StatisticalModel} object.
#' @param listOfVariables A list of character string giving the variables to be defined.
#' @return The \code{StatisticalModel} object with the new variables.
setGeneric( "defineVariables",
function(.Object, listOfVariables ) #
{
standardGeneric( "defineVariables" )
}
)
setMethod(f = "defineVariables",
signature = "StatisticalModel",
definition = function(.Object, listOfVariables )
{
for (i in 1:length( listOfVariables)) {
modelVariable = ModelVariable( listOfVariables[[i]] )
variableName = modelVariable@name
.Object@variables[[ variableName ]] = modelVariable
}
return(.Object)
}
)
##########################################################################################################
# END Class "StatisticalModel"
##########################################################################################################
Try the PFIM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.