переменной длины отличаются R (линейное моделирование с lme4)

Question

Мой входного файла:

 
Treat1 Treat2 Batch gene1 gene2 
High Low  1  92.73 4.00 
Low  Low  1  101.85 6.00 
High High 1  136.00 4.00 
Low  High 1  104.00 3.00 
High Low  2  308.32 10.00 
Low  Low  2  118.93 3.00 
High High 2  144.47 3.00 
Low  High 2  189.66 4.00 
High Low  3  95.12 2.00 
Low  Low  3  72.08 6.00 
High High 3  108.65 2.00 
Low  High 3  75.00 3.00 
High Low  4  111.39 5.00 
Low  Low  4  119.80 4.00 
High High 4  466.55 11.00 
Low  High 4  125.00 3.00 

Там десятки тысяч дополнительных столбцов, каждый из которых имеет заголовок и список номеров, такой же длины, как «gene1» колонка.

Мой код:

library(lme4) 
library(lmerTest) 

# Import the data. 
mydata <- read.table("input_file", header=TRUE, sep="\t") 

# Make batch into a factor 
mydata$Batch <- as.factor(mydata$Batch) 

# Check structure 
str(mydata) 

# Get file without the factors, so that names(df) gives gene names. 
genefile <- mydata[c(4:2524)] 

# Loop through all gene names and run the model once per gene and print to file. 
for (i in names(genefile)){ 
    lmer_results <- lmer(i ~ Treat1*Treat2 + (1|Batch), data=mydata) 
    lmer_summary <- summary(lmer_results) 
    write(lmer_summary,file="results_file",append=TRUE, sep="\t", quote=FALSE) 
} 

Состав:

'data.frame':  16 obs. of 2524 variables: 
$ Treat1   : Factor w/ 2 levels "High","Low": 1 2 1 2 1 2 1 2 1 2 ... 
$ Treat2   : Factor w/ 2 levels "High","Low": 2 2 1 1 2 2 1 1 2 2 ... 
$ Batch   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ... 
$ gene1   : num 92.7 101.8 136 104 308.3 ... 
$ gene2   : num 4 6 4 3 10 3 3 4 2 6 ... 

Мое сообщение об ошибке:

Ошибка в model.frame.default (данные = MYDATA, падение .unused.levels = TRUE, formula = i ~: переменной длины отличаются (найдено 'Treat1') Calls: lmer ... -> Eval -> Eval -> -> model.frame.default Выполнение остановлено

Я попытался исследовать все объекты, задействованные и не вижу различий в переменной длине, и я также удостоверился, что нет отсутствующих данных. Запуск с помощью na.exclude ничего не меняет.

Любая идея о том, что происходит?

Answer 1

@ диагноза Роланда (lmer ищет переменную с именем я, а не переменная, имя является i: обязательно Lewis Carroll reference) правильно, я считаю. Самым непосредственным образом, чтобы справиться с этим было бы с reformulate(), что-то вроде:

for (i in names(genefile)){ 
    form <- reformulate(c("Treat1*Treat2","(1|Batch)"),response=i) 
    lmer_results <- lmer(form, data=mydata) 
    lmer_summary <- summary(lmer_results) 
    write(lmer_summary,file="results_file", 
      append=TRUE, sep="\t", quote=FALSE) 
} 

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

wfun <- function(x) write(summary(x), 
     file="results_file", append=TRUE, sep="\t",quote=FALSE) 
mod0 <- lmer(geneAAA ~ Treat1*Treat2 + (1|Batch), data=mydata) 
wfun(mod0) 
for (i in names(genefile)[-1]) { 
    mod1 <- refit(mod0,mydata[[i]]) 
    wfun(mod1) 
} 

(Кстати, я не уверен, что ваша команда write() делает ничего толкового ...)