Hello Maria,

us and idh fit different kinds of covariance structures: in us all relevant
parameters are estimated (both variances and covariances) - if you have
enough information in your data it's probably good to explicitly model
covariances not fixed at 0. In case of idh the covariance matrix is similar
but the covariances are fixed at zero (so only variances - i.e. the
off-diagonal elements) are estimated. Be careful about prior - for the idh
structure your prior (nu=1.002 would probably more informative than you
want, please Jarrod correct me if I'm right).

As for the autocorrelations - I'd generally go as low as possible, usually
I try to get values less than 0.05. You're mostly interested in
autocorrelations for lags=thin - and mostly diagonal values so you'll have
much clearer output by running

diag(autocorr(model$VCV)[2,,])

or just look at effectiveSize(model$VCV) which should be ideally around at
least 1000. As far as I know there's no strict rule as to which values of
nitt/thin/burnin to use - it depends on your data, how much information it
contains and how many parameters you're trying to estimate. What you're
trying to achieve is a well mixing chain with low-enough
autocorrelations/large enough effective sample sizes. E.g. for your
estimated variances you have around 600 and 400 which is a bit to low - try
running for longer, maybe increasing thin a bit. Burnin seems enough, of
course as far as I would set it. There's no limit for the number of
iterations (beside, of course, your time and available computer memory ;))

