I am very interested on Model-Based Meta-Analysis and therefore, Non Linear Mixed Effect models. I have been using the library "nlmixr2" on R but I have been struggling to identify what the error is in the model.
The code developed is:
library(nlmixr2)
library(xpose)
library(ggplot2)
library(xpose.nlmixr2)
library(rxode2)
library(symengine)
### Auxiliar function to aid on the synthetic data creation
Fevcreator<-function() {
FEVw0 = c(rnorm(102,mean = 1.4,sd = 0.42),rnorm(105,mean = 1.4,sd = 0.44))
FEVw14 = c((rnorm(102,mean = 0.073 ,sd = 0.0816 )-rnorm(102,mean = 0.066 ,sd = 0.0930 )),
(rnorm(105,mean = -0.027 ,sd = 0.08 )-rnorm(105,mean = -0.041 ,sd = 0.0911 )))
FEVw4 = FEVw0+FEVw14
FEVw18 = c((rnorm(102,mean = 0.043 ,sd = 0.0728 )-rnorm(102,mean = 0.033 ,sd = 0.0847 )),
(rnorm(105,mean = -0.055 ,sd = 0.0714 )-rnorm(105,mean = -0.074 ,sd = 0.0830 )))
FEVw8 = FEVw0+FEVw18
FEVw12 = c(rnorm(102,mean = 1.5,sd = 0.43),rnorm(105,mean = 1.6,sd = 0.45))
FEVwx = rep(NA,207)
FEV = c(FEVw0,FEVwx,FEVwx,FEVwx,FEVw4,FEVwx,FEVwx,FEVwx,FEVw8,FEVwx,FEVwx,FEVwx,FEVw12)
FEV <- FEV*100
return(FEV)
}
##If there is a Synthetic Patient data file in the working directory, load it,
##If not, create it.
if (file.exists("ModellingData.RData")) {
load("ModellingData.RData")
}else{
ID<-1:207
ID<-rep(ID,13)
time<-c(rep("0.00",207),rep("1.00",207),rep("2.00",207),rep("3.00",207),rep("4.00",207),rep("5.00",207),
rep("6.00",207),rep("7.00",207),rep("8.00",207),rep("9.00",207),rep("10.00",207),
rep("11.00",207),rep("12.00",207))
drug = c(rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105),
rep("Roflumilast", 102), rep("Placebo", 105))
dose=c(rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105),
rep(3.5, 102), rep(0, 105))
age<-c(rnorm(102, mean = 66.9, sd = 6.95), rnorm(105, mean = 66.9, sd = 7.71))
age<-round(age)
age<-c(age,age,age,age,age,age,age,age,age,age,age,age,age)
body_mass<- c(rnorm(102, mean = 22.9, sd = 2.90), rnorm(105, mean = 22.4, sd = 2.68))
body_mass<-c(body_mass,body_mass,body_mass,body_mass,body_mass,body_mass,
body_mass,body_mass,body_mass,body_mass,body_mass,body_mass,body_mass)
smoking_status <- c(rep("Current Smoker", 36), rep("Former Smoker", 66),
rep("Current Smoker", 38), rep("Former Smoker", 67))
smoking_status<-rep(smoking_status,13)
COPDsev <- c(rep("Very Severe", 4), rep("Severe", 44), rep("Moderate", 53), rep("Mild", 1),
rep("Very Severe", 6), rep("Severe", 44), rep("Moderate", 53), rep("Mild", 2))
COPDsev <- rep(COPDsev,13)
DV<- Fevcreator()
DV2<- Fevcreator()
DV3<- Fevcreator()
DV4<- Fevcreator()
DV5<- Fevcreator()
dataf <- data.frame(ID,time,drug,age,body_mass,smoking_status,COPDsev,dose,DV)
dataf <- dataf[order(dataf$ID),]
datatoaproximate<-data.frame(ID,time,DV,DV2,DV3,DV4,DV5)
datatoaproximate <- datatoaproximate[order(datatoaproximate$ID),]
pc<-pca(datatoaproximate,method = "bpca",nPcs = 5,center = TRUE)
cObs<-completeObs(pc)
cObs<-data.frame(cObs)
cObs$mean <- rowMeans(cObs[,c('DV', 'DV2', 'DV3', 'DV4', 'DV5')], na.rm=TRUE)
dataf$DV<-cObs$mean
rm(ID,time,drug,age,body_mass,smoking_status,COPDsev,dose,datatoaproximate,pc,cObs,DV,DV2,DV3,DV4,DV5)
save(dataf, file = "ModellingData.RData")
}
###Model creation
jpkpdmodel <- function() {
ini({
tcl <- -log(0.8) # typical value of clearance
tv <- log(0.6)
thl<- 0.03571429 ##Half life in weeks
tEmax <- 104.9
eta.cl + eta.v~c(1,0.1, 1)
eta.bsln~10
eta.hl~0.1
add.err <- 0.1
})
model({
cl <- exp(tcl-eta.cl)
v <- exp(tv-eta.v)
Ke<-cl/v
hl <- exp(thl-eta.hl)
Ka <- (0.69314718056)/hl
baseline_FEV<-DV - exp(eta.bsln)
Emax <- tEmax*0.1
eff(0) <- baseline_FEV
conc <- Central/v
PD = Emax*conc
d/dt(Central) = Ka - cl / v * Central
d/dt(eff) = Ka*PD -Ke*eff
eff ~ add(add.err)
})
}
###Model execution
fit.s <- nlmixr2(jpkpdmodel, dataf,est="saem",control=list(print=0),table=list(cwres=TRUE, npde=TRUE))
When run, the execution reads:
ℹ parameter labels from comments will be replaced by 'label()'
→ loading into symengine environment...
→ pruning branches (`if`/`else`) of saem model...
✔ done
→ finding duplicate expressions in saem model...
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00
→ optimizing duplicate expressions in saem model...
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00
✔ done
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
Calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ loading into symengine environment...
→ pruning branches (`if`/`else`) of saem model...
✔ done
→ finding duplicate expressions in saem predOnly model 0...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ finding duplicate expressions in saem predOnly model 1...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in saem predOnly model 1...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ finding duplicate expressions in saem predOnly model 2...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in saem predOnly model 2...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
✔ done
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
→ Calculating residuals/tables
→ loading into symengine environment...
→ pruning branches (`if`/`else`) of full model...
✔ done
→ calculate jacobian
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ calculate sensitivities
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ calculate ∂(f)/∂(η)
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ calculate ∂(R²)/∂(η)
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ finding duplicate expressions in inner model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in inner model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ finding duplicate expressions in EBE model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in EBE model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ compiling inner model...
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
✔ done
→ finding duplicate expressions in FD model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in FD model...
[====|====|====|====|====|====|====|====|====|====] 0:00:00
→ compiling EBE model...
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
✔ done
→ compiling events FD model...
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
✔ done
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 70000 are needed
IDID=-2, larger nmax is needed
using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
using SDK: ‘MacOSX14.0.sdk’
Error : cannot figure out numeric time
Error: cannot figure out numeric time
I have tried changing the differential equations in the model, as well as changing the initial values.
The expectation is for the model to predict the effect of a bronchodilator on patients with COPD by measuring their FEV (Forced Expiratory Volume) at different stages in the trial.
I believe that it requires numeric time, here you have time as a character vector.
It was also answered on the (more active) github discussion group: