Issue with uneven number of samples in groups

[lindauer1980]

Issue with uneven number of samples in groups

[andrewhooker]

Hi Andreas,

I have a hole in the poped code where model_switch, when entered as a list, is not transformed to the needed matrix format. I have added a fix to the latest version of the development version of PopED. If you don't want to upgrade, then a simple fix is the following:

library(PopED)

Nfull <- 20
sampling.times.full <- c(0,671.5,2015.5,2024,2028,2032,2040,2064,2088,2112,2184,4416,6432,8732) # full sampling group
N.full <- Nfull# number of subjects in the full sampling group
sampling.times.sparse <- c(0,671.5,2064,4416,6432,8732) # sparse sampling group
N.sparse <- 100-Nfull # number of subjects in the sparse group

# fix by making model_switch a matrix before adding to function argument
model_switch = list(sort(rep(c(1,2,3),length(sampling.times.full))),
                    sort(rep(c(1,2,3),length(sampling.times.sparse))))

model_switch <- t(sapply(model_switch,'[',seq(max(sapply(model_switch,length)))))

create_design(
  groupsize=c(N.full,N.sparse),
  xt=list(rep(sampling.times.full,3),
          rep(sampling.times.sparse,3)),
  model_switch = model_switch
) 

Best regards,
Andy

[lindauer1980]

Thank you Andy for your swift reply,

While the code to convert the model_switch from list to matrix works, running create.poped.database gives the following error:

Error in if (mat[i, j] == 0) { : missing value where TRUE/FALSE needed

Both the xt and the model_switch objects contain NAs for the 'missing' timepoints in the second group, may this cause the error ?

Regards, Andreas.

[andrewhooker]

Hi

This will not cause an error. Can you send code? How are you using create.poped.database?

Best regards,
Andy

On 21 Jun 2016, at 14:41 , lindauer1980 [email protected] wrote:

Thank you Andy for your swift reply,

While the code to convert the model_switch from list to matrix works, running create.poped.database gives the following error:

Error in if (mat[i, j] == 0) { : missing value where TRUE/FALSE needed

Both the xt and the model_switch objects contain NAs for the 'missing' timepoints in the second group, may this cuase the error ?

Regards, Andreas.

—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub #3 (comment), or mute the thread https://github.com/notifications/unsubscribe/AFLmZPsBqeUr25cza6zoFsew7PIAiESyks5qN9wJgaJpZM4I5xi5.

[lindauer1980]

here is the code:

create.poped.database(ff_file="ff.model.ode.compiled",
                                           fError_file="feps.compiled",
                                           fg_file="fg",
                                           groupsize=c(N.full,N.sparse),
                                           m=2,      #number of groups
                                           sigma=c(propRUV.tcon^2,addRUV.tcon^2,
                                                   propRUV.hgh^2,addRUV.hgh^2,
                                                   propRUV.igf^2,addRUV.igf^2),
                                           bpop=thetas,
                                           d=omegas,
                                           notfixed_bpop=not.fixed.th,
                                           notfixed_d = not.fixed.om,
                                           xt=list(rep(sampling.times.full,3),
                                                   rep(sampling.times.sparse,3)),
                                           model_switch=model_switch,
                                           minxt=0,
                                           maxxt=max.time,#max(sampling.times),
                                           a=c(doseBW,tau,Ndoses,BW), 
                                           maxa=c(10,300,60,100),
                                           mina=c(0,8,1,1))

All works well if the xt and model_switch are of the same size without NAs

Thanks.

[lindauer1980]

Since I can't share more details on my model, I have adapted one of your PKPD examples to reproduce the error. Hope it helps:

library(PopED)

# This option is used to make this script run fast but without convergence 
# (fast means a few seconds for each argument at the most).
# This allows you to "source" this file and easily see how things work
# without waiting for more than 10-30 seconds.
# Change to FALSE if you want to run each function so that
# the solutions have converged (can take many minutes).
fast <- TRUE 

EAStepSize <- ifelse(fast,40,1)


ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption + inhibitory imax
  ## -- works for both mutiple and single dosing  
  with(as.list(parameters),{

    y=xt
    MS <- model_switch

    # PK model
    N = floor(xt/TAU)+1
    CONC=(DOSE*Favail/V)*(KA/(KA - CL/V)) * 
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - 
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))  

    # PD model
    EFF = E0*(1 - CONC*IMAX/(IC50 + CONC))

    y[MS==1] = CONC[MS==1]
    y[MS==2] = EFF[MS==2]

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

sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c( V=bpop[1]*exp(b[1]),
                KA=bpop[2]*exp(b[2]),
                CL=bpop[3]*exp(b[3]),
                Favail=bpop[4],
                DOSE=a[1],
                TAU = a[2],
                E0=bpop[5]*exp(b[4]),
                IMAX=bpop[6],
                IC50=bpop[7])
  return( parameters ) 
}



feps <- function(model_switch,xt,parameters,epsi,poped.db){
  ## -- Residual Error function
  returnArgs <- ff(model_switch,xt,parameters,poped.db) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]

  MS <- model_switch

  pk.dv <- y*(1+epsi[,1])+epsi[,2]
  pd.dv <-  y*(1+epsi[,3])+epsi[,4]

  y[MS==1] = pk.dv[MS==1]
  y[MS==2] = pd.dv[MS==2]

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


N.full <- 20
sampling.times.full   <- c(0,671.5,2015.5,2024,2028,2032,2040,2064,2088,2112,2184,4416,6432,8732) # full sampling group
N.full                <- Nfull# number of subjects in the full sampling group
sampling.times.sparse <- c(0,671.5,2064,4416,6432,8732) # sparse sampling group
N.sparse              <- 100-Nfull # number of subjects in the sparse group


# fix by making model_switch a matrix before adding to function argument
model_switch = list(sort(rep(c(1,2),length(sampling.times.full))),
                    sort(rep(c(1,2),length(sampling.times.sparse))))

model_switch <- t(sapply(model_switch,'[',seq(max(sapply(model_switch,length)))))


poped.db <- create.poped.database(ff_file="ff",
                                  fError_file="feps",
                                  fg_file="sfg",
                                  groupsize=c(N.full,N.sparse),
                                  m=2,
                                  bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,E0=1120,IMAX=0.807,IC50=0.0993),  
                                  notfixed_bpop=c(1,1,1,0,1,1,1),
                                  d=c(V=0.09,KA=0.09,CL=0.25^2,E0=0.09), 
                                  sigma=c(0.04,5e-6,0.09,100),
                                  notfixed_sigma=c(0,0,0,0),
                                  #xt=c( 1,2,8,240,240,1,2,8,240,240),
                                  xt=list(rep(sampling.times.full,2),
                                          rep(sampling.times.sparse,2)),
                                  #minxt=c(0,0,0,240,240,0,0,0,240,240),
                                  #maxxt=c(10,10,10,248,248,10,10,10,248,248),
                                  #G_xt=c(1,2,3,4,5,1,2,3,4,5),
                                  #model_switch=c(1,1,1,1,1,2,2,2,2,2),
                                  model_switch = model_switch,
                                  a=cbind(c(20,40),c(24,24)),
                                  bUseGrouped_xt=1,
                                  ourzero=0,
                                  maxa=c(200,40),
                                  mina=c(0,2))

