[R-sig-ME] lme interaction result strange
Ista Zahn
istazahn at gmail.com
Wed May 2 16:41:56 CEST 2012
I still think the difference is in the default contrasts. Try this
example, and notice that the default output differs for Age, and that
the same results are obtained when the contrasts in R are set to
contr.sum. I would be very surprised if this is not the source of your
confusion.
* SAS
input Person Gender $ y1 y2 y3 y4;
y=y1; Age=8; output;
y=y2; Age=10; output;
y=y3; Age=12; output;
y=y4; Age=14; output;
drop y1-y4;
datalines;
1 F 21.0 20.0 21.5 23.0
2 F 21.0 21.5 24.0 25.5
3 F 20.5 24.0 24.5 26.0
4 F 23.5 24.5 25.0 26.5
5 F 21.5 23.0 22.5 23.5
6 F 20.0 21.0 21.0 22.5
7 F 21.5 22.5 23.0 25.0
8 F 23.0 23.0 23.5 24.0
9 F 20.0 21.0 22.0 21.5
10 F 16.5 19.0 19.0 19.5
11 F 24.5 25.0 28.0 28.0
12 M 26.0 25.0 29.0 31.0
13 M 21.5 22.5 23.0 26.5
14 M 23.0 22.5 24.0 27.5
15 M 25.5 27.5 26.5 27.0
16 M 20.0 23.5 22.5 26.0
17 M 24.5 25.5 27.0 28.5
18 M 22.0 22.0 24.5 26.5
19 M 24.0 21.5 24.5 25.5
20 M 23.0 20.5 31.0 26.0
21 M 27.5 28.0 31.0 31.5
22 M 23.0 23.0 23.5 25.0
23 M 21.5 23.5 24.0 28.0
24 M 17.0 24.5 26.0 29.5
25 M 22.5 25.5 25.5 26.0
26 M 23.0 24.5 26.0 30.0
27 M 22.0 21.5 23.5 25.0
;
proc mixed data=pr method=ml covtest;
class Person Gender;
model y = Age Gender Gender*Age/solution;
random intercept/type=un subject=Person;
run;
Type 3 Tests of Fixed Effects
Effect Num DF Den DF F Value Pr > F
Age 1 79 111.10 <.0001
Gender 1 79 0.47 0.4960
Age*Gender 1 79 6.46 0.0130
# R
dat <- read.table(text = 'Person,Gender,y,Age
1,F,21,8
1,F,20,10
1,F,21.5,12
1,F,23,14
2,F,21,8
2,F,21.5,10
2,F,24,12
2,F,25.5,14
3,F,20.5,8
3,F,24,10
3,F,24.5,12
3,F,26,14
4,F,23.5,8
4,F,24.5,10
4,F,25,12
4,F,26.5,14
5,F,21.5,8
5,F,23,10
5,F,22.5,12
5,F,23.5,14
6,F,20,8
6,F,21,10
6,F,21,12
6,F,22.5,14
7,F,21.5,8
7,F,22.5,10
7,F,23,12
7,F,25,14
8,F,23,8
8,F,23,10
8,F,23.5,12
8,F,24,14
9,F,20,8
9,F,21,10
9,F,22,12
9,F,21.5,14
10,F,16.5,8
10,F,19,10
10,F,19,12
10,F,19.5,14
11,F,24.5,8
11,F,25,10
11,F,28,12
11,F,28,14
12,M,26,8
12,M,25,10
12,M,29,12
12,M,31,14
13,M,21.5,8
13,M,22.5,10
13,M,23,12
13,M,26.5,14
14,M,23,8
14,M,22.5,10
14,M,24,12
14,M,27.5,14
15,M,25.5,8
15,M,27.5,10
15,M,26.5,12
15,M,27,14
16,M,20,8
16,M,23.5,10
16,M,22.5,12
16,M,26,14
17,M,24.5,8
17,M,25.5,10
17,M,27,12
17,M,28.5,14
18,M,22,8
18,M,22,10
18,M,24.5,12
18,M,26.5,14
19,M,24,8
19,M,21.5,10
19,M,24.5,12
19,M,25.5,14
20,M,23,8
20,M,20.5,10
20,M,31,12
20,M,26,14
21,M,27.5,8
21,M,28,10
21,M,31,12
21,M,31.5,14
22,M,23,8
22,M,23,10
22,M,23.5,12
22,M,25,14
23,M,21.5,8
23,M,23.5,10
23,M,24,12
23,M,28,14
24,M,17,8
24,M,24.5,10
24,M,26,12
24,M,29.5,14
25,M,22.5,8
25,M,25.5,10
25,M,25.5,12
25,M,26,14
26,M,23,8
26,M,24.5,10
26,M,26,12
26,M,30,14
27,M,22,8
27,M,21.5,10
27,M,23.5,12
27,M,25,14',
header = TRUE,
sep = ",")
dat$Person <- factor(dat$Person)
library(nlme)
fit1 <- lme(y~Gender*Age,
data=dat,
method = "ML",
random=~1|Person)
anova(fit1, type="marginal", adjustSigma=F)
contrasts(dat$Gender) <- contr.sum(n = 2)
fit2 <- lme(y~Gender*Age,
data=dat,
method = "ML",
random=~1|Person)
anova(fit2, type="marginal", adjustSigma=F)
On Wed, May 2, 2012 at 9:31 AM, Charles Determan Jr <deter088 at umn.edu> wrote:
> Thank you Thierry,
>
> To answer your questions, I am sure it applies the same test as the simple
> model matches perfectly. Furthermore, two of the interaction model
> components match but the two 'middle' effects are very different. As such,
> I was wondering if there was something in the R syntax that would throw this
> off and haven't preoccupied myself with the calculation of the F-value as
> this is a multivariate dataset. The relevance here is that this is a
> practice set to see if I can replicate the mixed model analysis in R.
>
> Regards,
> Charles
>
>
> On Wed, May 2, 2012 at 8:24 AM, ONKELINX, Thierry <Thierry.ONKELINX at inbo.be>
> wrote:
>>
>> Dear Charles,
>>
>> Are you sure that SAS and lme apply the same tests? So you'll need to know
>> how the F-values are calculated. Futhermore: what is the relevance of
>> testing the effect of a main effect when an interaction with that main
>> effect is present?
>>
>> Best regards,
>>
>> Thierry
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
>> Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>> + 32 2 525 02 51
>> + 32 54 43 61 85
>> Thierry.Onkelinx at inbo.be
>> www.inbo.be
>>
>> To call in the statistician after the experiment is done may be no more
>> than asking him to perform a post-mortem examination: he may be able to say
>> what the experiment died of.
>> ~ Sir Ronald Aylmer Fisher
>>
>> The plural of anecdote is not data.
>> ~ Roger Brinner
>>
>> The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of data.
>> ~ John Tukey
>>
>>
>> -----Oorspronkelijk bericht-----
>> Van: r-sig-mixed-models-bounces at r-project.org
>> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens Charles Determan Jr
>> Verzonden: woensdag 2 mei 2012 15:14
>> Aan: Ista Zahn
>> CC: r-sig-mixed-models at r-project.org
>> Onderwerp: Re: [R-sig-ME] lme interaction result strange
>>
>> I have checked the contrasts and 'unfortunately' they are identical. Is
>> that a statistically correct avenue to just take the interaction term from
>> the interaction model and use the first ones? That seems odd to me as the
>> numbers aren't the same even in the SAS output between the models. I
>> certainly appreciate your prior suggestion but perhaps you or someone else
>> has another idea? Something else seems to be strange, perhaps the syntax
>> needs to be different when doing an interaction term?
>>
>> Thanks,
>> Charles
>>
>> On Tue, May 1, 2012 at 2:02 PM, Ista Zahn <istazahn at gmail.com> wrote:
>>
>> > Hi Charles,
>> >
>> > Ask yourself what the Event_name and Died terms represent in this
>> > model. When you understand that, you'll understand why you need to
>> > know what contrasts were used if you hope to correctly interpret these
>> > terms.
>> >
>> > Alternatively, you can interpret the Event_name and Died terms from
>> > the first model (without the interaction term), and interpret just the
>> > interaction term from this model (so-called type II sums of squares by
>> > some).
>> >
>> > Best,
>> > Ista
>> >
>> > On Tue, May 1, 2012 at 2:54 PM, Charles Determan Jr <deter088 at umn.edu>
>> > wrote:
>> > > I have yet to attempt the contrast suggestion yet but here is the
>> > > SAS
>> > output
>> > > to complement the R.
>> > >
>> > > SAS output for the simple and interaction models. I have presented
>> > > the
>> > Type
>> > > III tests for simplicity. As you can see from my prior R output,
>> > > the DFs match exactly. The only variation is with the F and p value
>> > > of the Event_name and Died in the interaction model.
>> > >
>> > > Type 3 Tests of Fixed Effects
>> > >
>> > > Effect NumDF DenDF F Value Pr > F
>> > >
>> > > group 1 23 0.65 0.4293
>> > >
>> > > Event_name 5 96 2.09 0.0738
>> > >
>> > > Died 1 23 1.83 0.1889
>> > >
>> > >
>> > >
>> > > Type 3 Tests of Fixed Effects
>> > >
>> > > Effect NumDF DenDF F Value Pr > F
>> > >
>> > > group 1 23 0.47
>> > > 0.5015
>> > >
>> > > Event_name 5 91 2.62 0.0290
>> > >
>> > > Died 1 23 2.44
>> > > 0.1318
>> > >
>> > > Event_name*Died 5 91 1.10 0.3637
>> > >
>> > >
>> > > R output of interaction model
>> > >
>> > > numDF denDF F-value p-value
>> > > (Intercept) 1 91 111.20483 <.0001
>> > >
>> > > group 1 23 0.46632 0.5015
>> > > Event_name 5 91 1.14042 *0.3449*
>> > > Died 1 23 0.50989 * 0.4824*
>> > > Event_name:Died 5 91 1.10436 0.3637
>> > > Done.
>> > >
>> > >
>> > > Regards,
>> > >
>> > > Charles
>> > >
>> > >
>> > >
>> > > On Tue, May 1, 2012 at 1:39 PM, Thompson,Paul
>> > > <Paul.Thompson at sanfordhealth.org> wrote:
>> > >>
>> > >> There is another issue. From the error df, it seems like this is a
>> > >> multi-level/RM/multiple obs study, and SAS and R do not always
>> > >> agree on
>> > the
>> > >> computation of the df, as well as the type of SS that is being
>> > computed. You
>> > >> need to present both outputs, so that we can see both.
>> > >>
>> > >> I know almost nothing about R, and so my comments may not be
>> > >> relevant.
>> > >>
>> > >> -----Original Message-----
>> > >> From: r-sig-mixed-models-bounces at r-project.org
>> > >> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ista
>> > Zahn
>> > >> Sent: Tuesday, May 01, 2012 12:51 PM
>> > >> To: Charles Determan Jr
>> > >> Cc: r-sig-mixed-models at r-project.org
>> > >> Subject: Re: [R-sig-ME] lme interaction result strange
>> > >>
>> > >> Hi Charles,
>> > >>
>> > >> My first guess is that you have (a) categorical variable(s) in your
>> > >> predictors, and that the contrasts in SAS are different than those
>> > >> in R.
>> > >>
>> > >> Best,
>> > >> Ista
>> > >>
>> > >> On Tue, May 1, 2012 at 1:39 PM, Charles Determan Jr
>> > >> <deter088 at umn.edu>
>> > >> wrote:
>> > >> > Dear R users,
>> > >> >
>> > >> > I have been working on replicating some linear mixed models from
>> > >> > SAS.
>> > >> > The
>> > >> > first one matches perfectly when the SAS model is simple with the
>> > three
>> > >> > separate factors.
>> > >> >
>> > >> > fit=lme(var~group+Event_name+Died,
>> > >> > data=liv34,
>> > >> > random=~1|ID)
>> > >> > anova.lme(fit, type="marginal", adjustSigma=F)
>> > >> >
>> > >> > However, once I put an interaction into the formula the values
>> > >> > don't match.
>> > >> >
>> > >> > fit=lme(var~group+Event_name+Died+Event_name*Died,
>> > >> > data=liv34,
>> > >> > random=~1|ID)
>> > >> > anova.lme(fit, type="marginal", adjustSigma=F)
>> > >> >
>> > >> > numDF denDF F-value p-value
>> > >> > (Intercept) 1 91 111.20483 <.0001
>> > >> > group 1 23 0.46632 0.5015
>> > >> > Event_name 5 91 1.14042 *0.3449*
>> > >> > Died 1 23 0.50989 * 0.4824*
>> > >> > Event_name:Died 5 91 1.10436 0.3637
>> > >> > Done.
>> > >> >
>> > >> > The numbers *bold* don't match up. They should be approximately
>> > >> > .0290 and
>> > >> > .1318 respectively. The other two are still exact matches. I
>> > >> > know looking for exact matches is ambitious but the numbers
>> > >> > should be at least similar that the conclusions don't change so
>> > >> > drastically.
>> > >> >
>> > >> > Any thoughts as to why this discrepancy is happening would be
>> > >> > most appreciated.
>> > >> >
>> > >> > 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
>> > >>
>> > >> _______________________________________________
>> > >> R-sig-mixed-models at r-project.org mailing list
>> > >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> > >>
>> > >> -------------------------------------------------------------------
>> > >> ---- Confidentiality Notice: This e-mail message, including any
>> > >> attachments, is for the sole use of the intended recipient(s) and
>> > >> may contain privileged and confidential information. Any
>> > >> unauthorized review, use, disclosure or distribution is prohibited.
>> > >> If you are not the intended recipient, please contact the sender by
>> > >> reply e-mail and destroy all copies of the original message.
>> > >>
>> > >
>> >
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> * * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
>> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver
>> weer en binden het INBO onder geen enkel beding, zolang dit bericht niet
>> bevestigd is door een geldig ondertekend document.
>> The views expressed in this message and any annex are purely those of the
>> writer and may not be regarded as stating an official position of INBO, as
>> long as the message is not confirmed by a duly signed document.
>
>
More information about the R-sig-mixed-models
mailing list