Я пытаюсь сопоставить модель с повторными измерениями с авторегрессией, которая была изначально подобрана с использованием SAS. Код SAS выглядит следующим образом ...
*** Repeated Measures ANOVA Analysis ***;
*** for cs covariance structure ***;
options pageno=1;
ods output dimensions=work._dim
(where=(ddescr ? 'Covariance') rename=(value=parms
descr=ddescr))
fitStatistics = work._fit;
proc mixed data=_proj_.Set_bisc1_set info;
class SET;
model ABS_PIN_DIFF = DAY SET DAY*SET / ddfm=satterth htype=3 solution
outp=WORK._PRED;
repeated SET / sub=DAY type=cs;
run;
** Accumulate info criteria for multiple cov structures **;
data work._ic;
length label $ 128;
retain aic bic ll2 . label;
label = "Compound symmetry";
if _n_ = 1 then set work._dim;
set work._fit end=_last_;
select (scan(descr,1,' '));
when ("AIC") aic = value;
when ("BIC") bic = value;
when ("-2") ll2 = value;
otherwise delete;
end;
if _last_ then output;
drop descr value;
run;
data work._ll2; set work._ic; drop aic bic;
rename parms=parms_0 ll2 = ll2_0;
run;
proc append base=work._icall data=work._ic;
proc datasets nolist library=work;
delete _dim _fit _ic;
run; quit;
goptions reset=all device=WIN;
И вывод, который я пытаюсь сопоставить, ...
The Mixed Procedure
Model Information
Data Set _PROJ_.SET_BISC1_SET
Dependent Variable Abs_Pin_Diff
Covariance Structure Autoregressive
Subject Effect Day
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
SET 3 BISC1SET1 BISC1SET2 BISC1SET3
Dimensions
Covariance Parameters 2
Columns in X 8
Columns in Z 0
Subjects 9
Max Obs Per Subject 3
Number of Observations
Number of Observations Read 73
Number of Observations Used 27
Number of Observations Not Used 46
Iteration History
Iteration Evaluations -2 Res Log Like Criterion
0 1 218.92633483
1 2 213.69046030 0.00000316
2 1 213.69018242 0.00000000
Convergence criteria met.
Covariance Parameter Estimates
Cov Parm Subject Estimate
AR(1) Day 0.5767
Residual 183.76
The Mixed Procedure
Fit Statistics
-2 Res Log Likelihood 213.7
AIC (smaller is better) 217.7
AICC (smaller is better) 218.4
BIC (smaller is better) 218.1
Null Model Likelihood Ratio Test
DF Chi-Square Pr > ChiSq
1 5.24 0.0221
Solution for Fixed Effects
Standard
Effect SET Estimate Error DF t Value Pr > |t|
Intercept 11.6248 8.2068 13.7 1.42 0.1790
Day 0.006454 0.009346 13.7 0.69 0.5014
SET BISC1SET1 -1.7517 9.4816 18.5 -0.18 0.8554
SET BISC1SET2 -14.1656 7.5510 13.2 -1.88 0.0829
SET BISC1SET3 0 . . . .
Day*SET BISC1SET1 0.01310 0.01080 18.5 1.21 0.2403
Day*SET BISC1SET2 0.009653 0.008599 13.2 1.12 0.2816
Day*SET BISC1SET3 0 . . . .
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
Day 1 7.08 3.40 0.1072
SET 2 13.3 2.57 0.1135
Day*SET 2 13.3 0.84 0.4520
Похоже, что они используют "Анализ ANOVA с повторными измерениями", который, насколько я знаю, не существует в R. Модели, которые я пытался подогнать, это ...
trymodel3<-gls(Abs_Pin_Diff~Day+SET+Days*SET,
correlation = corCAR1(form=~Day|SET),
data = Andrea3, na.action = na.exclude, method="REML")
trymodel4<-lme(Abs_Pin_Diff~Day+SET+Day*SET, random =~Day|SET,
correlation = corCAR1(form=~Day|SET),
data = Andrea3, na.action = na.exclude, method="REML")
Вывод, который я получаю, немного отличается от того, который производит SAS. Следовательно, мне интересно, связано ли это с различиями в двух программах или если мне не удалось собрать всю модель SAS в R. Мои данные выглядят примерно так ...
dput(Andrea3)
structure(list(SET = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L), .Label = c("BISC1SET1", "BISC1SET2", "BISC1SET3"
), class = "factor"), Day = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2,
3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9), Abs_Pin_Diff = c(0,
-7.0175, 25.2675, 47.4775, 33.8325, 26.605, 37.7725, 27.1425,
26.785, 0, -8.0575, -6.06, 33.23, 2.4275, 4.7975, 23.6675, 13.5,
19.875, 0, -1.5775, 31.2025, 34.6925, 16.87, 14.125, 23.0575,
11.6675, 17.1475)), row.names = c(NA, -27L), vars = "SET", drop = TRUE, class
= c("grouped_df",
"tbl_df", "tbl", "data.frame"))