[R-sig-ME] 2nd attempt - conflict of dfs or f value in lme
Joshua Wiley
jwiley.psych at gmail.com
Tue May 22 17:42:46 CEST 2012
Hi Charles,
Unless you have a maniacal employer who demands that you replicate the
fixed effects from SAS, I do not think that adding DFs to the gls
output will really solve your problem. You are still talking about
rather different models from PROC MIXED to gls(). You may be able to
get equivalent tests of the fixed effects, but I can promise the
random effects are not identical. So you want to modify stable R
code, in order to get output from R to match a portion of output from
a different model altogether in SAS? I appreciate that many people
are only interested in the fixed effects and just consider the
nonindependence in their data a nuissance factor, so you may have no
strong preference for random effects versus gls, but we are not
talking about most people here. We are talking about R which is used
by all sorts of people.
If you need that in your work, you already showed in your original
email you can do it. You have the F-values, you just want to
calculate p-values based on the DF from SAS, and you can do so. If
you want a SAS independent solution, you can get the DF from lme() and
use those with gls().
Doug Bates has removed DF calculations from his more recent package,
lme4 which is sort of the new nlme, so the precedent is actually for
not providing any, leaving to the user to choose sensibly. Finally
because new development is happening in lme4 not nlme, I rather doubt
any of the authors are going to want to put effort into altering code
in nlme (particularly because it could potentially break other code).
If you are doing this a lot and want a niceish interface, my
suggestion would be to hack the anova.gls() code to include
denominator DF, call it anovadf.gls() or something, and just put that
code in your .Rprofile so that function is always available to you.
Cheers,
Josh
On Tue, May 22, 2012 at 8:01 AM, Charles Determan Jr <deter088 at umn.edu> wrote:
> 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
>> >
>
>
--
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