r - 使用 apply(,2,) 对巨大数据集的列执行混合效果 (lme())

Question

我是这方面的新手，所以我无法确定这是否愚蠢。
基本上，我想在一个巨大的数据集中的所有连续变量之间做成对的混合效应模型。显而易见的替代方案是简单的 spearman 相关，但我有我的理由，要解释我为什么要使用混合效应模型需要很长时间。

数据看起来像这样:

0     X1507.07  XAB1524.33  XAB1624.21  XAB1808.09...(~4000 columns)
1       12         19           12         45
2       15         35           2          25
3       22         23           65         33
4       0          55           23         67
5       12         10           90         94
6       34         22           11         2
...
90      13         8            14         45

目标是所有列的成对模型。
这是脚本有问题的部分:

for(i in 1:ncol(dat))
{
ni<-names(dat)[i]
pvalue <- apply(dat, 2, function(x)
    {
formula<-as.formula(paste(ni,"~", x," + Location",sep=""))
model<-do.call("lme", args = list(formula, random=~1|Subject, data=dat))    
summary(model)$tTable[2,5]
    })

错误:

invalid model formula in ExtractVars

对于那些感到困惑的人:我使用 as.formula 因为如果您尝试:

model<-lme(X1507.07~x+Region,random=~1|Subject, data=dat)

错误:

Error in eval(expr, envir, enclos) : object 'x' not found

(“位置”和“主题”是数据框 dat 中的因素)。我只关心一个 p 值(我知道它与混合效应有争议)。我尝试在 as.formula() 中传递 x as.matrix(x) 和 colnames(x) 但似乎没有任何效果。重点是:有谁知道这是否可能？如果我必须循环大约 10 ^ 7 次，它不值得花时间(年)，所以 apply() 是我能想到的唯一合理的选择。

Answer 1

我仍然认为这可能是愚蠢的，并且您可以通过手动计算出瞬时方法的答案来获得更快的解决方案，但是(除非我在某处犯了错误)这样做的时机并不像灾难性的那样我想可能是这样。

tl;博士如果您没有遇到任何其他缩放问题，整个过程应该需要大约 3-4 天的蛮力。 lmer实际上更慢(虽然它被设计为更快地解决大问题，但由于设置成本，单个小问题的时间实际上可能更慢)。因为lme是一个不平凡的问题，我认为循环实际上是总计算成本的一小部分。

补一些数据:

set.seed(101)
n <- 10
nobs <- 90
dat <- as.data.frame(matrix(rpois(nobs*n,20),nrow=nobs))
Subject <- rep(1:n,each=nobs/n)
Location <- runif(n)

使用 nlme 进行操作:

library(nlme)
fun.lme <- function() {
    r <- numeric(n*(n-1)/2)
    k <- 1
    for (i in 2:n) {
        for (j in 1:(i-1)) {
            m <- lme(y~x+Location,random=~1|Subject,
                     data=data.frame(x=dat[,i],y=dat[,j],
                     Location,Subject))
            tt <- summary(m)$tTable[2,5]
            r[k] <- tt
            k <- k+1
        }
    }
    r
}
t1 <- system.time(r1 <- fun.lme())
detach("package:nlme")

(在使用 nlme 之前分离 lme4 是个好主意)

fun.lmer <- function(...) {
    r <- numeric(n*(n-1)/2)
    k <- 1
    for (i in 2:n) {
        for (j in 1:(i-1)) {
            m <- lmer(y~x+Location+(1|Subject),
                     data=data.frame(x=dat[,i],y=dat[,j],
                     Location,Subject),...)
            tt <- coef(summary(m))[2,2]  
            r[k] <- tt
            k <- k+1
        }
    }
    r
}

稳定的测试时序 lme4 :

library(lme4.0) ## 'stable' version (the same as you get by installing
                ## lme4 from CRAN
t2 <- system.time(r2 <- fun.lmer())
detach("package:lme4.0")

现在正在开发(r-forge) lme4 ，具有标准和非标准优化器选择:

library(lme4)
t3 <- system.time(r3 <- fun.lmer())
t4 <- system.time(r4 <- fun.lmer(optim="bobyqa"))
detach("package:lme4")

检查时间:

tvals <- c(lme=t1["elapsed"],lme4.0=t2["elapsed"],
           lme4=t3["elapsed"],lme4_bobyqa=t4["elapsed"])

将我们花费的时间缩放到完整工作的时间:

totsecs <- (3789*3788/2)*tvals/(n*(n-1)/2)
totdays <- totsecs/(60*60*24)

round(totdays,1)
##   lme.elapsed   lme4.0.elapsed   lme4.elapsed lme4_bobyqa.elapsed 
##           3.0              4.4            4.1                 3.8 

您不妨使用 t 统计量作为 p 值进行比较；据我所知，除了作为关联强度的指标之外，p 值将完全无关紧要。