[R-sig-ME] MCMCglmm: priors for ordinal regression

m.fenati at libero.it m.fenati at libero.it
Fri Jul 20 15:07:02 CEST 2012 

I will try as soon as possible.
I am very interested to include MCMCglmm as possible tool for my analyses 
where ordinal data are very usual.
Thanks again. 
Massimo

>----Messaggio originale----
>Da: j.hadfield a ed.ac.uk
>Data: 20/07/2012 11.32
>A: "m.fenati a libero.it"<m.fenati a libero.it>
>Cc: <r-sig-mixed-models a r-project.org>
>Ogg: Re: [R-sig-ME] MCMCglmm: priors for ordinal regression
>
>Dear Massimo,
>
>The function below (prior.scale) may be useful. It returns a  
>covariance matrix that can be passed to prior$B$V.  If the predictors  
>had been scaled according to the procedure outlined in Gelman et al.  
>(2006) then this would induce an identity prior covariance matrix on  
>the new regression coefficients (i.e. they're iid a priori). Gelman  
>recommends a scaled Cauchy or scaled t prior on these new regression  
>ceofficients. I don't think this is possible in MCMCglmm* but a normal  
>with variance equal to the sum of the variance components + 1 (probit)  
>  pi^2/3 (logit) I think is not too unreasonable - but it does put  
>less weight on very extreme probabilities than the Cauchy and t. To  
>achieve this just multiply the output of prior.scale by the variance  
>you want to use.
>
>Bear in mind that the inputs are scaled, not the columns of the design  
>matrix - this means that interactions are penalised more than main  
>effects. In my field just-so stories involving 'significant'  
>interactions  are often used to resuscitate a `failed' experiment, so  
>penalising them by default may be no bad thing.
>
>Cheers,
>
>Jarrod
>
>*setting a weaker prior on the non-identified residual variance rather  
>than fixing it may achieve a t-prior - but really not sure
>
>Gelman et al. A WEAKLY INFORMATIVE DEFAULT PRIOR DISTRIBUTION FOR
>LOGISTIC AND OTHER REGRESSION MODELS The Annals of Applied Statistics
>2008, Vol. 2, No. 4, 1360–1383
>
># formula = the fixed formula in the model
># data = the data set
>
>prior.scale<-function(formula, data){
>   X1<-model.matrix(formula, data)
>   X2<-get_all_vars(formula, data)
>   X2<-as.data.frame(lapply(X2, function(x){if(is.numeric(x)){scale(x,  
>scale=sd(x)*2)}else{x}}))
>   X2<-model.matrix(formula, data=X2)
>   X2[,-1]<-apply(X2[,-1], 2,  
>function(x){if(any(!x%in%c(0,1))){x}else{scale(x,  
>center=sum(x)/length(x), scale=1)}})
>   crossprod(t(solve(t(X1)%*%X1, t(X1)%*%X2)))
>}
>
>
>Quoting "m.fenati a libero.it" <m.fenati a libero.it> on Mon, 16 Jul 2012  
>11:04:45 +0200 (CEST):
>
>>
>> Thank you again for your great suggestions.
>>
>> Regards
>>
>> Massimo
>>
>>
>>
>>> ----Messaggio originale----
>>> Da: j.hadfield a ed.ac.uk
>>> Data: 12/07/2012 15.20
>>> A: "m.fenati a libero.it"<m.fenati a libero.it>
>>> Cc: <r-sig-mixed-models a r-project.org>
>>> Ogg: Re: [R-sig-ME] MCMCglmm: priors for ordinal regression
>>>
>>> Hi,
>>>
>>> If the prior variance on your fixed effects is V+1 where V is the sum
>>> of the variance components (including the residual) then the marginal
>>> prior on the fixed effects is as flat as possible on the probability
>>> interval (0,1). However, you have to set up the contrasts correctly.
>>>
>>> If you still get numerical problems I'm afraid you will have to find
>>> another way of doing the analysis. I have no solution, and no one has
>>> suggested any:
>>>
>>> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q1/017976.html
>>>
>>> Cheers,
>>>
>>> Jarrod
>>>
>>>
>>>
>>> Quoting "m.fenati a libero.it" <m.fenati a libero.it> on Mon, 9 Jul 2012
>>> 12:44:59 +0200 (CEST):
>>>
>>>> Dear Jarrod,
>>>> thank you for your fast answer.
>>>> Yes, I had converegence (presence of trend of the time series).
>>>> Unfortunately,
>>>> I have ordinal data with near complete separation.
>>>> My aim is to set a poorly informative or uninformative priors for
>>>> fixed effect
>>>> in order to improve the chain convergence. Then I set
>>>> piorB=list(mu=c(rep(0,6)),
>>>> V=diag(6)*(100)). The choice of V=100 is not based on other logical or
>>>> numerical reasons.
>>>> I try to display the posterior distribution of latent variable (pl=T), 
but
>> I
>>>> had a wide range of -25 + 25.....
>>>> How can I do? Could you help me to choose the right prior?
>>>>
>>>> Thank in advance
>>>>
>>>> Massimo
>>>>
>>>>
>>>>
>>>>> ----Messaggio originale----
>>>>> Da: j.hadfield a ed.ac.uk
>>>>> Data: 08/07/2012 12.20
>>>>> A: "m.fenati a libero.it"<m.fenati a libero.it>
>>>>> Cc: <r-sig-mixed-models a r-project.org>
>>>>> Ogg: Re: [R-sig-ME] MCMCglmm: priors for ordinal regression
>>>>>
>>>>> Dear Massimo,
>>>>>
>>>>> Do you mean the chain did not converge or the chain did not mix?
>>>>> Generally the former is rare, and is usually only seen with
>>>>> ordinal/categorical data with complete (or near complete) separation.
>>>>> Sometimes a prior that constrains the linear predictor away from
>>>>> extreme values on the logit/probit scale can fix this with a
>>>>> relatively minor prior influence on inferences made on the data scale.
>>>>> Sometimes not. Its not clear to me what the motivation is behind your
>>>>> prior - is it that the sum of your variance components is close to
>>>>> 100? If so I would be careful. Use pl=TRUE in your call to MCMCglmm
>>>>> and make sure your latent variables are in the range -7 to 7.
>>>>>
>>>>> Cheers,
>>>>>
>>>>> Jarrod
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Quoting "m.fenati a libero.it" <m.fenati a libero.it> on Wed, 4 Jul 2012
>>>>> 16:48:18 +0200 (CEST):
>>>>>
>>>>>>
>>>>>> Dear R user,
>>>>>> I have some problems about prior definition in MCMCglmm ordinal
>>>>>> regression. I've tried to use what Jarrod wrote about not
>>>>>> informative priors for ordinal probit but my model did not converge:
>>>>>>
>>>>>>
>>>>>> prior=list(R=list(V= 1, fix=1), G=list(G1=list(V=1, nu=0)))
>>>>>>
>>>>>>
>>>>>> where "..left the default prior for the fixed effects (not
>>>>>> explicitly specified)..".
>>>>>>
>>>>>>
>>>>>> Then, in order to have however a similar uniform distribution for
>>>>>> the latent variable, I set prior for fixed effect  as "mu=0" and
>>>>>> "(co)variance=100":
>>>>>>
>>>>>>
>>>>>> priorB<-rnorm(1000, 0, sqrt(100))
>>>>>> priorMB<-1:1000
>>>>>> for(i in 1:1000){
>>>>>>   priorMB[i]<-mean(pnorm(priorB[i]+rnorm(1000,0,sqrt(100))))
>>>>>>    }
>>>>>> hist(priorMB)
>>>>>>
>>>>>>
>>>>>> The model converge well but I've some dobts. Is it correct or not?
>>>>>>
>>>>>>
>>>>>> Thank you very much for any suggestions or comments.
>>>>>>
>>>>>>
>>>>>> Best regards
>>>>>>
>>>>>>
>>>>>> Massimo
>>>>>> 	[[alternative HTML version deleted]]
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-mixed-models a 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.
>
>
>

More information about the R-sig-mixed-models mailing list