----Messaggio originale----
Da: j.hadfield at ed.ac.uk
Data: 20/07/2012 11.32
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,
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),