Storing model estimates after running several models in a for loop in R-CodePudding

I want to get the model estimates that I am creating in my for loop and save them all in a data frame in r. First part of the code just to simulate a similar dataset to mine

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

This the code to simulate a dataset then I just BRMS to get posterior estimates and save to a dataframe

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)

Now I want to use each of these posterior estimates to create a dataset and run a model


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

my issue is I want to save each model estimate into a dataframe. But when I run my loop I just get the results of the last iteration. I am unsure on how to save each one

the above code shows what I tired and one the last model gets saved not all. Note when you run rand <-tidy(fit, effects = "ran_pars", scales = "vcov") fixed <- tidy(fit, effects = "fixed") you get data frames with 1 row 5 variables for fixed and a dataframe with 2 rows 5 variables for rand. This is for one model

CodePudding user response：

I second Paul's suggestion. You can store the results from tidy in a list, then transform the list in a data frame. But you will have to use double bracket to index the list (i.e. rands[[i]] <- tidy(fit)). Try something like:

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