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用来重复。结果是一个具有维度（# sims）（# 响应列）（# 方差分量）的 3D 数组。

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

这在我的笔记本电脑上花了大约 45 秒，所以大约需要 15 分钟来完成 1000 次重复......

于 2012-11-24T15:02:09.147 回答