R help - Feb 2009 - glmmBUGS: logistic regression on proportional data

Hello,

I am trying to run a logistic regression with random effects on 
proportional data in glmmBUGS.  I am a newcomer to this package, and 
wondered if anyone could help me specify the model correctly. 

I am trying to specify the response variable, /yseed/, as # of successes 
out of total observations... but I suspect that given the error below, 
that is not correct.  Also, Newsect should be a factor, whereas Newdist 
is continuous.

Thanks,
John


Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
each=10), Newdist=rep(1:5, 2),
                   y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))

yseed<-cbind(Newdat$y, Newdat$tot)

mod<-glmmBUGS(yseed~Newsect + Newdist, effects="Newtree", 
family="binomial", data=Newdat)

Error in `[.data.frame`(data, , response, drop = FALSE) :
  undefined columns selected

John Poulsen <jpoulsen <at> zoo.ufl.edu> writes:
> I am trying to run a logistic regression with random effects on 
> proportional data in glmmBUGS.  I am a newcomer to this package, and 
> wondered if anyone could help me specify the model correctly. 
> 
> I am trying to specify the response variable, /yseed/, as # of successes 
> out of total observations... but I suspect that given the error below, 
> that is not correct.  Also, Newsect should be a factor, whereas Newdist 
> is continuous.
> 
> Thanks,
> John
> 
> Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
> each=10), Newdist=rep(1:5, 2),
>                    y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))
> 
> yseed<-cbind(Newdat$y, Newdat$tot)
> 
> mod<-glmmBUGS(yseed~Newsect + Newdist, effects="Newtree", 
> family="binomial", data=Newdat)
> 
First, a typo, there is no yseed. Second, after the error message
"must be between 0 and 1", this looks more like poisson, because
you have the counts, not the events. 

This might come close

mod<-glmmBUGS(y~Newsect + Newdist, effects="Newtree", 
  family="poisson", data=Newdat)


Dieter

On Feb 8, 2009, at 10:46 AM, Dieter Menne wrote:
> John Poulsen <jpoulsen <at> zoo.ufl.edu> writes:
>
>> I am trying to run a logistic regression with random effects on
>> proportional data in glmmBUGS.  I am a newcomer to this package, and
>> wondered if anyone could help me specify the model correctly.
>>
>> I am trying to specify the response variable, /yseed/, as # of  
>> successes
>> out of total observations... but I suspect that given the error  
>> below,
>> that is not correct.  Also, Newsect should be a factor, whereas  
>> Newdist
>> is continuous.
>>
>> Thanks,
>> John
>>
>> Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
>> each=10), Newdist=rep(1:5, 2),
>>                   y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))
>>
>> yseed<-cbind(Newdat$y, Newdat$tot)
>>
>> mod<-glmmBUGS(yseed~Newsect + Newdist, effects="Newtree",
>> family="binomial", data=Newdat)
>>
>
> First, a typo, there is no yseed. Second, after the error message
> "must be between 0 and 1", this looks more like poisson, because
> you have the counts, not the events.
Puzzled. I see yseed defined above as a two column vector, as is  
sometimes used to handle grouped data input to the glm response side  
of a formula.
>
>
> This might come close
>
> mod<-glmmBUGS(y~Newsect + Newdist, effects="Newtree",
>  family="poisson", data=Newdat)
Reasoning only by analogy from the experience with ordinary glm()  
input to create a Poisson model and having no experience with glmmBUGS:
  How you are accounting for the tot (presumably totals) from which it  
appears the y variable is being considered as forming a proportion?  
Would have expected to see an offset=log(tot) or perhaps a weights=tot  
in that call.


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

Installing the glmmBUGS package and a bit of experimentation produces  
this minor modification of your code that seems to run without error.  
It's back to you to see if the output is sensible.

library(glmmBUGS)
Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
each=10), Newdist=rep(1:5, 2),
                  y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))

yseed<-cbind(Newdat$y, Newdat$tot)

mod<-glmmBUGS(y/tot~Newsect + Newdist, effects="Newtree",
family="binomial", data=Newdat)

I am guessing that the interpreter was looking for a variable yseed  
within Newdat and not finding it. It, however, is able to "see" y and
tot within Newdat. Also appears that using cbind(y, tot) on the LHS of  
hte formula will avoid the error you were getting.

-- 
David Winsemius
Heritage Labs

