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

[Joerg Luedicke]

Hi Ben and Fredrik,

Thanks for the comments and the proof! Next time I will be less hasty
in responding...

Cheers,

Joerg

On Mon, Mar 19, 2012 at 3:09 AM, Fredrik Nilsson <laf.nilsson at gmail.com> wrote:
> 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