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

[Lenth, Russell V]

For some reason, there were scads of "?"s in your posting; I removed them in the quote below so it's easier to read.

You did not fit the same model in R as you did in SAS. The SAS one has unstructured covariance (3 covariance parameters plus residual error), whereas the R one has just one, and corresponds to the SAS model with the REPEATED statement replaced by "RANDOM SBJ;" 

I don't think you can fit the unstructured model with lmer (but I'm sure somebody will correct me if I'm wrong). It is possible to do that with lme in the nlme package, but you can't get the K-R d.f. with that.

Russ

Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science   
The University of Iowa  -  Iowa City, IA 52242  USA   
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017



-----Original Message-----

Message: 1
Date: Wed, 10 Jun 2015 04:18:00 +0000 (UTC)
From: knouri <nouri4 at yahoo.com>
To: "r-sig-mixed-models at r-project.org"
	<r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Cross-over Data with Kenward-Roger correction
Message-ID:
	<180812510.9495.1433909880558.JavaMail.yahoo at mail.yahoo.com>
Content-Type: text/plain; charset="UTF-8"

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