Эквивалентный скрипт R для повторных мер SAS-скрипт, вложенный дизайн

Question

У меня возникли проблемы с преобразованием скрипта SAS в соответствующий R-скрипт.

Модель представляет собой анализ повторных измерений отклика (resp), основанный на обработке (trt) с участком (plot) вложен в лечении.

SAS код:

data data_set; 
input trt $ plot time resp; 
datalines; 
Burn 1 1 27 
Burn 1 9 25 
Burn 1 12 18 
Burn 1 15 21 
Burn 2 1 5 
Burn 2 9 15 
Burn 2 12 10 
Burn 2 15 12 
... 
Unburn 1 1 57 
Unburn 1 9 46 
Unburn 1 12 49 
Unburn 1 15 51 
Unburn 2 1 43 
Unburn 2 9 59 
Unburn 2 12 59 
Unburn 2 15 60 

proc mixed data = data_set; 
class trt plot time; 
model resp = trt time trt*time/ddfm = kr; 
repeated time/subject = trt(plot) type = vc rcorr; 
run; 

код R попытка (загрузка набор данных из файла CSV):

library(nlme) 

data.set <- read.csv("data_set.csv") 
data.set$plot <- factor(data.set$plot) 
data.set$time <- factor(data.set$time) 

model1 <- lme(resp ~ trt + time + trt:time, data = data.set, random = ~1 | plot) 

Это работает, но не нужная модель. Другие попытки, которые я пробовал, как правило, приводит к ошибке:

Error in getGroups.data.frame(dataMix, groups) : 
    invalid formula for groups 

В основном я уезжаю в сорняках здесь ...

Вопрос 1: как указать ту же модель в R как то, что является уже указано в SAS?

Вопрос 2: Я хочу иметь возможность изменять матрицу ковариации для тиражирования другой работы, выполненной в SAS. Я считаю, что знаю, как это сделать с параметром корреляции для функции lme. Но, пожалуйста, поправьте меня, если я ошибаюсь.

Заранее спасибо.

Answer 1

Спецификации модели в R бы быть логически:

model1 <- lme(resp ~ trt + time + trt:time, data = data.set, random = ~1 | trt:plot) 

Это при условии, что участок вложен в лечении в соответствии с кодированием, или в качестве альтернативы, существует взаимодействие между участком и лечением. Однако если он указан как таковой, то она выдает предупреждение упомянутое:

Error in getGroups.data.frame(dataMix, groups) : invalid formula for groups

проблема встречается имеет дело с уровнями введенных (я думаю) путем использования такого взаимодействия. Независимо от точного вопроса, то проблема может быть решена путем создания комбинированного лечения переменного участок предсказатель:

data.set$trtplot <- with(data.set, factor(paste(trt, plot, sep = "."))) 

А затем выполняет анализ следующим образом:

model1 <- lme(resp ~ trt + time + trt:time, data = data.set, random = ~ 1 | trtplot) 

Для полноты это может так же легко быть следующее, где каждая переменная предиктор добавляется плюс взаимодействия:

model1 <- lme(resp ~ trt * time, data = data.set, random = ~ 1 | trtplot) 

Это соответствует затем результаты, достигнутые в SAS, когда Compound Symmetry (CS) определена ковариационная структура (хотя критерий АИК отличается от другого - не совсем понятно). Так что немного отличается от кода SAS выше, где указана структура ковариации Variance Components (VC), но это всего лишь вопрос изменения типа структуры в коде SAS.

Что касается сравнения различных структур ковариации, это, по-видимому, является более сложной задачей.Ковариационная структура, которые я хотел бы исследовать, является:

• Соединения Symmetry (CS) - сделано
• Компоненты дисперсии (VC)
• неструктурированной (ООН)
• Пространственной мощности (SP)

Любые мысли были бы очень желанными!