在R中的for循环中运行多个模型后存储模型估计



我想得到我在for循环中创建的模型估计,并将它们全部保存在r中的数据帧中。代码的第一部分只是模拟一个类似的数据集来挖掘

library(readxl)
library(dplyr)
library(drc)
library(purrr)
library(readr)
library(ggplot2)
library(broom)
library(broom.mixed)
library(brms)
library(lme4)
########## Simulation #########################################################
#number of assays 
nassay= 12 
#number of plates
nplate= 3
#fixed mean from above model 
mu=0.94587
### mean SD 
musd=0.04943
#standard deviation of assay from intial model
sda=0.06260
#standard deviation of residual between plates 
sd= 0.07793

set.seed(16)

(assay = rep(LETTERS[1:nassay], each = nplate))
( plate =1:(nassay*nplate) )
plate2 =rep(c(1, 2, 3), times = 12)
### simulate assay level effect
( assayeff = rnorm(nassay, 0, sda) )
#### each assay has 3 plates the assay must be repeated for each plate
( assayeff = rep(assayeff, each = nplate) )

#### every plate measurement has an effect on potency so we draw for every observation based on sd of residual 
plateeff= (rnorm(nassay*nplate, 0, sd))
###### simulate different means 
(musims= rnorm(nassay*nplate, mu, musd))
( dat = data.frame(assay, assayeff, plate, plateeff,musims) )
sim_dat <- dat
#### now combine all estimates to get rel potency 
( dat$relpot = with(dat, mu + assayeff + plateeff ) )
sim_dat <- dat
fit1 = lmer(relpot ~ 1 + (1|assay), data = dat)
fit1 

这是模拟数据集的代码,然后我只需要BRMS来获得后验估计并保存到数据帧

fit<-brms::brm(relpot ~ 1 + (1 | Filename), data = dat,iter=100,warmup=10,seed=355545)
post_dat <- posterior_samples(fit,fixed=TRUE)
summary(fit)
plot(fit)
post_fit_use <- post_dat %>% dplyr::select(b_Intercept, sd_Filename__Intercept, sigma)
post_fit_use <- post_fit_use %>% mutate(assay_var=(sd_Filename__Intercept)^2) %>% mutate(platevar=(sigma)^2)

现在我想使用这些后验估计中的每一个来创建一个数据集并运行一个模型


for (i in 1:nrow(post_fit_use)) { #fixed mean from above model 
mu=post_fit_use$b_Intercept[i]

#standard deviation of assay from intial model
sda=post_fit_use$sd_Filename__Intercept[i]

#standard deviation of residual between plates 
sd= post_fit_use$sigma[i]

(assay = rep(LETTERS[1:nassay], each = nplate))

( plate =1:(nassay*nplate) )

plate2 =rep(c(1, 2, 3), times = 12)

### simulate assay level effect

( assayeff = rnorm(nassay, 0, sda) )

#### each assay has 3 plates the assay must be repeated for each plate

( assayeff = rep(assayeff, each = nplate) )


#### every plate measurement has an effect on potency so we draw for every observation based on sd of residual 

plateeff= (rnorm(nassay*nplate, 0, sd))

###### simulate different means 


( dat = data.frame(assay, assayeff, plate, plateeff) )
sim_dat <- dat

#### now combine all estimates to get rel potency 
( dat$relpot = with(dat, mu + assayeff + plateeff ) )

sim_dat <- dat

fit = lmer(relpot ~ 1 + (1|assay), data = dat)
rand <-tidy(fit, effects = "ran_pars", scales = "vcov")
fixed <- tidy(fit, effects = "fixed")
}

我的问题是我想把每个模型估计保存到一个数据帧中。但当我运行循环时,我只得到最后一次迭代的结果。我不确定如何保存每个

上面的代码显示了我厌倦了什么,最后一个模型保存了一个——不是全部。请注意,当您运行rand <-tidy(fit, effects = "ran_pars", scales = "vcov") fixed <- tidy(fit, effects = "fixed")时,您会得到一个包含1行5个变量的固定数据帧和一个包含2行5个变数的rand数据帧。这是一个型号的

我赞同保罗的建议。您可以将tidy的结果存储在列表中,然后在数据帧中转换列表。但是您必须使用双括号来索引列表(即rands[[i]] <- tidy(fit)(。试试类似的东西:
library(broom)
rands <- list()

for(i in 2:ncol(mtcars)){   
mod <- lm(mtcars[,1] ~ mtcars[,i])   
rands[[i]] <- tidy(mod) 
}
df <- do.call(rbind, rands) 
df

最新更新