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