[R-sig-ME] 2nd attempt - conflict of dfs or f value in lme
Baldwin, Jim -FS
jbaldwin at fs.fed.us
Tue May 22 17:20:52 CEST 2012
Have you considered the other direction? SAS offers 5 different options for constructing the denominator degrees of freedom: BETWITHIN, CONTAIN, KENWARDROGER, RESIDUAL, SATTERTHWAITE. Also reading the associated references for those methods will help with understanding why there isn't a single method that is "correct" for all situations. (My apologies for showing part of a SAS manual on this list.)
Jim
Jim Baldwin
Station Statistician
Pacific Southwest Research Station
USDA Forest Service
Albany, California
-----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 Charles Determan Jr
Sent: Tuesday, May 22, 2012 8:01 AM
To: Douglas Bates
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] 2nd attempt - conflict of dfs or f value in lme
Dr. Bates,
After another thought, would it be possible to modify the gls source code to incorporate the same denDF calculation used in the lme function from this package? I have been looking through the source code of the nlme package and by adding the denDF to the gls function I would solve my initial problem of having consistent F values and dfs. However, I would prefer to work with one of the authors regarding any modification. Perhaps I should contact R core mailing list?
Regards,
Charles
On Mon, May 21, 2012 at 12:22 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> I don't really think there is a way of defining the degrees of freedom
> for all lme and gls model fits so that they correspond to the
> "correct" values, which often seems to mean the values returned by
> SAS. If you feel that SAS provides the correct analysis in all cases
> then you should use SAS. We can't reproduce the SAS results in R
> because SAS is proprietary code and we don't know how those results
> are calculated.
>
> The whole question of degrees of freedom for certain tests is very
> complicated. Often there is no answer because the distribution of the
> test statistic is not a t or an F distribution.
>
> On Mon, May 21, 2012 at 11:16 AM, Charles Determan Jr
> <deter088 at umn.edu>
> wrote:
> > Thank you for giving this some thought Joshua. I will be sure to
> > keep troubleshooting this as well. Would it be related to how gls
> > calculates its values that is somehow reports the 'correct' F value perhaps?
> >
> > Regards,
> > Charles
> >
> > On Mon, May 21, 2012 at 11:08 AM, Joshua Wiley
> ><jwiley.psych at gmail.com
> >wrote:
> >
> >> That suggests they are fitting different models somehow (not just
> >> that the specifc effect estimates are different). No brilliant
> >> ideas on in what way though...
> >>
> >> I checked REML versus ML but that was not it. I'll let you know if
> >> I have any ideas.
> >>
> >> Cheers,
> >>
> >> Josh
> >>
> >> On Mon, May 21, 2012 at 9:02 AM, Charles Determan Jr
> >> <deter088 at umn.edu>
> >> wrote:
> >> > Hi Joshua,
> >> >
> >> > The -2 Log-likelihood from SAS is -58.6 The log-likelihood from R
> >> > fit.13 is 21.9. I assume this would be
> >> multiplied
> >> > by -2 which would give -43.8.
> >> >
> >> > Regards,
> >> > Charles
> >> >
> >> >
> >> > On Mon, May 21, 2012 at 10:02 AM, Joshua Wiley <
> jwiley.psych at gmail.com>
> >> > wrote:
> >> >>
> >> >> Hi Charles,
> >> >>
> >> >> Out of curiosity, what is the log likelihood from SAS? Does it
> >> >> match that from lme()?
> >> >>
> >> >> logLik(fit.13)
> >> >>
> >> >> will give it to you.
> >> >>
> >> >> Cheers,
> >> >>
> >> >> Josh
> >> >>
> >> >> On Mon, May 21, 2012 at 5:45 AM, Charles Determan Jr <
> deter088 at umn.edu>
> >> >> wrote:
> >> >> > Greetings R users,
> >> >> >
> >> >> >
> >> >> > I am trying to match some SAS output of a mixed model. After
> running
> >> >> > PROC
> >> >> > MIXED with the covariance structure as AR1 the output below is
> given.
> >> >> > Now
> >> >> > when I try to replicate this with lme, I get the correct
> >> >> > degrees of freedom and somewhat close values. If I try with
> >> >> > gls, I get the exact F
> value
> >> >> > but
> >> >> > the incorrect denominator degrees of freedom. Is there some
> syntax or
> >> >> > parameter I can adjust to get lme to have the same F values as gls?
> >> >> > That
> >> >> > way the correct degrees of freedom would be applied and I
> >> >> > would get
> >> the
> >> >> > replicated model. The data is also included below. Thanks to
> >> >> > any
> who
> >> >> > can
> >> >> > help with this issue.
> >> >> >
> >> >> >
> >> >> > Regards,
> >> >> >
> >> >> > Charles
> >> >> >
> >> >> >
> >> >> > Row ID Group Died Event_name var
> >> >> > 1 12510 3 YES B -1.05257
> >> >> > 2 12510 3 YES S45 -1.00000
> >> >> > 3 12510 3 YES FR2 -1.14630
> >> >> > 4 12510 3 YES FR8 -1.08831
> >> >> > 5 12510 3 YES FR20 -1.03339
> >> >> > 6 21510 3 NO B -0.87290
> >> >> > 7 21510 3 NO S45 -1.22185
> >> >> > 8 21510 3 NO FR2 -1.01592
> >> >> > 9 21510 3 NO FR8 -1.06550
> >> >> > 10 21510 3 NO PR48 -0.65326
> >> >> > 11 22210 3 YES B -1.37059
> >> >> > 12 30810 3 NO B -0.92959
> >> >> > 13 30810 3 NO S45 -1.26680
> >> >> > 14 30810 3 NO FR2 -0.99396
> >> >> > 15 30810 3 NO FR8 -1.04769
> >> >> > 16 30810 3 NO FR20 -1.04769
> >> >> > 17 31510 4 YES B -1.50585
> >> >> > 18 31510 4 YES S45 -1.35262
> >> >> > 19 31510 4 YES FR2 -1.26520
> >> >> > 20 31510 4 YES FR8 -1.11407
> >> >> > 21 31510 4 YES FR20 -0.97757
> >> >> > 22 40510 2 YES B -1.04721
> >> >> > 23 40510 2 YES S45 -1.13371
> >> >> > 24 51010 3 NO B -1.21467
> >> >> > 25 51010 3 NO S45 -1.18177
> >> >> > 26 51010 3 NO FR2 -1.23582
> >> >> > 27 51010 3 NO FR8 -1.35262
> >> >> > 28 51010 3 NO FR20 -0.69897
> >> >> > 29 51010 3 NO PR48 -1.14086
> >> >> > 30 51710 4 NO B -1.13371
> >> >> > 31 51710 4 NO S45 -1.24336
> >> >> > 32 51710 4 NO FR2 -0.85201
> >> >> > 33 51710 4 NO FR8 -0.72700
> >> >> > 34 51710 4 NO FR20 -0.71220
> >> >> > 35 51710 4 NO PR48 -1.10735
> >> >> > 36 51809 1 NO B -0.70752
> >> >> > 37 51809 1 NO FR2 -0.58120
> >> >> > 38 51809 1 NO FR8 -0.71332
> >> >> > 39 51809 1 NO FR20 -0.37212
> >> >> > 40 61410 4 NO B -1.25493
> >> >> > 41 61410 4 NO S45 -1.24949
> >> >> > 42 61410 4 NO FR2 -1.18310
> >> >> > 43 61410 4 NO FR8 -0.87517
> >> >> > 44 61410 4 NO FR20 -0.34679
> >> >> > 45 61410 4 NO PR48 -0.76472
> >> >> > 46 62110 3 NO B -1.34199
> >> >> > 47 62110 3 NO S45 -1.23657
> >> >> > 48 62110 3 NO FR2 -1.20412
> >> >> > 49 62110 3 NO FR8 -0.72011
> >> >> > 50 62110 3 NO FR20 -0.72308
> >> >> > 51 62110 3 NO PR48 -0.89688
> >> >> > 52 62209 3 NO B -0.74041
> >> >> > 53 62209 3 NO S45 -0.99183
> >> >> > 54 62209 3 NO FR2 -1.06854
> >> >> > 55 62209 3 NO FR8 -0.94386
> >> >> > 56 62209 3 NO FR20 -0.38817
> >> >> > 57 71210 4 YES B -1.23657
> >> >> > 58 71210 4 YES S45 -0.94808
> >> >> > 59 71309 4 YES B -1.14630
> >> >> > 60 71309 4 YES S45 -1.10568
> >> >> > 61 71309 4 YES FR2 -0.94462
> >> >> > 62 71309 4 YES FR8 -0.89211
> >> >> > 63 71309 4 YES FR20 -0.98970
> >> >> > 64 71309 4 YES PR48 -0.83387
> >> >> > 65 81009 4 YES B -1.15989
> >> >> > 66 81009 4 YES S45 -0.93779
> >> >> > 67 81009 4 YES FR2 -0.86710
> >> >> > 68 81009 4 YES FR8 -0.69379
> >> >> > 69 81610 1 NO B -1.04769
> >> >> > 70 81610 1 NO S45 -1.02780
> >> >> > 71 81610 1 NO FR2 -0.96098
> >> >> > 72 81610 1 NO FR8 -0.91829
> >> >> > 73 81610 1 NO FR20 -0.68550
> >> >> > 74 81610 1 NO PR48 -0.66374
> >> >> > 75 81709 4 YES B -0.65521
> >> >> > 76 82410 3 NO B -1.06702
> >> >> > 77 82410 3 NO S45 -1.30103
> >> >> > 78 82410 3 NO FR2 -0.96617
> >> >> > 79 82410 3 NO FR8 -0.84497
> >> >> > 80 82410 3 NO FR20 -0.47716
> >> >> > 81 82410 3 NO PR48 -0.68131
> >> >> > 82 91310 4 NO B -1.30803
> >> >> > 83 91310 4 NO S45 -1.07212
> >> >> > 84 91310 4 NO FR2 -0.94615
> >> >> > 85 91310 4 NO FR8 -1.29414
> >> >> > 86 91310 4 NO FR20 -0.81304
> >> >> > 87 91310 4 NO PR48 -0.86519
> >> >> > 88 91409 4 NO B -1.12321
> >> >> > 89 91409 4 NO S45 -1.08831
> >> >> > 90 91409 4 NO FR2 -1.11862
> >> >> > 91 91409 4 NO FR8 -0.87976
> >> >> > 92 91409 4 NO FR20 -0.87517
> >> >> > 93 91409 4 NO PR48 -0.68867
> >> >> > 94 92109 3 NO B -1.15120
> >> >> > 95 92109 3 NO S45 -0.99525
> >> >> > 96 92109 3 NO FR2 -0.95429
> >> >> > 97 92109 3 NO FR8 -0.54914
> >> >> > 98 92109 3 NO FR20 -0.67264
> >> >> > 99 92109 3 NO PR48 -0.74256
> >> >> > 100 92110 4 YES B -1.05012
> >> >> > 101 92110 4 YES S45 -1.14086
> >> >> > 102 92110 4 YES FR2 -0.91901
> >> >> > 103 92110 4 YES FR8 -0.80521
> >> >> > 104 92110 4 YES FR20 -0.88807
> >> >> > 105 92710 4 NO B -1.20412
> >> >> > 106 92710 4 NO S45 -1.21325
> >> >> > 107 92710 4 NO FR2 -1.19518
> >> >> > 108 92710 4 NO FR8 -0.92959
> >> >> > 109 92710 4 NO FR20 -0.69315
> >> >> > 110 92710 4 NO PR48 -0.73779
> >> >> > 111 100410 1 NO B -0.92300
> >> >> > 112 100410 1 NO S45 -1.07314
> >> >> > 113 100410 1 NO FR2 -0.80162
> >> >> > 114 100410 1 NO FR8 -0.63283
> >> >> > 115 100410 1 NO FR20 -0.73119
> >> >> > 116 100410 1 NO PR48 0.10102
> >> >> > 117 101209 3 NO B -1.24718
> >> >> > 118 101209 3 NO S45 -1.00524
> >> >> > 119 101209 3 NO FR2 -1.01592
> >> >> > 120 101209 3 NO FR8 -0.85201
> >> >> > 121 101209 3 NO FR20 -0.51670
> >> >> > 122 101209 3 NO PR48 -0.72239
> >> >> > 123 101909 3 NO B -1.13549
> >> >> > 124 101909 3 NO S45 -0.92665
> >> >> > 125 101909 3 NO FR2 -0.98548
> >> >> > 126 101909 3 NO FR8 -1.03152
> >> >> > 127 101909 3 NO FR20 -0.71715
> >> >> > 128 101909 3 NO PR48 -1.07935
> >> >> > 129 102510 2 YES B -1.30103
> >> >> > 130 102510 2 YES S45 -0.83387
> >> >> > 131 110810 3 YES B -1.21968
> >> >> > 132 110810 3 YES S45 -0.68048
> >> >> > 133 110810 3 YES FR2 -1.09474
> >> >> > 134 110909 2 YES B -1.09259
> >> >> > 135 110909 2 YES S45 -0.91293
> >> >> > 136 111510 3 NO B -1.39794
> >> >> > 137 111510 3 NO S45 -1.42366
> >> >> > 138 111510 3 NO FR2 -1.22621
> >> >> > 139 111510 3 NO FR8 -0.98885
> >> >> > 140 111510 3 NO FR20 -0.66055
> >> >> > 141 111510 3 NO PR48 -0.84497
> >> >> > 142 111609 1 NO B -0.79942
> >> >> > 143 111609 1 NO S45 -0.87517
> >> >> > 144 111609 1 NO FR2 -0.74041
> >> >> > 145 111609 1 NO FR8 -0.64092
> >> >> > 146 111609 1 NO FR20 -0.53835
> >> >> > 147 111609 1 NO PR48 -0.68719
> >> >> > 148 120610 4 NO B -1.05110
> >> >> > 149 120610 4 NO S45 -1.45100
> >> >> > 150 120610 4 NO FR2 -1.21183
> >> >> > 151 120610 4 NO FR8 -1.04144
> >> >> > 152 120610 4 NO FR20 -0.61101
> >> >> > 153 120610 4 NO PR48 -0.81702
> >> >> > 154 120709 1 NO B -0.84497
> >> >> > 155 120709 1 NO S45 -0.99439
> >> >> > 156 120709 1 NO FR2 -1.13430
> >> >> > 157 120709 1 NO FR20 -0.57187
> >> >> > 158 120709 1 NO PR48 -0.80632
> >> >> > 159 121310 3 YES B -1.16115
> >> >> > 160 121310 3 YES S45 -1.18310
> >> >> >
> >> >> > 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)
> >> >> >
> >> >> >
> >> >> > # What the results should be from SAS AR1 var #Type 3 Tests
> >> >> > of Fixed Effects
> >> >> > #Effect NumDF DenDF F Value Pr > F
> >> >> > #Event_name 5 91 7.43 <.0001
> >> >> > #Died 1 24 0.08 0.7756
> >> >> > #Event_name*Died 5 91 3.09 0.0128
> >> >> >
> >> >> > fit.13=lme(var~Event_name*Died,
> >> >> > data=liver34,
> >> >> > random=~1|ID,
> >> >> > corr=corAR1())
> >> >> > anova(fit.13, type="marginal", adjustSigma=F)
> >> >> > # numDF denDF F-value p-value
> >> >> > #(Intercept) 1 91 1342.7364 <.0001
> >> >> > #Event_name 5 91 8.7143 <.0001
> >> >> > #Died 1 24 0.0527 0.8204
> >> >> > #Event_name:Died 5 91 3.0909 0.0127
> >> >> >
> >> >> > ######################################################
> >> >> > ###THIS DOES IT but not lme function (need dfs)
> >> >> > fit.16=gls(var~Event_name*Died,
> >> >> > data=liver34,
> >> >> > corr=corAR1(, ~1|ID))
> >> >> > anova(fit.16, type="marginal", adjustSigma=F) #Denom. DF: 115
> >> >> > # numDF F-value p-value
> >> >> > #(Intercept) 1 1496.9032 <.0001
> >> >> > #Event_name 5 7.4304 <.0001
> >> >> > #Died 1 0.0831 0.7736
> >> >> > #Event_name:Died 5 3.0871 0.0119
> >> >> >
> >> >> > #Give exact replicated answers 1-pf(7.4304, 5, 91)
> >> >> > 1-pf(.0831, 1, 24) 1-pf(3.0871, 5, 91)
> >> >> > #######################################################
> >> >> >
> >> >> > [[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/
> >>
> >
> > [[alternative HTML version deleted]]
> >
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
[[alternative HTML version deleted]]
This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.
More information about the R-sig-mixed-models
mailing list