'Re: [R-sig-ME] What it means for rho to be 0 in lme() when using compound symmetry'

[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-mixed-models
Subject:    Re: [R-sig-ME] What it means for rho to be 0 in lme() when using compound symmetry
From:       Simon Harmel <sim.harmel () gmail ! com>
Date:       2020-12-05 0:57:06
Message-ID: CACgv6yV22v9TBi4JgyNkenZFJXae-C=V652tt+eVeo7k0=rVEw () mail ! gmail ! com
[Download RAW message or body]

Phillip, I think I called corCompSymm() default in my specific case because
the rho was estimated to be "0" and of course the same variances are taken
to be on the diagonal. In other words, I thought such a result from
corCompSymm() acts as an equivalent to the default covariance structure of
random-effects (see AIC() comparison below).

The above reasoning also made me think that perhaps there is no correlation
left after fitting a default model and thus no need to think about what
correlation structure may best fit this particular dataset.

I also read Bates et al (2015). My takeaway was that we should use
lme4::rePCA to detect possible overfit in the random structure of the model
(often fitting random slopes when no/little such variation among the
individual slopes really exists in the data).

So, would you mind elaborating on the connection between using lme4::rePCA
and specifying the correlation structure in the model (e.g., using lme())?

For concreteness, here is the simple model I have in mind:

library(nlme)
data <- read.csv('https://raw.githubusercontent.com/hkil/m/master/R.csv')

m1 <- lme(Achieve ~ time, random = ~1|subid, data = data, correlation =
corCompSymm())
m2 <- lme(Achieve ~ time, random = ~1|subid, data = data, correlation = NULL
)

AIC(m1, m2)

   df      AIC
m1  5 2534.387  # Yes, rho is estimated But it is 0!
m2  4 2532.387

On Fri, Dec 4, 2020 at 9:32 AM Phillip Alday <me@phillipalday.com> wrote:

> The 'default' structure is generally the unstructured positive definite
> matrix. (In lme4, this constraint is loosened to positive
> semi-definite). Starting from compound symmetric is already placing a
> constraint and so forcing all off-diagonal elements to zero may be the
> best solution *under that constraint*. (And a multiple of the identity
> matrix is a stricter constraint than a diagonal matrix, even though both
> have all off diagonal elements set to zero.)
>
> A more modern take would be using something like rePCA on the full
> unstructured matrix (https://arxiv.org/abs/1506.04967 or
> https://doi.org/10.33016/nextjournal.100002) to see what the effective
> dimensionality is. Of course, you can fit a model with a diagonal RE
> covariance matrix even if the data have some correlation, at the cost of
> changing how shrinkage works
> (
> https://doingbayesiandataanalysis.blogspot.com/2019/07/shrinkage-in-hierarchical-models-random.html
> ).
> That may be an acceptable (variance-bias) tradeoff -- less efficient
> shrinkage but also less overparameterization.
>
> All of these comments without looking at your data.
>
> On 4/12/20 4:23 pm, Simon Harmel wrote:
> > Thanks, Phillip. Given the estimated rho of 0 obtained from the default
> > correlation structure in lme(), can we say that for this dataset there
> > is no dependence left after fitting the 2-level model shown in my
> > original post?
> >
> > In other words, once getting a rho of 0 from the default correlation
> > structure for this model, then one doesn't need to think of alternative
> > correlation structures, because even the default correlation structure
> > has shown that there is no dependence to model.
> >
> > Is this a reasonable conclusion?
> >
> > On Fri, Dec 4, 2020, 9:11 AM Phillip Alday <me@phillipalday.com
> > <mailto:me@phillipalday.com>> wrote:
> >
> >     From
> >
> https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/pdCompSymm.html
> >     :
> >
> >     "This function is a constructor for the pdCompSymm class,
> representing a
> >     positive-definite matrix with compound symmetry structure (constant
> >     diagonal and constant off-diagonal elements)."
> >
> >     Any multiple of the identity matrix is technically compound
> symmetric,
> >     because all the off-diagonal elements are the same (0).
> >
> >     Phillip
> >
> >     On 28/11/20 2:30 am, Simon Harmel wrote:
> >     > Hello All,
> >     >
> >     > Below, I'm using corCompSymm() (compound symmetry) for my simple
> >     model.
> >     >
> >     > The rho is estimated to be 0. I was wondering what it means for
> >     rho in the
> >     > var-covariance matrix to be "0"? Is my var-covariance matrix below
> >     valid?
> >     > -- Thank you all, Simon
> >     > #----------------------------------------------------------------
> >     > library(nlme)
> >     > data <-
> >     read.csv('https://raw.githubusercontent.com/hkil/m/master/R.csv')
> >     >
> >     > m <- lme(Achieve ~ time, random = ~1|subid, data = data,
> correlation =
> >     > corCompSymm())
> >     >
> >     >   aa <- corMatrix(m$modelStruct$corStruct)[[1]]
> >     >   aa * sigma(m)^2
> >     >
> >     >          [,1]     [,2]     [,3]     [,4]
> >     > [1,] 112.5003   0.0000   0.0000   0.0000
> >     > [2,]   0.0000 112.5003   0.0000   0.0000
> >     > [3,]   0.0000   0.0000 112.5003   0.0000
> >     > [4,]   0.0000   0.0000   0.0000 112.5003
> >     >
> >     >       [[alternative HTML version deleted]]
> >     >
> >     > _______________________________________________
> >     > R-sig-mixed-models@r-project.org
> >     <mailto: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
[prev in list] [next in list] [prev in thread] [next in thread] 