[R-sig-ME] my first random effects logistic regression

[Jarrod Hadfield]

Hi Dieter,

I was writing this as Thierry posted. My recommendation is similar,  
although you wont be able to fit  farm:flock in glmer because it will  
complain about the number of random effects being equal to the number  
of data. Anyhow....

I'm not sure if I have understood the sampling design correctly, but  
it seems there are 9 flocks over 3 farms, 3 in each farm? If this is  
the case you are going to run into problems with (flock|farm). This  
specification is trying to estimate a 9X9 matrix which is the between  
farm variance for each flock along the diagonals, and the between farm  
covariances between flocks in the off-diagonals.  This is 45  
(co)variance parameters from 9 data points. I think I would try  
something simpler:

model2<-glm(cbind(positives, broilers-positives)~farm,  
data=broilers.dat, family=quasibinomial)

The quasi bit is capturing the flock effects, and the farm effects are  
fitted as fixed. Another way to do it is to expand the binomial  
response into a series of binary data and fit flock as a random effect:


broilers.dat2<-broilers.dat[rep(1:9, broilers.dat$broilers),]
broilers.dat2$positives<- 
(broilers.dat2$positives>=unlist(sapply(broilers.dat$broilers,  
function(x){1:x})))

tapply(broilers.dat2$positives, broilers.dat2$flock, sum) # check its  
OK - it is

model3<-glmer(positives~farm+(1|flock), data=broilers.dat2,  
family=binomial)


you could try farm as well, but with only 3 I would be careful

model4<-glmer(positives~(1|farm)+(1|flock), data=broilers.dat2,  
family=binomial)


This model:

model5<-glmer(positives~(1|flock), data=broilers.dat2, family=binomial)

is equivalent to the one posted by Thierry , but it will run.

Cheers,

Jarrod


On 18 Aug 2010, at 08:28, dieter.anseeuw wrote:

> Hi all,
>
>
>
> I am new to using lme4 and only just discovered the mailing list (I  
> already shortly introduced my problem in the epi-mailinglist, sorry  
> for cross-posting, but my questions seem more appropriate for the  
> mixed models discussion group).
>
>
>
> A friend has inspected three randomly chosen farms (random factor  
> 'farm'). At each farm three randomly chosen series of chickens  
> (random factor 'flock' nested within 'farm') were each inspected for  
> the presence of a certain bacteria. The contaminated chickens were  
> counted (response variable 'positives'). The sample sizes per flock  
> are given by the variable 'broilers'. We want to have a look at  
> within-broilers, within-farm and between-farm variability.
>
>
>
> This is the closest I get to analysing her data:
>
>
>
>> broilers.dat<- 
>> data.frame(farm=c("FA","FA","FA","FB","FB","FB","FC","FC","FC"),  
>> flock=c("a","b","c","d","e","f","g","h","i"), broilers=c(50,  
>> rep(25,8)), positives=c(7,2,0,7,2,0,0,0,2))
>
>> library(lme4)
>
>> model1<-glmer(positives~1 + (flock|farm), data=broilers.dat,  
>> family=binomial(link="logit"), weights=broilers)
>
> Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1
>
>
>
> Hence, my code doesn't work. Could anybody help me out where and why  
> I go wrong?
>
>
>
> Many thanks in advance,
>
> Dieter
>
>
>
> --
>
> Dr. Ir. Dieter Anseeuw
>
> Katho Campus Roeselare
>
> Wilgenstraat 32
>
> 8800 Roeselare Belgium
>
>
>
> Direct phone: +32 51 23 29 68
>
> http://www.katho.be/hivb
>
> http://www.linkedin.com/in/dieteranseeuw
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.