使用應用（2）超過一個巨大的數據集

我這個新手的列進行混合效應（LME（）），所以我不能確定這是否是愚蠢與否。基本上，我想在一個巨大的數據集中的所有連續變量之間做成對的混合效應模型。明顯的選擇是簡單的spearman相關性，但我有我的理由，解釋爲什麼我要使用混合效果模型需要很長時間。使用應用（2）超過一個巨大的數據集

數據看起來是這樣的：

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（）是我能想到的唯一合理的選擇。

來源

2012-05-03 Ian

'apply'不會比循環快得多。 'lmer'（來自'lme4'）*可能會更快;你可以使用任何你喜歡的自由度（包括lme會猜到的值）來獲得t統計量並將它們轉換爲p值。你的數據集有多大（行*列）？你可以請一個小**重現**（http://tinyurl.com/reproducible-000）的例子？（PS：它聽起來有點像愚蠢......你能給出一個簡短的解釋，爲什麼你更喜歡混合模式的相關性嗎？你認爲http://stats.stackexchange.com？） –

感謝您的答覆！我也嘗試了lmer，但遇到了類似的問題。該數據集特別是92行x 3789列。使用混合效應的簡單解釋：從每個受試者收集2個相似但不同的生物樣品（來自胃腸道中的兩個位置），並且我想要獲得每個樣品中發現的每個測試參數之間的關係的p值而不用將採集的位置分割數據集。我已經進行了相關性分析，並希望使用ME作爲補充，不需要將數據按樣本網站劃分 – Ian

因此，您真的想要做所有3789 * 3788/2> 700萬兩兩比較嗎？我認爲，在混合效應模型的合理時間內，這將非常困難。你能爲你的特定情況制定出矩量法（即方差分解）嗎？ –

我仍然認爲這可能是愚蠢的，那你可能是能夠獲得通過鍛鍊方法 - 矩量更快的解決方案手回答，但（除非我什麼地方搞錯了）的這個時機並不像我想的那樣災難性的。

tl; dr整個事情應該需要大約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")

（與lme4工作之前分離nlme是一個好主意）

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-鍛）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值將是完全不相關的，除非作爲聯合力量的指標。

來源

2012-05-05 01:18:56

哇。你是個不好的屁股。我甚至沒有機會重現我的問題。非常感謝！！！這真太了不起了。 – Ian

使用應用（2）超過一個巨大的數據集

回答

相關問題