As always, complete code for this article is available on my GitHub page<\/span><\/a><\/strong>.<\/p>\n Side Note: <\/em><\/strong>This is a very code heavy article as the goal is to show not only simulations of data and models but how to extract the model parameters from the list of simulations. As such, this a rather long article with minimal interpretation of the models.<\/p>\n First a linear model<\/strong><\/span><\/p>\n Our data will look at training load for 10 players in two different positions, Forwards (F) and Mids (M). The Forwards will be simulated to have less training load than the Mids and we will add some random error around the difference in these group means. We will simulate this as a linaer model where b0 represents the model intercept and is the mean value for Forwards and b1 represents the coefficient for when the group is Mids.<\/p>\n Calculate some summary statistics<\/p>\n Build a linear model<\/p>\n Now, we wrap all the steps in a function so that we can use replicate()<\/strong> to create thousands of simulations or to change parameters of our simulation. The end result of the function is going to be the linear regression model.<\/p>\n Try the function out with the default parameters<\/p>\n Use replicate()<\/strong> to create many simulations.<\/p>\n Doing this simulation once doesn’t help us. We want to be able to do this thousands of times. All of the articles in this series have used for()<\/strong> loops up until this point. But, if you recall the first article in the series where I laid out several helpful functions for coding simulations, I showed an example of the replicate()<\/strong> function, which will take a function and run it’s result for as many times as you specify. I found this function while I was working through the simulation chapter of Gelman et al’s book, Regression and Other Stories<\/a><\/span><\/strong>. I think in cases like a mixed model simulation, where you can have many layers and complexities to the data, writing a simple function and then replicating it thousands of times is much easier to debug and much cleaner for others to read than having a bunch of nested for()<\/strong> loops.<\/p>\n Technical Note:<\/strong><\/em> We specify the argument simplify = FALSE<\/strong> so that the results are returned in a list format. This makes more sense since the results are the regression summary results and not a data frame.<\/p>\n Coefficient Results<\/p>\n Compare the simulated results to the results of the original model fit.<\/p>\n Model fit parameters<\/p>\n Mixed Model 1<\/strong><\/span><\/p>\n Now that we have the general frame work for building a simulation function and using replicate()<\/strong> we will to build a mixed model simulation.<\/p>\n Above we had a team with two positions groups and individual players nested within those position groups. In this mixed model, we will add a second team so that we can explore hierarchical data.<\/p>\n We will simulate data from 3 teams, each with 2 positions (Forward & Mid). This is a pretty simple mixed model. We will build a more complex one after we get a handle on the code below.<\/p>\n Construct the mixed model and evaluate the outputs<\/p>\n Write a mixed model function<\/p>\n Try out the function<\/p>\n Now use replicate()<\/strong> and create 1000 simulations of the model and look at the first model in the list<\/p>\n Store the coefficient results of the 1000 simulations in a data frame, create plots of the model coefficients, and compare the results of the simulation to the original mixed model.<\/p>\n Extract the random effects for the intercept and the residual for each of the simulated models<\/p>\n Mixed Model 2<\/strong><\/span><\/p>\n Above was a pretty simple model, just to get our feet wet. Let’s create a more complicated model. Usually in sport and exercise science we have repeated measures of individuals. Often, researchers will set the individual players as random effects with the fixed effects being the component that the researcher is attempting to make an inference about.<\/p>\n In this example, we will set up a team of 12 players with three positions (4 players per position): Forward, Mid, Defender. The aim is to explore the training load differences between position groups while accounting for repeated observations of individuals (in this case, each player will have 20 training sessions). Similar to our first regression model, we will build a data frame of everything we need and then calculate the outcome variable (training load) with a regression model using parameters that we specify. To make this work, we will need to specific an intercept and slope for the position group and an intercept and slope for the individual players as well as a model sigma value. Once we’ve done that, we will fit a mixed model, write a function, and then create 1000 simulations.<\/p>\n Construct the mixed model and evaluate the output<\/p>\n Now use `replicate()` and create 1000 simulations of the model and look at the first model in the list.<\/p>\n Store the coefficient results from the simulations, summarize them, and compare them to the original mixed model.<\/p>\n\r\nlibrary(tidyverse)\r\nlibrary(broom)\r\nlibrary(lme4)\r\n\r\ntheme_set(theme_light())\r\n## Model Parameters\r\nn_positions <- 2\r\nn_players <- 10\r\nb0 <- 80\r\nb1 <- 20\r\nsigma <- 10\r\n\r\n## Simulate\r\nset.seed(300)\r\npos <- rep(c("M", "F"), each = n_players)\r\nplayer_id <- rep(1:10, times = 2)\r\nmodel_error <- rnorm(n = n_positions*n_players, mean = 0, sd = sigma)\r\ntraining_load <- b0 + b1 * (pos == "M") + model_error\r\n\r\nd <- data.frame(pos, player_id, training_load)\r\nd<\/pre>\n\r\nd %>%\r\n ggplot(aes(x = pos, y = training_load)) +\r\n geom_boxplot()\r\n\r\nd %>%\r\n group_by(pos) %>%\r\n summarize(N = n(),\r\n avg = mean(training_load),\r\n sd = sd(training_load)) %>%\r\n mutate(se = sd \/ sqrt(N)) %>%\r\n knitr::kable()\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.40.32\u202fPM-975x1024.png\")
<\/a> ![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.40.37\u202fPM.png\")
<\/a><\/p>\n\r\nfit <- lm(training_load ~ pos, data = d)\r\nsummary(fit)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.42.04\u202fPM.png\")
<\/a><\/p>\n\r\nsim_func <- function(n_positions = 2, n_players = 10, b0 = 80, b1 = 20, sigma = 10){\r\n\r\n ## simulate data\r\n pos <- rep(c("M", "F"), each = n_players)\r\n player_id <- rep(1:10, times = 2)\r\n model_error <- rnorm(n = n_positions*n_players, mean = 0, sd = sigma)\r\n training_load <- b0 + b1 * (pos == "M") + model_error\r\n\r\n ## store in data frame\r\n d <- data.frame(pos, player_id, training_load)\r\n \r\n ## construct linear model\r\n fit <- lm(training_load ~ pos, data = d)\r\n summary(fit)\r\n}\r\n<\/pre>\n\r\nsim_func()\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.44.38\u202fPM.png\")
<\/a><\/p>\n\r\nteam_training <- replicate(n = 1000,\r\n sim_func(),\r\n simplify = FALSE)\r\n<\/pre>\n
\r\nteam_training %>%\r\n map_df(tidy) %>%\r\n select(term, estimate) %>%\r\n ggplot(aes(x = estimate, fill = term)) +\r\n geom_density() +\r\n facet_wrap(~term, scales = "free_x")\r\n\r\nteam_training %>%\r\n map_df(tidy) %>%\r\n select(term, estimate) %>%\r\n group_by(term) %>%\r\n summarize(avg = mean(estimate),\r\n SD = sd(estimate))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.48.50\u202fPM.png\")
<\/a> ![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.49.02\u202fPM-975x1024.png\")
<\/a><\/p>\n\r\ntidy(fit)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.49.48\u202fPM.png\")
<\/a><\/p>\n\r\nteam_training %>%\r\n map_df(glance) %>%\r\n select(adj.r.squared, sigma) %>%\r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.51.54\u202fPM.png\")
<\/a>
\nCompare these results to the original fit<\/p>\n\r\nfit %>% \r\n glance()\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-7.52.36\u202fPM-1024x139.png\")
<\/a><\/p>\n\r\n## Model Parameters\r\nn_teams <- 3\r\nn_positions <- 2\r\nn_players <- 10\r\n\r\nteam1_fwd_avg <- 130\r\nteam1_fwd_sd <- 15\r\nteam1_mid_avg <- 100\r\nteam1_mid_sd <- 5\r\n\r\nteam2_fwd_avg <- 150\r\nteam2_fwd_sd <- 20\r\nteam2_mid_avg <- 90\r\nteam2_mid_sd <- 10\r\n\r\nteam3_fwd_avg <- 180\r\nteam3_fwd_sd <- 15\r\nteam3_mid_avg <- 150\r\nteam3_mid_sd <- 15\r\n\r\n\r\n## Simulated data frame\r\nteam <- rep(c("Team1","Team2", "Team3"), each = n_players * n_positions)\r\npos <- c(rep(c("M", "F"), each = n_players), rep(c("M", "F"), each = n_players), rep(c("M", "F"), each = n_players))\r\nplayer_id <- as.factor(round(seq(from = 100, to = 300, length = length(team)), 0))\r\n\r\nd <- data.frame(team, pos, player_id) d %>%\r\n head()\r\n\r\n# simulate training loads\r\nset.seed(555)\r\ntraining_load <- c(rnorm(n = n_players, mean = team1_mid_avg, sd = team1_mid_sd), rnorm(n = n_players, mean = team1_fwd_avg, sd = team1_fwd_sd), rnorm(n = n_players, mean = team2_mid_avg, sd = team2_mid_sd), rnorm(n = n_players, mean = team2_fwd_avg, sd = team2_fwd_sd), rnorm(n = n_players, mean = team3_mid_avg, sd = team3_mid_sd), rnorm(n = n_players, mean = team3_fwd_avg, sd = team3_fwd_sd)) d ,- d %>%\r\n bind_cols(training_load = training_load) \r\n\r\nd %>%\r\n head()\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.04.17\u202fPM.png\")
<\/a>
\nCalculate summary statistics<\/p>\n\r\n## Average training load by team\r\nd %>%\r\n group_by(team) %>%\r\n summarize(avg = mean(training_load),\r\n SD = sd(training_load))\r\n\r\n## Average training load by pos\r\nd %>%\r\n group_by(pos) %>%\r\n summarize(avg = mean(training_load),\r\n SD = sd(training_load))\r\n\r\n## Average training load by team & position\r\nd %>%\r\n group_by(team, pos) %>%\r\n summarize(avg = mean(training_load),\r\n SD = sd(training_load)) %>%\r\n arrange(pos)\r\n\r\n## Plot\r\nd %>%\r\n ggplot(aes(x = training_load, fill = team)) +\r\n geom_density(alpha = 0.5) +\r\n facet_wrap(~pos)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.05.15\u202fPM-975x1024.png\")
<\/a> ![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.05.25\u202fPM.png\")
<\/a><\/p>\n\r\n## Mixed Model\r\nfit_lmer <- lmer(training_load ~ pos + (1 |team), data = d)\r\nsummary(fit_lmer)\r\ncoef(fit_lmer)\r\nfixef(fit_lmer)\r\nranef(fit_lmer)\r\nsigma(fit_lmer)\r\nhist(residuals(fit_lmer))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.09.06\u202fPM-1024x947.png\")
<\/a><\/p>\n\r\nsim_func_lmer <- function(n_teams = 3, \r\n n_positions = 2, \r\n n_players = 10, \r\n team1_fwd_avg = 130,\r\n team1_fwd_sd = 15,\r\n team1_mid_avg = 100,\r\n team1_mid_sd = 5,\r\n team2_fwd_avg = 150,\r\n team2_fwd_sd = 20,\r\n team2_mid_avg = 90,\r\n team2_mid_sd = 10,\r\n team3_fwd_avg = 180,\r\n team3_fwd_sd = 15,\r\n team3_mid_avg = 150,\r\n team3_mid_sd = 15){\r\n\r\n ## Simulated data frame\r\n team <- rep(c("Team1","Team2", "Team3"), each = n_players * n_positions)\r\n pos <- c(rep(c("M", "F"), each = n_players), rep(c("M", "F"), each = n_players), rep(c("M", "F"), each = n_players))\r\n player_id <- as.factor(round(seq(from = 100, to = 300, length = length(team)), 0))\r\n \r\n d <- data.frame(team, pos, player_id)\r\n\r\n # simulate training loads\r\n training_load <- c(rnorm(n = n_players, mean = team1_mid_avg, sd = team1_mid_sd),\r\n rnorm(n = n_players, mean = team1_fwd_avg, sd = team1_fwd_sd),\r\n rnorm(n = n_players, mean = team2_mid_avg, sd = team2_mid_sd),\r\n rnorm(n = n_players, mean = team2_fwd_avg, sd = team2_fwd_sd),\r\n rnorm(n = n_players, mean = team3_mid_avg, sd = team3_mid_sd),\r\n rnorm(n = n_players, mean = team3_fwd_avg, sd = team3_fwd_sd))\r\n \r\n ## construct the mixed model\r\n fit_lmer <- lmer(training_load ~ pos + (1 |team), data = d)\r\n summary(fit_lmer)\r\n}\r\n<\/pre>\n\r\nsim_func_lmer()\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.10.40\u202fPM.png\")
<\/a><\/p>\n\r\nteam_training_lmer <- replicate(n = 1000,\r\n sim_func_lmer(),\r\n simplify = FALSE)\r\n\r\n## look at the first model in the list\r\nteam_training_lmer[[1]]$coefficient\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.13.01\u202fPM.png\")
<\/a><\/p>\n\r\n\r\nlmer_coef <- matrix(NA, ncol = 5, nrow = length(team_training_lmer))\r\ncolnames(lmer_coef) <- c("intercept", "intercept_se", "posM", "posM_se", 'model_sigma')\r\n\r\nfor(i in 1:length(team_training_lmer)){\r\n \r\n lmer_coef[i, 1] <- team_training_lmer[[i]]$coefficients[1]\r\n lmer_coef[i, 2] <- team_training_lmer[[i]]$coefficients[3]\r\n lmer_coef[i, 3] <- team_training_lmer[[i]]$coefficients[2]\r\n lmer_coef[i, 4] <- team_training_lmer[[i]]$coefficients[4]\r\n lmer_coef[i, 5] <- team_training_lmer[[i]]$sigma\r\n \r\n}\r\n\r\nlmer_coef <- as.data.frame(lmer_coef) head(lmer_coef) ## Plot the coefficient for position lmer_coef %>%\r\n ggplot(aes(x = posM)) +\r\n geom_density(fill = "palegreen")\r\n\r\n## Summarize the coefficients and their standard errors for the simulations\r\nlmer_coef %>% \r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n\r\n## compare to the original model\r\nbroom.mixed::tidy(fit_lmer)\r\nsigma(fit_lmer)\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.23.58\u202fPM-1024x541.png\")
<\/a><\/p>\n\r\nranef_sim <- matrix(NA, ncol = 2, nrow = length(team_training_lmer))\r\ncolnames(ranef_sim) <- c("intercept_sd", "residual_sd")\r\n\r\nfor(i in 1:length(team_training_lmer)){\r\n \r\n ranef_sim[i, 1] <- team_training_lmer[[i]]$varcor %>% as.data.frame() %>% select(sdcor) %>% slice(1) %>% pull(sdcor)\r\n ranef_sim[i, 2] <- team_training_lmer[[i]]$varcor %>% as.data.frame() %>% select(sdcor) %>% slice(2) %>% pull(sdcor)\r\n \r\n}\r\n\r\nranef_sim <- as.data.frame(ranef_sim) head(ranef_sim) ## Summarize the results ranef_sim %>%\r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n\r\n## Compare with the original model\r\nVarCorr(fit_lmer)\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.25.57\u202fPM.png\")
<\/a><\/p>\n\r\n## Set up the data frame\r\nn_pos <- 3\r\nn_players <- 12\r\nn_obs <- 20\r\nplayers <- as.factor(round(seq(from = 100, to = 300, length = n_players), 0))\r\n\r\n\r\ndat <- data.frame( player_id = rep(players, each = n_obs), pos = rep(c("Fwd", "Mid", "Def"), each = n_players\/n_pos * n_obs), training_day = rep(1:n_obs, times = n_players) ) dat %>%\r\n head()\r\n\r\n## Create model parameters\r\n# NOTE: Defender will be the intercept\r\nset.seed(6687)\r\npos_intercept <- 150\r\nposF_coef <- 170\r\nposM_coef <- -70\r\nindividual_intercept <- 50\r\nindividual_slope <- 10\r\nsigma <- 10\r\nmodel_error <- rnorm(n = nrow(dat), mean = 0, sd = sigma)\r\n\r\n\r\n## we will also create some individual player variance\r\nindividual_player_variance <- c()\r\n\r\nfor(i in players){\r\n\r\n individual_player_variance[i] <- rnorm(n = 1, \r\n mean = runif(min = 2, max = 10, n = 1), \r\n sd = runif(min = 2, max = 5, n = 1))\r\n \r\n}\r\n\r\nindividual_player_variance <- rep(individual_player_variance, each = n_obs)\r\n\r\ndat$training_load <- pos_intercept + posF_coef * (dat$pos == "Fwd") + posM_coef * (dat$pos == "Mid") + individual_intercept + individual_slope * individual_player_variance + model_error dat %>%\r\n head()\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.27.32\u202fPM.png\")
<\/a>
\nCalculate summary stats<\/p>\n\r\n## Average training load by pos\r\ndat %>%\r\n group_by(pos) %>%\r\n summarize(avg = mean(training_load),\r\n SD = sd(training_load))\r\n\r\n## Plot\r\ndat %>%\r\n ggplot(aes(x = training_load, fill = pos)) +\r\n geom_density(alpha = 0.5)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.28.36\u202fPM.png\")
<\/a> ![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.28.50\u202fPM-967x1024.png\")
<\/a><\/p>\n\r\n## Mixed Model\r\nfit_lmer_pos <- lmer(training_load ~ pos + (1 | player_id), data = dat)\r\nsummary(fit_lmer_pos)\r\ncoef(fit_lmer_pos)\r\nfixef(fit_lmer_pos)\r\nranef(fit_lmer_pos)\r\nsigma(fit_lmer_pos)\r\nhist(residuals(fit_lmer_pos))\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.32.39\u202fPM-1024x566.png\")
<\/a>
\nCreate a function for the simulation<\/p>\n\r\nsim_func_lmer2 <- function(n_pos = 3,\r\n n_players = 12,\r\n n_obs = 20,\r\n pos_intercept = 150,\r\n posF_coef = 170,\r\n posM_coef = -70,\r\n individual_intercept = 50,\r\n individual_slope = 10,\r\n sigma = 10){\r\n \r\n players <- as.factor(round(seq(from = 100, to = 300, length = n_players), 0))\r\n\r\n dat <- data.frame(\r\n player_id = rep(players, each = n_obs),\r\n pos = rep(c("Fwd", "Mid", "Def"), each = n_players\/n_pos * n_obs),\r\n training_day = rep(1:n_obs, times = n_players)\r\n )\r\n \r\n model_error <- rnorm(n = nrow(dat), mean = 0, sd = sigma)\r\n \r\n individual_player_variance <- c()\r\n\r\n for(i in players){\r\n\r\n individual_player_variance[i] <- rnorm(n = 1, \r\n mean = runif(min = 2, max = 10, n = 1), \r\n sd = runif(min = 2, max = 5, n = 1))\r\n }\r\n\r\n individual_player_variance <- rep(individual_player_variance, each = n_obs)\r\n\r\n dat$training_load <- pos_intercept + posF_coef * (dat$pos == "Fwd") + posM_coef * (dat$pos == "Mid") + individual_intercept + individual_slope * individual_player_variance + model_error\r\n\r\n fit_lmer_pos <- lmer(training_load ~ pos + (1 | player_id), data = dat)\r\n summary(fit_lmer_pos)\r\n}\r\n<\/pre>\n\r\nplayer_training_lmer <- replicate(n = 1000,\r\n sim_func_lmer2(),\r\n simplify = FALSE)\r\n\r\n## look at the first model in the list\r\nplayer_training_lmer[[1]]$coefficient\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.35.56\u202fPM.png\")
<\/a><\/p>\n\r\nlmer_player_coef <- matrix(NA, ncol = 7, nrow = length(player_training_lmer))\r\ncolnames(lmer_player_coef) <- c("intercept", "intercept_se","posFwd", "posFwd_se", "posMid", "posMid_se", 'model_sigma')\r\n\r\nfor(i in 1:length(player_training_lmer)){\r\n \r\n lmer_player_coef[i, 1] <- player_training_lmer[[i]]$coefficients[1]\r\n lmer_player_coef[i, 2] <- player_training_lmer[[i]]$coefficients[4]\r\n lmer_player_coef[i, 3] <- player_training_lmer[[i]]$coefficients[2]\r\n lmer_player_coef[i, 4] <- player_training_lmer[[i]]$coefficients[5]\r\n lmer_player_coef[i, 5] <- player_training_lmer[[i]]$coefficients[3]\r\n lmer_player_coef[i, 6] <- player_training_lmer[[i]]$coefficients[6]\r\n lmer_player_coef[i, 7] <- player_training_lmer[[i]]$sigma\r\n \r\n}\r\n\r\nlmer_player_coef <- as.data.frame(lmer_player_coef) head(lmer_player_coef) ## Plot the coefficient for position lmer_player_coef %>%\r\n ggplot(aes(x = posFwd)) +\r\n geom_density(fill = "palegreen")\r\n\r\nlmer_player_coef %>%\r\n ggplot(aes(x = posMid)) +\r\n geom_density(fill = "palegreen")\r\n\r\n\r\n## Summarize the coefficients and their standard errors for the simulations\r\nlmer_player_coef %>% \r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n\r\n## compare to the original model\r\nbroom.mixed::tidy(fit_lmer_pos)\r\nsigma(fit_lmer_pos)\r\n<\/pre>\n![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.39.43\u202fPM-1024x476.png\")
<\/a>
\nExtract the random effects for the intercept and the residual for each of the simulated models.<\/p>\n\r\nranef_sim_player <- matrix(NA, ncol = 2, nrow = length(player_training_lmer))\r\ncolnames(ranef_sim_player) <- c("player_sd", "residual_sd")\r\n\r\nfor(i in 1:length(player_training_lmer)){\r\n \r\n ranef_sim_player[i, 1] <- player_training_lmer[[i]]$varcor %>% as.data.frame() %>% select(sdcor) %>% slice(1) %>% pull(sdcor)\r\n ranef_sim_player[i, 2] <- player_training_lmer[[i]]$varcor %>% as.data.frame() %>% select(sdcor) %>% slice(2) %>% pull(sdcor)\r\n \r\n}\r\n\r\nranef_sim_player <- as.data.frame(ranef_sim_player) head(ranef_sim_player) ## Summarize the results ranef_sim_player %>%\r\n summarize(across(.cols = everything(),\r\n ~mean(.x)))\r\n\r\n## Compare with the original model\r\nVarCorr(fit_lmer_pos)\r\n<\/pre>\n
![\"\"](\"https:\/\/optimumsportsperformance.com\/blog\/wp-content\/uploads\/2023\/11\/Screenshot-2023-11-05-at-8.42.08\u202fPM-1024x380.png\")
<\/a><\/p>\n