Преобразование кода повторяющихся мер в SAS в эквивалент в R

Question

Я пытаюсь преобразовать существующий код SAS, который у меня есть для исследовательского проекта в R. Я, к сожалению, совершенно не понимаю, как подойти к этому для повторных мер ANOVA, несмотря на несколько часов рассмотрения вопросов других людей как в StackExchange, так и в Интернете в целом. Я подозреваю, что это может быть, по крайней мере, если не в целом, из-за того, что я не знаю правильных вопросов для запроса и ограниченности статистики.

Во-первых, я приведу некоторые примеры данных (с разделителями табуляции, которые я не уверен, будет сохранен на SE), а затем объясните, что я пытаюсь сделать, а затем код, который я написал по этому вопросу момент.

Образец данных:

Full data frame at: http://grandprairiefriends.org/document/data.df 

Obs SbjctID Sex Treatment Measured BirthDate DateStarted DateAssayed SubjectAge_Start_days SubjectAgeAssay.d. PreMass_mg PostMass_mg DiffMass_mg PerCentMassDiff Length_mm Width_mm PO1_abs_min PO1_r2 PO2_abs_min PO2_r2 ProteinConc_ul Protein1_net_abs Protein1_mg_ml Protein1_adjusted_mg_ml Protein2_net_abs Protein2_mg_ml Protein2_adjusted_mg_ml zPO_avg_abs_min z_Protein_avg_adjusted_mg_ml POPer_ug_Protein POPer_ug_Protein_x1000 ImgDarkness1 ImgDarkness2 ImgDarkness3 ImgDarkness4 DarknessAvg AGV_1_1 AGV_1_2 AGV_2_1 AGV_2_2 AGV_12_1 AGV_12_2 z_AGV predicted_premass resid_premass predicted_premass_calculated resid_premass_calculated predicted_postmass_calculated resid_postmass_calculated predicted_postmass resid_postmass ln_premass_mg ln_postmass_mg ln_length ln_melanization ln_po sqrt_p 
1 aF001 Female a PO_P 08/05/09 09/06/09 09/13/09 32 39 282.7 309.4 26.66 9.43 10.1 5.3 0.0175 0.996 0.0201 0.996 40 0.227 0.960 0.960 0.234 1.030 1.030 0.0188 0.995 0.00031 0.31491 33.7045 35.9165 28.8383 30.3763 32.2089 NA NA NA NA NA NA NA 5.660963 -0.016576413 4.077123 1.567263 4.077123 1.657382 5.660963 0.0735429694 8.143128 8.273329 3.336283 NA -5.733124 -0.007231569 
2 aF002 Female a PO_P 08/02/09 09/06/09 09/13/09 35 42 298.9 313.1 14.23 4.76 10.0 5.9 0..999 0.0134 0.996 40 0.213 0.840 0.840 0.219 0.860 0.860 0.0129 0.850 0.00025 0.25196 31.8700 31.8800 32.4680 32.3020 32.1300 NA NA NA NA NA NA NA 5.640012 0.059996453 4.056173 1.643836 4.056173 1.690350 5.640012 0.1065103847 8.223519 8.290480 3.321928 NA -6.276485 -0.234465254 
3 aF003 Female a PO_P 08/03/09 09/06/09 09/13/09 34 41 237.1 270.6 33.53 14.14 9.4 5.3 0.0227 0.992 0.0248 0.994 40 0.245 1.120 1.120 0.235 1.030 1.030 0.0238 1.075 0.00037 0.36822 36.0565 41.9355 41.6260 40.0180 39.9090 NA NA NA NA NA NA NA 5.509734 -0.041209334 3.925894 1.542630 3.925894 1.674895 5.509734 0.0910560222 7.889352 8.080018 3.232661 NA -5.392895 0.104336660 
82 bM001 Male b PO_P 08/02/09 08/31/09 09/07/09 29 36 468.1 371.7 -96.38 -20.59 10.7 6.8 0.0049 0.999 0.0056 1.000 40 0.228 0.350 0.350 0.222 0.330 0.330 0.0053 0.340 0.00026 0.25735 NA NA NA NA NA NA NA NA NA NA NA NA 5.782468 0.366214334 4.198628 1.950054 4.198628 1.719513 5.640012 -0.0844204671 8.870673 8.537995 3.419539 NA -7.559792 -1.556393349 
157 cM022 Male c PO_P 08/03/09 10/31/09 11/07/09 89 96 451.1 402.4 -48.71 -10.80 11.3 6.9 0.0024 0.995 0.0026 0.995 10 0.091 0.110 0.028 NA NA NA 0.0025 0.028 0.00152 1.51515 NA NA NA NA NA NA NA NA NA NA NA NA 5.897342 0.214325251 4.313502 1.798165 4.313502 1.683895 5.897342 0.1000552907 8.817303 8.652486 3.498251 NA -8.643856 -5.158429363 

