1 year ago
#226051
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