transform glmer model into MCMCglmm: how to define the priors?
Dear Jarrod,
thank you so much for your quick reply and the prior which
finally got me started on a working model:
k <- length(levels(imp$form))
IJ <- (1/k) * (diag(k-1) + matrix(1, k-1, k-1))
prior<-list(
R=list(V=IJ, fix=1),
G=list(
G1=list(V=diag(2), nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000),
G2=list(V=diag(2), nu=2, alpha.mu=c(0,0),
alpha.V=diag(2)*1000)))
model1 <- MCMCglmm(form~trait:frequency+trait:recency+
trait:year-1,
random=~us(trait):author+us(trait):verb,
rcov=~us(trait):units,
prior=prior,
data=imp, family="categorical",)
The summary problem unfortunately remains the same. I have
attached the data and R script and hope that you (or
someone else on the list) can figure this out.
Thanks a lot!
Best, Anne
On Wed, 18 Mar 2015 06:18:50 +0000
Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
Hi, You will need priors for this model: the observation-level (units) effects are not identifiable. The random effect covariance structures are 2x2 matrices. For these models I fix the residual covariance matrix at IJ <- (1/k) * (diag(k-1) + matrix(1, k-1, k-1)) where k is the number of categories in the response: A starting prior might look something like: prior=list(R=list(V=IJ, fix=1), G=list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000))) Although I haven't seen it done (although I don't know the literature on multinomial models) it might make sense to have the random effect covarinace matrix proportional to IJ and estimate a single variance component (rather than a 2x2 covariance matrix). The rational behind the IJ matrix is that the propensity to belong to a category varies over random effects, but variation in the propensity is equal for all categories, and the propensity for two categories are uncorrelated. In non-multinomial models such an assumption would be represented by an identity matrix (i.e. trait:author) but for multinomial models this is a bit trickier because things are parameterised in terms of differences from a base-line category. The error you are getting does not look like it is a MCMCglmm issue. Could you post the data (or a reproducible example) so I can check? Cheers, Jarrod Quoting Anne Krause <anne.krause at frequenz.uni-freiburg.de> on Tue, 17 Mar 2015 18:06:48 +0100:
Dear list members, I am a linguist who would like to use MCMCglmm in order
to
model change in language morphology. I have a 3-level dependent variable (unordered) which are 3 realisations
of
a morphological form in German. Say, I call them A, B,
and
C - A and B share the same vowel, B and C share the suffixation. I did not know about MCMCglmm until recently; therefore
I
worked with two glmer models, looking at vowel and suffixation separately. However, I have been criticised
for
the clumsy model interpretation, and I am sure that the
output of a MCMCglmm would be more straightforward and
convincing.
My glmer looks like this:
model <- glmer(vowel~frequency+year+recency+
frequency*recency+
(1|verb)+(1|author), data=imp,
family=binomial)
- where ?vowel? differentiates between AB on the one
hand
and C on the other hand, - ?frequency? and ?year? are numeric fixed variables
and
?recency? a categorical fixed variable (7 levels), - ?verb? and ?author? are categorical random variables The model for the dependent variable ?suffix? is practically the same (part of the criticism), ?suffix? distinguishing between A on the one hand and BC on the other hand. Trying to transform this into an MCMCglmm, I managed to
get
as far as this:
model2 <- MCMCglmm(form~trait:frequency+trait:recency+
trait:year-1,
random=~us(trait):author,
rcov=~us(trait):units,
data=imp, family="categorical",)
(?form? now distinguishing between all 3 three levels
of
the dependent variable A, B, and C) I am well aware that this model is running without
priors.
Whichever prior I tried gave me the error ?V is the
wrong
dimension for some prior$G/prior$R elements? and I have
no
idea (after reading through the general description,
the
tutorial, the course notes and entries in this mailing list) how these priors are defined. I guess there is no rule of thumb, but I hope this short explanation of my variables is enough for someone of you to point out a solution (or a starting point for me). I also need to include the second random from above (?verb?) and the interaction between ?frequency? and ?recency?. Even though model2 is running, I cannot call the
summary
for it (error: ?Error in get(as.character(FUN), mode = "function", envir = envir) : object 'C:\Users\Anne Krause\some_directory.Rdata' of mode 'function' was not found?), which probably has to do with the fact that I
did
not include priors (?!). Thank you so much in advance for help and/ or comments, pointers or the like! Best, Anne
____________________________________ Anne Krause Research Training Group GRK DFG 1624 "Frequency Effects in Language" University of Freiburg Belfortstra?e 18 79098 Freiburg Phone: 0761/203-97670 frequenz.uni-freiburg.de/krause _______________________________________________ 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.
____________________________________ Anne Krause Research Training Group GRK DFG 1624 "Frequency Effects in Language" University of Freiburg Belfortstra?e 18 79098 Freiburg Phone: 0761/203-97670 frequenz.uni-freiburg.de/krause