if (FALSE) {
#----------------------------------------------------------------------------
# Mplus Example 3.18: Moderated Mediation with a Plot of the Indirect Effect
#..........
# Trace Plots
# Example 1a: Default setting
mplus.plot("ex3.18.gh5")
# Example 1b: Exclude first half of each chain
mplus.plot("ex3.18.gh5", burnin = FALSE)
# Example 1c: Print all parameters
mplus.plot("ex3.18.gh5", param = "all")
# Example 1d: Print user-specified parameters
mplus.plot("ex3.18.gh5", param = "param")
# Example 1e: Arrange panels in three columns
mplus.plot("ex3.18.gh5", ncol = 3)
# Example 1f: Specify "Pastel 1" palette for the hcl.colors function
mplus.plot("ex3.18.gh5", palette = "Pastel 1")
#..........
# Posterior Distribution Plots
# Example 2a: Default setting, i.e., posterior median and equal-tailed interval
mplus.plot("ex3.18.gh5", plot = "post")
# Example 2b: Display posterior mean and maximum a posteriori
mplus.plot("ex3.18.gh5", plot = "post", point = c("m", "map"))
# Example 2c: Display maximum a posteriori and highest density interval
mplus.plot("ex3.18.gh5", plot = "post", point = "map", ci = "hdi")
# Example 2d: Do not display any point estimates and credible interval
mplus.plot("ex3.18.gh5", plot = "post", point = "none", ci = "none")
# Example 2d: Do not display histograms
mplus.plot("ex3.18.gh5", plot = "post", hist = FALSE)
#..........
# Autocorrelation Plots
# Example 3a: Default setting, i.e., first chain
mplus.plot("ex3.18.gh5", plot = "auto")
# Example 3b: Use second chain
mplus.plot("ex3.18.gh5", plot = "auto", chain = 2)
# Example 3b: Modify limits and breaks of the y-axis
mplus.plot("ex3.18.gh5", plot = "auto",
ylim = c(-0.05, 0.05), ybreaks = seq(-0.1, 0.1, by = 0.025))
#..........
# Posterior Predictive Check Plots
# Example 4a: Default setting, i.e., 95% Interval
mplus.plot("ex3.18.gh5", plot = "ppc")
# Example 4b: Default setting, i.e., 99% Interval
mplus.plot("ex3.18.gh5", plot = "ppc", conf.level = 0.99)
#..........
# Loop Plot
# Example 5a: Default setting
mplus.plot("ex3.18.gh5", plot = "loop")
# Example 5b: Do not fill area and draw vertical lines
mplus.plot("ex3.18.gh5", plot = "loop", area = FALSE)
#..........
# Save Plots
# Example 6a: Save all plots in pdf format
mplus.plot("ex3.18.gh5", saveplot = "all")
# Example 6b: Save all plots in png format with 300 dpi
mplus.plot("ex3.18.gh5", saveplot = "all", file = "Mplus_Plot.png", dpi = 300)
# Example 6a: Save loop plot, specify width and height of the plot
mplus.plot("ex3.18.gh5", plot = "none", saveplot = "loop",
width = 7.5, height = 7)
#----------------------------------------------------------------------------
# Plot from misty.object
# Create misty.object
object <- mplus.plot("ex3.18.gh5", plot = "none")
# Trace plot
mplus.plot(object, plot = "trace")
# Posterior distribution plot
mplus.plot(object, plot = "post")
# Autocorrelation plot
mplus.plot(object, plot = "auto")
# Posterior predictive check plot
mplus.plot(object, plot = "ppc")
# Loop plot
mplus.plot(object, plot = "loop")
#----------------------------------------------------------------------------
# Create Plots Manually
# Load ggplot2 package
library(ggplot2)
# Create misty object
object <- mplus.plot("ex3.18.gh5", plot = "none")
#..........
# Example 7: Trace Plots
# Extract data in long format
data.post <- object$data$post$long
# Extract ON parameters
data.trace <- data.post[grep(" ON ", data.post$param), ]
# Plot
ggplot(data.trace, aes(x = iter, y = value, color = chain)) +
annotate("rect", xmin = 0, xmax = 15000, ymin = -Inf, ymax = Inf,
alpha = 0.4, fill = "gray85") +
geom_line() +
facet_wrap(~ param, ncol = 2, scales = "free") +
scale_x_continuous(name = "", expand = c(0.02, 0)) +
scale_y_continuous(name = "", expand = c(0.02, 0)) +
scale_colour_manual(name = "Chain",
values = hcl.colors(n = 2, palette = "Set 2")) +
theme_bw() +
guides(color = guide_legend(nrow = 1, byrow = TRUE)) +
theme(plot.margin = margin(c(4, 15, -10, 0)),
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.box.margin = margin(c(-16, 6, 6, 6)),
legend.background = element_rect(fill = "transparent"))
#..........
# Example 8: Posterior Distribution Plots
# Extract data in long format
data.post <- object$data$post$long
# Extract ON parameters
data.post <- data.post[grep(" ON ", data.post$param), ]
# Discard burn-in iterations
data.post <- data.post[data.post$iter > 15000, ]
# Drop factor levels
data.post$param <- droplevels(data.post$param,
exclude = c("[Y]", "[M]", "Y", "M", "INDIRECT", "MOD"))
# Plot
ggplot(data.post, aes(x = value)) +
geom_histogram(aes(y = after_stat(density)), color = "black", alpha = 0.4,
fill = "gray85") +
geom_density(color = "#0072B2") +
geom_vline(data = data.frame(param = unique(data.post$param),
stat = tapply(data.post$value, data.post$param, median)),
aes(xintercept = stat, color = "Median"), linewidth = 0.6) +
geom_vline(data = data.frame(param = unique(data.post$param),
low = tapply(data.post$value, data.post$param,
function(y) quantile(y, probs = 0.025))),
aes(xintercept = low), linetype = "dashed", linewidth = 0.6) +
geom_vline(data = data.frame(param = unique(data.post$param),
upp = tapply(data.post$value, data.post$param,
function(y) quantile(y, probs = 0.975))),
aes(xintercept = upp), linetype = "dashed", linewidth = 0.6) +
facet_wrap(~ param, ncol = 2, scales = "free") +
scale_x_continuous(name = "", expand = c(0.02, 0)) +
scale_y_continuous(name = "Probability Density, f(x)",
expand = expansion(mult = c(0L, 0.05))) +
scale_color_manual(name = "Point Estimate", values = c(Median = "#D55E00")) +
labs(caption = "95% Equal-Tailed Interval") +
theme_bw() +
theme(plot.margin = margin(c(4, 15, -8, 4)),
plot.caption = element_text(hjust = 0.5, vjust = 7),
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.box.margin = margin(c(-30, 6, 6, 6)),
legend.background = element_rect(fill = "transparent"))
#..........
# Example 9: Autocorrelation Plots
# Extract data in long format
data.auto <- object$data$auto$long
# Select first chain
data.auto <- data.auto[data.auto$chain == 1, ]
# Extract ON parameters
data.auto <- data.auto[grep(" ON ", data.auto$param), ]
# Plot
ggplot(data.auto, aes(x = lag, y = cor)) +
geom_bar(stat = "identity", alpha = 0.4, color = "black", fill = "gray85",
width = 0.8) +
facet_wrap(~ param, ncol = 2) +
scale_x_continuous(name = "Lag", breaks = seq(1, 30, by = 2), expand = c(0.02, 0)) +
scale_y_continuous(name = "Autocorrelation", limits = c(-0.1, 0.1),
breaks = seq(-0.1, 1., by = 0.05), expand = c(0.02, 0)) +
theme_bw() +
theme(plot.margin = margin(c(4, 15, 4, 4)))
#..........
# Example 10: Posterior Predictive Check (PPC) Plots
# Extract data
data.ppc <- object$data$ppc
# Scatter plot
ppc.scatter <- ggplot(data.ppc, aes(x = obs, y = rep)) +
geom_point(shape = 21, fill = "gray85") +
geom_abline(slope = 1) +
scale_x_continuous("Observed", limits = c(0, 45), breaks = seq(0, 45, by = 5),
expand = c(0.02, 0)) +
scale_y_continuous("Recpliated", limits = c(0, 45), breaks = seq(0, 45, by = 5),
expand = c(0.02, 0)) +
theme_bw() +
theme(plot.margin = margin(c(2, 15, 4, 4)))
# Histogram
ppc.hist <- ggplot(data.ppc, aes(x = diff)) +
geom_histogram(color = "black", alpha = 0.4, fill = "gray85") +
geom_vline(xintercept = mean(data.ppc$diff), color = "#CC79A7") +
geom_vline(xintercept = quantile(data.ppc$diff, probs = 0.025),
linetype = "dashed", color = "#CC79A7") +
geom_vline(xintercept = quantile(data.ppc$diff, probs = 0.975),
linetype = "dashed", color = "#CC79A7") +
scale_x_continuous("Observed - Replicated", expand = c(0.02, 0)) +
scale_y_continuous("Count", expand = expansion(mult = c(0L, 0.05))) +
theme_bw() +
theme(plot.margin = margin(c(2, 15, 4, 4)))
# Combine plots using the patchwork package
patchwork::wrap_plots(ppc.scatter, ppc.hist)
#..........
# Example 11: Loop Plot
# Extract data
data.loop <- object$data$loop
# Plot
plot.loop <- ggplot(data.loop, aes(x = xval, y = estimate)) +
geom_line(linewidth = 0.6, show.legend = FALSE) +
geom_line(aes(xval, low)) +
geom_line(aes(xval, upp)) +
scale_x_continuous("MOD", expand = c(0.02, 0)) +
scale_y_continuous("INDIRECT", expand = c(0.02, 0)) +
scale_fill_manual("Statistical Significance",
values = hcl.colors(n = 2, palette = "Set 2")) +
theme_bw() +
theme(plot.margin = margin(c(4, 15, -6, 4)),
legend.position = "bottom",
legend.title = element_text(size = 10),
legend.text = element_text(size = 10),
legend.box.margin = margin(-10, 6, 6, 6),
legend.background = element_rect(fill = "transparent"))
# Significance area
for (i in unique(data.loop$group)) {
plot.loop <- plot.loop + geom_ribbon(data = data.loop[data.loop$group == i, ],
aes(ymin = low, ymax = upp, fill = sig), alpha = 0.4)
}
# Vertical lines
plot.loop + geom_vline(data = data.loop[data.loop$change == 1, ],
aes(xintercept = xval, color = sig), linewidth = 0.6,
linetype = "dashed", show.legend = FALSE)
}
Run the code above in your browser using DataLab