[andrewhooker]

Hi Andreas,

Sorry for your inconvenience.

This is a bug that has to do with PopED checking if the range of possible sample times includes zero. If so, then the default is to move that range slightly away from zero, since in PK studies, zero observations give no information. I have fixed the bug in the development version. To fix in the current version I suggest using the ourzero=NULL option in create.poped.database.

Best Regards
Andy

library(PopED)


sfg <- function(x,a,bpop,b,bocc){
  parameters=c(CL=bpop[1]*exp(b[1]),
               V=bpop[2]*exp(b[2]),
               KA=bpop[3]*exp(b[3]),
               Favail=bpop[4],
               DOSE=a[1])
  return(parameters) 
}

Nfull <- 20
sampling.times.full <- c(0,671.5,2015.5,2024,2028,2032,2040,2064,2088,2112,2184,4416,6432,8732) # full sampling group
N.full <- Nfull# number of subjects in the full sampling group
sampling.times.sparse <- c(0,671.5,2064,4416,6432,8732) # sparse sampling group
N.sparse <- 100-Nfull # number of subjects in the sparse group

# fix by making model_switch a matrix before adding to function argument
model_switch = list(sort(rep(c(1,2,3),length(sampling.times.full))),
                    sort(rep(c(1,2,3),length(sampling.times.sparse))))

model_switch <- t(sapply(model_switch,'[',seq(max(sapply(model_switch,length)))))

design <- create_design(
  groupsize=c(N.full,N.sparse),
  xt=list(rep(sampling.times.full,3),
          rep(sampling.times.sparse,3)),
  model_switch = model_switch
) 


## -- Define initial design  and design space
poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL",
                                  fg_file="sfg",
                                  fError_file="feps.add.prop",
                                  bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), 
                                  notfixed_bpop=c(1,1,1,0),
                                  d=c(CL=0.07, V=0.02, KA=0.6), 
                                  sigma=c(0.01,0.25),
                                  ourzero = NULL,
                                  groupsize=c(N.full,N.sparse),
                                  xt=list(rep(sampling.times.full,3),
                                          rep(sampling.times.sparse,3)),
                                  model_switch = model_switch,
                                  a=70)

[andrewhooker]

By the way... if you use the ourzero=NULL option then you need to be more careful about the choice of the maximum and minimum allowed values of your design parameters (make sure they do not produce no information). In the next release version I will have a more transparent way of dealing with this issue.

[lindauer1980]

Thank you Andy, and sorry for continuing to bother you... when I try to generate prediction from the above code I get the following:

