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

[Anks2030]

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:

1. How to correctly transfer the combined error model code to popED?
2. 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(ERRPROP2*IPRED2 + ERRADD**2)
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]

Hi

If you estimate the residual variability on the SD scale then you should square that value when assuming that you estimate the variances in PopED. The expected SE estimates will be approximately twice as large as when you estimate the SD values, but should not affect the design optimization.

Best regards,
Andy