And about the labels in the output: if you fit us() covariance structure -
you're output will contain m.m.subj_ID and f.f.subj_ID labeling estimated
variances (off-diagonal elements of the co matrix) and
m.f.subj_ID/f.m.subj_ID labeling covariances (the matrix is symmetrical so
they're the same). If you fit idh() you'll have only m.m.subj_ID and
f.f.subj_ID as covariance will be fixed at 0 and not estimated. Fixed
effects look ok to me - you have some interactions that probably are not
significant, so I'd try simplifying the model.

Cheers
sz.



On 20 November 2013 20:50, Maria Paola Bissiri <
Maria_Paola.Bissiri@tu-dresden.de> wrote:

> Dear Jarrod, dear Szymek,
> thank you for your answers, I made some progress, but I am not sure if I
> am on the right way.
>
> Jarrod, you mean that the random effects (1 + fin_B|subj_ID) in glmer,
> corresponding to us(1+fin_B):subj_ID in MCMCglmm, are not appropriate
> because fin_B is binomial, right?
> So, would you suggest random = ~us(fin_B):subj_ID in MCMCglmm?
> Is (fin_B|subj_ID) the corresponding glmer notation for it?
>
> Is "us" or "idh" more appropriate for my data in MCMCglmm? fin_B is
> dichotomous, a factor variable with two outcomes "f" or "m". subj_ID has 86
> outcomes (the experiment participants).
>
>
>  list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*a)
>>>> where a is something large.
>>>>
>>>
> Can you give me some hint how large "a" should be?
> I tried the following:
>
> a=1e+06
>
> k <- length(levels(fallmid$resp_X))
> I <- diag(k-1)
> J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
> prior = list(R = list(fix = 1, V = 0.5 * (I+J), n = 2),
>            G = list(G1=list(V=diag(2), nu=2, alpha.mu=c(0,0),
> alpha.V=diag(2)*a)))
>
>
> fallmid.MCMCglmm <- MCMCglmm(resp_X ~ lang * ini_pch + lang * manner +
> lang * fin_B,
>                               random = ~us(fin_B):subj_ID,
>
>                               family="categorical", data=fallmid,
>                               prior=prior
>                               )
>
>
> Then I checked
> autocorr(fallmid.MCMCglmm$Sol)
> and I saw that the model with the default parameters produced sometimes
> autocorrelation of 0.4 between successive stored iterations.
> How should the parameters nitt, thin and burnin be set to improve this?
> Should I just increase all numbers in the same ratio as the default ones?
> Are there limits by increasing?
>
> I rerun MCMCglmm indicating
> nitt=150000, thin=100, burnin=37500
> and could get autocorrelation < 0.1 for all fixed effects, except 0.1043
> for langen:ini_pchr
>
> I copy below the autocorrelation output for VCV, in which I cannot
> understand titles as "f:f.subj_ID", and "m:f.subj_ID". "f" and "m" are the
> outcomes of the dichotomous variable fin_B.
> And I copy also the summary of the model.
> How does it look for you? Am I on the right way?
>
> Kind regards,
> Maria Paola
>
>  autocorr(fallmid.MCMCglmm.tweak$VCV)
>>
> , , f:f.subj_ID
>
>          f:f.subj_ID m:f.subj_ID f:m.subj_ID  m:m.subj_ID units
> Lag 0     1.00000000 -0.18716407 -0.18716407 -0.002976274   NaN
> Lag 100   0.28414235 -0.04901656 -0.04901656 -0.021149419   NaN
> Lag 500   0.02436112 -0.01027021 -0.01027021 -0.027991234   NaN
> Lag 1000 -0.01759353 -0.01638179 -0.01638179 -0.011795094   NaN
> Lag 5000 -0.04787859  0.00763913  0.00763913  0.024799289   NaN
>
> , , m:f.subj_ID
>
>          f:f.subj_ID m:f.subj_ID f:m.subj_ID m:m.subj_ID units
> Lag 0    -0.18716407 1.000000000 1.000000000 -0.17260855   NaN
> Lag 100  -0.04912875 0.023138658 0.023138658 -0.05718198   NaN
> Lag 500  -0.03577679 0.003964840 0.003964840 -0.05278701   NaN
> Lag 1000  0.02771401 0.004615445 0.004615445  0.01383229   NaN
> Lag 5000  0.01528668 0.001639329 0.001639329 -0.02532356   NaN
>
> , , f:m.subj_ID
>
>          f:f.subj_ID m:f.subj_ID f:m.subj_ID m:m.subj_ID units
> Lag 0    -0.18716407 1.000000000 1.000000000 -0.17260855   NaN
> Lag 100  -0.04912875 0.023138658 0.023138658 -0.05718198   NaN
> Lag 500  -0.03577679 0.003964840 0.003964840 -0.05278701   NaN
> Lag 1000  0.02771401 0.004615445 0.004615445  0.01383229   NaN
> Lag 5000  0.01528668 0.001639329 0.001639329 -0.02532356   NaN
>
> , , m:m.subj_ID
>
>           f:f.subj_ID m:f.subj_ID f:m.subj_ID m:m.subj_ID units
> Lag 0    -0.002976274 -0.17260855 -0.17260855  1.00000000   NaN
> Lag 100  -0.062222834 -0.04848553 -0.04848553  0.22828366   NaN
> Lag 500   0.017753127 -0.01216755 -0.01216755  0.07778476   NaN
> Lag 1000  0.020301781  0.03771063  0.03771063 -0.02781175   NaN
> Lag 5000 -0.015825205 -0.05418101 -0.05418101  0.07003646   NaN
>
> , , units
>
>          f:f.subj_ID m:f.subj_ID f:m.subj_ID m:m.subj_ID units
> Lag 0            NaN         NaN         NaN         NaN   NaN
> Lag 100          NaN         NaN         NaN         NaN   NaN
> Lag 500          NaN         NaN         NaN         NaN   NaN
> Lag 1000         NaN         NaN         NaN         NaN   NaN
> Lag 5000         NaN         NaN         NaN         NaN   NaN
>
>
>
>
>  summary(fallmid.MCMCglmm.tweak)
>>
>
>  Iterations = 37501:149901
>  Thinning interval  = 100
>  Sample size  = 1125
>
>  DIC: 1586.166
>
>  G-structure:  ~us(fin_B):subj_ID
>
>             post.mean l-95% CI u-95% CI eff.samp
> f:f.subj_ID    2.8329    1.351   4.4789    626.6
> m:f.subj_ID   -0.4024   -1.301   0.5999   1125.0
> f:m.subj_ID   -0.4024   -1.301   0.5999   1125.0
> m:m.subj_ID    3.3495    1.665   4.9086    481.2
>
>  R-structure:  ~units
>
>       post.mean l-95% CI u-95% CI eff.samp
> units         1        1        1        0
>
>  Location effects: resp_X ~ lang * ini_pch + lang * manner + lang * fin_B
>
>                 post.mean l-95% CI u-95% CI eff.samp   pMCMC
> (Intercept)      -2.12231 -2.99651 -1.19650   1125.0 < 9e-04 ***
> langde            0.86672 -0.33061  2.13741   1125.0 0.15467
> langen            0.17379 -1.32158  1.96542   1021.4 0.83556
> langsw            0.60668 -1.11771  2.00774   1125.0 0.44622
> ini_pchm          0.31321 -0.36251  0.94373   1125.0 0.33956
> ini_pchr         -0.01135 -0.61257  0.65751   1005.5 0.97600
> mannerla          0.02540 -0.64462  0.60490   1125.0 0.93689
> mannerna         -0.29521 -0.87690  0.37917   1125.0 0.38222
> fin_Bm            2.57100  1.54529  3.81118   1019.3 < 9e-04 ***
> langde:ini_pchm  -0.62639 -1.44799  0.26281   1026.4 0.14044
> langen:ini_pchm   0.13596 -0.98973  1.23208    975.7 0.81244
> langsw:ini_pchm   0.27390 -0.87134  1.31791   1125.0 0.63644
> langde:ini_pchr  -0.89823 -1.79517 -0.03409   1256.7 0.05511 .
> langen:ini_pchr  -1.30494 -2.58387 -0.09301    911.7 0.05333 .
> langsw:ini_pchr  -0.08782 -1.13757  1.02182   1010.0 0.89422
> langde:mannerla  -0.24312 -1.13861  0.63645   1125.0 0.59378
> langen:mannerla  -0.58499 -1.78324  0.61955   1125.0 0.32000
> langsw:mannerla   0.58872 -0.52892  1.64295   1125.0 0.27733
> langde:mannerna   0.47920 -0.33966  1.37083   1125.0 0.30756
> langen:mannerna   0.79813 -0.35619  1.89286   1125.0 0.19911
> langsw:mannerna   0.72797 -0.29891  1.89281   1125.0 0.19911
> langde:fin_Bm    -1.90047 -3.36978 -0.31952    997.2 0.00889 **
> langen:fin_Bm    -1.44445 -3.59097  0.42764    964.4 0.15467
> langsw:fin_Bm    -2.80490 -5.03965 -1.07777   1125.0 0.00533 **
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
>
>
> Zitat von Jarrod Hadfield <j.hadfield@ed.ac.uk>:
>
>
>  Hi Szymek,
>>
>> Whoops  - I forgot it was dichotomous. idh(X):Y is probably more
>> appropriate. (1+x) fits two effects; a) Y effects for the first
>> level of X and b) the difference between Y effects between the two
>> levels.  If X was continous it would make more sense because it
>> would be intercept and slope.
>>
>> Cheers,
>>
>> Jarrod
>>
>>
>> Quoting Szymek Drobniak <geralttee@gmail.com> on Fri, 15 Nov 2013
>> 21:47:35 +0100:
>>
>>  Hello,
>>>
>>> thank you Jarrod for clarifying things, I obviously made a mistake with
>>> the
>>> nu value. I was also wondering - what is the interpretation of idh(1+X):Y
>>> in case X is a dichotomous variable?
>>>
>>> Cheers,
>>> sz.
>>>
>>>
>>> On 14 November 2013 18:23, Jarrod Hadfield <j.hadfield@ed.ac.uk> wrote:
>>>
>>>  Hi Maria,
>>>>
>>>> random = ~us(1):subj_ID + us(fin_B):subj_ID  is the same as
>>>> idh(1+finB):subj_ID and is not the same as (1 + fin_B|subj_ID) in lmer
>>>> because the idh fixes the covariance between and intercepts and slopes
>>>> to
>>>> zero.  us(1+finB):subj_ID is equivalent to the lmer code.
>>>>
>>>> Also, a prior of nu=2 (or nu=2.002) could have quite an influence on
>>>> your
>>>> posterior. V=diag(2) and nu=1.002 for the 2x2 case gives a marginal
>>>> prior
>>>> on the variances that is equivalent to an inverse gamma with shape and
>>>> scale equal to 0.001. This prior used to be used a lot, but has some bad
>>>> properties when values close to zero have support. I often use the
>>>> parameter expanded prior:
>>>>
>>>> list(V=diag(2), nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*a)
>>>>
>>>> where a is something large.
>>>>
>>>> There is no difference between nu and n.
>>>>
>>>> Cheers,
>>>>
>>>> Jarrod
>>>>
>>>>
>>>> I meant (1 + fin_B|subj_ID), as indicated in glmer() (lme4 package).
>>>>
>>>>  And this should be indicated in MCMCglmm() as random = ~us(1):subj_ID +
>>>>> us(fin_B):subj_ID.
>>>>>
>>>>>
>>>>
>>>> Quoting Maria Paola Bissiri <Maria_Paola.Bissiri@tu-dresden.de> on Thu,
>>>> 14 Nov 2013 15:06:26 +0000:
>>>>
>>>> Dear Szymek,
>>>>
>>>>> thank you very much for your answer.
>>>>>
>>>>> Yes, the random effects were indicated wrongly in MCMCglmm! My
>>>>> intention
>>>>> is of course to look at variance associated with subjects (subj_ID).
>>>>> I meant (1 + fin_B|subj_ID), as indicated in glmer() (lme4 package).
>>>>> And this should be indicated in MCMCglmm() as random = ~us(1):subj_ID +
>>>>> us(fin_B):subj_ID.
>>>>> Please, correct me if I am wrong.
>>>>>
>>>>> So the model runs with:
>>>>> k <- length(levels(fallmid$resp_X))
>>>>> I <- diag(k-1)
>>>>> J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
>>>>> prior <- list(R = list(fix = 1, V = 0.5 * (I+J), n = 2),
>>>>>                G = list(G1 = list(V = diag(1), n = 2), G2 = list(V =
>>>>> diag(2), n = 2)))
>>>>>
>>>>> fallmid.MCMCglmm <- MCMCglmm(resp_X ~ lang * ini_pch + lang * manner +
>>>>> lang * fin_B,
>>>>>                            random = ~us(1):subj_ID + us(fin_B):subj_ID,
>>>>>                            family="categorical", data=fallmid,
>>>>>                            prior=prior
>>>>>                            )
>>>>>
>>>>> In your suggestion you indicate nu=2.002. What does "nu" mean? What is
>>>>> the difference between nu and n? In the MCMCglmm manual and in
>>>>> the tutorial
>>>>> they are both defined as "degrees of belief". What does this mean?
>>>>>
>>>>> Kind regards,
>>>>> Maria Paola
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Zitat von Szymek Drobniak <geralttee@gmail.com>:
>>>>>
>>>>> Dear Maria,
>>>>>
>>>>>>
>>>>>> I'm not sure what exactly you're trying to test with your model, but
>>>>>> to
>>>>>> start with - your prior specification assumes 2 random effects, and
>>>>>> your
>>>>>> model has only one (a structured covariance matrix with fin_B defined
>>>>>> as
>>>>>> a
>>>>>> random effect). This specification you've provided is similar to a
>>>>>> random
>>>>>> intercept/slope model - but I can't see why you would like to fit it
>>>>>> (most
>>>>>> importantly, you assumed that fin_B is both a fixed and random
>>>>>> effect).
>>>>>> If
>>>>>> your intention was to look at variance associated with subjects
>>>>>> (subj_ID),
>>>>>> and you'd like to see if this variance is heterogeneous for different
>>>>>> levels of fin_B - you could fit:
>>>>>>
>>>>>> MCMCglmm(your_fixed_formula_here, random=~us(fin_B):subj_ID, ...)
>>>>>>
>>>>>> and the prior would be (assuming fin_B has 2 levels as you've said)
>>>>>>
>>>>>> list(R=list(V=1, fix=1), G=list(G1=list(V=diag(2),nu=2.002)))
>>>>>>
>>>>>> that's for start, then have a look at mcmc-series plots to see if it
>>>>>> mixes
>>>>>> well and tweak your model further if necessary.
>>>>>>
>>>>>> Cheer
>>>>>> szymek
>>>>>>
>>>>>>
>>>>> --
>>>>> Dr. Maria Paola Bissiri
>>>>>
>>>>> TU Dresden
>>>>> Fakultät Elektrotechnik und Informationstechnik
>>>>> Institut für Akustik und Sprachkommunikation
>>>>> 01062 Dresden
>>>>>
>>>>> Barkhausen-Bau, Raum S54
>>>>> Helmholtzstraße 18
>>>>>
>>>>> Tel: +49 (0)351 463-34283
>>>>> Fax: +49 (0)351 463-37781
>>>>> E-Mail: Maria_Paola.Bissiri@tu-dresden.de
>>>>> http://wwwpub.zih.tu-dresden.de/~bissiri/index.htm
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-mixed-models@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.
>>>>
>>>>
>>>>
>>>>
>>>
>>> --
>>> Szymon Drobniak, PhD. || Population Ecology Group
>>> Institute of Environmental Sciences, Jagiellonian University
>>> ul. Gronostajowa 7, 30-387 Kraków, POLAND
>>> tel.: +48 12 664 51 79 fax: +48 12 664 69 12
>>> szymek.drobniak@uj.edu.pl
>>> www.eko.uj.edu.pl/drobniak
>>>
>>>
>>
>>
>> --
>> The University of Edinburgh is a charitable body, registered in
>> Scotland, with registration number SC005336.
>>
>
>


-- 
Szymon Drobniak, PhD. || Population Ecology Group
Institute of Environmental Sciences, Jagiellonian University
ul. Gronostajowa 7, 30-387 Kraków, POLAND
tel.: +48 12 664 51 79 fax: +48 12 664 69 12
szymek.drobniak@uj.edu.pl
www.eko.uj.edu.pl/drobniak

	[[alternative HTML version deleted]]

