[R-sig-ME] gls vs lme covariance structures
Joshua Wiley
jwiley.psych at gmail.com
Thu May 3 22:59:15 CEST 2012
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/
More information about the R-sig-mixed-models
mailing list