Объяснение того, что я ищу, чтобы выполнить:

Этот эксперимент пытается определить, является ли конкретный режим кормления (лечение) оказало влияние на после-экспериментальной массы субъекта (ln_postmass_mg). Масса каждого человека измерялась дважды, один раз в начале (ln_premass_mg) и один раз в конце режима питания. Секс, Лечение и Измерение являются категориальными переменными.

Я создал некоторый код R, но результат не соответствует коду SAS, которого он не должен, поскольку я не верю, что он закодирован для повторных мер. Мне непонятно, нужно ли мне переставлять или иным образом манипулировать моей файловой рамкой в ​​R, чтобы выполнять дополнительные анализы, или что. Кажется, я читаю несколько разных подходов к проблемам повторных мер, и я не уверен, какие из них применимы к моей конкретной проблеме. Если кто-то может поставить меня на правильный путь, чтобы узнать, как писать дополнительные строки кода, необходимые для эквивалента R, или есть предложения, я бы очень его оценил.

SAS Код:

/* test for effect of diet regime */ 
/* repeated measures ANOVA for mass */ 
proc glm data=No_diet_lab; 
class measured sex Treatment; 
model ln_premass ln_postmass=Measured Sex Treatment Measured*Sex Measured*Treatment Sex*Treatment Measured*Sex*Treatment /nouni; 
repeated time 2; 

R Код:

options(contrasts=c("contr.sum","contr.poly")) 
model <- lm(cbind(ln_premass_mg, ln_postmass_mg) ~ Sex + Treatment + Measured + Sex:Treatment + Sex:Measured + Measured:Treatment + Sex:Treatment:Measured, data = diet_lab_data, na.action=na.omit) 

Answer 1

Это должно , надеюсь реплицировать выход SAS:

Сначала мы поставить данные i n длинной формы:

df <- subset(diet_lab_data, select = c("SubjectID", "Sex", "Treatment", "Measured", 
             "ln_premass_mg", "ln_postmass_mg")) 

dfL <- reshape(df, varying = list(5:6), idvar = "SubjectID", direction = "long", 
       v.names = "ln_mass_mg") 
dfL$time <- factor(dfL$time, levels = 1:2, labels = c("pre", "post")) 
head(dfL); tail(dfL) 

     SubjectID Sex Treatment Measured time ln_mass_mg 
aF001.1  aF001 Female   a  PO_P pre 8.143128 
aF002.1  aF002 Female   a  PO_P pre 8.223519 
aF003.1  aF003 Female   a  PO_P pre 7.889352 
aF004.1  aF004 Female   a  PO_P pre 8.521993 
aF005.1  aF005 Female   a  PO_P pre 8.335390 
aF006.1  aF006 Female   a  PO_P pre 8.259743 
     SubjectID Sex Treatment Measured time ln_mass_mg 
cM033.2  cM033 Male   c Melaniz post 8.163398 
bF037.2  bF037 Female   b Melaniz post 8.222070 
cM032.2  cM032 Male   c Melaniz post 8.422485 
cF030.2  cF030 Female   c Melaniz post 8.580447 
cM039.2  cM039 Male   c Melaniz post 8.710118 
cM036.2  cM036 Male   c Melaniz post 8.049849 

Это лучше. Теперь мы подходим к модели с использованием aov и укажем time как фактор внутри предмета.

aovMod <- aov(ln_mass_mg ~ Sex * Treatment * Measured * time + 
       Error(SubjectID/time), data = dfL) 

Все, что было сказано, я не уверен, что это соответствующий анализ, так как ваша конструкция является несбалансированной. Рассмотрим модель смешанных эффектов.