[R-sig-ME] Fitting different error structures with lme()

Mon Feb 18 10:17:38 CET 2013

Dear Cengiz,

IMHO you are overfitting the model. You have only 3 occasions per person? Note that a model with a random intercept and correlated random slope requires 3 parameters to be fit.

The first argument of corAR1() and corCompSymm() is value. It is used a starting values for the parameters. Unless use add fixed = TRUE, then it fixed at the those values. This is described in ?corAR1

You might want to look at the dataset shipped with nlme. Have a look at the examples in the helpfiles.

Best regards,

Thierry

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
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

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

-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces op r-project.org [mailto:r-sig-mixed-models-bounces op r-project.org] Namens Cengiz Zopluoglu
Verzonden: zondag 17 februari 2013 7:09
Aan: r-sig-mixed-models op r-project.org
Onderwerp: [R-sig-ME] Fitting different error structures with lme()

Hi,

I have a question about fitting a random coefficients model with different error structures in lme() function. The data I have in the long format as the following:

     id gender Occasion Strength
1.0   1   Male        0      161
1.1   1   Male        1      210
1.2   1   Male        2      230
2.0   2   Male        0      215
2.1   2   Male        1      245
2.2   2   Male        2      265
3.0   3   Male        0      134
3.1   3   Male        1      215
................................................................
................................................................

When I fit the following model by imposing a compound symmetry structure on the error var-cov matrix:

fit <- lme(Strength ~ 1 + Occasion,
           random = ~ 1 + Occasion | id,
           data = grip.L,
           correlation=corCompSymm(),
           control=list(maxIter=1000,msMaxIter=1000, niterEM=1000)
           )

Linear mixed-effects model fit by REML
  Data: grip.L
  Log-restricted-likelihood: -873.9892
  Fixed: Strength ~ 1 + Occasion
(Intercept)    Occasion
144.8304598  -0.6637931

Random effects:
 Formula: ~1 + Occasion | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr
(Intercept) 72.69465 (Intr)
Occasion    16.14353 -0.366
Residual    16.45378

Correlation Structure: Compound symmetry
 Formula: ~1 | id
 Parameter estimate(s):
Rho
  0
Number of Observations: 174
Number of Groups: 58

The "Rho" is estimated as zero. When I specify corCompSymm(.3), the output
is:

Linear mixed-effects model fit by REML
  Data: grip.L
  Log-restricted-likelihood: -873.9892
  Fixed: Strength ~ 1 + Occasion
(Intercept)    Occasion
144.8304598  -0.6637931

Random effects:
 Formula: ~1 + Occasion | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr
(Intercept) 72.02091 (Intr)
Occasion    16.14343 -0.369
Residual    19.18853

Correlation Structure: Compound symmetry
 Formula: ~1 | id
 Parameter estimate(s):
      Rho
0.2647256
Number of Observations: 174
Number of Groups: 58

The rho is estimated .264. In this case, I am wondering what really .3 within
corCompSymm(.3) represents. Or, when I say corCompSymm(.3) , what does it me an?

Also, when I specify corAR1(), I am getting an error message about convergence.

fit <- lme(Strength ~ 1 + Occasion,
           random = ~ 1 + Occasion | id,
           data = grip.L,
           correlation=corAR1(),
           control=list(maxIter=1000,msMaxIter=1000, niterEM=1000)
           )

Error in lme.formula(Strength ~ 1 + Occasion, random = ~1 + Occasion |  :
  nlminb problem, convergence error code = 1
  message = singular convergence (7)

After a trial-error process, I was able to fit the model by specifying
corAR1(.05):

Linear mixed-effects model fit by REML
 Data: grip.L
       AIC      BIC    logLik
  1761.052 1783.084 -873.5258

Random effects:
 Formula: ~1 + Occasion | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr
(Intercept) 39.491786 (Intr)
Occasion     5.554041 -0.999
Residual    62.765655

Correlation Structure: AR(1)
 Formula: ~1 | id
 Parameter estimate(s):
      Phi
0.9021773
Fixed effects: Strength ~ 1 + Occasion
                Value Std.Error  DF   t-value p-value
(Intercept) 143.52300  9.736396 115 14.740875  0.0000
Occasion     -0.66379  2.617492 115 -0.253599  0.8003
 Correlation:
         (Intr)
Occasion -0.396

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-1.4632691 -0.6942647 -0.1709512  0.4677042  2.0371593

Number of Observations: 174
Number of Groups: 58

Here, the estimated Phi is .9021. I have the same question again. What (.05)represents within corAR1(.05)? How is it related to the estimated phi .9021? Should I always do a trial-error process to figure out what value I should put in paranthesis to be able fit a model?

I also tried other error structures for the same dataset, but none of them converged. I used a different larger dataset, but I always get convergence error when I try to fit different error structures. Any suggestions?

This is for didactic purposes. We would like to teach students how to fit different error structures for mixed effects models using lme() function, and I can't get it worked until now for none of the datasets I have.

Thanks in advance.

--

Cengiz Zopluoglu

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models op r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.