1 year ago

#226051

test-img

Constantin

Fitting a NLME model in jags

Generate some data

##### Load packages
library(nlme); library(lme4)
library(dclone)

### GENERATE SOME DATA
## Example data with a weibull curve 
k       <- 120;
nindivs <- 30;
pwr     <- 0.2
a       <- 500

set.seed(123)
sim.data <- data.frame()
for(i in 1:nindivs){
  time  <- 1:20
  theta <- runif(1,0.5,0.6)
  #y     <- (k/theta)* ( (time/theta)^(k-1) ) * (exp(-(time/theta))^k )
  y     <- k - a*exp( -exp(theta)*time^pwr) + rnorm(length(time))
  ID    <- i

  sim.data <- rbind(data.frame(ID, time, y), sim.data)
}

sim.data    <- na.omit(sim.data)
sim.data$ID <- as.factor(sim.data$ID)
plot(sim.data)

I am trying to fit a nonlinear mixed effects model in JAGS. Using the nlme package, I would simply write the following:

NLME CODE

#### NLME 
nlme.data <- groupedData(y ~ time | ID, data=sim.data)
nlme_mod <- nlme(y ~ SSweibull(time, Asym, Drop, lrc, pwr), 
                 data  = nlme.data, 
                  fixed  = Asym + Drop + pwr ~ 1, 
                  random = lrc ~ 1, 
                  start  = c(Asym=10, Drop=10, pwr=0.5))
summary(nlme_mod)
plot(augPred(nlme_mod, level = 0:1),layout = c(5,1))

How would I do this in JAGS? Here is my attempt:

## Parameters to load
parms  <- c("a", "pwr", "k", "theta")

### Data to run the Bayesian model
len  <- dim(sim.data)[1]; nind=length(unique(sim.data$ID))
bayes.list <- list(X = sim.data$time, Y= sim.data$y, ID = sim.data$ID,
                 len = len, nind=nind)
#### JAGS
weibullBayes <- function(){
    
    ## --------------------------------------
    ## Priors
    shape ~ dunif(0, 100)
    rate  ~ dunif(0, 0.5)
    a     ~ dunif(100, 1000)
    pwr   ~ dunif(0, 0.5)
    k     ~ dunif(0, 0.5)
    
    ## --------------------------------------   
    ## Random effects
    for(r in 1:nind){
        theta[r] ~ dgamma(shape, rate)
    }
    
    ## --------------------------------------
    ## Likelihood
    for(i in 1:len){
        lambda[i] <- k-a*exp( -exp(theta[ID[i]])*X[i]^pwr)
        Y[i] ~ dpois(lambda[i])
    }
}
bayes.fit <- jags.fit(bayes.list, parms, weibullBayes, 
                      n.chains = 50, n.adapt=100, n.update=500, n.iter=50000, thin=50)
summary(bayes.fit)[1]; summary(nlme_mod)

The output for the JAGS model seem wrong, especially the theta and k parameters. Where am I going wrong?

r

nlme

rjags

0 Answers

Your Answer

Accepted video resources