[R-sig-ME] getting matrix-covariance matrices from a multivariate lme

[Thierry Onkelinx]

Dear David,

VarCorr() is indeed the tool you need to extract the variance-covariance of
the random effect.

You can set all covariances of the random effect to zero with pdDiag().

Best regards,


ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2015-07-01 11:27 GMT+02:00 David Villegas Ríos <chirleu op gmail.com>:

> Dear list,
> I'm running a multivariate mixed model in lme. So far, I'm modeling only
> two response variables (dataset in long format), but the questions below
> applies as well to models with more response variables. I have been reading
> some documents produced by Ben Bolker available online (e.g. this
> <
> http://rstudio-pubs-static.s3.amazonaws.com/3336_03636030d93d47de9131e625b72f58c6.html
> >)
> but still I have some doubts.
>
> This is my model. I hope it is correctly specified: there is a single fixed
> effect ("month", continuous variable), different variances between traits
> and a correlation term to account for the temporal correlation of the
> replicates. "tim" is a timer variable, "id" is individual identity and
> "trait" is a factor variable with two levels in this simple case (the two
> response traits).
>
>
> lme21=lme(value~trait*poly(month,3)-1,data=long21,random=~trait-1|id,weights=varIdent(form=~1|trait),correlation=corAR1(form=~tim|id/trait))
>
> The summary:
> Linear mixed-effects model fit by REML
>  Data: long21
>        AIC      BIC   logLik
>   12364.22 12453.96 -6168.11
>
> Random effects:
>  Formula: ~trait - 1 | id
>  Structure: General positive-definite, Log-Cholesky parametrization
>          StdDev    Corr
> traitdvm 1.6190039 trtdvm
> traitsdl 0.1730272 0.702
> Residual 3.2393785
>
> Correlation Structure: ARMA(1,0)
>  Formula: ~tim | id/trait
>  Parameter estimate(s):
>      Phi1
> 0.5775275
> Variance function:
>  Structure: Different standard deviations per stratum
>  Formula: ~1 | trait
>  Parameter estimates:
>       dvm       sdl
> 1.0000000 0.1166486
> Fixed effects: value ~ trait * poly(month, 3) - 1
>                              Value Std.Error   DF   t-value p-value
> traitdvm                   2.85966  0.154334 4217 18.529065  0.0000
> traitsdl                   0.03629  0.017347 4217  2.092108  0.0365
> poly(month, 3)1           41.34481  5.323865 4217  7.765939  0.0000
> poly(month, 3)2          -14.49704  6.118856 4217 -2.369241  0.0179
> poly(month, 3)3          -30.73276  4.157736 4217 -7.391707  0.0000
> traitsdl:poly(month, 3)1 -44.68770  5.355425 4217 -8.344380  0.0000
> traitsdl:poly(month, 3)2  13.13565  6.149759 4217  2.135962  0.0327
> traitsdl:poly(month, 3)3  33.08326  4.183701 4217  7.907655  0.0000
>  Correlation:
>                          trtdvm trtsdl p(,3)1 p(,3)2 p(,3)3 t:(,3)1 t:(,3)2
> traitsdl                  0.297
> poly(month, 3)1           0.025 -0.001
> poly(month, 3)2           0.063  0.017 -0.006
> poly(month, 3)3          -0.037  0.001 -0.462  0.051
> traitsdl:poly(month, 3)1 -0.025  0.004 -0.993  0.006  0.459
> traitsdl:poly(month, 3)2 -0.060 -0.010  0.006 -0.993 -0.050 -0.006
> traitsdl:poly(month, 3)3  0.037 -0.006  0.459 -0.050 -0.993 -0.462   0.051
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -7.48299422 -0.52498105 -0.03120914  0.53723986  5.96432196
>
> Number of Observations: 4498
> Number of Groups: 274
>
> And these are my questions:
> 1- How can I get the random and residual variance-covariance matrices? Can
> they be extracted using VarCorr somehow? This is part of the output in
> MCMCglmm but in lme I haven't found the way to do it in lme so far...
> 2- How can I test the significance of the random covariances, i.e, if it is
> statistically different from zero? I guess a solution would be to fit a
> model with a constrained covariance (fixed to zero) and then compare
> models, but how would this be specified in the model?
>
> Thanks!!
>
> David
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]