[R-sig-ME] 2nd attempt - conflict of dfs or f value in lme
Douglas Bates
bates at stat.wisc.edu
Mon May 21 19:22:25 CEST 2012
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
>
More information about the R-sig-mixed-models
mailing list