[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