[R] help with winbugs glm

Fri Feb 24 18:02:19 CET 2012

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.
>>> >
>>> >
>>
>>
>