[R-sig-ME] Bayesian Hierarchical Model - Estimating Rho (correlation between intercepts and slopes)

[Jarrod Hadfield]

Hi, 

There's probably little (no) variation in the intercepts and/or slopes.
Is this the case? If so, I think it is a good thing that the correlation
has such wide credible intervals!

Jarrod






On Mon, 2011-07-11 at 11:44 +0100, Xavier Harrison wrote:
> Dear List
> 
> I'm trying to run a hierarchical model to look at the effect of body condition and environmental variables on breeding probability (1/0). The model I'm using is adapted from Gelman and Hill (2006) pg 376, with the modifications to add a group level predictor, for both intercepts and slopes, listed on pg 379. I'm modelling body condition at the individual level, and NAO at the group (year level). These models are being run in JAGS using Denwood's 'runjags' interface for R.
> 
> The issue I'm having is that the 95% credible intervals for rho, the correlation among intercepts and slopes, is -.95 to 0.93, so essentially is spanning the entire available range. I would expect that if the model were 'behaving', then this estimate would take some kind of sensible value i.e. perhaps with  a posterior mode of zero, but much more certainty that it's zero. This may be telling me simply that there is no relationship between intercepts and slopes and that the correlated model is unjustified, but I would appreciate any advice on the model (listed below) which I quite likely have coded incorrectly. If anyone would like simulated data to go with this code I can happily provide it.
> 
> As an aside, what are the likely effects of fiddling with the the degrees of freedom (set to 3 here)? Is this akin to the 'nu' parameter in MCMCglmm?
> 
> Any help is greatly appreciated
> 
> Regards
> 
> Xavier
> 
> #########################################
> #Separate Intercepts and Slopes - CORRELATED Scaled inverse-Wishart Distribution AND Group-Level Predictor
> #########################################
> 
> cfgroupmod<-"model{
> 
> #Individual model
> 
> for (i in 1:n.samp){
> 	logit.p[i]<-a[yearcode[i]] + b[yearcode[i]]*mass[i] 
> 	p[i]<-exp(logit.p[i])/(1+exp(logit.p[i]))
> 	breed[i]~dbern(p[i])
> 	}
> 	
> 	
> #Group Model 
> 
> for (j in 1:5){
> 	a[j]<-xi.a*B.raw[j,1]
> 	b[j]<-xi.b*B.raw[j,2]
> 	B.raw[j,1:2]~dmnorm(B.raw.hat[j,],Tau.B.raw[,])
> 	B.raw.hat[j,1]<-g.a.0.raw+g.a.1.raw*naoj[j]
> 	B.raw.hat[j,2]<-g.b.0.raw+g.b.1.raw*naoj[j]
> 		}
> 		
> 	#Recover Parameters on Original Scale
> 	g.a.0<-xi.a*g.a.0.raw
> 	g.a.1<-xi.a*g.a.1.raw
> 	g.b.0<-xi.b*g.b.0.raw
> 	g.b.1<-xi.b*g.b.1.raw
> 	
> 	#Priors for scaled Parameters
> 	g.a.0.raw~dnorm(0,0.0001)
> 	g.a.1.raw~dnorm(0,0.0001)
> 	g.b.0.raw~dnorm(0,0.0001)
> 	g.b.1.raw~dnorm(0,0.0001)
> 	
> 	xi.a~dunif(0,100)	
> 	xi.b~dunif(0,100)
> 	
> 	Tau.B.raw[1:2,1:2]~dwish(W[,],df)
> 	df<-3
> 	Sigma.B.raw[1:2,1:2]<-inverse(Tau.B.raw[,])
> 	sigma.a<-xi.a*sqrt(Sigma.B.raw[1,1])
> 	sigma.b<-xi.b*sqrt(Sigma.B.raw[2,2])
> 	rho<-Sigma.B.raw[1,2]/sqrt(Sigma.B.raw[1,1]*Sigma.B.raw[2,2])
> 
> }"
> 	
> W<-diag(2)
> 
> cfgroupdat<-dump.format(list(n.samp=nrow(fem),n.year=5,breed=fem$futurebreed,mass=fem$seasonbci,yearcode=fem$yearcode,naoj=naoframe$junnaoresids,W=W))
> 
> cfgroupparams<-c("a","b","g.a.0","g.a.1","g.b.0","g.b.1","sigma.a","sigma.b","rho")
> burncfgroup=20000
> sampcfgroup=200000
> adaptcfgroup=10000
> cfgroupthin=200
> 
> 
> 
> ---------------------------------------------------------
> Dr Xavier Harrison
> BBSRC Postdoctoral Research Associate
> Centre for Ecology and Conservation
> University of Exeter
> Tremough Campus
> Penryn
> Cornwall
> TR10 9EZ
> Tel: (01326) 371872
> xav.harrison at gmail.com
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.