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

（「場所」と「件名」はデータフレームデータの要素です）。私は1つのp値のみを気にします（混合効果で物議を醸していることを知っています）。as.formula（）でx as.matrix（x）とcolnames（x）を渡そうとしましたが、実際には何も機能していないようです。ポイントは：これが可能かどうか誰かが知っていますか？私がそれを〜10 ^ 7回ループする必要がある場合、それは時間（年）の価値がないので、apply（）は私が考えることができる唯一の合理的な代替手段です。

4

1 に答える 1

4

これはおそらく愚かなことだと思います。モーメント法の答えを手作業で計算することで、より迅速な解決策を得ることができるかもしれませんが、（どこかで間違えない限り）これのタイミングはそれほど悲惨ではありません。そうかもしれないと思いました。

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")

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

比較のためのp値としてt統計量を使用することもできます。私が知る限り、p値は、関連性の強さの指標を除いて、まったく無関係です。

于 2012-05-05T01:18:56.233 に答える