[R-sig-ME] 2nd attempt - conflict of dfs or f value in lme
Joshua Wiley
jwiley.psych at gmail.com
Mon May 21 17:02:54 CEST 2012
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/
More information about the R-sig-mixed-models
mailing list