R help - Jun 2007 - lme vs. SAS proc mixed. Point estimates and SEs are the same, DFs are different

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.

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}}

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.
> 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.
>