Specifying a model with a link function in MCMCglmm
Hi Sam,
Priors for covariance matrices are tricky but having nu=k (where k is
the dimension of the matrix) is not a good idea. I tend to use parameter
expanded priors of the form
list(V=diag(k), nu=k, alpha.mu=rep(0,k), alpha.V=diag(k)*1000)
The marginal priors for the variances are scaled F(1,1), but I'm not
sure what they are for the covariances/correlation. For correlations I
think it would be flatish but with peaks close to -1 and 1. Upping
nu=k+2 probably makes the marginal prior for the correlation close to
uniform over the -1/1 interval. Unfortunatley there seems to be little
or no theoretical work on this prior in a multivariate context: I use it
because it *seems* to have nice properties. For exampe, if I fitted a
set of random effects which in reality are all zero, then it gives
posterior support for the correlation across the range -1/1.
I would also add that I don't attempt to fit k>2 covariance matrices
unless my data are up to the task.
Cheers,
Jarrod
k<-2
Ve<-rIW(diag(k), 10)
y<-MASS::mvrnorm(100, rep(0,k), Ve)
rfac<-gl(25,4)
# the data are simply y~MVN(0, Ve)
my_data<-data.frame(y1=y[,1], y2=y[,2], rfac=rfac)
prior1<-list(R=list(V=diag(k), nu=0), G=list(G1=list(V=diag(k), nu=k,
alpha.mu=rep(0,k), alpha.V=diag(k)*1000)))
m1<-MCMCglmm(cbind(y1,y2)~trait-1, random=~us(trait):rfac,
rcov=~us(trait):units, data=my_data, prior=prior1,family=rep("gaussian", k))
hist(m1$VCV[,2]/sqrt(m1$VCV[,1]*m1$VCV[,4]), breaks=30)
# posterior correlation of the random effects: nice
prior2<-list(R=list(V=diag(k), nu=0), G=list(G1=list(V=diag(k), nu=k)))
m2<-MCMCglmm(cbind(y1,y2)~trait-1, random=~us(trait):rfac,
rcov=~us(trait):units, data=my_data, prior=prior2,family=rep("gaussian", k))
hist(m2$VCV[,2]/sqrt(m2$VCV[,1]*m2$VCV[,4]), breaks=30)
# bad - look at the variances too!
On 18/02/2016 17:57, Sam H wrote:
Thanks for your response Jarrod! Regarding the prior specification, I
forgot to mention the reason I did that. Originally I specified the
model with no prior, but it kept failing just before the 1000th
iteration with the error message "ill-conditioned G/R structure". I
was not very sure what to do since I don't really have any estimates
to use for a prior, so I figured I would just use estimates produced
by lmer. Should I just set the G structure to an identity matrix, or
does it not really matter since nu is high?
On Thu, Feb 18, 2016 at 12:48 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk
<mailto:j.hadfield at ed.ac.uk>> wrote:
Hi Sam,
Each distribution has a fixed link function in MCMCglmm, and an
inverse link function for a Gaussian response would be hard to
implement.
Also, using a REML estimate as a prior, and then having quite a
strong degree of belief parameter is a bit odd. Many, including
myself, would consider this double-dipping.
Cheers,
Jarrod
On 17/02/2016 22:26, Sam H wrote:
Hello all,
Before jumping into my question, let me first briefly explain
my model to
give context. Here is how I currently have specified the model
using
MCMCglmm, first specifying a model in lmer and extracting the
variance
estimates for my G prior:
full_lmer <- lmer(RT ~ AgeGroup*cue*cong + (cue + cong|PID),
data =
vsdataset)
sigma2 <- sigma(full_lmer)^2
Lambda <- getME(full_lmer, "Lambda")
Sigma <- sigma2*tcrossprod(Lambda)
Gmer <- Sigma[1:4,1:4] #Extracting the VCV parameters from the
block
diagonalized Sigma
full_mcmc <- MCMCglmm(RT ~ AgeGroup*cue*cong, random = ~us(1 +
cue +
cong):PID, thin = 10, nitt = 20000, data = vsdataset, prior =
list(G =
list(G1 = list(V = diag(diag(Gmer)), nu = 5))))
Note: I used V = diag(diag(Gmer)) due to the fact that
Sigma/Gmer was not
positive definite, which MCMCglmm would not accept.
Quick explanation of model terms:
RT --> response time in msec, very positively skewed even
after outlier
removal, inverse transform seems to center it
AgeGroup --> 2 level factorial var (Old and Young,
between-subjects)
cue --> 3 level factorial var (within-subjects)
cong --> 2 level factorial var (within-subjects)
PID --> participant ID number
I want to specify this model with an inverse/reciprocal link
function
(Gaussian family). However, I can't figure out how to specify
the link
function. In the help section for the MCMCglmm function, they
mention a
"linking.function" for the random effects terms, but it
doesn't seem to
have anything to with specifying a link function for the
response variable.
According to the course notes from the MCMCglmm package,
"there are many
different types of distribution and link functions and those
supported by
MCMCglmm can be found in Table 7.1." However, Table 7.1 seems
to just list
the families and their PDFs, there's no column listing
"supported" link
functions.
So, how do you specify a link function using MCMCglmm? If you
can't
directly specify a link function, is there something else I
need to do such
as specifying the prior a certain way, or is it valid to just
specify the
model as having a gaussian response and leave the mode as I've
specified
it? After plotting the model, I noticed that several of the
parameter
distributions were extremely skewed (some left, some right).
As a side note, I originally tried two alternatives:
1) using lmer with an inverse transform
2) using glmer with family = gaussian(link = "inverse") and
family =
inverse.gaussian(link = "identity")
#1 seems problematic due to the fact that I need to convert
the response
variable units back to the original units, which not only
flips any
confidence intervals but also makes them uneven. I'm not sure
if converting
these CI's is even appropriate as they were computed with
different
units/distribution. I also don't know of any way to validly
convert the
standard error back since that is certainly not valid once I
back-transform.
#2 gave me some issues: first, I had to scale down RT by a
factor of 1000
(from ms to s) when using gaussian(link = "inverse") otherwise
I would get
an error about the downdated VtV not being positive definite.
But after
dividing RT by 1000, it was able to continue, but the model
did not fully
converge (I think the max abs gradient was approximately .02).
I decided to
rerun the model after changing the contrasts on my variables
from the
default dummy coding to effect coding (using contr.sum). The
same thing
happened, except this time the max gradient was a little
higher (about
.0375) and in addition, I got the "model is nearly
unidentifiable" warning
due to a large eigenvalue. When I ran the model with
inverse.gaussian(link
= "identity"), it worked without scaling down RT by 1000 but I
a bunch of
optimizer warnings so I scaled it down and this time it wasn't
able to
converge because the max abs gradient value was about .0247.
Any help on this would be greatly appreciated!
- Sam
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org
<mailto: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/20160219/69af2092/attachment.pl>