[R-sig-ME] Specifying a model with a link function in MCMCglmm

[Sam H]

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]]