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

Charles Determan Jr deter088 at umn.edu
Fri May 4 15:00:25 CEST 2012


Joshua,

My apologies, I never intended to make it more difficult for people to
help.  I only wish to be very careful when I every work with data.  I
have created a subset as you requested in the attached csv file.  The
two covariance structures I wish to do with a mixed model are (UN) for
var, and (AR(1)) for var2.

dat=read.table("C:/…/subset.csv",sep=",",header=TRUE, na.strings=".")
attach(dat)

dat34=dat[Pig_group %in% c("3", "4"),]
attach(dat34)
dat34=within(dat34, {
	group=factor(group)
	Event_name=factor(Event_name)
	Died=factor(Died)
	ID=factor(ID)
})
attach(dat34)

contrasts(dat34$Event_name)=contr.sum(n=6)
contrasts(dat34$group)=contr.sum(n=2)
contrasts(dat34$Died)=contr.sum(n=2)


The models I wish to fit with a mixed model are:
UN
var~group+Event_name+Died,
random=~1|ID
corr=corSymm()?
anova(fit, type=”marginal”, adjustSigma=F)

#Output should be: (SAS mixed procedure with 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

AR(1)
var2~group+Event_name+Died+group*Event_name
random=~1|ID
corr=corAR1()?
anova(fit, type=”marginal”, adjustSigma=F)

Output should be:
Type 3 Tests of Fixed Effects (SAS mixed procedure with AR(1))
Effect 		         	NumDF     DenDF F Value Pr > F
Pig_group                              1  	     23    0.99      0.3293
Event_name                           5  	     91    16.23    <.0001
Died                                      1	             23    1.70      0.2047
Pig_group*Event_name           5  	     91    3.04     0.0140


Many thanks, I truly appreciate the guidance,
Charles

On Thu, May 3, 2012 at 3:59 PM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
>
> 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