using lme for cov structures?

Dear Thierry and Wolfgang,

Thanks for responding so quickly. Here is the reproducible example of
the two models that I'm interested in.

# Read data.
opposites <-
read.table("https://stats.idre.ucla.edu/stat/r/examples/
alda/dat/opposites_pp.txt",header=TRUE,sep=",")

# Open library.
library(nlme)

# Define group "one" with value 1 for all cases.
one <- rep(1, 140)

#----- Model estimated with gls.
comsym.gls <- gls(opp~1,opposites,
                   correlation=corCompSymm(form = ~ 1 |id), method="REML")
summary(comsym.gls)


#----- Same model estimated with lme.
comsym.lme <- lme(opp~1,opposites,
                   random= ~1|one,
                   correlation=corCompSymm(form = ~ 1 |one/id),
method="REML")
summary(comsym.lme)


#----- Wolfgang's gls suggestion for heterogeneous CS.

summary(gls(opp ~ 1, opposites, correlation = corCompSymm(form = ~ 1 |
id), weights = varIdent(form = ~ 1 | wave)))

# Does not work with lme.
summary(hetercom <- lme(opp ~ 1,opposites,
                 correlation=corCompSymm(form = ~ 1 |id),
                 weights=varIdent(form = ~1|wave), method="REML"))

# But does work with the "one" trick.
summary(lme(opp ~ 1,opposites,
     random = ~1|one,
     correlation=corCompSymm(form = ~ 1 |one/id),
     weights=varIdent(form = ~1|wave), method="REML"))



The main reason of my mailing to the list is this. I was wondering
whether with lme it is possible to also estimate models with the
individuals nested in higher levels like schools or hospitals etc and at
the same time letting the residuals correlate within individuals over
time with one of the nice covar structures. However, I do NOT have a
reproducible example of such more complex models at the moment. I only
hoped that if the lme version of the model presented above has no
further problems than a (slightly) different AIC, BIC and std. error of
the fixed intercept, it would be meaningful to proceed in this way, i.e.
using lme instead of gls, since lme provides the important possibility
of additional random effects in the model.

As to Wolfgang's suggestion for heretogeneous CS using:

gls(opp ~ 1, correlation = corCompSymm(form = ~ 1 | id), weights =
varIdent(form = ~ 1 | timepoint))

I didn't find a way to estimate such a model with lme, except for the
method with the "trick". Using:

lme(opp ~ 1,opposites,
             correlation=corCompSymm(form = ~ 1 |id),
             weights=varIdent(form = ~1|wave), method="REML")

leads to error-message:

Error in lme.formula(opp ~ 1, opposites, correlation = corCompSymm(form
= ~1 | : incompatible formulas for groups in 'random' and 'correlation'


whereas

lme(opp ~ 1,opposites,
             random = ~1|one,
             correlation=corCompSymm(form = ~ 1 |one/id),
             weights=varIdent(form = ~1|wave), method="REML")

leads to similar results as your gls suggestion.


Best regards, Ben.


On 26-7-2017 14:44, Thierry Onkelinx wrote:

Dear Ben,

Please be more specific on the kind of model you want to fit. That
would lead to a more relevant discussion that your potential misuse of
lme. A reproducible example is always useful...

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

2017-07-26 13:36 GMT+02:00 Ben Pelzer <b.pelzer at maw.ru.nl
<mailto:b.pelzer at maw.ru.nl>>:

    Dear list,

    With longitudinal data, the package nlme offers the possibility to
    specify particular covariance structures for the residuals. In the
    examples below I used data from 35 individuals measured at 4 time
    points, so we have 4 observations nested in each of 35 individuals.
    These data are discussed in Singer and Willett, Applied Longitudinal
    Data Analysis, chapter 7.

    In several sources I found that e.g. compound symmetry can easily be
    obtained with gls from package nlme, by using the correlation
    structure
    corCompSymm, as in:

         comsym <- gls(opp~1,opposites,
                       correlation=corCompSymm(form = ~ 1 |id),
    method="REML")

    where id is subject-identifier for the individual. With gls
    however one
    cannot specify random effects, as opposed to lme. To have lme

estimate

    compound symmetry is simple, one would not even need have to specify

a

    particular correlation structure, it would suffice to say
    "random=~1|id". However, for more complex covariance structures, e.g.
    heterogeneous compound symmetry, I was only able to find syntax
    for gls,
    but not for lme.

    Then I thought that tricking lme might be an option by having a
    kind of
    "fake" random effect. That is, I constructed a variable "one" which
    takes value 1 for all cases, and let the intercept be random
    across "all
    one groups". This led to the following lme model:


         comsym.lme <- lme(opp~1,opposites, random= ~1|one,
                           correlation=corCompSymm(form = ~ 1 |one/id),
    method="REML")

    And actually to my surprise, the results of this lme are highly
    similar
    to those of the gls above.
    The output of both is attached below. The loglik's are identical, the
    AIC and BIC are not, which I can understand, as the lme has an extra
    variance to be estimated, compared to the gls. Also, the fixed
    intercept's standard error is different, the one of the lme being
    about
    twice as large.

    I added some independent variables (not shown below) but the
    similarities between gls and lme remain, with only the AIC and BIC

and

    the standard error of the fixed intercept being different for the two
    models.

    My question is threefold.
    1) Suppose the individuals (say pupils) would be nested in
    schools, then
    I assume that with lme I could add school as a random effect, and
    run a
    "usual" model, with pupils nested in schools and observations in
    pupils,
    and then use any of the possible residual covariance structures
    for the
    observations in pupils. Would you agree with this? (with gls one
    cannot
    use an additional random effect, like e.g. school)
    2) Are the lme results indeed to be trusted when using this "fake"
    random effect, apart from the differences with gls mentioned?
    Could you
    imagine a situation where lme with this trick would produce wrong
    results?
    3) I don't understand the variance of the intercept in the lme

