[R-sig-ME] MCMCglmm help- information about 'units' term

[Jarrod Hadfield]

Dear Karen,

I think the model is specified correctly. The units term is the  
variance of observation-level random effects, or a residual variance  
if you like. MCMCglmm always fits this term because a) I think over- 
dispersed models should be the default and b) the algorithm is  
designed so that the mixing properties of the chain are a function of  
the residual variance. This poses a problem for ordinal, binary and  
single-outcome multinomial responses because observation-level  
heterogeneity cannot be estimated from the data. Given this, I suggest  
that the residual (units) variance is fixed at some value. The default  
in most (all?) other packages is to fix the residual variance at zero,  
but MCMCglmm will not mix under this assumption.  In fact the mixing  
properties improve as the units variance is increased (see van Dyck  
and Meng's The art of data augmentation) although at some point  
underflow/overflow problems start to occur. For this reason I suggest  
fixing the residual variance to one in such models, as you have done.   
The interpretation of the parameters are identical to those under a  
probit model without a residual variance although the inverse link is  
now pnorm(x, 0, sqrt(2)) rather than  pnorm(x, 0, sqrt(1)).

Very small values for the ID variance seem to have some support.  The  
mixing of the ID variance can be poor when the value is close to zero  
(it can get stuck at zero).  For reasons very similar to those above  
(again see van Dyck and Meng's paper) mixing can be improved by adding  
a redundant non-identified parameter. Parameter expansion can be  
specified in MCMCglmm through the prior. For example, you may find  
that a change from

prior=list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=0)))

to

prior=list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=1, alpha.mu=0,  
alpha.V=100)))

hardly changes the posterior distribution but the autocorrelation in  
the chain is vastly improved. Most of this is documented in the  
CourseNotes and I recommend Gelman's paper "Prior distributions for  
variance parameters in hierarchical models" in Bayesian Analysis for a  
good introduction into these prior distributions.

Cheers,

Jarrod

On 30 Jun 2010, at 12:12, Karen Lamb wrote:

> Hi all,
>
> I have been working with mixed effects models for a couple of years  
> but I am relatively new to MCMCglmm() and MCMC techniques in general  
> so I hope someone may be able to shed some light on an issue I have.
>
> I am currently trying to fit a 3 level ordinal multinomial mixed  
> model. To begin, I fitted a very simplistic model to try out the  
> approach with only an intercept term and the random effect of ID  
> using the following code:
>
> prior=list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=0)))
>
> m1<-MCMCglmm(newbmi~1, random=~ID, family="ordinal", data=data1,  
> prior=prior)
>
>
>
>> From assessment of  autocorr(m1$Sol)it appears that convergence is  
>> ok.
>
>
>
> The issue I have concerns the random effects.  If I assess anything  
> involving m1$VCV l obtain really strange results. For example,
>
>
>
>> HPDinterval(m1$VCV[, "ID"]/(m1$VCV[,"ID"]+m1$VCV[,"units"]))
>            lower     upper
> var1 1.269593e-11 0.1352029
> attr(,"Probability")
> [1] 0.95
>
>
>> cor(m1$VCV)
>      ID units
> ID     1    NA
> units NA     1
> Warning message:
> In cor(m1$VCV) : the standard deviation is zero
>
>
>> diag(autocorr(m1$VCV)[2,,])
>           ID         units
> -0.0006257816           NaN
>
>
> Can anyone explain what this output means? My ID effect is tiny and  
> I am not sure that it is necessary to have the random effect in the  
> model after all. However, I intend to fit ID level explanatory  
> variables in the model so wish to retain this random effect. I just  
> don't understand what the problem is with the units term. Am I  
> specifying the model incorrectly?
>
>
>
> Any help/suggestions would be greatly appreciated!
>
>
>
> Cheers,
> Karen
>
>
>
> 	[[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.