[R] glmmBUGS: logistic regression on proportional data

David Winsemius dwinsemius at comcast.net
Sun Feb 8 21:36:09 CET 2009


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




More information about the R-help mailing list