Преобразование формулы смешанной модели из SAS в R

Question

Я хочу разместить смешанную модель, используя пакет nlme в R, что эквивалентно следующим кодам SAS:

proc mixed data = one;
class var1 var2  year loc rep;
model yld = var1 * var2;
random loc year(loc) rep*year(loc);

РЕДАКТИРОВАТЬ: объяснение, что такое эксперимент о

одна и та же комбинация var1 и var2 была протестирована в повторностях (повторения пронумерованы 1: 3). Копия (повтор) считается случайной. Этот набор экспериментов повторяется по местоположениям (loc) и годам (year). Хотя реплики пронумерованы 1: 3 для каждого местоположения и года для удобства, поскольку они не имеют никакого имени, репликация 1 внутри местоположения и года не имеет корреляции репликации 1 в другом местоположении и другом году

Я пробовал следующие коды:

 require(nlme) 
    fm1 <- lme(yld ~ var1*var2, data = one, random = loc + year / loc + rep * year / loc)  

Верны ли мои коды?

РЕДАКТИРОВАТЬ: данные и модель на основе предложений Вы можете скачать файл данных примера по следующей ссылке: https://sites.google.com/site/johndatastuff/mydata1.csv

data$var1 <- as.factor(data$var1)
data$var2 <- as.factor(data$var2)
data$year <- as.factor(data$year)
data$loc <- as.factor(data$loc)
data$rep <- as.factor(data$rep)

following suggestions from the comments below:
fm1 <- lme(yld ~ var1*var2, data = data, random = ~ loc + year / loc + rep * year / loc)

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

ОЖИДАЕТСЯ НА ОСНОВЕ ВЫХОДА SAS

Type 3 tests of fixed effects 
var1*var2         14         238       F value 16.12 Pr >F = < 0.0001

Covariance parameters:
loc = 0, year(loc) = 922161, year*rep(loc) = 2077492, residual = 1109238 

Я попробовал следующую модель, но все еще получаю ошибки:

Edits: Just for information I tried the following model
require(lme4)  
 fm1 <- lmer(yld ~ var1*var2 + (1|loc) +  (1|year / loc) + (1|rep : (year / loc)),  
            data = data)  
Error in rep:`:` : NA/NaN argument 
In addition: Warning message: 
In rep:`:` : numerical expression has 270 elements: only the first used

Answer 1

Спасибо за более подробную информацию. Я сохранил данные в d, чтобы избежать путаницы с функцией и параметром data; Команды работают в любом случае, но избегание data обычно считается хорошей практикой.

Обратите внимание, что взаимодействие трудно приспособить из-за отсутствия баланса между var и var2; для справки вот кросс-таблицы:

> xtabs(~var1 + var2, data=d)
    var2
var1  1  2  3  4  5
   1 18 18 18 18 18
   2  0 18 18 18 18
   3  0  0 18 18 18
   4  0  0  0 18 18
   5  0  0  0  0 18

Обычно, чтобы просто соответствовать взаимодействию (а не основным эффектам), вы бы использовали : вместо *, но здесь лучше всего работает для создания одного фактора, например:

d$var12 <- factor(paste(d$var1, d$var2, sep=""))

Затем с nlme, попробуйте

fm1 <- lme(yld ~ var12, random = ~ 1 | loc/year/rep, data = d)
anova(fm1)

и с lme4, попробуйте

fm1 <- lmer(yld ~ var12 + (1 | loc/year/rep), data=d)
anova(fm1)

Также обратите внимание, что, поскольку nlme и lme4 имеют перекрытия в именах своих функций, вам нужно загружать только одну за раз в R-сеанс; для переключения нужно закрыть R и перезапустить. (Существуют и другие способы, но это проще всего объяснить.)