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@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@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.



	[[alternative HTML version deleted]]


