[R-sig-ME] gls vs lme covariance structures

Fri May 4 15:18:39 CEST 2012

I just wanted to point out that for simple repeated measures, longitudinal data, it --is-- possible to mimic some simple reasonable lme models with gls.

If lme uses the random effects specification

	random = ~1 | Subject

Then gls can mimic that by using 

	correlation = corCompSymm(form = ~1 | Subject)

The resulting parameterizations are different, but they give the same predictions and SEs for the "fixed"/model effects.

-----

I often take this gls paradigm a step further and make it a habit to use

	correlation = corrExp(form = LongitudinalTimeVariable | Subject, nugget = TRUE)

In the limit of the range parameter going to +infinity, the model reduces to the above corCompSymm model equivalent to an lme model.  In the limit as the nugget parameter going to zero, the model reduces to a CAR1 model.  (see ?corrExp and other corXXX)  I see cases of both of these.  But what astonishes me is that how often this two parameter correlation model fits my data so well: it mirrors nearly exactly what a full unrestricted correlation structure, using corSymm, gives.  (of course, that is my data, YMMV.)  Furthermore, the product of the nugget parameter and the overall variance gives an estimate of cross sectional replication variance, a nicety since I don't believe one can use replicated data with the corXXX functions.  

John
John  Szumiloski,  Ph.D.

Senior Biometrician
Biometrics Research
WP53B-120
Merck Research Laboratories
P.O. Box 0004
West Point, PA 19486-0004

>(215) 652-7346 (PH)
>(215) 993-1835 (FAX)
>

-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of r-sig-mixed-models-request at r-project.org
Sent: Thursday, May 03, 2012 5:17 PM
To: r-sig-mixed-models at r-project.org
Subject: R-sig-mixed-models Digest, Vol 65, Issue 12

Send R-sig-mixed-models mailing list submissions to
	r-sig-mixed-models at r-project.org

To subscribe or unsubscribe via the World Wide Web, visit
	https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
	r-sig-mixed-models-request at r-project.org

You can reach the person managing the list at
	r-sig-mixed-models-owner at r-project.org

When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-mixed-models digest..."

Today's Topics:

   1. Re: gls vs lme covariance structures (Joshua Wiley)
   2. Re: Nested subject-longitudinal logit design (arun)
   3. Re: Extracting variances of the estimated variance components
      in lme4 (Joshua Wiley)

----------------------------------------------------------------------

Message: 1
Date: Thu, 3 May 2012 13:59:15 -0700
From: Joshua Wiley <jwiley.psych at gmail.com>
To: Charles Determan Jr <deter088 at umn.edu>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] gls vs lme covariance structures
Message-ID:
	<CANz9Z_JYcSput9=kB8ibkvooKcABZWfTRDm4QZU9tmYYcxnOHQ at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8

Hi Charles,

Well you could post a subset of it, or make up some data that is sharable (whether the data make any sense is not important to us, just nice to have runable code, for example your previous thread about contrasts could have been solved in one email if we could have shown you how to set the contrasts on your data and then it matched your SAS output).  In any case, whether you use lme or gls really depends on your question and goals, I think.  Generalized least squares is not the same as a random effects model.  If you want a random effect, you cannot use gls.  If you just want correlated errors, gls is fine.

This part of your code strikes me as atypical though I cannot promise it is wrong/not what you want: corr=corAR1(ID)

Cheers,

Josh