On Feb 8, 2009, at 8:29 AM, John Poulsen wrote:
> Hello,
>
> I am trying to run a logistic regression with random effects on  
> proportional data in glmmBUGS.  I am a newcomer to this package, and  
> wondered if anyone could help me specify the model correctly.
> I am trying to specify the response variable, /yseed/, as # of  
> successes out of total observations... but I suspect that given the  
> error below, that is not correct.  Also, Newsect should be a factor,  
> whereas Newdist is continuous.
>
> Thanks,
> John
>
>
> Newdat<-data.frame(Newtree=rep(1:3, each=20),  
> Newsect=rep(c("a","b"), each=10), Newdist=rep(1:5, 2),
>                  y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))
>
> yseed<-cbind(Newdat$y, Newdat$tot)
>
> mod<-glmmBUGS(yseed~Newsect + Newdist, effects="Newtree",  
> family="binomial", data=Newdat)
>
> Error in `[.data.frame`(data, , response, drop = FALSE) :
> undefined columns selected
>
> ______________________________________________
> 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.

David Winsemius <dwinsemius <at> comcast.net> writes:
> >>
> >> Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
> >> each=10), Newdist=rep(1:5, 2),
> >>                   y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))
> >>
> >> yseed<-cbind(Newdat$y, Newdat$tot)
> >>
> >> mod<-glmmBUGS(yseed~Newsect + Newdist,
effects="Newtree",
> >> family="binomial", data=Newdat)
> >>
> >
> > First, a typo, there is no yseed. Second, after the error message
> > "must be between 0 and 1", this looks more like poisson,
because
> > you have the counts, not the events.
> 
> Puzzled. I see yseed defined above as a two column vector, as is  
> sometimes used to handle grouped data input to the glm response side  
> of a formula.
Sorry, my short-circuit. You example is in analogy to the budworm example
in MASS, giving successes and failures. I have no solution, but a quick
look shows that the formula interface of glmmPQL is called, so you
might cross check if glmmPQL support the success/failure form. I am not
on an R machine currently to check it.

Dieter

Just make sure it IS help. Following Dieter's advice I checked to see  
what format V&R use to pass data as success/failure and now wonder if  
the the proper invocation following their examples on pg 49, 59 and  
134 in <http://cran.r-project.org/web/packages/VR/VR.pdf> would be:

mod<-glmmBUGS( cbind(y, not.y = tot - y) ~ Newsect + Newdist,  
effects="Newtree", family="binomial", data=Newdat)

That invocation runs without complaint, but this two step one does not:

 > Newdat$not.y <- with(Newdat, tot - y)
 > mod<-glmmBUGS(cbind(y, not.y) ~ Newsect + Newdist,  
effects="Newtree", family="binomial", data=Newdat)
Error in family$linkfun(mustart) : Value 0 out of range (0, 1)

It worries me when two submission that I would expect to be equivalent  
produce different results. It's not that the failure vector has zeros,  
but I cannot immediately see the source of the problem:

 > Newdat$not.y
  [1] 13 11  7  5  4 13 10  6  6  5  9 12  9  5  5 10 10  8  6  5 12  
10  8  6  5 14  7  8  5  4 11 11  9  8  3
[36] 10  7  9  8  1 13 11  9  6  6 11 10  7  7  2 13 10  9  5  6 10   
8  8  7  2

You probably need either to consult a reference that is more  
authoritative or to construct an example for which you are certain of  
the correct result with your chosen package.

-- 
David Winsemius
Heritage Labs


On Feb 8, 2009, at 12:48 PM, John Poulsen wrote:
> Thanks for your help! -- John
>
> David Winsemius wrote:
>> Installing the glmmBUGS package and a bit of experimentation  
>> produces this minor modification of your code that seems to run  
>> without error. It's back to you to see if the output is sensible.
>>
>> library(glmmBUGS)
>> Newdat<-data.frame(Newtree=rep(1:3, each=20),
Newsect=rep(c("a","b"),
>> each=10), Newdist=rep(1:5, 2),
>>                 y=rpois(60,2), tot=rep(c(14,12,10,8,6), 12))
>>
>> yseed<-cbind(Newdat$y, Newdat$tot)
>>
>> mod<-glmmBUGS(y/tot~Newsect + Newdist, effects="Newtree",
>> family="binomial", data=Newdat)
>>
>> I am guessing that the interpreter was looking for a variable yseed  
>> within Newdat and not finding it. It, however, is able to
"see" y
>> and tot within Newdat. Also appears that using cbind(y, tot) on the  
>> LHS of hte formula will avoid the error you were getting.
>>
>