[R-sig-ME] Cross-over Data with Kenward-Roger correction

[knouri]

 Dear all:I have a SAS code but need to write it in R. Any advice will be greatly appreciated.
Here is the data ( a two-period) cross-over design comparing Trt 1 vs. Trt 2:

Sbj  Seq  Per  Trt  PEF
1      1       1     1   310
1      1       2     2   270
4      1       1     1   310
 4      1       2     2   260
6      1       1     1   370
6      1       2     2   300
7      1       1     1   410
7      1       2     2   390
10    1       1    1    250
10    1       2    2   210
11    1       1    1   380
11    1       2    2   350
14    1       1    1   330
14    1       2    2   365
2      2       1    2   370
2      2       2   1   385
3      2      1    2   310
3      2      2    1   400
5      2      1    2   380
5      2      2    1   410
9      2      1   2   290
9      2      2   1   320
12    2      1   2   260
12    2      2   1   340
13    2      1   2    90
13    2      2   1   220

Here is the SAS code:
#########################################################
proc mixed data=one method=reml;
class Sbj Per Trt;
model PEF = Per Trt /ddfm=kr;
repeated Trt / sub=Sbj type=un r;
lsmeans Trt / cl alpha=0.05;
estimate 'B vs. A' Trt -1  1 / alpha=0.1 cl;
run;
###################################################
where kr option in model statement stands for "Kenward-Roger" correction. The above SAS code produces (selected portions):

       
|  Covariance Parameter Estimates |
|  Cov Parm |  Subject |  Estimate |
|  UN(1,1) |  Sbj |  3567.58 |
|  UN(2,1) |  Sbj |  4437.00 |
|  UN(2,2) |  Sbj |  6760.65 |





| Type 3 Tests of Fixed Effects |
| Effect | Num DF | Den DF | F Value | Pr > F |
| Per | 1 | 11 | 2.59 | 0.1359 |
| Trt | 1 | 11.6 | 19.26 | 0.0010 |




     
|  Estimates |
|  Label |  Estimate |  Standard
 Error |  DF |  t Value |  Pr > |t| |  Alpha |  Lower |  Upper |
|  B vs. A |  -46.5150 |  10.5999 |  11.6 |  -4.39 |  0.0010 |  0.1 |  -65.4629 |  -27.5671 |



      
|  Least Squares Means |
|  Effect |  Trt |  Estimate |  Standard
 Error |  DF |  t Value |  Pr > |t| |  Alpha |  Lower |  Upper |
|  Trt |  1 |  341.72 |  16.5696 |  12 |  20.62 |  <.0001 |  0.05 |  305.62 |  377.82 |
|  Trt |  2 |  295.20 |  22.8073 |  12 |  12.94 |  <.0001 |  0.05 |  245.51 |  344.90 |






I used several R functions including lme, lmer. For example,  the following R code:
require(lme4)
require(pbkrtest)
require(lmerTest)
mydat$Trt = as.factor(mydat$Trt)
mydat$Per = as.factor(mydat$Per)fit.lmer <- lmer(PEF ~ Trt + Per + (1| Sbj), REML = TRUE, data=mydat)
anova(fit.lmer, ddf = "Kenward-Roger", type=3)( lsmeans(fit.lmer, test.effs="Trt") )
difflsmeans(fit.lmer, test.effs="Trt")

produces rather different results (compared to SAS), especially for df's, standard errors, etc.
Analysis of Variance Table of type 3  with  Kenward-Roger approximation for degrees of freedom
       Sum Sq       Mean Sq       NumDF       DenDF        F.value        Pr(>F)   
Trt    14036         14036               1                11              18.70         0.0012 **
Per   1632           1632                 1                11                2.17          0.1683   


Least Squares Means table:
        Trt    Per    Estimate    Standard Error   DF   t-value     Lower CI Upper CI        p-value    
Trt  1  1.0  NA    341.8           20.0               13.9    17.1           299      385                <2e-16 ***
Trt  2  2.0  NA    295.2           20.0               13.9    14.8           252      338                <2e-16 ***

Differences of LSMEANS:
                Estimate      Standard Error   DF       t-value      Lower CI      Upper CI     p-value   
Trt 1 - 2     46.6           10.8                    11.0    4.32              22.9            70.3           0.001 **


I am not quite sure if I am using the right R functions. Any advice will be greatly appreciated,
Thanks again,
Keramat

  
  
	[[alternative HTML version deleted]]