WinBUGS: Nonlinear regression model

Revisiting Example 3.7 (pages 102-3) on non-linear regression with exponential decay

y_i ~ N(b₁- b₁exp(-b₂x_i),t ^-2)

where y are oxygen concentrations and x are the number of days the experiments are carried out.

Data: x = (1,2,3,4,5,7) and y = (8.3,10.3,19,16,15.6,19.8).

Prior: p( b₁,b₂,t²) = p( b₁) p( b₂) p(t²) and

β₁ ~ U(-20,50)

β₂ ~ U(-2,6)

t² ~ G(0.001,0.001)

WINBUGS CODE FOR THE NONLINEAR MODEL

# This model is stored in file “exponentialdecaymodel.bug”

# NONLINEAR MODEL

model{

for (i in 1:N){

mu[i] <- beta1*(1-exp(-beta2*x[i]))

y[i]~dnorm(mu[i],tau2)

}

beta1~dunif(-20,50)

beta2~dunif(-2,6)

tau2 ~ dgamma(0.001,0.001)

}

# DATA

list(x=c(1,2,3,4,5,7),y=c(8.3,10.3,19,16,15.6,19.8),N=6)

# INITIAL VALUES

list(beta1=18,beta2=1.78,tau2=1)

Evoking WinBUGS from R

library(R2WinBUGS)

quant025 = function(x){quantile(x,0.025)}

quant975 = function(x){quantile(x,0.975)}

x = c(1,2,3,4,5,7)

y = c(8.3,10.3,19,16,15.6,19.8)

N = 6

data = list("x","y","N")

inits = function(){

list(beta1=18,beta2=1.78,tau2=1)

}

decay = bugs(data,inits,

model.file="exponentialdecaymodel.bug",

parameters=c("beta1","beta2","tau2"),

n.chains=1,n.iter=100000,n.burnin=80000,n.thin=1,

bugs.directory="c:/Program Files/WinBUGS14/",

codaPkg=FALSE)

# Replicating Figure 3.4, page 104

beta1 = decay$sims.array[,1,1]

beta2 = decay$sims.array[,1,2]

xs = seq(1,7,length=50)

f = matrix(0,20000,50)

for (i in 1:50)

f[,i] = beta1*(1-exp(-beta2*xs[i]))

q025 = apply(f,2,quant025)

q050 = apply(f,2,median)

q975 = apply(f,2,quant975)

par(mfrow=c(1,2))

plot(decay$sims.array[,1,1],decay$sims.array[,1,2],xlim=c(0,50),ylim=c(-1,6),

xlab=expression(beta[1]),ylab=expression(beta[2]),pch=”.”)

plot(x,y,xlab="time (days)",ylab="BOD (mg/liter)",ylim=c(0,30))

lines(xs,q025,col=2,lty=2)

lines(xs,q050,col=2)

lines(xs,q975,col=2,lty=2)