Re: [R] lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Tue, 05 Jun 2007 11:14:17 +0200

Peter Dalgaard wrote:
> John Sorkin wrote:
>
>> R 2.3
>> Windows XP
>>
>> I am trying to understand lme. My aim is to run a random effects regression in which the intercept and jweek are random effects. I am comparing output from SAS PROC MIXED with output from R. The point estimates and the SEs are the same, however the DFs and the p values are different. I am clearly doing something wrong in my R code. I would appreciate any suggestions of how I can change the R code to get the same DFs as are provided by SAS.
>>
>>
> This has been hashed over a number of times before. In short:
>
> 1) You're not necessarily doing anything wrong
> 2) SAS PROC MIXED is not necessarily doing it right
> 3) lme() is _definitely_ not doing it right in some cases
> 4) both work reasonably in large sample cases (but beware that this is
> not equivalent to having many observation points)
>
> SAS has an implementation of the method by Kenward and Rogers, which
> could be the most reliable general DF-calculation method around (I don't
> trust their Satterthwaite option, though). Getting this or equivalent
> into lme() has been on the wish list for a while, but it is not a
> trivial thing to do.
>

Forgot to say: All DF-based corrections are wrong if you have non-normally distributed data (they depend on the 3rd and 4th moment of the error distribution(s)), although they can be useful as warning signs even in those cases. I also forgot to point to the simulate.lme() function which can simulate the LR statistics directly.
>
>> SAS code:
>> proc mixed data=lipids2;
>> model ldl=jweek/solution;
>> random int jweek/type=un subject=patient;
>> where lastvisit ge 4;
>> run;
>>
>> SAS output:
>> Solution for Fixed Effects
>>
>> Standard
>> Effect Estimate Error DF t Value Pr > |t|
>>
>> Intercept 113.48 7.4539 25 15.22 <.0001
>> jweek -1.7164 0.5153 24 -3.33 0.0028
>>
>> Type 3 Tests of Fixed Effects
>>
>> Num Den
>> Effect DF DF F Value Pr > F
>> jweek 1 24 11.09 0.0028
>>
>>
>> R code:
>> LesNew3 <- groupedData(LDL~jweek | Patient, data=as.data.frame(LesData3), FUN=mean)
>> fit3 <- lme(LDL~jweek, data=LesNew3[LesNew3[,"lastvisit"]>=4,], random=~1+jweek)
>> summary(fit3)
>>
>> R output:
>> Random effects:
>> Formula: ~1 + jweek | Patient
>> Structure: General positive-definite, Log-Cholesky parametrization
>>
>>
>> Fixed effects: LDL ~ jweek
>> Value Std.Error DF t-value p-value
>> (Intercept) 113.47957 7.453921 65 15.224144 0.0000
>> jweek -1.71643 0.515361 65 -3.330535 0.0014
>>
>> John Sorkin M.D., Ph.D.
>> Chief, Biostatistics and Informatics
>> University of Maryland School of Medicine Division of Gerontology
>> Baltimore VA Medical Center
>> 10 North Greene Street
>> GRECC (BT/18/GR)
>> Baltimore, MD 21201-1524
>> (Phone) 410-605-7119
>> (Fax) 410-605-7913 (Please call phone number above prior to faxing)
>>
>> Confidentiality Statement:
>> This email message, including any attachments, is for the so...{{dropped}}
>>
>> ______________________________________________
>> R-help_at_stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
> ______________________________________________
> R-help_at_stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue 05 Jun 2007 - 09:27:01 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 05 Jun 2007 - 09:31:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.