Hi @andrewhooker, I have a NM code in which the error model is described as fixed effect. I wish to correctly transfer it to the popED. I have two questions pertaining to the code below:

How to correctly transfer the combined error model code to popED?
In NM the proportional and additive effects are modeled as fixed effects. Hence I have SD estimated as 0.268 and 20.8. How I should be transfering these values to the popED under sigma=c(0.268,20.8), # decides the variances of RV? Thanks, -Ankur

My NM Code as follows

; $SUBROUTINE ADVAN3 TRANS4 ; $PK CLALLO = (BWT/93.0)THETA(7) VALLO = (BWT/93.0)THETA(8) TVCL = THETA(3)CLALLO TVV1 = THETA(4)VALLO TVQ = THETA(5) TVV2 = THETA(6) CL = TVCL EXP(ETA(1)) V1 = TVV1 EXP(ETA(2)) Q = TVQ EXP(ETA(3)) V2 = TVV2 EXP(ETA(4)) S1 = V1/1000 S2 = V2/1000

$ERROR IPRED = F IRES = DV-IPRED ERRADD = THETA(1) ERRPROP = THETA(2) W = SQRT(ERRPROP*2IPRED2 + ERRADD2) DEL = 0 IF(W.EQ.0) DEL = 1 IWRES = IRES/(W + DEL) Y = IPRED + W*EPS(1)

$THETA
(0,100) ; TH1_ERRADD (0,1) ; TH2_ERRPROP (0,0.1) ; TH3_CL (0,2) ; TH4_V1 (0,1) ; TH5_Q (0,2) ; TH6_V2 0.75 ; TH7_CLALLO 1.0 ; TH8_VALLO

$OMEGA
0.1 ; n_CL 0.1 ; n_V1 0 FIX ; n_Q 0.1 ; n_V2

$SIGMA 1 FIX $ESTIMATION METHOD=1 INTER POSTHOC MAXEVAL=9999 NOABORT PRINT=5

My popED code is as below:

PK_2comp_ode_inf <- function(Time, State, Pars){ with(as.list(c(State, Pars)), { CL=CL(WT/93)^(WT_CL) V1=V1(WT/93)^(WT_V1) if( Time %% TAU < TINF ){In <- DOSE/TINF} else { In <- 0 } dA1 <- In - CL/V1A1- Q/V1A1+ Q/V2A2 dA2 <- Q/V1A1-Q/V2*A2

return(list(c(dA1, dA2)))

})

} PK_2comp_ode_inf_ff <- function(model_switch, xt, parameters, poped.db){ with(as.list(parameters),{

initial conditions of ODE system

A_ini <- c(A1=0, A2=0)

#Set up time points to get ODE solutions

times_xt <- drop(xt)# sample times event_times <- seq(from=0,to=max(times_xt),by=TAU)# dose times times <- c(0,times_xt,event_times) ## add extra time for start of integration

times<-seq(0,1344,10)

times <- sort(times) 
times <- unique(times) # remove duplicates

# Dosing
eventdat <- data.frame(var = c("A1"), 
                       time = event_times,
                       value = c(DOSE), method = c("add"))
out <- ode(A_ini, times, PK_2comp_ode_inf, parameters,events = list(data = eventdat),atol=1e-8,rtol=1e-8)
y = out[, "A1"]/(V1)
y=y[match(times_xt,out[,"time"])]
y=cbind(y)
return(list(y=y,poped.db=poped.db))

}) }

-- parameter definition function

-- names match parameters in function ff

sfg <- function(x,a,bpop,b,bocc){

parameters=c( CL=bpop[1]exp(b[1]), V1=bpop[2]exp(b[2]), Q= bpop[3]exp(b[3]), V2=bpop[4]exp(b[4]),

BL=bpop[5]*exp(b[5]),

            WT_CL=bpop[5],
            WT_V1=bpop[6],
            DOSE=a[1],
            TINF=a[2],
            TAU= a[3],
            WT=  a[4])

return( parameters )

}

-- Residual unexplained variablity (RUV) function

-- Proportional + additive

feps.add.prop <- function(model_switch,xt,parameters,epsi,poped.db){ returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) returnArgs <- ff(model_switch,xt,parameters,poped.db) y <- returnArgs[[1]] poped.db <- returnArgs[[2]]

y = y*(1+epsi[,1])+epsi[,2]

return(list( y= y,poped.db =poped.db )) }

Original Design

poped.db_original <- create.poped.database(ff_file = "PK_2comp_ode_inf_ff", fg_file="sfg", fError_file="feps.add.prop",

                                       bpop=c(CL=0.221,V1=4.07,Q=0.425,V2=3.63, WT_CL=1.15,WT_V1=0.795), # run557.lst  # units L/day and L for CL and V resp
                                       notfixed_bpop =  c(1,1,0,1,1,1), #use 1 to denote "not fixed" (estimated) and 0 to denote "not not fixed" (not estimated)

                                       d=c(CL=0.241,V1=0.12,Q=0,V2=0.456), # decides the variances of BSV
                                       notfixed_d = c(1,1,0,1), #use 1 to denote "not fixed" (estimated) and 0 to denote "not not fixed" (not estimated)

                                       sigma=c(0.268,20.8), # decides the variances of RV
                                       notfixed_sigma = c(0,0), #use 1 to denote "not fixed" (estimated) and 0 to denote "not not fixed" (not estimated)

                                       groupsize=1095, #proposed study, see slde deck mentioned above ref
                                       m=2,      #number of groups

                                       #        ± 1 ± 1 ± 1 ± 1 ± 2 ± 2 ± 2 ± 2 ± 2 ± 2 ± 2 interval
                                       xt=   c( -0.5,168, 2016, 4032, 8736),

                                       a=list(c(DOSE=5,TINF=0.016,TAU=672, WT=95),c(DOSE=0,TINF=0.016,TAU=672, WT=95)))# Unit of dose is mg, # added a very fast IV infusion of 1 min

andrewhooker / PopED

How to transfer PROP and ADD error to popED if estimated as fixed effects in NM? #75