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

[Joshua Wiley]

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/