output.

    Even when I specify a very simple model: lme(opp ~ 1, random =

~1|one,

    opposites, method="REML")
    lme is able to estimate the intercept variance...but what is this
    variance estimate expressing?

    Thanks a lot for any help!!!

    Ben.





    *----------------- gls
    ---------------------------------------------------

     > comsym <- gls(opp~1,opposites,

                     correlation=corCompSymm(form = ~ 1 |id),
    method="REML")

     > summary(comsym)

    Generalized least squares fit by REML
    Model: opp ~ 1
    Data: opposites
          AIC      BIC    logLik
    1460.954 1469.757 -727.4769

    Correlation Structure: Compound symmetry
    Formula: ~1 | id
    Parameter estimate(s):
       Rho
    0.2757052

    Coefficients:
                    Value Std.Error  t-value p-value
    (Intercept) 204.8143  5.341965 38.34063       0

    Standardized residuals:
       Min          Q1         Med          Q3         Max
    -2.75474868 -0.71244027  0.00397158  0.56533908  2.24944156

    Residual standard error: 46.76081
    Degrees of freedom: 140 total; 139 residual


    #------------------ lme
    ---------------------------------------------------

     > comsym.lme <- lme(opp~1,opposites, random= ~1|one,

                     correlation=corCompSymm(form = ~ 1 |one/id),
    method="REML")

     > summary(comsym.lme)

    Linear mixed-effects model fit by REML
    Data: opposites
          AIC      BIC    logLik
    1462.954 1474.692 -727.4769

    Random effects:
       Formula: ~1 | one
             (Intercept) Residual
    StdDev:    10.53875 46.76081

    Correlation Structure: Compound symmetry
    Formula: ~1 | one/id
    Parameter estimate(s):
       Rho
    0.2757052
    Fixed effects: opp ~ 1
                    Value Std.Error  DF  t-value p-value
    (Intercept) 204.8143  11.81532 139 17.33463       0

    Standardized Within-Group Residuals:
       Min          Q1         Med          Q3         Max
    -2.75474868 -0.71244027  0.00397158  0.56533908  2.24944156

    Number of Observations: 140
    Number of Groups: 1


            [[alternative HTML version deleted]]

        [[alternative HTML version deleted]]