On Thu, May 3, 2012 at 1:44 PM, Charles Determan Jr <deter088 at umn.edu> wrote:
> Hi Joshua,
>
> Thanks for your response.? It is probably best that I don't post the 
> data as some of it is not yet published.? My main question is whether 
> UN and AR(1) should be done with gls or if I have done the syntax incorrectly with lme.
> Since AR(1) is replicated perfectly if I put the correct dendf, I can 
> work with it.? And UN is close, so I just want to be sure my use and 
> syntax are correct, not necessarily modifying the data.
>
> Regards,
> Charles
>
>
> On Thu, May 3, 2012 at 3:21 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
>>
>> Hi Charles,
>>
>> Could you upload the dataset you are using somewhere and post the 
>> link? ?Something like:
>>
>> ##########
>> dat34 <- read.csv("http:/wherever/you/uploaded/yourdata.csv", header 
>> =
>> TRUE)
>> ## code to convert to factors anything that needs to be etc.
>> ##########
>>
>> Then it is easier for us to try things that way.
>>
>> corAR1 and corSymm seem appropriate. ?Have you checked the examples 
>> in their documents? ?I found them helpful.
>>
>> Cheers,
>>
>> Josh
>>
>> On Thu, May 3, 2012 at 12:58 PM, Charles Determan Jr 
>> <deter088 at umn.edu>
>> wrote:
>> > Greetings R users,
>> >
>> > I have been attempting to replicate various covariance structures 
>> > from SAS's analysis of mixed models. ?I have successfully been able 
>> > to replicate compound symmetry, however it becomes confusing with 
>> > autoregression and unstructured. ?As such, there are two questions 
>> > regarding this issue.
>> >
>> > Autoregression
>> > SAS output (Type III fixed effects) for covariance structure AR(1)
>> >
>> > *Type 3 Tests of Fixed Effects*
>> >
>> > *Effect ? ? ? ? ? ? ? ? ? ?NumDF DenDF F Value Pr > F*
>> >
>> > *group ? ? ? ? ? ? ? ? ? ? ? ? ?*1 ? ? ? ? 23 ? ? ? ? ?0.99 ? ? 
>> > ?0.3293
>> >
>> > *Event_name ? ? ? ? ? ? ? *5 ? ? ? ? 91 ? ? ? ? 16.23 ? ?<.0001
>> >
>> > *Died ? ? ? ? ? ? ? ? ? ? ? ? ? ?*1 ? ? ? ? 23 ? ? ? ? ?1.70 ? ? 
>> > ?0.2047
>> >
>> > *group*Event_name ?*5 ? ? ? ?91 ? ? ? ? ?3.04 ? ? 0.0140
>> >
>> > R output (corAR1=AR(1)?)
>> > I can replicate these results if I run the following:
>> >
>> > fit.18=gls(var~group+Event_name+Died+group*Event_name,
>> > ? ?data=dat34,
>> > ? ?corr=corAR1(, ~1|ID),
>> > ? ?weight=varIdent(~1|Event_name))
>> > anova(fit.18, type="marginal", adjustSigma=F)
>> >
>> > #the DenDF are off with gls, so use the 'correct' ones 1-pf(.9935, 
>> > 1, 23) 1-pf(16.2323, 5, 91) 1-pf(1.7041, 1, 23) 1-pf(3.0367, 5, 91) 
>> > #and the output matches exactly
>> >
>> > However, I can not get the lme function to run the autoregression. 
>> > ?The output is very different:
>> >
>> > fit.11=lme(var~group+Event_name+Died+group*Event_name,
>> > ? ?data=dat34,
>> > ? ?random=~1|ID,
>> > ? ?corr=corAR1(ID),
>> > ? ?weight=varIdent(~1|Event_name))
>> > anova.lme(fit.11, type="marginal", adjustSigma=F)
>> >
>> > ? ? ? ? ? ? ? ? ? ? ? numDF denDF ?F-value ? ? ?p-value
>> > (Intercept) ? ? ? ? ? ? ? 1 ? ?96 ? ? ? ?9.816419 ?0.0023 group ? ? 
>> > ? ? ? ? ? ? ? ? ?1 ? ?23 ? ? ? 0.131950 ?0.7197 Event_name ? ? ? ? 
>> > ? ?5 ? ?96 ? ? ? ?1.081785 ?0.3754 Died ? ? ? ? ? ? ? ? ? ? ? 1 ? 
>> > ?23 ? ? ? ?0.074428 ?0.7874
>> >
>> > Is this type of covariance structure only done with gls and I 
>> > should continue with the analysis as such or am I doing something 
>> > silly with lme?
>> >
>> > My second question is with regards to the unstructured covariance.
>> > SAS output (UN)
>> >
>> > *Type 3 Tests of Fixed Effects*
>> >
>> > *Effect ? ? ? ?NumDF DenDF F Value Pr > F*
>> >
>> > *Pig_group ? ?*1 ? ? ? ? ? ?23 ? ? ? ? 2.73 ? ?0.1120
>> >
>> > *Event_name *5 ? ? ? ? ? ?23 ? ? ? ? 1.11 ? ?0.3806
>> >
>> > *Died ? ? ? ? ? ? ?*1 ? ? ? ? ? ?23 ? ? ? ? 0.51 ? ?0.4833
>> >
>> > R output (corSymm = UN?)
>> > fit.11=lme(var2~group+Event_name+Died,
>> > ? ?data=dat34,
>> > ? ?random=~1|ID,
>> > ? ?corr=corSymm(, ~1|ID),
>> > ? ?weight=varIdent(~1|Event_name))
>> > anova.lme(fit.11, type="marginal", adjustSigma=F)
>> >
>> > #same as corAR1???
>> > ? ? ? ? ? ? ? ? ? ? ? numDF denDF ?F-value ? ? ?p-value
>> > (Intercept) ? ? ? ? ? ? ? 1 ? ?96 ? ? ? ?9.816419 ?0.0023 group ? ? 
>> > ? ? ? ? ? ? ? ? ?1 ? ?23 ? ? ? 0.131950 ?0.7197 Event_name ? ? ? ? 
>> > ? ?5 ? ?96 ? ? ? ?1.081785 ?0.3754 Died ? ? ? ? ? ? ? ? ? ? ? 1 ? 
>> > ?23 ? ? ? ?0.074428 ?0.7874
>> >
>> > but with gls
>> > fit.18=gls(var~group+Event_name+Died,
>> > ? ?data=dat34,
>> > ? ?corr=corSymm(~1|ID),
>> > ? ?weight=varIdent(~1|Event_name))
>> > anova(fit.18, type="marginal", adjustSigma=F)
>> >
>> > 1-pf(2.869837, 1, 23)
>> > 1-pf(1.126747, 5, 23)
>> > 1-pf(.514726, 1, 23)
>> >
>> > [1] 0.1037549
>> > [1] 0.3742309
>> > [1] 0.4803239
>> >
>> > #close but not exact (however I can work with this if it is indeed
>> > correct)
>> >
>> > Overall, I want to clarify the difference between gls and lme and 
>> > if I am simply making some weird syntax error with lme that I can't 
>> > seem to get the covariance structures to match.
>> >
>> >
>> > Apologies for lots of information in one go, but hopefully this 
>> > provides necessary information to point me in the correct direction.
>> >
>> > Thanks to any and all who give their time to these questions,
>> >
>> > Regards,
>> > Charles
>> >
>> > ? ? ? ?[[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > R-sig-mixed-models at r-project.org mailing list 
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>> --
>> Joshua Wiley
>> Ph.D. Student, Health Psychology
>> Programmer Analyst II, Statistical Consulting Group University of 
>> California, Los Angeles https://joshuawiley.com/
>
>

-- 
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/

[...snip...]

------------------------------

_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

End of R-sig-mixed-models Digest, Vol 65, Issue 12
**************************************************
Notice:  This e-mail message, together with any attachme...{{dropped:11}}