help with winbugs glm
My apologies to R-list members, this discussion diverged and is no longer R related. Probably more fitting in one of the BUGS forums, but I keep cc ing for any future interested reader that stumbled upon this post in r-help. Now to the point. Thank you so much for the reference, I was not aware of this book. I will check it out as it seems can be useful for me too. However, compare McArthy's model to your own: He is using Refuge[i] i=1,2,3 BUT ONLY AS A PLACE HOLDER for b[1:3] !!! "... b[Refuge[i]]*Density[i] " where density is numeric - this is as you said ANCOVA - fitting different slopes to Density BY levels of Refuge. Not the same thing as your "... b[Depth[i]] * Depth[i] " Which is senseless as already discussed. Side note, he can enter an explicit intercept "a" because it is the estimate for the Med level, in the data he has dummy vars only for Small and Large (e.g. in the first two columns of his design matrix, lines 6 and 7 are 0 0 ). Thus, HIS model is identifiable. Brad Efron (1986) said "Bayesian theory requires a great deal of thought about the given situation to apply sensibly". And maybe this is a good time to revisit another statement I think attributed to him: "Recommending that scientists use Bayes? theorem is like giving the neighborhood kids the key to your F-16. " Present company excluded :) as you seem to be more careful in your application than some. Cheers, Elai On Fri, Feb 24, 2012 at 8:46 AM, Adan Jordan-Garza
<ajordangarza2009 at my.fit.edu> wrote:
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
<ajordangarza2009 at 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 <keren at math.montana.edu> wrote:
On Thu, Feb 23, 2012 at 8:32 PM, Adan Jordan-Garza <ajordangarza2009 at 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 <keren at 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 <ajordangarza2009 at 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 at 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.