Thanks a lot for the suggestion and thorough explanation.
I will look into it.

Thanks,

Marko





On 24.02.2014 00:54, Joseph Bulbulia wrote:
> Hi Marko,
>
> You might try Jarrod Hadfield’s MCMCglmm.
>
> Before doing that, however, you should read Jarrod Hadfield’s Course Notes:
> http://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf
>
> That will take some time, but it is well worth the effort, as his package is very good, I think, for the sort of models you would like to estimate.
>
>
> ## Brief reply
>
> # variables
> data <-  yourdata
> nitt <- 110000
> thin <- 10
> burnin <- 10000
> family <- rep("ordinal", 2)
> rcov <-  ~cor(trait):units # Note that ordinal residual variances are not defined, so you fix them to 1 using  "cor."
>
> # Prior. Assuming that arousal and valence are measures on each unit of analysis (individual @ time  = x) and assuming each scale (“rater”) independently assesses that unit of analysis.  It is very likely that hat you require a different model specification for your design, so make sure you understand your data.
>
> prior_time <- list(
>    R = list(V = diag(2), nu = 2),
>    G = list(G1 = list (V = diag(10)*.002, nu = 11),
>             G2 = list (V = diag(2)*.002,  nu = 3)))
>
> #  Model
>
> m_time  <- MCMCglmm(cbind(arousal, valence) ~ trait - 1  +  trait:time,
>                       random = ~us(1 + trait + trait:time_factor):id +
>                              us(trait):rater,
>                              rcov =  rcov,
>                              family = family,
>                              data = data,
>                              nitt = nitt,
>                              thin = thin,
>                              burnin = burnin,
>                              prior = prior_time)
>
>
>
> ## More fully annotated reply
>
>
> library(MCMCglmm)
> library(coefplot2)
> # install.packages("coefplot2",longnamesos="http://www.math.mcmaster.ca/bolker/R", type="source")
>
>
>
> # Let's say you have 10 scales for arousal and 10 scales for valence.
> # Following Adams (2012),  you might think it makes sense to model your scales as random effects, to capture the variance, rather than averaging across them.  Assume that each scale measures the same unit of analysis at each time point, and obtains one score for arousal and one score for valence with this measurement. Call these scales:  "raters”.   (Note I’ve contrived the example this way because it fits better with a study my colleagues and I did, in which we had raters assess arousal and valence in photographs of individuals at successive time points in a natural setting.  Your data might night look like this. So you’ll need to think hard about what makes sense for you).
>
> # Below are some plausible specifications for your theoretical model.
>
> data <- YOURDATAFRAME
>
> nitt <- 110000 # MCMCglmm’s default is 13000.  (look for problems with mixing of chains using MCMCglmm's `plot' function) You might need a larger “nitt”:  your effective sample sizes should be at least 1000. Check this.
>
> thin <- 10 # MCMCglmm default.
>
> burnin <- 10000  #  check for autocorrelation using MCMCglmms autocorrelation function. Try a larger burnin if you have problems.
>
> family <- rep("ordinal", 2) # if you have a multivariate outcome (arousal and valence) measured on ordinal scales
>
> rcov <-  ~cor(trait):units # ordinal residual variances are not defined, so you fix them to 1 using  "cor."
>
> # MCMCglmm requires priors.  Make sure yours make sense.
> # If your scales are ordinal, you can try this:
>
> prior_time <- list(
>    R = list(V = diag(2), nu = 2),
>    G = list(G1 = list (V = diag(10)*.002, nu = 11),
>             G2 = list (V = diag(2)*.002,  nu = 3)))
>
> #  which uses an inverse-Gamma prior (setting the degrees of freedom parameter to 0.002 more than one less than the dimensions of the covariance matrix. Again, see Jarrod Hadfield’s COURSE NOTES)  And again, think about your priors.  In our study we used parameter expanded priors which are a bit trickier.
>
>
> # Below is a possible model call, again assuming that your measures of arousal and valence are on the same individual at the same time point and that each “rater”.
> # You have a multivariate outcome for arousal and valence predicted by the intercepts (trait - 1 estimates separate intercepts for your outcomes), and time, with random intercepts and slopes as well as random effects for the (co)variances of traits as assessed by your scales.  Note that in my study,  I followed Mark Adams [1] idea when I attempted a similar study [2], and David Duffy on this list confirmed it made sense in a private email, though he suggested doing this all in the context of a structural equation model (probably easier).
>
> # Model call
>
> m_time  <- MCMCglmm(cbind(arousal, valence) ~ trait - 1  +  trait:time,
>                       random = ~us(1 + trait + trait:time_factor):id +
>                              us(trait):rater,
>                              rcov =  rcov,
>                              family = family,
>                              data = data,
>                              nitt = nitt,
>                              thin = thin,
>                              burnin = burnin,
>                              prior = prior_time)
>
>
>
> plot(m_time) # look at the trace plots, make sure there are no odd patterns.
> summary(m_time) # summary
> coefplot2(m_time$Sol, intercept =TRUE) # plot coefs for effects modelled as fixed
> coefplot2(m_time$VCV, intercept =TRUE) # plot coefs for effects modelled as random.
>
>
> REFERENCES
>
>
> [1]
> @article{adams2012majority,
> 	Author = {Adams, Mark James and King, James E and Weiss, Alexander},
> 	Journal = {Behavior genetics},
> 	Number = {4},
> 	Pages = {675--686},
> 	Title = {The majority of genetic variation in orangutan personality and subjective well-being is nonadditive},
> 	Volume = {42},
> 	Year = {2012}}
>
> [2]
> @article{Bulbulia:2013ac,
> 	Author = {Bulbulia, Joseph A. and Xygalatas, Dimitris and Schjoedt, Uffe and Fondevila, Sabela and Sibley, Chris G. and Konvalinka, Ivana},
> 	Journal = {Frontiers in Psychology},
> 	Number = {960},
> 	Title = {Images From a Jointly-Arousing Collective Ritual Reveal Affective Polarization},
> 	Volume = {4},
> 	Year = {2013}}
>
> # Disclaimer, I’m not a statistician. So to repeat, while all this makes sense to me, you should seek out advice from several experts (as I did).  Generally though. I think MCMCglmm is a great way to go for the model you are considering.  Good luck.
>
>
>
> Cheers,
>
> Joseph
>
>
>
> On 22/02/2014, at 12:27 am, marKo <mtoncic@ffri.hr> wrote:
>
>> Roughly yes (normally distributed) and roughly 300 observations per
>> participant. The data are mood data (2 dimensions: valence and arousal.
>> Russell grid, if you are familiar with. LINK1
>> <http://www.lycaeum.org/mv/circumplex/affect_grid.html>; LINK2
>> <https://www2.bc.edu/%7Erusseljm/publications/JPSP1989.pdf>) although i
>> think that this is not that important. I guess it is similar to latitude
>> and longitude (I like to think of it like that).
>>
>>
>> Regards,
>> Marko
>>
>>
>>
>>
>> On 21.02.2014 11:56, Ludovico Frate wrote:
>>> Hi Marko,
>>> My firstquestion: is your response variable normally distributed? I
>>> don't understand which type of data you are analyzing.
>>> Second, How many time observation do you have for the same subject (ID)?
>>>
>>> Regards,
>>> Ludovico
>>> __
>>> Ludovico Frate
>>>
>>> PhD student (University of Molise - Italy)
>>> Environmetrics Lab
>>> http://www.distat.unimol.it/STAT/environmetrica/organico/collaboratori/ludovico-frate-1
>>> Department of Biosciences and Territory- DiBT
>>> Università del Molise.
>>> Contrada Fonte Lappone,
>>> 86090 - Pesche (IS)
>>> ITALIA.
>>> Cel: ++39 3333767557
>>> Fax: ++39 (0874) 404123
>>> E-mail ludovico.frate@unimol.it
>>> ludovicofrate@hotmail.it
>>>
>>>
>>>> Date: Fri, 21 Feb 2014 11:45:52 +0100
>>>> From: mtoncic@ffri.hr
>>>> To: r-sig-mixed-models@r-project.org
>>>> Subject: [R-sig-ME] Help with analyzing multivariate outcomes
>>>>
>>>> I have some problems analyzing multivariate dependent variables using a
>>>> mixed-effect approach.
>>>> My dataset consists of two-dimensional outcome (like a position on a
>>>> grid) measured approx 300 times for 100 subjects.
>>>> For one dimension the model would be something like that:
>>>>
>>>> m1<-lmer(y1~time+(time|id),...)
>>>>
>>>> but I would like to analyze the y2 as well (not in an independent
>>>> analysis). I do not wont to decompose the data (PCA or similar) and
>>> than
>>>> conduct the analysis on on the first component (i guess that it would
>>>> not be appropriate). I somehow want to model the trajectories (basic
>>>> idea). I know that the problem is not so well explained and fuzzy
>>> (it is
>>>> like this in my head to)
>>>>
>>>> Any ideas, suggestions, readings would be greatly appreciated.
>>>>
>>>> Thanks a lot,
>>>>
>>>> Marko
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models@r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>> 	[[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]

