Code

setwd("D:/OneDrive - HKUST Connect/practice/mixed models")

require("pacman")

p_load(sjPlot, DT, gridExtra, ggthemes, ggsci, tidyverse, stargazer, lme4, gt, gtExtras, skimr, kableExtra, merTools, tidyr, tidytext, broom, jtools, ggstance, huxtable, tidyfst)

Code

fread("shrimp.csv", sep = ",", header = T, stringsAsFactors = T) -> shrimp

datatable(shrimp, 
          filter = "top", 
          options = list(pageLength = 5, autoWidth = TRUE))

Code

shrimp %>% 
  ggplot(aes(SexID, M2BW, color = SexID)) + 
  geom_boxplot() + 
  geom_dotplot(binaxis = "y", 
               stackdir = "center", 
               position = "dodge", 
               binwidth = 0.17) +
  scale_color_d3() + 
  theme_bw() +
  theme(legend.position = "none") + 
  labs(x = NULL,
       y = "weight") -> dp_sex

shrimp %>% 
  ggplot(aes(TankID, M2BW, color = TankID)) + 
  geom_boxplot() + 
  geom_dotplot(binaxis = "y", 
               stackdir = "center", 
               position = "dodge", 
               binwidth = 0.17) +
  scale_color_d3() + 
  theme_bw() +
  theme(legend.position = "none") + 
  labs(x = NULL,
       y = NULL) -> dp_tank

grid.arrange(dp_sex, dp_tank, ncol = 2)

Code

shrimp %>% 
  lm(M2BW ~ SexID, data = .) -> sex

shrimp %>% 
  lm(M2BW ~ SexID + TankID, data = .) -> sex_tank

shrimp %>% 
  lm(M2BW ~ SexID + TankID + M1BW, data = .) -> sex_tank_before

shrimp %>% 
  lm(M2BW ~ SexID + TankID + M1BW + TankID:M1BW, data = .) -> sex_tank_before_inter

shrimp %>% 
  lm(M2BW ~ SexID + TankID + M1BW + TankID:M1BW + SexID:M1BW , data = .) -> sex_tank_before_inter_two

export_summs(sex, sex_tank, sex_tank_before, sex_tank_before_inter, sex_tank_before_inter_two, scale = TRUE, results = "asis")

	Model 1	Model 2	Model 3	Model 4	Model 5
(Intercept)	34.38 ***	35.86 ***	35.86 ***	35.86 ***	35.86 ***
	(0.11)	(0.13)	(0.12)	(0.12)	(0.12)
SexIDMale	-6.21 ***	-6.21 ***	-6.23 ***	-6.23 ***	-6.23 ***
	(0.16)	(0.15)	(0.14)	(0.14)	(0.14)
TankIDT2		-2.96 ***	-2.95 ***	-2.95 ***	-2.95 ***
		(0.15)	(0.14)	(0.14)	(0.14)
M1BW			1.32 ***	1.38 ***	1.45 ***
			(0.07)	(0.10)	(0.12)
TankIDT2:M1BW				-0.12	-0.12
				(0.14)	(0.14)
SexIDMale:M1BW					-0.15
					(0.15)
N	4282	4282	4241	4241	4241
R2	0.27	0.33	0.38	0.38	0.38
All continuous predictors are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

shrimp %>% 
  group_by(FamilyID) %>% 
  summarise(M1BW = mean(M1BW, na.rm = T),
            M2BW = mean(M2BW, na.rm = T)) %>% 
  arrange(-desc(M1BW)) -> per_family

per_family %>% 
  ggplot(aes(M1BW, M2BW)) +
  geom_point(shape = 20) + 
  geom_smooth(method = "lm", color = "#FF0000") + 
  theme_bw() +
  theme(legend.position = "none") + 
  labs(x = "before",
       y = "after")

`geom_smooth()` using formula 'y ~ x'

Interaction effects represent the combined effects of factors on the dependent measure. When an interaction effect is present, the impact of one factor depends on the level of the other factor.

Code

shrimp %>% 
  lm(M2BW ~ SexID + SexID:M1BW, data = .) -> model_1

shrimp %>% 
  lm(M2BW ~ SexID + TankID + SexID:M1BW, data = .) -> model_2

shrimp %>% 
  lm(M2BW ~ SexID + TankID + SexID:M1BW + SexID:TankID, data = .) -> model_3

export_summs(model_1, model_2, model_3, scale = TRUE, results = "asis")

	Model 1	Model 2	Model 3
(Intercept)	34.38 ***	35.86 ***	35.92 ***
	(0.11)	(0.12)	(0.14)
SexIDMale	-6.23 ***	-6.23 ***	-6.35 ***
	(0.15)	(0.14)	(0.20)
SexIDFemale:M1BW	1.36 ***	1.39 ***	1.39 ***
	(0.10)	(0.10)	(0.10)
SexIDMale:M1BW	1.30 ***	1.24 ***	1.24 ***
	(0.11)	(0.11)	(0.11)
TankIDT2		-2.95 ***	-3.07 ***
		(0.14)	(0.20)
SexIDMale:TankIDT2			0.25
			(0.29)
N	4241	4241	4241
R2	0.32	0.38	0.38
All continuous predictors are mean-centered and scaled by 1 standard deviation. * p < 0.001; p < 0.01; * p < 0.05.

Code

plot_model(model_3, show.p = T, show.values = T, title = "") + 
  geom_hline(yintercept = 0, linetype = "longdash", color = "#ffa200") + 
  theme_bw()

1 Mixed Model

Code

shrimp %>% 
  ggplot(aes(PopID, M2BW, color = PopID)) + 
  geom_boxplot() + 
  geom_jitter(alpha = 0.2) +
  scale_color_d3() + 
  theme_bw() +
  theme(legend.position = "none") + 
  labs(x = NULL,
       y = NULL)

Code

shrimp %>% 
  ggplot(aes(PopID, M2BW, fill = FamilyID)) + 
  geom_boxplot(outlier.size = 0) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_color_d3() + 
  labs(x = NULL)