[R-sig-ME] gls vs lme covariance structures
Joshua Wiley
jwiley.psych at gmail.com
Sat May 5 17:05:26 CEST 2012
Hi Charles,
I know you were not trying to make it difficult, no worries. Most of
the results seem to match fairly closely. Here is some stuff you can
try if you want (also added some residual plots and examination of the
distribution of random effects which is assumed gaussian).
Cheers,
Josh
####################################################
dat=read.table("subset.csv",sep=",",header=TRUE, na.strings=".")
dat34=dat[dat$group %in% c("3", "4"),]
dat34=within(dat34, {
group=factor(group)
Event_name=factor(Event_name)
Died=factor(Died)
ID=factor(ID)
})
contrasts(dat34$Event_name)=contr.sum(n=6)
contrasts(dat34$group)=contr.sum(n=2)
contrasts(dat34$Died)=contr.sum(n=2)
dat34$rep <- unlist(sapply(rle(as.numeric(dat34$ID))$lengths, seq.int))
require(nlme)
require(lattice)
m.run <- lme(var ~ group + Event_name + Died,
random = ~ 1 | ID, corr = corSymm(), data = dat34)
## store residuals and plot
dat34$r.mrun <- resid(m.run)
xyplot(r.mrun ~ rep | ID, data = dat34)
## examine distribution of random effects
qqmath(ranef(m.run)[[1]])
summary(m.run)
## Linear mixed-effects model fit by REML
## Data: dat34
## AIC BIC logLik
## 47.06365 116.5417 1.468174
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.3934213 0.2635925
## Correlation Structure: General
## Formula: ~1 | ID
## Parameter estimate(s):
## Correlation: <-- you can see the unstructure correlation matrix
## 1 2 3 4 5
## 2 0.289
## 3 0.553 0.774
## 4 0.225 0.532 0.752
## 5 0.532 0.635 0.844 0.662
## 6 0.717 -0.052 0.361 -0.053 0.535
anova(m.run, type = "marginal", adjustSigma = FALSE)
## numDF denDF F-value p-value
## (Intercept) 1 96 9.816408 0.0023
## group 1 23 0.131949 0.7197
## Event_name 5 96 1.081781 0.3754
## Died 1 23 0.074427 0.7874
#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
## group*Event_name expands to: group + Event_name + group:Event_name
## i.e., main effects + interaction, like "group | Event_name in proc mixed
m.run2 <- lme(var2 ~ group*Event_name + Died,
random = ~ 1 | ID, corr = corSymm(), data = dat34)
dat34$r.run2 <- resid(m.run2)
xyplot(r.run2 ~ rep | ID, data = dat34)
## examine distribution of random effects
qqmath(ranef(m.run2)[[1]])
summary(m.run2)
## auto regressive model
m.rar1 <- lme(var2 ~ group*Event_name + Died,
random = ~ 1 | ID, corr = corAR1(), data = dat34)
## store residuals and plot
dat34$r.rar1 <- resid(m.rar1)
xyplot(r.rar1 ~ rep | ID, data = dat34)
## examine distribution of random effects
qqmath(ranef(m.rar1)[[1]])
summary(m.rar1)
## Linear mixed-effects model fit by REML
## Data: dat34
## AIC BIC logLik
## 413.8188 457.598 -190.9094
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.383367 1.004379
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.1092874 <-- phi (rho) autocorrelation estimate
anova(m.rar1, type = "marginal", adjustSigma = FALSE)
## numDF denDF F-value p-value
## (Intercept) 1 91 337.2176 <.0001
## group 1 23 0.9015 0.3523
## Event_name 5 91 15.9816 <.0001
## Died 1 23 1.4741 0.2370
## group:Event_name 5 91 3.4790 0.0064
## 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
## compare auto regressive with unstructured
anova(m.run2, m.rar1)
####################################################
On Fri, May 4, 2012 at 6:00 AM, Charles Determan Jr <deter088 at umn.edu> wrote:
> 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/
--
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