Closed andreasmeid closed 1 year ago
mattfidler
commented
2 years ago Hi @andreasmeid
I am unclear what is wrong here. Could you provide some sample dataset so I can try to reproduce what you are experiencing?
In the new nlmixr2 mu referencing is no longer needed, so this may not be taken care of on this repository.
andreasmeid
commented
2 years ago To be honest, I only saw it now, the random effects were assigned with the wrong symbol ("<-" instead of "~"). With the correct notation the model runs. Maybe it could be pointed out more explicitly, obviously I didn't recognize it in the tutorials.
The data is available for download as csv-file on the ddmore repository (see link above). There are also the original NONMEM files. In the NONMEM code there is, among other things, also the NONMEM-specific assignment IF(NEWIND.NE.2) OCB=CB. This I have now picked up in advance in the R data management, hopefully correct. If something like this doesn't exist yet, could some kind of dictionary between NONMEM code and nlmixr code be thought of for the future?
But for now the issue seems to be solved, thanks a lot
mattfidler
commented
2 years ago In theory the newind flag exists
mattfidler
commented
2 years ago The meaning may be slightly different. Details about functions and variables like newind are here:
https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-syntax.html
andreasmeid
commented
2 years ago Thank you very much for the link to the ODE syntx. This was indeed very helpful. When I tackle the second part of the tumor growth model (http://repository.ddmore.foundation/model/DDMODEL00000218#Files), I get an error message that is hard to interpret: "Error in RxODE::rxState(rxode) : Evaluation error: syntax errors (see above)." I have no idea which syntax error might be meant here.
OS <- read.table(here("data/Simulated_OS.csv"), header=T, sep=",")
OS2 <- OS %>% group_by(ID) %>% mutate(group_IND = case_when(any(AUC1 > 0) ~ "CB + G", # and group
all(AUC1 == 0) ~ "CB")) %>%
rename(NWLS=NWLSCOV) %>% select(-CCOM) %>%
mutate(dvid = case_when(EVID == 0 ~ "Y", # indicator for dependent var (dv)
EVID != 0 ~ "other")) %>% rename(dv=DV, id=ID) %>%
mutate(dvid=as.factor(dvid)) %>% # adjust variable types
mutate(AUC0 = case_when(TIME == 0 ~ lead(AUC0), # assign first value to baseline
TIME != 0 ~ AUC0),
AUC1 = case_when(TIME == 0 ~ lead(AUC1),
TIME != 0 ~ AUC1),
NWLS = case_when(TIME == 0 ~ lead(NWLS),
TIME != 0 ~ NWLS))
myModel_OS <- function(){ # TTE model for OS. Weibull baseline hazard.
ini({
tLAM <- log(0.001) # scale parameter
tSHP <- log(2) # shape parameter
tBSLD0 <- log(0.1) # parameter for SLD (sum of longest diameter) at enrolment
tBTSR <- log(0.1) # parameter for TSR(t) (tumour size ratio)
tBNWLS <- log(0.1) # parameter for NewLesion(t)
tBECOG <- log(0.1) # parameter for ECOG at enrolment
eta.LAM ~ 0.01 # inter-individual variability in scale parameter
add.err <- 0
})
model({
LAM = exp(tLAM + eta.LAM)
SHP = exp(tSHP)
BSLD0 = exp(tBSLD0)
BTSR = exp(tBTSR)
BNWLS = exp(tBNWLS)
BECOG = exp(tBECOG)
#Time constant covariates
TVSLD0 = 70 # average SLD at enrolment
NSLD0 = SLD0/TVSLD0 # normalised SLD at enrolment
IECOG=ECOG # ECOG at enrolment
# Model for dSLD(t)
A(0) = IBASE*1000 # SLD baseline
MMBAS = IBASE*1000 # Parameter of the SLD(t) model were estimated previously (IPP approach)
d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + + KD1/100 * AUC1) * A
TUM = A
TSR = (TUM-MMBAS)/MMBAS
if(time==0){ #TSR(t) for t<= week 12 and TSR(week12) for t> week 12
WTS = 0
TM12 = 0
}
if(time<=84){
WTS = TSR
TM12 = WTS
}
else {
WTS = WTS
}
# survival model
DEL = 1E-6
d/dt(A2) = LAM*SHP*(LAM*(T+DEL))**(SHP-1) *EXP(BSLD0*NSLD0+BTSR*WTS+BNWLS*NWLS+BECOG*IECOG)
# Death hazard
CHZ = A2
SUR=EXP(-CHZ)
DELX = 1E-6
xTUM = A
xTSR = (xTUM-MMBAS)/MMBAS
if(time==0){
XWTS=0
XTM12=0
}
if(time<=84){
XWTS = XTSR
XTM12 = XWTS
}
else {
XWTS=XWTS
}
HAZN = LAM*SHP*(LAM*(TIME+DELx))**(SHP-1)*EXP(BSLD0*NSLD0+BTSR*XWTS+BNWLS*NWLS+BECOG*IECOG)
# prediction
if(FLAG == 9 & EVID == 0 & OSCENS == 1){
IPRED=SUR # probability of survival (censored event)
Y = IPRED # Y is probability for TTE data
Y ~ add(add.err)
}
if(FLAG == 9 & EVID == 0 & OSCENS == 0){
IPRED=SUR*HAZN # probability of event (death) at time=TIME
Y = IPRED # Y is probability for TTE data
Y ~ add(add.err)
}
})
}
fit.OSModel <- nlmixr(myModel_OS,OS2,est="saem",control=saemControl(n.burn=50,n.em=100,print=50),
RxODE.syntax.assign.state = TRUE);
What could be the reason for this or where could one specifically look in the code?
mattfidler
commented
2 years ago The current version of nlmixr stores the errors, but not the messages; put a rxClean() at the beginning should see the errors:
myModel_OS <- function(){ # TTE model for OS. Weibull baseline hazard.
ini({
tLAM <- log(0.001) # scale parameter
tSHP <- log(2) # shape parameter
tBSLD0 <- log(0.1) # parameter for SLD (sum of longest diameter) at enrolment
tBTSR <- log(0.1) # parameter for TSR(t) (tumour size ratio)
tBNWLS <- log(0.1) # parameter for NewLesion(t)
tBECOG <- log(0.1) # parameter for ECOG at enrolment
eta.LAM ~ 0.01 # inter-individual variability in scale parameter
add.err <- 0
})
model({
LAM = exp(tLAM + eta.LAM)
SHP = exp(tSHP)
BSLD0 = exp(tBSLD0)
BTSR = exp(tBTSR)
BNWLS = exp(tBNWLS)
BECOG = exp(tBECOG)
#Time constant covariates
TVSLD0 = 70 # average SLD at enrolment
NSLD0 = SLD0/TVSLD0 # normalised SLD at enrolment
IECOG=ECOG # ECOG at enrolment
# Model for dSLD(t)
A(0) = IBASE*1000 # SLD baseline
MMBAS = IBASE*1000 # Parameter of the SLD(t) model were estimated previously (IPP approach)
d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + + KD1/100 * AUC1) * A
TUM = A
TSR = (TUM-MMBAS)/MMBAS
if(time==0){ #TSR(t) for t<= week 12 and TSR(week12) for t> week 12
WTS = 0
TM12 = 0
}
if(time<=84){
WTS = TSR
TM12 = WTS
}
else {
WTS = WTS
}
# survival model
DEL = 1E-6
d/dt(A2) = LAM*SHP*(LAM*(T+DEL))**(SHP-1) *EXP(BSLD0*NSLD0+BTSR*WTS+BNWLS*NWLS+BECOG*IECOG)
# Death hazard
CHZ = A2
SUR=EXP(-CHZ)
DELX = 1E-6
xTUM = A
xTSR = (xTUM-MMBAS)/MMBAS
if(time==0){
XWTS=0
XTM12=0
}
if(time<=84){
XWTS = XTSR
XTM12 = XWTS
}
else {
XWTS=XWTS
}
HAZN = LAM*SHP*(LAM*(TIME+DELx))**(SHP-1)*EXP(BSLD0*NSLD0+BTSR*XWTS+BNWLS*NWLS+BECOG*IECOG)
# prediction
if(FLAG == 9 & EVID == 0 & OSCENS == 1){
IPRED=SUR # probability of survival (censored event)
Y = IPRED # Y is probability for TTE data
Y ~ add(add.err)
}
if(FLAG == 9 & EVID == 0 & OSCENS == 0){
IPRED=SUR*HAZN # probability of event (death) at time=TIME
Y = IPRED # Y is probability for TTE data
Y ~ add(add.err)
}
})
}
library(nlmixr)
rxClean()
nlmixr(myModel_OS)
#> i parameter labels from comments will be replaced by 'label()'
#> RxODE model syntax error:
#> ================================================================================
#> :001: TVSLD0 = 70
#> :002: NSLD0 = SLD0/TVSLD0
#> :003: IECOG = ECOG
#> :004: A(0) = IBASE * 1000
#> :005: MMBAS = IBASE * 1000
#> :006: d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + +KD1/100 * AUC1) * A
#> :007: TUM = A
#> :008: TSR = (TUM - MMBAS)/MMBAS
#> :009: if (time == 0) {
#> :010: WTS = 0
#> :011: TM12 = 0
#> :012: }
#> :013: if (time <= 84) {
#> :014: WTS = TSR
#> :015: TM12 = WTS
#> :016: }
#> :017: else {
#> :018: WTS = WTS
#> :019: }
#> :020: DEL = 1e-06
#> :021: function 'EXP' is not supported in RxODE:
#> d/dt(A2) = LAM * SHP * (LAM * (T + DEL))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * WTS + BNWLS * NWLS + BECOG * IECOG)
#> ^~~
#> :021: d/dt(A2) = LAM * SHP * (LAM * (T + DEL))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * WTS + BNWLS * NWLS + BECOG * IECOG)
#> :022: CHZ = A2
#> :023: function 'EXP' is not supported in RxODE:
#> SUR = EXP(-CHZ)
#> ^~~
#> :023: SUR = EXP(-CHZ)
#> :024: DELX = 1e-06
#> :025: xTUM = A
#> :026: xTSR = (xTUM - MMBAS)/MMBAS
#> :027: if (time == 0) {
#> :028: XWTS = 0
#> :029: XTM12 = 0
#> :030: }
#> :031: if (time <= 84) {
#> :032: XWTS = XTSR
#> :033: XTM12 = XWTS
#> :034: }
#> :035: else {
#> :036: XWTS = XWTS
#> :037: }
#> :038: function 'EXP' is not supported in RxODE:
#> HAZN = LAM * SHP * (LAM * (TIME + DELx))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * XWTS + BNWLS * NWLS + BECOG * IECOG)
#> ^~~
#> :038: if (FLAG == 9 & EVID == 0 & OSCENS == 1) {
#> :039: IPRED = SUR
#> :040: Y = IPRED
#> :041: }
#> :042: if (FLAG == 9 & EVID == 0 & OSCENS == 0) {
#> :043: IPRED = SUR * HAZN
#> :044: Y = IPRED
#> :045: }
#> ================================================================================
#> RxODE model syntax error:
#> ================================================================================
#> :001: TVSLD0 = 70
#> :002: NSLD0 = SLD0/TVSLD0
#> :003: IECOG = ECOG
#> :004: A(0) = IBASE * 1000
#> :005: MMBAS = IBASE * 1000
#> :006: d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + +KD1/100 * AUC1) * A
#> :007: TUM = A
#> :008: TSR = (TUM - MMBAS)/MMBAS
#> :009: if (time == 0) {
#> :010: WTS = 0
#> :011: TM12 = 0
#> :012: }
#> :013: if (time <= 84) {
#> :014: WTS = TSR
#> :015: TM12 = WTS
#> :016: }
#> :017: else {
#> :018: WTS = WTS
#> :019: }
#> :020: DEL = 1e-06
#> :021: function 'EXP' is not supported in RxODE:
#> d/dt(A2) = LAM * SHP * (LAM * (T + DEL))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * WTS + BNWLS * NWLS + BECOG * IECOG)
#> ^~~
#> :021: d/dt(A2) = LAM * SHP * (LAM * (T + DEL))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * WTS + BNWLS * NWLS + BECOG * IECOG)
#> :022: CHZ = A2
#> :023: function 'EXP' is not supported in RxODE:
#> SUR = EXP(-CHZ)
#> ^~~
#> :023: SUR = EXP(-CHZ)
#> :024: DELX = 1e-06
#> :025: xTUM = A
#> :026: xTSR = (xTUM - MMBAS)/MMBAS
#> :027: if (time == 0) {
#> :028: XWTS = 0
#> :029: XTM12 = 0
#> :030: }
#> :031: if (time <= 84) {
#> :032: XWTS = XTSR
#> :033: XTM12 = XWTS
#> :034: }
#> :035: else {
#> :036: XWTS = XWTS
#> :037: }
#> :038: function 'EXP' is not supported in RxODE:
#> HAZN = LAM * SHP * (LAM * (TIME + DELx))^(SHP - 1) * EXP(BSLD0 * NSLD0 + BTSR * XWTS + BNWLS * NWLS + BECOG * IECOG)
#> ^~~
#> :038: if (FLAG == 9 & EVID == 0 & OSCENS == 1) {
#> :039: IPRED = SUR
#> :040: Y = IPRED
#> :041: }
#> :042: if (FLAG == 9 & EVID == 0 & OSCENS == 0) {
#> :043: IPRED = SUR * HAZN
#> :044: Y = IPRED
#> :045: }
#> ================================================================================
#> Error in RxODE::rxState(rxode): Evaluation error: syntax errors (see above).Created on 2022-03-30 by the reprex package (v2.0.1)
mattfidler
commented
2 years ago This parses correctly:
myModel_OS <- function(){ # TTE model for OS. Weibull baseline hazard.
ini({
tLAM <- log(0.001) # scale parameter
tSHP <- log(2) # shape parameter
tBSLD0 <- log(0.1) # parameter for SLD (sum of longest diameter) at enrolment
tBTSR <- log(0.1) # parameter for TSR(t) (tumour size ratio)
tBNWLS <- log(0.1) # parameter for NewLesion(t)
tBECOG <- log(0.1) # parameter for ECOG at enrolment
eta.LAM ~ 0.01 # inter-individual variability in scale parameter
add.err <- 0
})
model({
LAM = exp(tLAM + eta.LAM)
SHP = exp(tSHP)
BSLD0 = exp(tBSLD0)
BTSR = exp(tBTSR)
BNWLS = exp(tBNWLS)
BECOG = exp(tBECOG)
#Time constant covariates
TVSLD0 = 70 # average SLD at enrolment
NSLD0 = SLD0/TVSLD0 # normalised SLD at enrolment
IECOG=ECOG # ECOG at enrolment
# Model for dSLD(t)
A(0) = IBASE*1000 # SLD baseline
MMBAS = IBASE*1000 # Parameter of the SLD(t) model were estimated previously (IPP approach)
d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + + KD1/100 * AUC1) * A
TUM = A
TSR = (TUM-MMBAS)/MMBAS
if(time==0){ #TSR(t) for t<= week 12 and TSR(week12) for t> week 12
WTS = 0
TM12 = 0
}
if(time<=84){
WTS = TSR
TM12 = WTS
}
else {
WTS = WTS
}
# survival model
DEL = 1E-6
d/dt(A2) = LAM*SHP*(LAM*(T+DEL))**(SHP-1) *exp(BSLD0*NSLD0+BTSR*WTS+BNWLS*NWLS+BECOG*IECOG)
# Death hazard
CHZ = A2
SUR=exp(-CHZ)
DELX = 1E-6
xTUM = A
xTSR = (xTUM-MMBAS)/MMBAS
if(time==0){
XWTS=0
XTM12=0
}
if(time<=84){
XWTS = XTSR
XTM12 = XWTS
}
else {
XWTS=XWTS
}
HAZN = LAM*SHP*(LAM*(TIME+DELx))**(SHP-1)*exp(BSLD0*NSLD0+BTSR*XWTS+BNWLS*NWLS+BECOG*IECOG)
# prediction
if(FLAG == 9 & EVID == 0 & OSCENS == 1){
IPRED=SUR # probability of survival (censored event)
Y = IPRED # Y is probability for TTE data
}
if(FLAG == 9 & EVID == 0 & OSCENS == 0){
IPRED=SUR*HAZN # probability of event (death) at time=TIME
Y = IPRED # Y is probability for TTE data
}
Y ~ add(add.err)
})
}
library(nlmixr)
rxClean()
nlmixr(myModel_OS)
#> i parameter labels from comments will be replaced by 'label()'
#> __ RxODE-based ODE model _______________________________________________________
#> -- Initialization: -------------------------------------------------------------
#> Fixed Effects ($theta):
#> tLAM tSHP tBSLD0 tBTSR tBNWLS tBECOG
#> -6.9077553 0.6931472 -2.3025851 -2.3025851 -2.3025851 -2.3025851
#>
#> Omega ($omega):
#> eta.LAM
#> eta.LAM 0.01
#>
#> Covariates or Uninitialized Parameters ($all.covs)
#> [1] "SLD0" "ECOG" "IBASE" "KG" "KD0" "AUC0" "KD1" "AUC1"
#> [9] "T" "NWLS" "XTSR" "TIME" "DELx" "FLAG" "EVID" "OSCENS"
#> -- mu-referencing ($muRefTable): -----------------------------------------------
#> +-------------------+
#> ¦ theta ¦ eta ¦
#> +---------+---------¦
#> ¦ tLAM ¦ eta.LAM ¦
#> +-------------------+
#>
#> -- Model: ----------------------------------------------------------------------
#>
#> LAM = exp(tLAM + eta.LAM)
#> SHP = exp(tSHP)
#> BSLD0 = exp(tBSLD0)
#> BTSR = exp(tBTSR)
#> BNWLS = exp(tBNWLS)
#> BECOG = exp(tBECOG)
#>
#> #Time constant covariates
#> TVSLD0 = 70 # average SLD at enrolment
#> NSLD0 = SLD0/TVSLD0 # normalised SLD at enrolment
#> IECOG=ECOG # ECOG at enrolment
#>
#> # Model for dSLD(t)
#> A(0) = IBASE*1000 # SLD baseline
#> MMBAS = IBASE*1000 # Parameter of the SLD(t) model were estimated previously (IPP approach)
#> d/dt(A) = KG/1000 * A - (KD0/1000 * AUC0 + + KD1/100 * AUC1) * A
#> TUM = A
#> TSR = (TUM-MMBAS)/MMBAS
#>
#> if(time==0){ #TSR(t) for t<= week 12 and TSR(week12) for t> week 12
#> WTS = 0
#> TM12 = 0
#> }
#> if(time<=84){
#> WTS = TSR
#> TM12 = WTS
#> }
#> else {
#> WTS = WTS
#> }
#>
#> # survival model
#> DEL = 1E-6
#> d/dt(A2) = LAM*SHP*(LAM*(T+DEL))**(SHP-1) *exp(BSLD0*NSLD0+BTSR*WTS+BNWLS*NWLS+BECOG*IECOG)
#>
#> # Death hazard
#> CHZ = A2
#> SUR=exp(-CHZ)
#> DELX = 1E-6
#> xTUM = A
#> xTSR = (xTUM-MMBAS)/MMBAS
#> if(time==0){
#> XWTS=0
#> XTM12=0
#> }
#>
#> if(time<=84){
#> XWTS = XTSR
#> XTM12 = XWTS
#> }
#> else {
#> XWTS=XWTS
#> }
#>
#> HAZN = LAM*SHP*(LAM*(TIME+DELx))**(SHP-1)*exp(BSLD0*NSLD0+BTSR*XWTS+BNWLS*NWLS+BECOG*IECOG)
#> # prediction
#> if(FLAG == 9 & EVID == 0 & OSCENS == 1){
#> IPRED=SUR # probability of survival (censored event)
#> Y = IPRED # Y is probability for TTE data
#> }
#> if(FLAG == 9 & EVID == 0 & OSCENS == 0){
#> IPRED=SUR*HAZN # probability of event (death) at time=TIME
#> Y = IPRED # Y is probability for TTE data
#> }
#> Y ~ add(add.err)
#>
#> ________________________________________________________________________________Created on 2022-03-30 by the reprex package (v2.0.1)
mattfidler
commented
2 years ago Also additive error initial estimate should not start at 0...
mattfidler
commented
2 years ago I looked through the ddmore files and they use phi(dum) and no error; This means they are specifying the likelihood. Unfortunately this isn't supported yet in nlmixr model functions
andreasmeid
commented
2 years ago Many thanks for these very helpful notes. Is there a way to fit the parametric survival model with another R package, for example? I could calculate the tumour size (ODE for A in the code above) with RxODE, for example. The subsequent weibull survival model can be seen separate and thus could be fit in a different step, I guess.
Hi,
I just tried to reproduce a ddmore-NONMEM analysis with nlmixr (data and NONMEM code are available from http://repository.ddmore.eu/model/DDMODEL00000217). However, I receive the error message "SAEM requires all ETAS to be mu-referenced" which I do not understand (as all etas appear to be defined correctly). Can you figure our what is wrong (with the model formula or data specification)?
Many thanks and best wishes Andreas
library(tidyverse) TS <- read.table("Simulated SLD.csv", header=T, sep=",") %>% select(-CCOM)
TS2 <- TS %>% group_by(ID) %>% mutate(BQL_IND = case_when(any(BQL == 1) ~ 1, all(BQL != 1) ~ 0)) %>% filter(BQL_IND != 1) %>% select(ID, TIME, DV, AUC0, AUC1, EVID, FLAG, CMT, ECOG) %>% mutate(dvid = case_when(EVID == 0 ~ "cp", EVID != 0 ~ "other")) %>% rename(cmt=CMT, dv=DV, id=ID) %>% mutate(dv=as.numeric(dv), cmt=as.numeric(cmt)) %>% filter(dvid=="cp") # simplifications for first run (e.g., remove subjects with values below LLOQ)
library(nlmixr)
myModel <- function() { # Claret-like model for describing CTS. ini({ tKG <- log(0.3) tKD0 <- log(0.03) tKD1 <- log(0.01) tIBASE <- log(0.01) eta.KG <- 0.08 eta.KD <- 0.1 eta.IBASE <- 0.1 add.err <- 20 # SD of additive unexpained variability, SD of additive error (mm) }) model({ KG = exp(tKG + eta.KG) # tumour growth rate constant (1/day) KD0 = exp(tKD0 + eta.KD) # carboplatin related death rate constant (1/day/AUC0) KD1 = exp(tKD1 + eta.KD) # gemcitabine related death rate constant (1/day/AUC1) IBASE = exp(tIBASE + eta.IBASE) # baseline SLD (m) E0 = AUC0 E1 = AUC1 A(0) = IBASE1000 # SLD baseline d/dt(A) = KG/1000 A - (KD0/1000 E0 + KD1/100 E1) * A # Model for dSLD(t) cp = A cp ~ add(add.err) }) }
fit.myModel <- nlmixr(myModel,TS2,est="saem",control=saemControl(n.burn=50,n.em=100,print=50));