[R-sig-dyn-mod] Fw: nonlinear mixed effect models with deSolve + brms

[Andras Farkas]

sorry here it is in plain 


Hello all,

wonder if you could help with setting up code for running ODE based models with brms package... ODE setup

parms<-list(k = 0.09,v=18)
tofit<-data.frame(ID=c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5),
                  t=c(c(1,3,5,7,10),c(1,3,5,7,10)-0.4,c(1,3,5,7,10)+0.6,c(1,3,5,7,10)-0.3,c(1,3,5,7,10)+1.2),
                  y=c(c(25,21,NA,15,12),c(20,NA,16,14,12),c(31,29,NA,NA,22),c(32,NA,22,NA,14),c(18,15,13,NA,8.2)))

intimes <- c(0,0.5,max(tofit$t)+0.1)
input   <- c(800,0,0)
forc <- approxfun(intimes, input, method="constant", rule=2)
derivs <- function(t, state, parms) {
  inp <- forc(t)
  dy1 <- - parms$k * state[1] + inp/parms$v
  return(list(c(dy1)))

}
model <- function(parms, times=seq(0, max(tofit$t), by = 0.1)) {
  state <- c(y = 0)
  return(ode(y = state, times = times, func = derivs, parms=parms,
             method="lsoda"))
}

#next, trying to follow example on https://rpubs.com/aforren1/orange-nonlinear and model nonlinear mixed effects on both #parameters with data to be fitted in data frame tofit for each 'ID' at time points 't' is named 'y' 

library(brms)
prior_1 <- c(set_prior("normal(0.07, 0.017)", nlpar = "k"),
             set_prior("normal(15, 3)", nlpar = "v")
             )

prior_2 <- rbind(prior_1, set_prior("normal(0.5,0.2)", class = "sd"))

n3_b <- brm(y~model(parms), 
            data = tofit,
            nonlinear = list(k ~ (1|ID),
                             v ~ (1|ID) 
                            ),
            prior = prior_2,                
            chains = 3,
            iter = 500)

is as far as I have got... important for parameters to be estimated is not to generate values zero or below, so I will have to further work to understand lognormal in STAN, but goal is to get results for fixed and and random effects... I am aware of Stan having an ODE solver but goal is to stay with deSolve for ODE as ... I am also aware of other algorithms/packages (nlmeODE, saemix, BUGS, JAGS, etc) available for parameter estimation but I would like to stay with a bayesian system for nonlinear mixed effect models that would allow continued use of deSolve... As an alternative to solving above dilemma I would be interested in adapting a working example of a bayesian nonlinear mixed effect models using deSolve if you are aware of one as my research turned out basically nothing out there that would include a working code for possible adaptation into my own practice... My interest is in pharrmacokinetics/dynamics but I presume I should be able to adapt working examples from other fields (also didn't find much out there that describe implementation)...  


appreciate your help...


Andras