Ohh! about the reference code for the categorical predictor, from the McCarthy book, Bayesian Methods for Ecology (ISBN978-0-521-85057-5) the following ANCOVA model has a Refuge as a categorical predictor with 3 levels coded as "1", "2" and "3".
This is the cited model: model { a ~ dnorm(1.01, 24.826) # intercept b[1] ~ dnorm(0.628, 56.70) # effect of density at resource level 1 (low) b[2] ~ dnorm(0.488, 110.83) # effect of density at resource level 2 (medium) b[3] ~ dnorm(0.180, 70.144) # effect of density at resource level 3 (high) intnS[1] ~ dnorm(0, 1.0E-6) # interaction terms for small plots at 3 refuge densities intnS[2] ~ dnorm(0, 1.0E-6) intnS[3] ~ dnorm(0, 1.0E-6) intnL[1] ~ dnorm(0, 1.0E-6) # interaction terms for large plots at 3 refuge densities intnL[2] ~ dnorm(0, 1.0E-6) intnL[3] ~ dnorm(0, 1.0E-6) prec ~ dgamma(0.001, 0.001) # precision for (i in 1:25) # for each of the 25 plots { pred[i] <- a + b[Refuge[i]]*Density[i] + intnS[Refuge[i]]*Small[i]*Density[i] + intnL[Refuge[i]]*Large[i]*Density[i] Mortality[i] ~ dnorm(pred[i], prec) } } Initial values list(a=0, b=c(0,0,0), intnS=c(0,0,0), intnL=c(0,0,0), prec=1) Data Small[] Large[] Density[] Refuge[] Mortality[] 0 1 0.3125 1 1.6 1 0 0.75 1 0 0 1 0.75 1 2.22 1 0 1 1 1.19 1 0 1.25 1 5.19 0 0 1.5 1 1.73 0 0 1.9375 1 5.15 1 0 1 2 0.83 0 1 1.1875 2 1.44 0 0 1.375 2 1.5 0 0 1.4375 2 1.73 1 0 1.75 2 0 1 0 1.75 2 2.5 1 0 3.5 2 2.82 0 0 4 2 1.81 1 0 1.75 3 1.5 0 1 1.84375 3 1.33 1 0 2 3 1.66 1 0 2.25 3 0.73 0 0 2.5 3 1.6 0 1 3 3 1.99 1 0 3.25 3 1.09 1 0 4.25 3 1.54 1 0 5 3 1.6 1 0 5.5 3 1.73 END On Fri, Feb 24, 2012 at 9:28 AM, Adan Jordan-Garza < ajordangarza2...@my.fit.edu> wrote: > Ok this makes a lot of sense, thank you very much Ilai! > Cheers > Guillermo > > On Fri, Feb 24, 2012 at 12:16 AM, ilai <ke...@math.montana.edu> wrote: > >> On Thu, Feb 23, 2012 at 8:32 PM, Adan Jordan-Garza >> <ajordangarza2...@my.fit.edu> wrote: >> > Hello Ilai, >> > thank you very much for your response, >> > can I bother you a little further? >> > What do you mean my model is not identifiable? >> >> You should read up on model identifiably or consult your local >> statistician for a complete explanation, especially if you are to be >> doing more analysis. An incomplete explanation just regarding your >> bugs code: you are estimating "a" (an overall mean) and the three >> levels b[1:3]. That's too much, you need to remove "a". Think what you >> are asking of bugs to estimate (I'll make up some numbers): for bugs, >> the solution >> a = 10 >> b[1] = 0 >> b[2] = -2 >> b[3] = -5 >> is no different than >> a = 8 >> b[1] = 2 >> b[2] = 0 >> b[3] = -3 >> or >> a = 5 >> b[1] = 5 >> b[2] = 3 >> b[3] = 0 >> >> This is why you're getting large intervals - even with convergence >> etc. the MCMC is simply skipping between these options with no way to >> "identify" which is it since it knows the means for depths 1:3, but >> nothing about "a". a could be 100 and all others -99,-95 etc. >> Makes sense? >> >> It is, for e.g., the reason why summary(glm(y~Depth)) gave you >> intercept (in R the first level) and estimates for the DIFFERENCES of >> all other factor levels from the first. It was not just to make it >> harder on you, it's because the same idea holds for the frequentist >> approach. >> >> It runs ok (10000 >> > iteractions) no autocorrelation, and it does converge but the credible >> > intervals are too big. Yes Depth has 3 levels and I coded them as 1, 2, >> and >> > 3. In all the examples I have seen for winBugs they use this type of >> coding >> > for categorical predictors. >> >> Like I said, winBugs may have some feature where it will build the >> design matrix internally, so you can pass categorical variables as >> column vectors (like R) but I seriously doubt it can guess your >> variable coded 1,2,3 is categorical. If you can reference an example >> like that from winBugs examples (not some random code from the >> Internet) I'd really appreciate it. >> What I think is happening is you're fitting >> y1 = 1b[1] + e1 >> y2 = 3b[3] + e2 >> ... >> y232 = ... + e232 >> >> But why should b[3] be 3 times b[1] ? that makes no sense for >> categories. And the estimates are completely off (as you've noted in >> your post) not just large variance. >> >> I have read about the matrix of dummy variables >> > but I am not sure if I am supposed to code that explicitly in winBugs >> or how >> > to do it. I already got the dummy matrix from R. >> >> There are R2winbugs and other packages to let you run winbugs from R >> or pass the data. I don't use them so I don't know. At the end what >> you want is Depth to be a 3 column matrix, either the same as glm to >> get estimate of Depth1,Depth2-1,Depth3-1 : >> 1 0 0 >> 1 0 0 >> 1 1 0 >> 1 1 0 >> 1 0 1 >> 1 0 1 >> >> Or you could have an estimate for every depth: >> 1 0 0 >> 1 0 0 >> 0 1 0 >> 0 1 0 >> 0 0 1 >> 0 0 1 >> >> And your code: >> model >> { >> for (i in 1:232) >> { >> Recruitment[i]~dpois(lambda[i]) >> log(lambda[i])<- b[1]*Depth[,1]+b[2]*Depth[,2]+b[3]*Depth[,3] >> } >> b[1]~dnorm(0,0.01) >> b[2]~dnorm(0,0.01) >> b[3]~dnorm(0,0.01) >> } >> >> Now compare the result to glm. >> >> > from R. Do you recommend JAGS over winBugs? >> >> I don't use winBugs because I'm not on a Win OS. I don't think JAGS >> has an advantage for simple models like this. >> >> Good luck, and please cc the r-help list in your answer so the thread >> is complete. >> >> Elai. >> >> >> > Thank you very much for your time, I really appreciate it. >> > >> > This is my whole code in winBugs: >> > >> > model >> > { >> > for (i in 1:232) >> > { >> > Recruitment[i]~dpois(lambda[i]) >> > log(lambda[i])<-a+b[Depth[i]]*Depth[i] >> > } >> > a~dnorm(0,0.01) >> > b[1]~dnorm(0,0.01) >> > b[2]~dnorm(0,0.01) >> > b[3]~dnorm(0,0.01) >> > } >> > list(a=0, b=c(0,0,0)) >> > >> > Depth[] Recruitment[] >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 1 >> > >> > 1 0 >> > >> > 1 1 >> > >> > 1 1 >> > >> > 1 1 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 1 >> > >> > 1 1 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 8 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 2 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 1 >> > >> > 1 2 >> > >> > 1 1 >> > >> > 1 0 >> > >> > 1 1 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 0 >> > >> > 1 8 >> > >> > 2 2 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 4 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 2 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 2 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 0 >> > >> > 2 1 >> > >> > 2 1 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 5 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 2 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 1 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 2 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 2 >> > >> > 3 3 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 3 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 1 >> > >> > 3 2 >> > >> > 3 4 >> > >> > 3 0 >> > >> > 3 0 >> > >> > 3 3 >> > >> > 3 1 >> > >> > 3 0 >> > >> > 3 0 >> > END >> > >> > >> > >> > >> > On Thu, Feb 23, 2012 at 1:19 PM, ilai <ke...@math.montana.edu> wrote: >> >> >> >> Adan, >> >> How many levels does Depth have? my wild guess: 3 and your bugs model >> >> is not identifiable. >> >> >> >> Second, I think you may have a critical error in the way you formatted >> >> the data for the bugs model. From your code it looks like you are just >> >> using the factor Depth and not a design matrix of dummy variables. >> >> >> >> I may be wrong with respect to WinBugs (I use JAGS), but if Depth is >> >> denoted as, for e.g., "low","med","high" wouldn't your multiply >> >> operation "...*Depth[i] " on line 5 fail ? >> >> >> >> More likely Depth is denoted "1","2","3" and WinBugs thinks it's >> >> numerical. Well, in that case clearly coefficients for this model >> >> don't make any sense (you'll only need one b for the slope). You can >> >> use model.matrix(~Depth) to get the proper format for your data. >> >> >> >> Hope this solves it. Next time, knowing n.chains n.iter and if they >> >> achieved convergence (with different starting values) can help sort >> >> through these sort of issues. >> >> >> >> Cheers >> >> >> >> Elai >> >> >> >> On Thu, Feb 23, 2012 at 9:57 AM, Adan Jordan-Garza >> >> <ajordangarza2...@my.fit.edu> wrote: >> >> > Hi, >> >> > I am running a model with count data and one categorical predictor >> >> > (simple >> >> > model for me to understand it fully), I did in R a glm like this: >> >> > glm(Recruitment~Depth, family=poisson). I get the coefficientes and >> >> > confidence intervals and all is ok. But then I want to do the same >> model >> >> > with Bayesian stats, here is my code: >> >> > >> >> > model >> >> > { for (i in 1:232) >> >> > { >> >> > Recruitment[i]~dpois(lambda[i]) >> >> > log(lambda[i])<-a+b[Depth[i]]*Depth[i] >> >> > } >> >> > a~dnorm(0,0.000001) >> >> > b[1]~dnorm(0,0.000001) >> >> > b[2]~dnorm(0,0.000001) >> >> > b[3]~dnorm(0,0.0000001) >> >> > } >> >> > list(a=0, b=c(0,0,0)) >> >> > >> >> > I have two problems: 1) the resulting credible intervals for the >> >> > coefficients (a, b1, b2 and b3) are HUGE don t make any reasonable >> >> > sense; >> >> > 2) Using OpenBugs and Winbugs I get different results, >> >> > >> >> > if anyone can help me I appreciate a lot your time, >> >> > >> >> > thanks >> >> > >> >> > Guillermo >> >> > >> >> > [[alternative HTML version deleted]] >> >> > >> >> > ______________________________________________ >> >> > R-help@r-project.org mailing list >> >> > https://stat.ethz.ch/mailman/listinfo/r-help >> >> > PLEASE do read the posting guide >> >> > http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html> >> >> > and provide commented, minimal, self-contained, reproducible code. >> > >> > >> > > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.