[R-sig-ME] modelling proportions, with aggregated data, and the new/old lme4

Fredrik Nilsson laf.nilsson at gmail.com
Mon Mar 19 11:09:17 CET 2012


Joerg,

What is the point in the Poisson model in that way (assuming
log(total) as offset)? This would be a binomial model; if you observe
N events from a Poisson distribution of which X are successes (or
marked in some way) with probability p, then X will be binomial(p, N).

Best regards,
Fredrik Nilsson

PS. Here's a practical "proof", the real is sketched below.
invlogit<-function(x) exp(x)/(1+exp(x))
n<-10
nn<-ceiling(exp(rnorm(n)+3))
pp<-seq(0.25,.75,length=n)
ca<-rbinom(rep(1,n),nn,prob=pp)
co<-nn-ca
fa<-gl(n,1)

#binomial model
test.glm<-glm(cbind(ca,co)~fa-1, family=binomial)
#conditional poisson
ptest.glm<-glm(ca~fa-1+offset(log(nn)), family=poisson)
summary(ptest.glm)
exp(coef(ptest.glm))
invlogit(coef(test.glm))

all.equal(exp(coef(ptest.glm)), invlogit(coef(test.glm)))

# proof.
X~Poisson(a) and  Y~Poisson(b), X & Y independent -> X+Y~Poisson(a+b)
(e.g. by generating functions).

Prob(X=x| X+Y=n) = Prob(X=x, Y=n-x)/Prob(X+Y=n) =
Prob(X=x)*Prob(Y=n-x)/Prob(X+Y=n) = a^x/x! exp(a) b^(n-x)/(n-x)!
exp(b) /((a+b)^n/n! exp(a+b)) =
n!/x!/(n-x)! (a/(a+b))^x (b/(a+b))^(n-x) = n!/x!/(n-x)! (a/(a+b))^x
(1-a/(a+b))^(n-x) , i.e. binomial(x,n,p) with p=a/(a+b).


2012/3/19 Ben Bolker <bbolker at gmail.com>:
> Joerg Luedicke <joerg.luedicke at ...> writes:
>
>> I would certainly check out a Poisson model with the number of
>> successes as outcome and successes+failures as an offset.
>
>  That seems odd to me; the Poisson+offset model
> should be appropriate when p<<1 (i.e. for very small p,
> the Poisson variance mu=n*p is approximately the same as
> the binomial variance n*p*(1-p); in this case p is not small.
>
>  Of course, I could be wrong.
>
>>
>> On Sun, Mar 18, 2012 at 11:51 AM, Rolf Turner <r.turner at ...> wrote:
>> > On 19/03/12 07:10, Ben Bolker wrote:
>> >
>> > <SNIP>
>> >
>> >
>> >>> head(dat)
>> >>>   successes failures id sex subdist
>> >>> 1       560      726  1   F       4
>> >>> 2       844      510  1   M       4
>
>  [snip]
>
>> >>>
>> >>> In community #1, which is in subdistrict #4, there are 560 women
>> >>
>> >>   you mean 510, right?
>> >
>> > <SNIP>
>> >
>> > Sure looks like 560 to me.  Time for a trek to the optometrist, Ben?
>
>  Looked at the wrong line (failures in men rather than successes in women).
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models




More information about the R-sig-mixed-models mailing list