[R-sig-ME] MCMCglmm variance estimates Poisson distribution

Wed Mar 5 11:13:57 CET 2014

Hi Jarrod

Ah ok, thanks so if I run 

re.uc11.cf <- glmer(totflct ~ (1|iobs), family=poisson, data=uc11)
summary(re.uc11.cf)

with a term for each observation I get

Generalized linear mixed model fit by the Laplace approximation 
Formula: totflct ~ (1 | iobs) 
   Data: uc11 
  AIC  BIC logLik deviance
 3960 3971  -1978     3956
Random effects:
 Groups Name        Variance Std.Dev.
 iobs   (Intercept) 0.69221  0.83199 
Number of obs: 1607, groups: iobs, 1607

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.2133     0.0249   48.74   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

which corresponds to the estimates from the MCMC model. 

My problem is now that the estimate it doesn't match the data (mean 4.69), but at least the models are in agreement.

Thanks
John

-----Original Message-----
From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk] 
Sent: 05 March 2014 08:18
To: Hodsoll, John
Cc: r-sig-mixed-models at r-project.org
Subject: RE: [R-sig-ME] MCMCglmm variance estimates Poisson distribution

Hi John,

The two models are different. The obs variable in the lmer model does not have a unique level for each observation:

obs: 1607, groups: obs, 168

Cheers,

Jarrod

Quoting "Hodsoll, John" <john.hodsoll at kcl.ac.uk> on Wed, 5 Mar 2014
07:18:39 +0000:

> Hi Jarrod
>  >  is there a reason that the data frames differ in each (uc11 and uc12)?
> Yes. Two different baseline conditions, cut and paste error.
>
> With the same data frame
>
> summary(mcmc.c11.cf2)
>  Iterations = 10001:99911
>  Thinning interval  = 90
>  Sample size  = 1000
>
>  DIC: 7489.396
>
>  R-structure:  ~units
>
>       post.mean l-95% CI u-95% CI eff.samp
> units    0.7111   0.6276   0.7912     1000
>
>  Location effects: totflct ~ 1
>
>             post.mean l-95% CI u-95% CI eff.samp  pMCMC
> (Intercept)     1.211    1.160    1.263     1000 <0.001 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>> posterior.mode(exp(mcmc.c11.cf$Sol+mcmc.c11.cf$VCV/2))
> (Intercept)
>      4.7921
>
>> exp(mean(mcmc.c11.cf2$Sol)+mean(mcmc.c11.cf2$VCV/2))
> [1] 4.793908
>
> And for glmer...
> Generalized linear mixed model fit by the Laplace approximation
> Formula: totflct ~ (1 | obs)
>    Data: uc11
>   AIC  BIC logLik deviance
>  6205 6216  -3100     6201
> Random effects:
>  Groups Name        Variance Std.Dev.
>  obs    (Intercept) 0.083672 0.28926
> Number of obs: 1607, groups: obs, 168
>
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept)  1.50622    0.02535   59.43   <2e-16 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 __
>
> for which exp(1.50622 + 0.083672/2) = 4.70232
>
> I take your point re prior. Number of observations is 1607 so I 
> thought this should be sufficient to limit the influence of the prior?
>
> Cheers
> John
> ______________________________________
> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
> Sent: 04 March 2014 17:13
> To: Hodsoll, John
> Cc: 'r-sig-mixed-models at r-project.org'
> Subject: RE: [R-sig-ME] MCMCglmm variance estimates Poisson 
> distribution
>
> Hi John,
>
> Perhaps the output from:
>
> summary(mcmc.c11.cf2)
>
> and
>
> summary(re.uc12.cf)
>
> will shed some light?  Also is there a reason that the data frames 
> differ in each (uc11 and uc12)?
>
> Failing something `obvious' then it must be the prior. How many 
> observations is this based on?
>
> Cheers,
>
> Jarrod
>
>
> Quoting "Hodsoll, John" <john.hodsoll at kcl.ac.uk> on Tue, 4 Mar 2014
> 16:40:11 +0000:
>
>> Dear Jarrod
>>
>> Thanks for your reply. I was using 4, but all give a similar answer
>>
>> 1. 4.796014
>>
>> 2. 4.792395
>>
>> 3. 4.798754
>>
>> 4  4.790677
>>
>> As I said I'm not sure I'm missing something obvious?
>>
>> Cheers
>> John
>>
>>
>> -----Original Message-----
>> From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
>> Sent: 04 March 2014 12:07
>> To: Hodsoll, John
>> Cc: 'r-sig-mixed-models at r-project.org'
>> Subject: Re: [R-sig-ME] MCMCglmm variance estimates Poisson 
>> distribution
>>
>> Dear John,
>>
>> How are you calculating the posterior expectation:
>>
>> 1/
>> posterior.mode(exp(mcmc.c11.cf2$Sol+mcmc.c11.cf2$VCV/2))
>> 2/
>> mean(exp(mcmc.c11.cf2$Sol+mcmc.c11.cf2$VCV/2))
>> 3/
>> exp(posterior.mode(mcmc.c11.cf2$Sol)+posterior.mode(mcmc.c11.cf2$VCV/
>> 2))
>> 4/
>> exp(mean(mcmc.c11.cf2$Sol)+mean(mcmc.c11.cf2$VCV/2))
>>
>> If it is not by method 1/ try that and see if there is less of a 
>> discrepancy.
>>
>> Cheers,
>>
>> Jarrod
>>
>> Quoting "Hodsoll, John" <john.hodsoll at kcl.ac.uk> on Tue, 4 Mar 2014
>> 11:25:17 +0000:
>>
>>> Dear list
>>>
>>> I'm trying to get some unadjusted estimates and 95% CI for a set of 
>>> correlated count data (due to repeated measures on the same cluster) .
>>> To do this I was trying to run an over-dispersed poisson model using 
>>> a glmer and MCMCglmm.
>>>
>>> I want to use MCMCglmm as that's the package I wish to use for my 
>>> main analysis. However, it seems to over-estimate the variance 
>>> meaning that the mean value I get from the intercept only model y = 
>>> XB + Var/2 (ch2 jarrod hadfield's course notes) is slightly greater 
>>> than the actual mean. For example, if I fit the model
>>>
>>> priortr <- list(R=list(V=1, nu=0.001))
>>>
>>> mcmc.c11.cf2 <- MCMCglmm(totflct ~ 1, family="poisson", prior = 
>>> priortr, data=uc11,
>>>                          nitt = 100000, burnin = 10000, thin = 90)
>>> summary(mcmc.c11.cf2)
>>>
>>> I'm ignoring the random effect and assuming the additive 
>>> over-dispersion term will capture all the extra variance. For a 
>>> count rate of 4.69 in the data I get 4.79 and for a count of 5.2 I get 5.52.
>>> On the other hand, if I use glmer including a per observation random 
>>> effect I get the correct means
>>>
>>> re.uc12.cf <- glmer(totflct ~ (1|obs), family=poisson, data=uc12)
>>> summary(re.uc12.cf)
>>>
>>> Is there something I missing here?
>>>
>>> Regards
>>> John Hodsoll
>>>
>>> _______________________________________________
>>> 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.
>>
>>
>>
>>
>
>
>
> --
> The University of Edinburgh is a charitable body, registered in 
> Scotland, with registration number SC005336.
>
>
>

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