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