Istnieje kilka pytań i postów na temat mieszanych modeli dla bardziej złożonych projektów eksperymentalnych, więc pomyślałem, że ten prostszy model pomoże innym początkującym w tym procesie, jak również I.

Więc moje pytanie brzmi: chciałbym sformułować powtarzane miary ancova w R z procedury mieszanej sas proc:

proc mixed data=df1;
FitStatistics=akaike
class GROUP person day;
model Y = GROUP X1 / solution alpha=.1 cl;
repeated / type=cs subject=person group=GROUP;
lsmeans GROUP;
run;

Oto wyjście SAS wykorzystujące dane utworzone w R (poniżej):

.           Effect       panel    Estimate       Error      DF    t Value    Pr > |t|     Alpha       Lower       Upper
            Intercept              -9.8693      251.04       7      -0.04      0.9697       0.1     -485.49      465.75
            panel        1         -247.17      112.86       7      -2.19      0.0647       0.1     -460.99    -33.3510
            panel        2               0           .       .        .         .             .           .           .
            X1                     20.4125     10.0228       7       2.04      0.0811       0.1      1.4235     39.4016

Poniżej przedstawiam, jak sformułowałem model w R przy użyciu pakietu „nlme”, ale nie uzyskuję podobnych szacunków współczynnika:

## create reproducible example fake panel data set:
set.seed(94); subject.id = abs(round(rnorm(10)*10000,0))

set.seed(99); sds = rnorm(10,15,5);means = 1:10*runif(10,7,13);trends = runif(10,0.5,2.5)

this = NULL; set.seed(98)
for(i in 1:10) { this = c(this,rnorm(6, mean = means[i], sd = sds[i])*trends[i]*1:6)}
set.seed(97)
that = sort(rep(rnorm(10,mean = 20, sd = 3),6))

df1 = data.frame(day = rep(1:6,10), GROUP = c(rep('TEST',30),rep('CONTROL',30)),
                 Y = this,
                 X1 = that,
                 person = sort(rep(subject.id,6)))

## use package nlme
require(nlme)

## run repeated measures mixed model using compound symmetry covariance structure:
summary(lme(Y ~ GROUP + X1, random = ~ +1 | person,
            correlation=corCompSymm(form=~day|person), na.action = na.exclude,
            data = df1,method='REML'))

Teraz wyjście R, które teraz sobie uświadamiam, jest podobne do wyjścia zlm():

                Value Std.Error DF    t-value p-value
(Intercept) -626.1622  527.9890 50 -1.1859379  0.2413
GROUPTEST   -101.3647  156.2940  7 -0.6485518  0.5373
X1            47.0919   22.6698  7  2.0772934  0.0764

Uważam, że jestem bliski specyfikacji, ale nie jestem pewien, jakiego fragmentu brakuje, aby wyniki były zgodne (w granicach rozsądku ..). Każda pomoc byłaby mile widziana!

AKTUALIZACJA: Używając kodu w odpowiedzi poniżej, wyjście R staje się:

> summary(model2)

Przewiń do dołu, aby uzyskać szacunki parametrów - patrz! identyczny z SAS.

Linear mixed-effects model fit by REML
 Data: df1 
      AIC      BIC   logLik
  776.942 793.2864 -380.471

Random effects:
 Formula: ~GROUP - 1 | person
 Structure: Diagonal
        GROUPCONTROL GROUPTEST Residual
StdDev:      184.692  14.56864 93.28885

Correlation Structure: Compound symmetry
 Formula: ~day | person 
 Parameter estimate(s):
         Rho 
-0.009929987 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | GROUP 
 Parameter estimates:
    TEST  CONTROL 
1.000000 3.068837

Fixed effects: Y ~ GROUP + X1 

                Value Std.Error DF    t-value p-value
(Intercept)   -9.8706 251.04678 50 -0.0393178  0.9688
GROUPTEST   -247.1712 112.85945  7 -2.1900795  0.0647
X1            20.4126  10.02292  7  2.0365914  0.0811