> plot_model_prediction(poped.db)
Error in seq.default(minv, maxv, length.out = model_num_points[j]) : 
  'from' cannot be NA, NaN or infinite
> model_prediction(poped.db)
Error in if (!all(models_to_use == unique(as.vector(model_switch[used_times ==  : 
  missing value where TRUE/FALSE needed

How can I fix this ?

[andrewhooker]

Hi Andreas,

Sorry for all your problems (all my fault!). I think the easiest way to go is to define the xt parameter and the model switch parameter as matrices that have dimensions of rows=(# of groups) and cols=(max number of observations in any group). Make the points that were previously NA be any number you like. You can the set the number of observations per group explicitly with the niargument (see below). All this is fixed in the development release, so if you want to try you can install that using

devtools::install_github("andrewhooker/PopED")

Otherwise, see my code below:

library(PopED)

# This option is used to make this script run fast but without convergence 
# (fast means a few seconds for each argument at the most).
# This allows you to "source" this file and easily see how things work
# without waiting for more than 10-30 seconds.
# Change to FALSE if you want to run each function so that
# the solutions have converged (can take many minutes).
fast <- TRUE 

EAStepSize <- ifelse(fast,40,1)


ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption + inhibitory imax
  ## -- works for both mutiple and single dosing  
  with(as.list(parameters),{

    y=xt
    MS <- model_switch

    # PK model
    N = floor(xt/TAU)+1
    CONC=(DOSE*Favail/V)*(KA/(KA - CL/V)) * 
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - 
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))  

    # PD model
    EFF = E0*(1 - CONC*IMAX/(IC50 + CONC))

    y[MS==1] = CONC[MS==1]
    y[MS==2] = EFF[MS==2]

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

sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c( V=bpop[1]*exp(b[1]),
                KA=bpop[2]*exp(b[2]),
                CL=bpop[3]*exp(b[3]),
                Favail=bpop[4],
                DOSE=a[1],
                TAU = a[2],
                E0=bpop[5]*exp(b[4]),
                IMAX=bpop[6],
                IC50=bpop[7])
  return( parameters ) 
}



feps <- function(model_switch,xt,parameters,epsi,poped.db){
  ## -- Residual Error function
  returnArgs <- ff(model_switch,xt,parameters,poped.db) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]

  MS <- model_switch

  pk.dv <- y*(1+epsi[,1])+epsi[,2]
  pd.dv <-  y*(1+epsi[,3])+epsi[,4]

  y[MS==1] = pk.dv[MS==1]
  y[MS==2] = pd.dv[MS==2]

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


Nfull <- 20
sampling.times.full   <- c(0,671.5,2015.5,2024,2028,2032,2040,2064,2088,2112,2184,4416,6432,8732) # full sampling group
N.full                <- Nfull# number of subjects in the full sampling group
sampling.times.sparse <- c(0,671.5,2064,4416,6432,8732) # sparse sampling group
N.sparse              <- 100-Nfull # number of subjects in the sparse group


# fix by making model_switch a matrix before adding to function argument
model_switch = list(sort(rep(c(1,2),length(sampling.times.full))),
                    sort(rep(c(1,2),length(sampling.times.sparse))))

model_switch <- t(sapply(model_switch,'[',seq(max(sapply(model_switch,length)))))

model_switch[is.na(model_switch)] <- 0

xt=list(rep(sampling.times.full,2),
        rep(sampling.times.sparse,2))
xt <- t(sapply(xt,'[',seq(max(sapply(xt,length)))))

xt[is.na(xt)] <- 0


poped.db <- create.poped.database(ff_file="ff",
                                  fError_file="feps",
                                  fg_file="sfg",
                                  groupsize=c(N.full,N.sparse),
                                  m=2,
                                  bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,E0=1120,IMAX=0.807,IC50=0.0993),  
                                  notfixed_bpop=c(1,1,1,0,1,1,1),
                                  d=c(V=0.09,KA=0.09,CL=0.25^2,E0=0.09), 
                                  sigma=c(0.04,5e-6,0.09,100),
                                  notfixed_sigma=c(0,0,0,0),
                                  #xt=c( 1,2,8,240,240,1,2,8,240,240),
                                  xt=xt,
                                  ni=c(28,12),
                                  #minxt=c(0,0,0,240,240,0,0,0,240,240),
                                  #maxxt=c(10,10,10,248,248,10,10,10,248,248),
                                  #G_xt=c(1,2,3,4,5,1,2,3,4,5),
                                  #model_switch=c(1,1,1,1,1,2,2,2,2,2),
                                  model_switch = model_switch,
                                  a=cbind(c(20,40),c(24,24)),
                                  #bUseGrouped_xt=1,
                                  ourzero=NULL,
                                  maxa=c(200,40),
                                  mina=c(0,2))


plot_model_prediction(poped.db)
det(evaluate.fim(poped.db ))

[lindauer1980]

Thanks a lot Andy for your help. Now it works fine!