我对基于模型的荟萃分析非常感兴趣，因此对非线性混合效应模型非常感兴趣。 我一直在 R 上使用库“nlmixr2”，但我一直在努力确定模型中的错误是什么。
开发的代码是：
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))
运行时，执行内容为：
ℹ 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
我尝试过更改模型中的微分方程，以及更改初始值。
模型期望通过在试验的不同阶段测量患者的 FEV（用力呼气量）来预测支气管扩张剂对 COPD 患者的效果。
我相信它需要数字时间，这里你有时间作为字符向量。
在（更活跃的）github 讨论组上也得到了回答：