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.
