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

[Jarrod Hadfield]

Dear Karen,

Imagine a latent variable (l) that conforms to the  standard linear  
model.

l = Xb+Zu+e

The probabilities of falling into each of the three categories are:


pnorm(-l)

pnorm(cp-l)-pnorm(-l)

1-pnorm(cp-l)

where cp is the cut-point.  If you have just an intercept, and you  
fixed the residual variance to 1 this gives:

pnorm(-Intercept, sd=sqrt(2))  = 0.589807

pnorm(cp-Intercept, sd=sqrt(2))-pnorm(-Intercept, sd=sqrt(2)) =  
0.2968132

1-pnorm(cp-Intercept, sd=sqrt(2)) = 0.1133798

where sd=sqrt(2) rather than sd=1 because we're marginalising the  
residuals. The numbers should correspond to the frequency of the  
categories in your data.

Cheers,

Jarrod






On 5 Jul 2010, at 12:24, Karen Lamb wrote:

> Thanks very much for your assistance in this matter. I have  
> persevered with the modelling but am now having problems  
> interpreting the output for the fixed effects from the modelling.
>
> I have a three category model (0=low, 1=medium, 2=high)  which I  
> have been treating as ordinal. I have only previously worked with  
> nominal multinomial models and am used to fitting a different  
> intercept estimate for each comparison. I had assumed when I was  
> fitting the model that I would obtain an intercept for the model  
> comparison of medium to low and one for the comparison of high to  
> low. Do I need to add in some further model specification in my code  
> to obtain this or am I simply misunderstanding the model fitting for  
> an ordinal response? As it is currently specified, I only obtain one  
> estimated intercept parameter and I do not know how to interpret  
> this. The results are:
>
> > posterior.mode(m1$Sol)
> cutpoint.traitnewbmi.1            (Intercept)
>              1.3883328             -0.3210951
>
> In addition, I have a number of explanatory variables I would like  
> to include in the multinomial model and am trying to figure out the  
> best approach to use when deciding which of these variables should  
> be retained in the model. Is it best to examine DIC in the model  
> comparison or is there some other way of testing the significance of  
> the predictor? Is there any particular measure of the goodness of  
> fit of the model that is routinely adopted? I have been asked to  
> provide some measure of the proportion of variance in the data  
> explained by the fixed effects included in the model and do not know  
> if this type of thing is possible with glms. I understand that there  
> are some suggested R-square type measures for linear mixed models  
> but am not aware of anything similar for a multinomial.
>
> Once again, any suggestions would be greatly appreciated.
>
> Cheers,
> Karen
>
>
>
> ----- Original Message -----
> From: Jarrod Hadfield
> To: Karen Lamb
> Cc: r-sig-mixed-models at r-project.org
> Sent: Friday, July 02, 2010 10:37 AM
> Subject: Re: [R-sig-ME] MCMCglmm help- information about 'units' term
>
> 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.
>
>

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20100706/27527ecf/attachment.pl>