REMixed : Regularisation & Estimation for Mixed effect model

Description

Suppose that we have a differential system of equations containing variables (S_{p})_{p\leq P} and R, that depends on some parameters. We define a non-linear mixed effects model from this system depending on individual parameters (\psi_i)_{i\leq N} and (\phi_i)_{i\leq N} that defines the parameters from the structural model as individuals. The first part (\psi_i)_{i\leq N} is supposed to derived from a generalized linear model for each parameter l\leq m :

h_l(\psi_{li}) = h_l(\psi_{l pop})+X_i\beta_l + \eta_{li}

with the covariates of individual i\leq N, (X_i)_{i\leq N}, random effects \eta_i=(\eta_{li})_{l\leq m}\overset{iid}{\sim}\mathcal{N}(0,\Omega), the population parameters \psi_{pop}=(\psi_{lpop})_{l\leq m} and \beta=(\beta_l)_{l\leq m_{re}} is the vector of covariates effects on parameters. The rest of the population parameters of the structural model, that hasn't random effetcs, are denoted by (\phi_i)_{i\leq N}, and are defined, for each parameters l\leq m_{no re} as

f_l(\phi_{li})=\phi_{l pop} + X_i \gamma_l

To simplify formula, we write the individual process over the time, resulting from the differential system for a set of parameters \phi_i = (\phi_{li})_{l\leq m_{no re}}, \psi_i = (\psi_{li})_{l\leq m_{re}} for individual i\leq N, as S_{p}(\cdot,\phi_i,\psi_i)=S_{pi}(\cdot), p\leq P and R(\cdot,\phi_i,\psi_i)=R_i(\cdot). We assume that individual trajectories (S_{pi})_{p\leq P,i\leq N} are observed through a direct observation model, up to a transformation g_p, p\leq P, at differents times (t_{pij})_{i\leq N,p\leq P,j\leq n_{ip}} :

Y_{pij}=g_p(S_{pi}(t_{pij}))+\epsilon_{pij}

with error \epsilon_p=(\epsilon_{pij})\overset{iid}{\sim}\mathcal N(0,\varsigma_p^2) for p\leq P. The individual trajectories (R_{i})_{i\leq N} are observed through $K$ latent processes, up to a transformation s_k, k\leq K, observed in (t_{kij})_{k\leq K,i\leq N,j\leq n_{kij}} :

Z_{kij}=\alpha_{k0}+\alpha_{k1} s_k(R_i(t_{kij}))+\varepsilon_{kij}

where \varepsilon_k\overset{iid}{\sim} \mathcal N(0,\sigma_k^2). The parameters of the model are then \theta=(\phi_{pop},\psi_{pop},B,\beta,\Omega,(\sigma^2_k)_{k\leq K},(\varsigma_p^2)_{p\leq P},(\alpha_{0k})_{k\leq K}) and \alpha=(\alpha_{1k})_{k\leq K}.

Author(s)

Maintainer: Auriane Gabaut auriane.gabaut@inria.fr

Authors:

• Ariane Bercu

• Mélanie Prague

• Cécile Proust-Lima

AIC for remix object

Description

Computes akaike information criterion from the output of remix as

AIC = -2\mathcal{LL}_{y}(\hat\theta,\hat\alpha)+k\times P

where P is the total number of parameters estimated and \mathcal{LL}_{y}(\hat\theta,\hat\alpha) the log-likelihood of the model.

Usage

## S3 method for class 'remix'
AIC(object, ..., k)

Arguments

Value

AIC.

References

Akaike, H. 1998. Information theory and an extension of the maximum likelihood principle, Selected papers of hirotugu akaike, 199-213. New York: Springer.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

AIC(res)

## End(Not run)

BIC for remix object

Description

Computes bayesian information criterion from the output of remix as

BIC = -2\mathcal{LL}_{y}(\hat\theta,\hat\alpha)+\log(N)P

where P is the total number of parameters estimated, N the number of subject and \mathcal{LL}_{y}(\hat\theta,\hat\alpha) the log-likelihood of the model.

Usage

## S3 method for class 'remix'
BIC(object, ...)

Arguments

Value

BIC.

References

Schwarz, G. 1978. Estimating the dimension of a model. The annals of statistics 6 (2): 461-464

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

