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

还有数万个附加列，每个列都有一个 header 和一个数字列表，长度与 "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 ...

我的错误信息：

Error in model.frame.default(data = mydata, drop.unused.levels = TRUE, formula = i ~ : variable lengths differ (found for 'Treat1') Calls: lmer ... -> eval -> eval -> -> model.frame.default Execution halted

我已尝试检查所有涉及的 objects，但看不出任何可变长度差异，而且我还确保没有丢失数据。运行它与 na.exclude 没有任何改变。

知道发生了什么吗？

Answer 1

@Roland的诊断（lmer是在找一个叫i的变量，不是name是i: obligatory 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() 命令是否合理...）

R 中的可变长度不同（使用 lme4 进行线性建模）

Variable lengths differ in R (linear modelling with lme4)

r

lme4