0

私はlmerと混合モデリングを行っており、データのランダム化の効果を調査したいと思います。次のコードは速度の点で改善できますか？現状では、現在の仕様では実際のデータを実行するのに数日かかるでしょう...私が見つけたものから、sapplyが進むべき道です。私は間違っていると思っています。

library(lme4)

##Generate real data
        real.data=data.frame(cat1=factor(rep(c("A","B","C","D"),500)),cat2=factor(rep(c("E","F","G","H"),500)), matrix(runif(12000),ncol=6))

##Apply lmer model for each data columm and extract variance estimates with VarCorr for 1000 randomizations of the real.data frame.

random=sapply(1:1000,function(z){print(z)

##Generate a randomized data set by sampling first two factor columns

    sample=data.frame(cat1=factor(sample(real.data$cat1)),cat2=factor(sample(real.data$cat2)), real.data[,3:8])

sapply(3:dim(sample)[2],function(y){print(y)

##Apply REML to each column of data, with 'cat1' and 'cat2' as random effects, including cat1:cat2 interaction

model=lmer(sample[,y]~(1|cat1)+(1|cat2)+(1|cat1:cat2), data=sample)

##Extract the estimates of the random effect terms
                c(as.numeric(VarCorr(model)$cat1),as.numeric(VarCorr(model)$cat2),as.numeric(VarCorr(model)$'cat1:cat2'))
})
})
4

1 に答える 1

1

sapplyなどはあまり時間の節約にはなりませんが、(場合によっては) よりクリーンになります。一方、因子列をスクランブルしたらrefit、時間の節約で異なる応答データ (つまり、列) のモデルを再調整するために使用できます。

plyr以下のいくつかの目的でパッケージを使用しました。

sample異なる応答変数を当てはめた順序をごちゃまぜにしていた理由がよくわからないので、その部分は省略しました...

予選:

library(lme4)
library(plyr)
set.seed(101)
##Generate real data
real.data=data.frame(cat1=factor(rep(c("A","B","C","D"),500)),
cat2 <- factor(rep(c("E","F","G","H"),500)), matrix(runif(12000),ncol=6))

因子列をランダム化し、列 3 から 8 のそれぞれに基づいてモデルを応答に適合させる関数は次のとおりです ...

sfun <- function() {
###Generate a randomized data set by sampling first two factor columns
    sampledat <- transform(real.data,
                           cat1=factor(sample(cat1)),
                           cat2=factor(sample(cat2)))
    ## fit first column
    m1 <- lmer(X1 ~ (1|cat1)+(1|cat2)+(1|cat1:cat2), data=sampledat)
    ## refit using every other column
    m_rest <- apply(real.data[,-(1:3)],2,
                     refit,object=m1)
    ## note this is 'laply' (from plyr), not 'lapply'
    laply(c(list(m1),m_rest),function(m) unlist(VarCorr(m)))
}

今すぐraply繰り返すために使用します。結果は、次元 (シム数) (応答列数) (分散成分数) の 3D 配列です。

nsim <- 50
sres <- raply(nsim,sfun(),.progress="text")

これには私のラップトップで約 45 秒かかったので、1000 回の繰り返しを行うには約 15 分かかりました...

于 2012-11-24T15:02:09.147 に答える