BIC(res)

## End(Not run)

BICc

Description

Computes corrected bayesian information criterion as

BICc = -2\mathcal{LL}_{y}(\hat\theta,\hat\alpha)+P_R\log(N)+P_F\log(n_{tot})

where P_F is the total number of parameters linked to fixed effects, P_R to random effects, N the number of subject, n_tot the total number of observations and \mathcal{LL}_{y}(\hat\theta,\hat\alpha) the log-likelihood of the model.

Usage

BICc(object, ...)

Arguments

Value

BICc.

References

Delattre M, Lavielle M, Poursat M-A. A note on BIC in mixed-effects models. Elect J Stat. 2014; 8(1): 456-475.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

BICc(res)

## End(Not run)

Compute final estimation

Description

Computes a final saem and wald test if 'test' on the final model found by remix algorithm.

Usage

computeFinalTest(
  remix.output,
  dynFUN,
  y,
  ObsModel.transfo,
  final.project = NULL,
  pop.set = NULL,
  prune = NULL,
  n = NULL,
  parallel = TRUE,
  ncores = NULL,
  print = TRUE,
  digits = 3,
  trueValue = NULL,
  test = TRUE,
  p.max = 0.05
)

Arguments

Details

For population parameter estimation settings, see (<https://monolixsuite.slp-software.com/r-functions/2024R1/setpopulationparameterestimationsettings>).

Value

a remix object on which final SAEM and test, if test is TRUE, have been computed.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

res = cv.remix(project = project,
               dynFUN = dynFUN_demo,
               y = y,
               ObsModel.transfo = ObsModel.transfo,
               alpha = alpha,
               selfInit = TRUE,
               eps1=10**(-2),
               ncores=8,
               nlambda=8,
               eps2=1)

res_with_test = computeFinalTest(retrieveBest(res0,criterion=BICc),
                                 dynFUN_demo,
                                 y,
                                 ObsModel.transfo)

## End(Not run)

REMixed algorithm over a grid of \lambda

Description

Regularization and Estimation in MIXed effects model, over a regularization path.

Usage

cv.remix(
  project = NULL,
  final.project = NULL,
  dynFUN,
  y,
  ObsModel.transfo,
  alpha,
  lambda.grid = NULL,
  alambda = 0.001,
  nlambda = 50,
  lambda_max = NULL,
  eps1 = 10^(-2),
  eps2 = 10^(-1),
  selfInit = FALSE,
  pop.set1 = NULL,
  pop.set2 = NULL,
  prune = NULL,
  n = NULL,
  parallel = TRUE,
  ncores = NULL,
  print = TRUE,
  digits = 3,
  trueValue = NULL,
  unlinkBuildProject = TRUE,
  max.iter = +Inf
)

Arguments

Details

See REMixed-package for details on the model. For each \lambda\in\Lambda, the remix is launched. For population parameter estimation settings, see (<https://monolixsuite.slp-software.com/r-functions/2024R1/setpopulationparameterestimationsettings>).

Value

A list of outputs of the final project and of the iterative process over each value of lambda.grid:

info

Information about the parameters.

project

The project path if not unlinked.

lambda

The grid of \lambda.

BIC

Vector of BIC values for the model built over the grid of \lambda.

BICc

Vector of BICc values for the model built over the grid of \lambda.

LL

Vector of log-likelihoods for the model built over the grid of \lambda.

LL.pen

Vector of penalized log-likelihoods for the model built over the grid of \lambda.

res

List of all REMixed results for each \lambda (see remix).

outputs

List of all REMixed outputs for each \lambda (see remix).

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

res = cv.remix(project = project,
               dynFUN = dynFUN_demo,
               y = y,
               ObsModel.transfo = ObsModel.transfo,
               alpha = alpha,
               selfInit = TRUE,
               eps1=10**(-2),
               ncores=8,
               nlambda=8,
               eps2=1)

## End(Not run)

Dynamic functions demo

Description

Example of solver for remix and cv.remix algorithm. It is perfectly adapted for the Monolix demo project (see getMLXdir).

Usage

dynFUN_demo

Format

dynFUN_demo function of t, y, parms :

t

vector of timepoint.

y

initial condition, named vector of form c(AB=<...>,S=<...>).

parms

named vector of model parameter ; should contain phi_S,delta_AB,delta_S.

Details

Suppose you have antibodies secreting cells -S- that produces antibodies -AB- at rate \varphi_S. These two biological entities decay respectively at rate \delta_S and \delta_{AB}. The biological mechanism behind is :

\left\{\begin{matrix} \frac{d}{dt}S(t) &=& -\delta_S S(t) \\ \frac{d}{dt} AB(t) &=& \varphi_S S(t) - \delta_{AB} AB(t) \\ (S(0),AB(0)) &=& (S_0,AB_0) \end{matrix}\right.

References

Pasin C, Balelli I, Van Effelterre T, Bockstal V, Solforosi L, Prague M, Douoguih M, Thiébaut R, for the EBOVAC1 Consortium. 2019. Dynamics of the humoral immune response to a prime-boost Ebola vaccine: quantification and sources of variation. J Virol 93 : e00579-19. https://doi.org/10.1128/JVI.00579-19

See Also

model.pasin, getMLXdir.

Examples

t = seq(0,300,1)
y =c(AB=1000,S=5)
parms = c(phi_S = 611, delta_AB = 0.03, delta_S=0.01)

res <- dynFUN_demo(t,y,parms)

plot(res[,"time"],
     log10(res[,"AB"]),
     ylab="log10(AB(t))",
     xlab="time (days)",
     main="Antibody titer over the time",
     type="l")

plot(res[,"time"],
     res[,"S"],
     ylab="S(t)",
     xlab="time (days)",
     main="Antibody secreting cells quantity over time",
     type="l")

eBIC

Description

Computes extended bayesian information criterion as

eBIC = -2\mathcal{LL}_{y}(\hat\theta,\hat\alpha)+P\log(N)+2\gamma\log(\binom(k,K))

where P is the total number of parameters estimated, N the number of subject, \mathcal{LL}_{y}(\hat\theta,\hat\alpha) the log-likelihood of the model, K the number of submodel to explore (here the numbre of biomarkers tested) and k the numbre of biomarkers selected in the model.

Usage

eBIC(object, ...)

Arguments

Value

eBIC.

References

Chen, J. and Z. Chen. 2008. Extended Bayesian information criteria for model selection with large model spaces. Biometrika 95 (3): 759-771.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

eBIC(res)

## End(Not run)

extract remix results from cvRemix object

Description

Extracts a build from a cvRemix object.

Usage

extract(fit, n)

Arguments

Value

outputs from remix algorithm of rank 'n' computed by cv.remix.

See Also

cv.remix, remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

cv.outputs = cv.Remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            ncores=8,
            eps2=1)

res <- extract(cv.outputs,6)

plotConvergence(res)

trueValue = read.csv(paste0(dirname(project),"/demoSMLX/Simulation/populationParameters.txt"))


plotSAEM(res,paramToPlot = c("delta_S_pop","phi_S_pop","delta_AB_pop"),trueValue=trueValue)

## End(Not run)

Get monolix demo project path

Description

Get monolix demo project path

Usage

getMLXdir()

Value

path to the monolix demo from REMix package.

See Also

dynFUN_demo.

Examples

print(getMLXdir())

Adaptive Gauss-Hermite approximation of log-likelihood derivatives

Description

Computes Adaptive Gauss-Hermite approximation of the log-likelihood and its derivatives in NLMEM with latent observation processes, see REMixed-package for details on the model.

Usage

gh.LL(
  dynFUN,
  y,
  mu = NULL,
  Omega = NULL,
  theta = NULL,
  alpha1 = NULL,
  covariates = NULL,
  ParModel.transfo = NULL,
  ParModel.transfo.inv = NULL,
  Sobs = NULL,
  Robs = NULL,
  Serr = NULL,
  Rerr = NULL,
  ObsModel.transfo = NULL,
  data = NULL,
  n = NULL,
  prune = NULL,
  parallel = TRUE,
  ncores = NULL,
  onlyLL = FALSE,
  verbose = TRUE
)

Arguments

Details

Based on notation introduced REMixed-package. The log-likelihood of the model LL(\theta,\alpha_1) for a set of population parameters \theta and regulatization parameters \alpha_1 is estimated using Adaptative Gausse-Hermite quadrature, using conditional distribution estimation to locate the mass of the integrand. If the project has been initialized as a Monolix project, the user can use readMLX function to retrieve all the project information needed here.

Value

A list with the approximation by Gauss-Hermite quadrature of the likelihood L, the log-likelihood LL, the gradient of the log-likelihood dLL, and the Hessian of the log-likelihood ddLL at the point \theta, \alpha provided.

Examples

## Not run: 
project <- getMLXdir()


ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

data <- readMLX(project,ObsModel.transfo,alpha)

LL <- gh.LL(dynFUN = dynFUN_demo,
            y = c(S=5,AB=1000),
            ObsModel.transfo=ObsModel.transfo,
            data = data)

print(LL)

## End(Not run)

Generate individual parameters

Description

Generate the individual parameters of indivual whose covariates are covariates and random effects eta_i.

Usage

indParm(theta, covariates, eta_i, transfo, transfo.inv)

Arguments

Details

The models used for the parameters are :

h_l(\psi_{li}) = h_l(\psi_{lpop})+X_i\beta_l + \eta_{li}

with h_l the transformation, \beta_l the vector of covariates effect and with \eta_i the random effects associated \psi_l parameter ;

g_k(\phi_{ki}) = g_k(\phi_{kpop})+X_i \gamma_l

with g_k the transformation and \gamma_k the vector of covariates effect associated \phi_k parameter.

Value

a list with phi_i and psi_i parameters.

See Also

model.clairon, model.pasin.

Examples

phi_pop = c(delta_S = 0.231, delta_L = 0.000316)
psi_pop = c(delta_Ab = 0.025,phi_S = 3057, phi_L = 16.6)
gamma = NULL
covariates = data.frame(cAGE = runif(1,-15,15), G1 = rnorm(1), G2 = rnorm(1))
beta = list(delta_Ab=c(0,1.2,0),phi_S = c(0.93,0,0),phi_L=c(0,0,0.8))

theta=list(phi_pop = phi_pop,psi_pop = psi_pop,gamma = gamma, beta = beta)
eta_i = c(delta_Ab = rnorm(1,0,0.3),phi_S=rnorm(1,0,0.92),phi_L=rnorm(1,0,0.85))
transfo = list(delta_Ab=log,phi_S=log,phi_L=log)
transfo.inv = list(delta_Ab = exp,phi_S=exp,phi_L=exp)

indParm(theta,covariates,eta_i,transfo,transfo.inv)

Initialization strategy

Description

Selecting an initialization point by grouping biomarkers of project and running the SAEM. Initial condition is then selected using maximum log-likelihood.

Usage

initStrat(
  project,
  alpha,
  ObsModel.transfo,
  Nb_genes = 2,
  trueValue = NULL,
  pop.set = NULL,
  useSettingsInAPI = FALSE,
  conditionalDistributionSampling = FALSE,
  print = TRUE,
  digits = 2,
  unlinkTemporaryProject = TRUE,
  seed = NULL
)

Arguments

Details

For population parameter estimation settings, see (<https://monolixsuite.slp-software.com/r-functions/2024R1/setpopulationparameterestimationsettings>).

Value

a list of outputs for every group of genes tested with composition of the group, final parameter estimates, final scores estimates (OFV, AIC, BIC, BICc), temporary project directory. The final selected set is initialize in the project.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

initStrat(project,alpha,ObsModel.transfo,seed=1710)

## End(Not run)

Model from Clairon and al.,2023

Description

Generates the dynamics of antibodies secreting cells -S- that produces antibodies -AB- over time, with two injection of vaccine at time t_0=0 and t_{inj}, using Clairon and al., 2023, model.

Usage

model.clairon(t, y, parms, tinj = 21)

Arguments

Details

Model is defined as

\displaystyle\left\{\begin{matrix} \frac{d}{dt} S(t) &=& f_{\overline M_k} e^{-\delta_V(t-t_k)}-\delta_S S(t) \\ \frac{d}{dt} Ab(t) &=& \theta S(t) - \delta_{Ab} Ab(t)\end{matrix}\right.

on each interval I_1=[0;t_{inj}[ and I_2=[t_{inj};+\infty[. For each interval I_k, we have t_k corresponding to the last injection date of vaccine, such that t_1=0 and t_2=t_{inj}. By definition, f_{\overline M_1}=1 (Clairon and al., 2023).

Value

Matrix of time and observation of antibody secreting cells S and antibody titer Ab.

References

Quentin Clairon, Melanie Prague, Delphine Planas, Timothee Bruel, Laurent Hocqueloux, et al. Modeling the evolution of the neutralizing antibody response against SARS-CoV-2 variants after several administrations of Bnt162b2. 2023. hal-03946556

See Also

indParm

Examples

y = c(S=1,Ab=0)

parms = c(fM2 = 4.5,
          theta = 18.7,
          delta_S = 0.01,
          delta_Ab = 0.23,
          delta_V = 2.7)

t = seq(0,35,1)

res <- model.clairon(t,y,parms)

plot(res)

Model from Pasin and al.,2019

Description

Generate trajectory of the Humoral Immune Response to a Prime-Boost Ebola Vaccine.

Usage

model.pasin(t, y, parms)

Arguments

Details

The model correspond to the dynamics of the humoral response, from 7 days after the boost immunization with antibodies secreting cells -S and L, characterized by their half lives- that produces antibodies -AB- at rate \theta_S and \theta_L. All these biological entities decay at rate repectively \delta_S, \delta_L and \delta_{Ab}. Model is then defined as

\left\{\begin{matrix}\frac{d}{dt} Ab(t) &=& \theta_S S(t) + \theta_L L(t) - \delta_{Ab} Ab(t) \\ \frac{d}{dt}S(t) &=& -\delta_S S(t) \\ \frac{d}{dt} L(t) &=& -\delta_L L(t)\end{matrix}\right.

Value

Matrix of time and observation of antibody titer Ab, and ASCs S and L.

References

Pasin C, Balelli I, Van Effelterre T, Bockstal V, Solforosi L, Prague M, Douoguih M, Thiébaut R, for the EBOVAC1 Consortium. 2019. Dynamics of the humoral immune response to a prime-boost Ebola vaccine: quantification and sources of variation. J Virol 93: e00579-19. https://doi.org/10.1128/JVI.00579-19

See Also

indParm, model.clairon.

Examples

y = c(Ab=0,S=5,L=5)
parms = c(theta_S = 611,
          theta_L = 3.5,
          delta_Ab = 0.025,
          delta_S = 0.231,
          delta_L = 0.000152)

t = seq(0,100,5)
res <- model.pasin(t,y,parms)
plot(res)

Generate trajectory of PK model

Description

The administration is via a bolus. The PK model has one compartment (volume V) and a linear elimination (clearance Cl). The parameter ka is defined as ka=\frac{Cl}{V}.

Usage

model.pk(t, y, parms)

Arguments

Value

Matrix of time and observation of Concentration C.

See Also

indParm.

Examples

res <- model.pk(seq(0,30,1),c(C=100),parms=c(ka=1))

plot(res)

Plot of cv.remix object

Description

Calibration plot for cvRemix object.

Usage

## S3 method for class 'cvRemix'
plot(x, criterion = BICc, trueValue = NULL, ...)

Arguments

Value

A plot.

See Also

cv.remix

Calibration plot

Description

Calibration plot

Usage

plotCalibration(
  fit,
  legend.position = "none",
  trueValue = NULL,
  criterion = BICc,
  dismin = TRUE
)

Arguments

Value

Calibration plot, over the lambda.grid.

See Also

remix, cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

res = cv.remix(project = project,
               dynFUN = dynFUN_demo,
               y = y,
               ObsModel.transfo = ObsModel.transfo,
               alpha = alpha,
               selfInit = TRUE,
               eps1=10**(-2),
               ncores=8,
               nlambda=8,
               eps2=1)

plotCalibration(res)

plotIC(res)

## End(Not run)

Log-likelihood convergence

Description

Log-likelihood convergence

Usage

plotConvergence(fit, ...)

Arguments

Value

Log-Likelihood values throughout the algorithm iteration.

See Also

remix, cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

plotConvergence(res)

trueValue = read.csv(paste0(dirname(project),"/demoSMLX/Simulation/populationParameters.txt"))#'

plotSAEM(res,paramToPlot = c("delta_S_pop","phi_S_pop","delta_AB_pop"),trueValue=trueValue)

## End(Not run)

IC plot

Description

IC plot

Usage

plotIC(fit, criterion = BICc, dismin = TRUE)

Arguments

Value

IC trhoughout the lambda.grid.

See Also

remix, cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

res = cv.remix(project = project,
               dynFUN = dynFUN_demo,
               y = y,
               ObsModel.transfo = ObsModel.transfo,
               alpha = alpha,
               selfInit = TRUE,
               eps1=10**(-2),
               ncores=8,
               nlambda=8,
               eps2=1)

plotCalibration(res)

plotIC(res)

## End(Not run)

Plot initialization

Description

Plot initialization

Usage

plotInit(init, alpha = NULL, trueValue = NULL)

Arguments

Value

log-likelihood value for all groups of genes tested.

See Also

initStrat.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

init <- initStrat(project,alpha,ObsModel.transfo,seed=1710)

plotInit(init)

## End(Not run)

Display the value of parameters at each iteration

Description

Display the value of parameters at each iteration

Usage

plotSAEM(fit, paramToPlot = "all", trueValue = NULL)

Arguments

Value

For each parameters, the values at the end of each iteration of remix algorithm is drawn. Moreover, the SAEM steps of each iteration are displayed.

See Also

remix, cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 1440

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

plotConvergence(res)

trueValue = read.csv(paste0(dirname(project),"/demoSMLX/Simulation/populationParameters.txt"))

plotSAEM(res,paramToPlot = c("delta_S_pop","phi_S_pop","delta_AB_pop"),trueValue=trueValue)

## End(Not run)

Extract Data for REMixed Algorithm from a Monolix Project

Description

This function retrieves all necessary information from a Monolix project file to format the input for the REMixed package. It gathers all relevant data required for the REMix algorithm.

Usage

readMLX(project = NULL, ObsModel.transfo, alpha)

Arguments

Details

To simplify its use, functions remix, cv.remix, gh.LL can be used with arguments data rather than all necessary informations "theta", "alpha1", "covariates", "ParModel.transfo", "ParModel.transfo.inv", "Sobs", "Robs", "Serr", "Rerr", "ObsModel.transfo" that could be extract from a monolix project. If the SAEM task of the project hasn't been launched, it's the initial condition and not the estimated parameters that are returned. If the conditional distribution estimation task has been launched, parameters "mu" and "Omega" are returned too.

Value

A list containing parameters, transformations, and observations from the Monolix project in the format needed for the REMixed algorithm :

• mu list of individuals random effects estimation (vector of r.e. need to be named by the parameter names), use to locate the density mass (if conditional distribution estimation through Monolix has been launched);

• Omega list of individuals estimated standard deviation diagonal matrix (matrix need to have rows and columns named by the parameter names), use to locate the density mass (if conditional distribution estimation through Monolix has been launched);

• theta list of model parameters containing i

  • phi_pop : named vector with the population parameters with no r.e. (\phi_{l\ pop})_{l\leq L} (NULL if none) ;

  • psi_pop : named vector with the population parameters with r.e. (\psi_{l\ pop})_{l\leq m} ;

  • gamma : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with no r.e. ;

  • beta : named list (for each parameters) of named vector (for each covariates) of covariate effects from parameters with r.e..

  • alpha0 : named vector of (\alpha_{0k})_{k\leq K} parameters (names are identifier of the observation model, such as in a Monolix project);

  • omega : named vector of estimated r.e. standard deviation;

• alpha1 named vector of regulatization parameters (\alpha_{1k})_{k\leq K}, with identifier of observation model as names;

• covariates matrix of individual covariates (size N x n). Individuals must be sorted in the same order than in mu and Omega;

• ParModel.transfo named list of transformation functions (h_l)_{l\leq m} and (s_k)_{k\leq K} for the individual parameter model (names must be consistent with phi_pop and psi_pop, missing entries are set by default to the identity function ;

• ParModel.transfo.inv named list of inverse transformation functions for the individual parameter model (names must be consistent with phi_pop and psi_pop ;

• Sobs ist of individuals trajectories for the direct observation models (Y_{pi})_{p \leq P,i\leq N}. Each element i\leq N of the list, is a list of p\leq P data.frame with time (t_{pij})_{j\leq n_{ip}} and observations (Y_{pij})_{j\leq n_{ip}}. Each data.frame is named with the observation model identifiers ;

• Robs list of individuals trajectories for the latent observation models (Z_{ki})_{k \leq K,i\leq N}. Each element i\leq N of the list, is a list of k\leq K data.frame with time (t_{kij})_{j\leq n_{ik}} and observations (Z_{kij})_{j\leq n_{ik}}. Each data.frame is named with the observation model identifiers ;

• Serr named vector of the estimated error mocel constants (\varsigma_p)_{p\leq P} with observation model identifiers as names ;

• Rerr named vector of the estimated error mocel constants (\sigma_k)_{k\leq K} with observation model identifiers as names ;

• ObsModel.transfo same as inputObsModel.transfo list.

See Also

remix, cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

res <- readMLX(project,ObsModel.transfo,alpha)

## End(Not run)

REMixed algorithm

Description

Regularization and Estimation in Mixed effects model.

Usage

remix(
  project = NULL,
  final.project = NULL,
  dynFUN,
  y,
  ObsModel.transfo,
  alpha,
  lambda,
  eps1 = 10^(-2),
  eps2 = 10^(-1),
  selfInit = FALSE,
  pop.set1 = NULL,
  pop.set2 = NULL,
  pop.set3 = NULL,
  prune = NULL,
  n = NULL,
  parallel = TRUE,
  ncores = NULL,
  print = TRUE,
  verbose = FALSE,
  digits = 3,
  trueValue = NULL,
  finalSAEM = FALSE,
  test = TRUE,
  max.iter = +Inf,
  p.max = 0.05
)

Arguments

Details

See REMixed-package for details on the model. For population parameter estimation settings, see (<https://monolixsuite.slp-software.com/r-functions/2024R1/setpopulationparameterestimationsettings>).

Value

a list of outputs of final project and through the iteration :

info

informations about the parameters (project path, regulatization and population parameter names, alpha names, value of lambda used, if final SAEM and test has been computed, parameters p.max and N) ;

finalRes

containing loglikelihood LL and penalized loglikelihood LL.pen values, final population parameters param and final regularization parameters alpha values, number of iterations iter and time needed , if computed, the estimated standard errors standardError and if test computed, the final results before test saemBeforeTest ;

iterOutputs

the list of all remix outputs, i.e. parameters, lieklihood, SAEM estimates and convergence criterion value over the iteration.

See Also

cv.remix.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)
lambda = 382.22

res = remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            eps2=1,
            lambda=lambda)

summary(res)

trueValue = read.csv(paste0(dirname(project),"/demoSMLX/Simulation/populationParameters.txt"))

plotSAEM(res,paramToPlot = c("delta_S_pop","phi_S_pop","delta_AB_pop"),trueValue=trueValue)

## End(Not run)

REMixed results

Description

Extracts the build minimizing an information criterion over a grid of lambda.

Usage

retrieveBest(fit, criterion = BICc)

Arguments

Value

outputs from remix algorithm achieving the best IC among those computed by cv.remix.

See Also

cv.remix, remix, BIC.remix, eBIC, AIC.remix, BICc.

Examples

## Not run: 
project <- getMLXdir()

ObsModel.transfo = list(S=list(AB=log10),
                        linkS="yAB",
                        R=rep(list(S=function(x){x}),5),
                        linkR = paste0("yG",1:5))

alpha=list(alpha0=NULL,
           alpha1=setNames(paste0("alpha_1",1:5),paste0("yG",1:5)))

y = c(S=5,AB=1000)

cv.outputs = cv.Remix(project = project,
            dynFUN = dynFUN_demo,
            y = y,
            ObsModel.transfo = ObsModel.transfo,
            alpha = alpha,
            selfInit = TRUE,
            eps1=10**(-2),
            ncores=8,
            eps2=1)

res <- retrieveBest(cv.outputs)

plotConvergence(res)

trueValue = read.csv(paste0(dirname(project),"/demoSMLX/Simulation/populationParameters.txt"))#'

plotSAEM(res,paramToPlot = c("delta_S_pop","phi_S_pop","delta_AB_pop"),trueValue=trueValue)

## End(Not run)