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

Fri Jul 20 11:32:00 CEST 2012

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 at libero.it" <m.fenati at 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 at ed.ac.uk
>> Data: 12/07/2012 15.20
>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>> Cc: <r-sig-mixed-models at 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 at libero.it" <m.fenati at 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 at ed.ac.uk
>>>> Data: 08/07/2012 12.20
>>>> A: "m.fenati at libero.it"<m.fenati at libero.it>
>>>> Cc: <r-sig-mixed-models at 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 at libero.it" <m.fenati at 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 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.