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.
#'
#' @param object A \code{StatisticalModel} object.
#' @return A list containing the parameters used for computing the gradient of the model.

function(object)
{
}
)

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,

predictedResponses = evaluationModel$evaluationResponse responseGradient = evaluationModel$responseGradient

# gradient of the parameters in model equations
for ( response in modelResponses )
{
nameResponse = response@name

if( length(  modelParametersInEquation ) == 1 )
{
}else{
}
}

return( list( fim = fim,
predictedResponses = predictedResponses, sensitivityIndicesModel = responseGradient ) )
}
else
{
# sampling time of the designs
evaluationModel = EvaluationModel( object,
samplingTimesModel,
cond_init_ode,
inputsModel,

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
}

modelParameters = getModelParameters( object )
modelParameters = modelParameters[indexModelParametersInEquation]

## 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 ) } )

}

if ( length( valueGradient ) !=0 )
{
}

# 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,
sampling_times,

parameterNames <- c( muNames, omegaNames )

}
else
if ( is( fim, "IndividualFim" ) )
{
MF = IndividualFIMEvaluateVariance(response,
modelEquations,
modelParametersNoFixedMu,
sampling_times,

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 )
{
}

# 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,
}
# 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,
}
})

# -------------------------------------------------------------------------------------------------------------------
#' 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.