I'm using WinBUGS to estimate parameters for individual growth models. I wanted to compare two models, one where the growth parameters are constant and one where some of the parameters are a function of covariates (sex and a spatial covariate I've called "gila"). I'm confused because pD, which I thought represented the effective number of parameters, is higher for a model without covariates (~715) than a model with covariates (~695). Since DIC is a function of pD, I'm wary about using DIC to compare support for the two models.
I was wondering whether anyone could help me understand why pD might be higher for models with fewer covariates? Maybe it's related to the fact that I have random measurement error incorporated into the models?
I attached relevant code for the two models below.
Any suggestions or help would be much appreciated!
Thank you,
Erin
Covariate model:
#############################################################
#Individual growth model (Richards) based on capture-recapture data
#Assuming parameters (A,k) a function of covariates
#Assuming hatchling size = 45mm
#############################################################
library(R2WinBUGS)
#objects for the bugs function
#mcl=length at each capture; time=time between captures (yrs); n=number of captures
data <- list ("ind","n","mcl","time","gila","sex")
inits <- function(){
tFCi<-sample(10:50,length(n),replace=T)
list (A.mu=runif(1,250,300),A.gila=rnorm(1,0,1),A.sex=rnorm(1,0,1),logit.kmu=rnorm(1,0,1),k.gila=rnorm(1,0,1),k.sex=rnorm(1,0,1),d=runif(1,0,2),meas.sigma=runif(1,0,10),tFC=tFCi)
}
params <- c("A.mu","A.gila","A.sex","k.mu","k.gila","k.sex","logit.kmu","d","meas.sigma")
#number of iterations
ni<-300000
#number to burn in
nb<-100000
#number to thin
nthin<-50
#number of chains
nc<-3
modelFilename = 'Growth_test.txt'
cat("
model{
for (i in 1:ind){
for(j in 1:1){
mcl[i,j]~dnorm(mcl.exp[i,j],meas.tau) #measurement error
mcl.exp[i,j]<-A[i]*pow((1+(pow((45/A[i]),(1-d))-1)*exp(-k[i]*tFC[i])),(1/(1-d))) #mcl at first capture is related to A, k, d, and age (tFC)
tFC[i]~dunif(0,50) #age at first capture could be anywhere from 0-50 years
} #close j loop
for (j in 2:n[i]){
mcl[i,j]~dnorm(mcl.exp[i,j],meas.tau)
#relating mcl from recapture to mcl from prior capture
mcl.exp[i,j]<-A[i]*pow((1+(pow((mcl.exp[i,j-1]/A[i]),(1-d))-1)*exp(-k[i]*time[i,j-1])),(1/(1-d)))
} #close j loop
A[i] <- A.mu + A.gila*gila[i] + A.sex*sex[i]
logit(k[i]) <- logit.kmu + k.gila*gila[i] + k.sex*sex[i]
} #close i loop
#Priors
A.mu~dnorm(270,0.001)
A.gila~dnorm(0,0.001)
A.sex~dnorm(0,0.001)
logit.kmu~dnorm(0,0.001)
k.gila~dnorm(0,0.001)
k.sex~dnorm(0,0.001)
d~dunif(0,2)
k.mu<-exp(logit.kmu)/(1+exp(logit.kmu))
meas.tau<-pow(meas.sigma,-2)
meas.sigma~dunif(0,100)
} #close model
",file=modelFilename)
Constant model:
#############################################################
#Individual growth model (Richards) based on capture-recapture data
#Assuming parameters (A,k) constant
#Assuming hatchling size = 45mm
#############################################################
library(R2WinBUGS)
#objects for the bugs function
data <- list ("ind","n","mcl","time") #mcl=length at each capture; time=time between captures (yrs)
inits <- function(){
tFCi<-sample(10:50,length(n),replace=T)
list (A.mu=runif(1,250,300),logit.kmu=rnorm(1,0,1),d=runif(1,0,2),meas.sigma=runif(1,0,10),tFC=tFCi)
}
params <- c("A.mu","k.mu","d","meas.sigma")
#number of iterations
ni<-300000
#number to burn in
nb<-100000
#number to thin
nthin<-50
#number of chains
nc<-3
modelFilename = 'Growth_test.txt'
cat("
model{
for (i in 1:ind){
for(j in 1:1){
mcl[i,j]~dnorm(mcl.exp[i,j],meas.tau)
mcl.exp[i,j]<-A[i]*pow((1+(pow((45/A[i]),(1-d))-1)*exp(-k[i]*tFC[i])),(1/(1-d)))
tFC[i]~dunif(0,50)
} #close j loop
for (j in 2:n[i]){
mcl[i,j]~dnorm(mcl.exp[i,j],meas.tau)
#relating mcl from recapture to mcl from prior capture
mcl.exp[i,j]<-A[i]*pow((1+(pow((mcl.exp[i,j-1]/A[i]),(1-d))-1)*exp(-k[i]*time[i,j-1])),(1/(1-d)))
} #close j loop
A[i] <- A.mu
logit(k[i]) <- logit.kmu
} #close i loop
#Priors
A.mu~dnorm(270,0.001)
logit.kmu~dnorm(0,0.001)
d~dunif(0,2)
k.mu<-exp(logit.kmu)/(1+exp(logit.kmu))
meas.tau<-pow(meas.sigma,-2)
meas.sigma~dunif(0,100)
} #close model
",file=modelFilename)
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|