Chapter 18:
Phil Wood
2025-08-16
Probability Density functions for the Logistic and Normal Distributions (Figure 18.1)
library(ggplot2)
library(dplyr)
# ---- Parameters ----
mu <- 0
sd_norm <- 1
# Logistic (unstandardized: scale = 1)
loc <- 0
scale <- 1
sd_logi <- pi / sqrt(3) # ~1.814
# ---- Grid & densities ----
x <- seq(-8, 8, length.out = 2001)
df_norm <- data.frame(
x = x,
density = dnorm(x, mean = mu, sd = sd_norm),
dist = "Normal(0, 1)"
) |>
mutate(
mean = mu,
sd = sd_norm,
crit_lo = qnorm(0.025, mean = mu, sd = sd_norm),
crit_hi = qnorm(0.975, mean = mu, sd = sd_norm)
)
df_logi <- data.frame(
x = x,
density = dlogis(x, location = loc, scale = scale),
dist = "Logistic(0, 1)"
) |>
mutate(
mean = loc,
sd = sd_logi,
crit_lo = qlogis(0.025, location = loc, scale = scale),
crit_hi = qlogis(0.975, location = loc, scale = scale)
)
df <- bind_rows(df_norm, df_logi)
# For shading tails
df_shade <- df |>
mutate(in_tail = (x <= crit_lo) | (x >= crit_hi))
# Label positions
labels_df <- df |>
group_by(dist) |>
summarize(
mean = first(mean),
sd = first(sd),
crit_lo = first(crit_lo),
crit_hi = first(crit_hi),
y_text = max(density) * 0.85,
y_crit = max(density) * 0.08,
.groups = "drop"
)
# ---- Plot ----
p <- ggplot(df, aes(x = x, y = density)) +
geom_ribbon(
data = df_shade |> mutate(ymin = 0, ymax = ifelse(in_tail, density, NA_real_)),
aes(ymin = ymin, ymax = ymax),
alpha = 0.25
) +
geom_line(linewidth = 1) +
geom_vline(aes(xintercept = mean), linetype = "solid", linewidth = 0.7) +
geom_vline(aes(xintercept = mean - sd), linetype = "dashed", linewidth = 0.6) +
geom_vline(aes(xintercept = mean + sd), linetype = "dashed", linewidth = 0.6) +
geom_vline(aes(xintercept = crit_lo), linetype = "dotdash", linewidth = 0.6) +
geom_vline(aes(xintercept = crit_hi), linetype = "dotdash", linewidth = 0.6) +
geom_text(
data = labels_df,
aes(x = mean, y = y_text, label = paste0("Mean = ", format(mean, digits = 3))),
vjust = 2.5, fontface = "bold"
) +
geom_text(
data = labels_df,
aes(x = mean, y = y_text, label = paste0("SD = ", round(sd, 3))),
vjust = 3.5
) +
geom_text(
data = labels_df,
aes(x = crit_lo, y = y_crit, label = paste0("2.5% = ", round(crit_lo, 3))),
hjust = 1.1, vjust = -.5
) +
geom_text(
data = labels_df,
aes(x = crit_hi, y = y_crit, label = paste0("97.5% = ", round(crit_hi, 3))),
hjust = -0.1, vjust = -.5
) +
facet_wrap(~ dist, ncol = 1, scales = "free_y") +
labs(
title = "Normal(0, 1) vs Logistic(0, 1)\nPDFs with Mean, +/- 1 SD, and 95% Critical Regions",
x = "x",
y = "Density"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.minor = element_blank(),
strip.text = element_text(face = "bold"),
plot.title = element_text(margin = margin(b = 8), lineheight = 1.15) # add space + nicer line spacing
)
print(p)
## PDFs
summary — Normal(μ = 0, σ = 1) vs Logistic(loc = 0, scale = 1, σ ≈
1.814) — 95% cutoffs: N[-1.96; 1.96], L[-3.664; 3.664] {#sec-pdfs}
Cut-off score of 0.5 for logistic distribution (Figure 18.2)
# ---- Packages ----
library(ggplot2)
library(dplyr)
# ---- Parameters ----
loc <- 0
scale <- 1
cut <- 0.5
# Areas under the curve using the logistic CDF
p_left <- plogis(cut, location = loc, scale = scale) # P(X <= cut)
p_right <- 1 - p_left # P(X >= cut)
# ---- Grid & density ----
x <- seq(-8, 8, length.out = 4001)
df <- data.frame(
x = x,
dens = dlogis(x, location = loc, scale = scale)
)
# For shading the two regions
df_left <- df %>% filter(x <= cut)
df_right <- df %>% filter(x >= cut)
# Positions for annotation labels
y_max <- max(df$dens)
label_y <- y_max * 0.35
label_left_x <- mean(range(df_left$x))
label_right_x <- mean(range(df_right$x))
# ---- Plot ----
ggplot(df, aes(x = x, y = dens)) +
# Shaded areas
geom_ribbon(data = df_left, aes(ymin = 0, ymax = dens, fill = "x <= 0.5"), alpha = 0.3) +
geom_ribbon(data = df_right, aes(ymin = 0, ymax = dens, fill = "x >= 0.5"), alpha = 0.3) +
# PDF curve
geom_line(size = 1) +
# Vertical reference line at cutoff
geom_vline(xintercept = cut, linetype = "dashed", linewidth = 0.7) +
# Annotations for areas
annotate("text",
x = label_left_x, y = label_y,
label = paste0("P(X <= 0.3) = ", sprintf("%.1f%%", 100 * p_left)),
hjust = 0.7, vjust = 0) +
annotate("text",
x = label_right_x, y = label_y,
label = paste0("P(X >= 0.3) = ", sprintf("%.1f%%", 100 * p_right)),
hjust = 0.5, vjust = 0) +
# Labels and styling
scale_fill_discrete(name = "Region") +
labs(
title = sprintf("Figure 18.2 Logistic(0, %g) PDF with Areas at Cutoff = %.1f", scale, cut),
x = "x",
y = "Density"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.minor = element_blank(),
legend.position = "top",
plot.title = element_text(margin = margin(b = 8), lineheight = 1.15)
)
#
Cutoff for the Normal Distribution (Figure 18.3)
# ---- Packages ----
library(ggplot2)
library(dplyr)
# ---- Parameters ----
mu <- 0
sigma <- 1
cut <- 0.3
# Areas using the Normal CDF
p_left <- pnorm(cut, mean = mu, sd = sigma) # P(X <= cut)
p_right <- 1 - p_left # P(X >= cut)
# ---- Grid & density ----
x <- seq(-4, 4, length.out = 4001)
df <- data.frame(
x = x,
dens = dnorm(x, mean = mu, sd = sigma)
)
# For shading
df_left <- df %>% filter(x <= cut)
df_right <- df %>% filter(x >= cut)
# Positions for annotation
y_max <- max(df$dens)
label_y <- y_max * 0.35
label_left_x <- mean(range(df_left$x))
label_right_x <- mean(range(df_right$x))
# ---- Plot ----
ggplot(df, aes(x = x, y = dens)) +
# Shaded areas
geom_ribbon(data = df_left, aes(ymin = 0, ymax = dens, fill = "x <= 0.3"), alpha = 0.3) +
geom_ribbon(data = df_right, aes(ymin = 0, ymax = dens, fill = "x >= 0.3"), alpha = 0.3) +
# PDF curve
geom_line(size = 1) +
# Vertical cutoff line
geom_vline(xintercept = cut, linetype = "dashed", linewidth = 0.7) +
# Percent annotations
annotate("text",
x = label_left_x, y = label_y,
label = paste0("P(X <= 0.3) = ", sprintf("%.1f%%", 100 * p_left)),
hjust = 1, vjust = 0) +
annotate("text",
x = label_right_x, y = label_y,
label = paste0("P(X >= 0.3) = ", sprintf("%.1f%%", 100 * p_right)),
hjust = 0, vjust = 0) +
# Labels and styling
scale_fill_discrete(name = "Region") +
labs(
title = "Figure 18.3 Normal(0,1) PDF with Areas at Cutoff = 0.3",
x = "x",
y = "Density"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.minor = element_blank(),
legend.position = "top"
)
Cumulative Distribution functions for the Normal and Logistic curves for the cutoff-examples (Figure 18.4)
# ---- Packages ----
library(ggplot2)
library(dplyr)
library(tidyr)
# ---- Grid and CDFs ----
x <- seq(-3, 3, length.out = 2001)
df <- tibble(
x = x,
Normal = pnorm(x, mean = 0, sd = 1),
Logistic = plogis(x, location = 0, scale = 1) # Logistic(0,1)
) |>
pivot_longer(cols = c(Normal, Logistic), names_to = "Distribution", values_to = "CDF")
# ---- Cutoffs ----
cuts <- c(0.3, 0.5)
# Annotation dataframe for vertical lines
ann_v <- tibble(
x = cuts,
y = 0.05, # near the bottom of the plot
label = paste0("x = ", cuts)
)
# ---- Horizontal reference ----
hline_y <- 0.62
ann_h <- tibble(
x = -2.9, # place near left edge
y = hline_y,
label = paste0("CDF = ", hline_y)
)
# ---- Plot ----
ggplot(df, aes(x = x, y = CDF, color = Distribution)) +
geom_line(size = 1) +
# Vertical reference lines
geom_vline(xintercept = cuts, linetype = "dashed", linewidth = 0.7) +
geom_text(data = ann_v, aes(x = x, y = y, label = label),
inherit.aes = FALSE, angle = 90, vjust = -0.5, hjust = 0,
size = 4, color = "black") +
# Horizontal reference line
geom_hline(yintercept = hline_y, linetype = "dashed", linewidth = 0.7) +
geom_text(data = ann_h, aes(x = x, y = y, label = label),
inherit.aes = FALSE, vjust = -0.8, hjust = 0, size = 4, color = "black") +
# Axes
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
scale_x_continuous(limits = c(-3, 3), breaks = -3:3) +
# Labels and styling
labs(
title = "Figure 18.4 CDFs of Normal(0,1) and Logistic(0,1)",
x = "x",
y = "Cumulative probability"
) +
theme_minimal(base_size = 14) +
theme(
panel.grid.minor = element_blank(),
legend.position = "top"
)
Logistic Regression in Mplus
Logit Link function
Using MplusAutomation
<a https://drive.google.com/file/d/110skb1FzgsZn-UeGvbbrK6SowNqR5R9X/view?usp=drive_link
<a https://docs.google.com/document/d/12Mv7XMAkTZKIkhkHQ5EZHJed4QTAaulnw-m37WWYUSo/edit?usp=drive_link
<a https://drive.google.com/file/d/1h2rz4F9Go1akOfe_TkgH4r-UPGATFDCu/view?usp=drive_link
(You’ll need to have Mplus to look at this file)
<a https://drive.google.com/file/d/1QorMQoe8317sAielktX12KsoD6ZI15Ha/view?usp=drive_link
library(MplusAutomation)
library(knitr)
# ---- Make knitting behave like your Rmd folder ----
if (!is.null(knitr::current_input())) {
knitr::opts_knit$set(root.dir = dirname(knitr::current_input()))
}
# ---- Tell MplusAutomation exactly where Mplus is (edit if needed) ----
if (.Platform$OS.type == "windows") {
options(Mplus_command = "C:/Program Files/Mplus/Mplus.exe")
} else if (Sys.info()[["sysname"]] == "Darwin") {
# Common macOS install path (adjust if different)
options(Mplus_command = "/Applications/Mplus/mplus")
} else {
# Linux/other (adjust if different)
options(Mplus_command = "/usr/local/bin/mplus")
}
stopifnot(!is.null(getOption("Mplus_command")),
nzchar(getOption("Mplus_command")),
file.exists(getOption("Mplus_command")))
# ---- File names ----
inp_file <- "rfdendplus5Sim.inp"
out_file <- sub("\\.inp$", ".out", inp_file)
stopifnot(file.exists(inp_file))
# ---- Run & read model ----
runModels(target = inp_file, showOutput = FALSE, replaceOutfile = "always")
mod <- readModels(target = out_file)
# ---- Helpers ----
fmt <- function(x, digits = 3) ifelse(is.na(x), "", formatC(x, format = "f", digits = digits))
fmt_p <- function(p) {
ifelse(is.na(p), "",
ifelse(p < .001, "< .001", paste0("= ", sub("^0", ".", formatC(p, format="f", digits=3)))))
}
# ---- Table 1: Model fit ----
S <- if (is.data.frame(mod$summaries)) mod$summaries[1, ] else as.data.frame(mod$summaries)
getv <- function(nm) if (nm %in% names(S)) as.numeric(S[[nm]]) else NA_real_
gets <- function(nm) if (nm %in% names(S)) as.character(S[[nm]]) else ""
chisq <- getv("ChiSqM_Value")
df <- getv("ChiSqM_DF")
pp <- getv("ChiSqM_PValue")
cfi <- getv("CFI")
tli <- getv("TLI")
rmsea <- getv("RMSEA_Estimate")
rm_lb <- getv("RMSEA_90CI_LB")
rm_ub <- getv("RMSEA_90CI_UB")
srmr <- getv("SRMR")
aic <- getv("AIC")
bic <- getv("BIC")
estm <- gets("Estimator")
atype <- gets("AnalysisType")
fit_row <- data.frame(
`χ² (df)` = if (is.na(chisq) | is.na(df)) "" else paste0(fmt(chisq, 2), " (", df, ")"),
`p` = fmt_p(pp),
`CFI` = fmt(cfi, 3),
`TLI` = fmt(tli, 3),
`RMSEA [90% CI]` = if (is.na(rmsea)) "" else paste0(fmt(rmsea, 3), " [", fmt(rm_lb, 3), ", ", fmt(rm_ub, 3), "]"),
`SRMR` = fmt(srmr, 3),
`AIC` = fmt(aic, 2),
`BIC` = fmt(bic, 2),
`Estimator` = estm,
`Analysis` = atype,
check.names = FALSE
)
kable(fit_row, align = "l", caption = "Table 1. Overall model fit.")| χ² (df) | p | CFI | TLI | RMSEA [90% CI] | SRMR | AIC | BIC | Estimator | Analysis |
|---|---|---|---|---|---|---|---|---|---|
| 2150.49 | 2161.66 | MLR | MISSING |
# ---- Table 2: Unstandardized & Standardized (STDYX) estimates ----
# Safely grab parameter frames if present
par_unstd <- try(mod$parameters$unstandardized, silent = TRUE)
if (inherits(par_unstd, "try-error") || is.null(par_unstd)) par_unstd <- data.frame()
par_stdyx <- try(mod$parameters$stdyx.standardized, silent = TRUE)
if (inherits(par_stdyx, "try-error") || is.null(par_stdyx)) par_stdyx <- data.frame()
# Helper accessors for column names that vary across Mplus versions
getcol <- function(df, nm, default = NA_real_) {
if (length(df) == 0) return(rep(default, 0))
if (nm %in% names(df)) return(df[[nm]])
# common alternates
alts <- c(est = c("est", "Estimate"),
se = c("se", "S.E."),
p = c("pval","P-Value","PValue"),
z = c("est_se","Z", "Z-value"))
if (nm %in% names(alts)) {
for (cand in alts[[nm]]) if (cand %in% names(df)) return(df[[cand]])
}
rep(default, nrow(df))
}
# Build unstandardized table
if (nrow(par_unstd)) {
U <- data.frame(
Section = if ("paramHeader" %in% names(par_unstd)) par_unstd$paramHeader else "",
Parameter = if ("param" %in% names(par_unstd)) par_unstd$param else "",
Estimate = as.numeric(getcol(par_unstd, "est")),
SE = as.numeric(getcol(par_unstd, "se")),
z = suppressWarnings(as.numeric(getcol(par_unstd, "z"))),
p = as.numeric(getcol(par_unstd, "p")),
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(is.na(U$z)) && any(is.finite(U$Estimate) & is.finite(U$SE) & U$SE > 0)) {
U$z <- with(U, ifelse(is.finite(Estimate) & is.finite(SE) & SE > 0, Estimate/SE, NA_real_))
}
if (all(is.na(U$p)) && any(is.finite(U$z))) {
U$p <- 2 * pnorm(-abs(U$z))
}
} else {
U <- data.frame()
}
# Build STDYX table and merge
if (nrow(par_stdyx)) {
Sx <- data.frame(
Section = if ("paramHeader" %in% names(par_stdyx)) par_stdyx$paramHeader else "",
Parameter = if ("param" %in% names(par_stdyx)) par_stdyx$param else "",
STDYX = as.numeric(getcol(par_stdyx, "est")),
stringsAsFactors = FALSE, check.names = FALSE
)
} else {
Sx <- data.frame(Section=character(), Parameter=character(), STDYX=numeric(), check.names = FALSE)
}
PE <- if (nrow(U)) merge(U, Sx, by = c("Section","Parameter"), all.x = TRUE) else data.frame()
# Keep common, interpretable sections first
section_order <- c("LOADINGS", "BY", "Latent Variable Loadings", "REGRESSIONS", "WITH",
"Latent Variable Covariances", "Residual Variances", "Variances",
"Means", "Intercepts", "Thresholds")
if (nrow(PE)) {
PE$Section <- as.character(PE$Section)
ord <- match(PE$Section, section_order)
PE <- PE[order(ifelse(is.na(ord), length(section_order)+1, ord), PE$Section, PE$Parameter), ]
# Format numbers for APA
PE_out <- transform(PE,
Estimate = fmt(Estimate, 3),
SE = fmt(SE, 3),
z = ifelse(is.na(z), "", fmt(z, 2)),
p = sapply(p, fmt_p),
STDYX = fmt(STDYX, 3))
kable(PE_out[, c("Section","Parameter","Estimate","SE","z","p","STDYX")],
align = "l", caption = "Table 2. Unstandardized and standardized (STDYX) parameter estimates.")
} else {
cat("*(No parameter estimates were found in the output.)*")
}| Section | Parameter | Estimate | SE | z | p | STDYX | |
|---|---|---|---|---|---|---|---|
| 2 | Thresholds | QPLS5MW1$1 | -0.882 | 0.054 | -16.33 | < .001 | -0.438 |
| 1 | QPLS5MW1.ON | RFDENHW1 | 0.247 | 0.015 | 16.47 | < .001 | 0.434 |
Probit Link function
Using MplusAutomation
Mplus program to estimate probit
<a https://drive.google.com/file/d/1S_wNI0l6ULhOu24uHXciDJTggXJYuW6U/view?usp=drive_link
library(MplusAutomation)
library(knitr)
# ---- Make knitting behave like your Rmd folder ----
if (!is.null(knitr::current_input())) {
knitr::opts_knit$set(root.dir = dirname(knitr::current_input()))
}
# ---- Tell MplusAutomation exactly where Mplus is (edit if needed) ----
if (.Platform$OS.type == "windows") {
options(Mplus_command = "C:/Program Files/Mplus/Mplus.exe")
} else if (Sys.info()[["sysname"]] == "Darwin") {
# Common macOS install path (adjust if different)
options(Mplus_command = "/Applications/Mplus/mplus")
} else {
# Linux/other (adjust if different)
options(Mplus_command = "/usr/local/bin/mplus")
}
stopifnot(!is.null(getOption("Mplus_command")),
nzchar(getOption("Mplus_command")),
file.exists(getOption("Mplus_command")))
# ---- File names ----
inp_file <- "rfdendplus5SimPROBIT.inp"
out_file <- sub("\\.inp$", ".out", inp_file)
stopifnot(file.exists(inp_file))
# ---- Run & read model ----
runModels(target = inp_file, showOutput = FALSE, replaceOutfile = "always")
mod <- readModels(target = out_file)
# ---- Helpers ----
fmt <- function(x, digits = 3) ifelse(is.na(x), "", formatC(x, format = "f", digits = digits))
fmt_p <- function(p) {
ifelse(is.na(p), "",
ifelse(p < .001, "< .001", paste0("= ", sub("^0", ".", formatC(p, format="f", digits=3)))))
}
# ---- Table 1: Model fit ----
S <- if (is.data.frame(mod$summaries)) mod$summaries[1, ] else as.data.frame(mod$summaries)
getv <- function(nm) if (nm %in% names(S)) as.numeric(S[[nm]]) else NA_real_
gets <- function(nm) if (nm %in% names(S)) as.character(S[[nm]]) else ""
chisq <- getv("ChiSqM_Value")
df <- getv("ChiSqM_DF")
pp <- getv("ChiSqM_PValue")
cfi <- getv("CFI")
tli <- getv("TLI")
rmsea <- getv("RMSEA_Estimate")
rm_lb <- getv("RMSEA_90CI_LB")
rm_ub <- getv("RMSEA_90CI_UB")
srmr <- getv("SRMR")
aic <- getv("AIC")
bic <- getv("BIC")
estm <- gets("Estimator")
atype <- gets("AnalysisType")
fit_row <- data.frame(
`χ² (df)` = if (is.na(chisq) | is.na(df)) "" else paste0(fmt(chisq, 2), " (", df, ")"),
`p` = fmt_p(pp),
`CFI` = fmt(cfi, 3),
`TLI` = fmt(tli, 3),
`RMSEA [90% CI]` = if (is.na(rmsea)) "" else paste0(fmt(rmsea, 3), " [", fmt(rm_lb, 3), ", ", fmt(rm_ub, 3), "]"),
`SRMR` = fmt(srmr, 3),
`AIC` = fmt(aic, 2),
`BIC` = fmt(bic, 2),
`Estimator` = estm,
`Analysis` = atype,
check.names = FALSE
)
kable(fit_row, align = "l", caption = "Table 1. Overall model fit.")| χ² (df) | p | CFI | TLI | RMSEA [90% CI] | SRMR | AIC | BIC | Estimator | Analysis |
|---|---|---|---|---|---|---|---|---|---|
| 0.00 (0) | < .001 | 1.000 | 1.000 | 0.000 [0.000, 0.000] | 0.000 | WLSMV | MISSING |
# ---- Table 2: Unstandardized & Standardized (STDYX) estimates ----
# Safely grab parameter frames if present
par_unstd <- try(mod$parameters$unstandardized, silent = TRUE)
if (inherits(par_unstd, "try-error") || is.null(par_unstd)) par_unstd <- data.frame()
par_stdyx <- try(mod$parameters$stdyx.standardized, silent = TRUE)
if (inherits(par_stdyx, "try-error") || is.null(par_stdyx)) par_stdyx <- data.frame()
# Helper accessors for column names that vary across Mplus versions
getcol <- function(df, nm, default = NA_real_) {
if (length(df) == 0) return(rep(default, 0))
if (nm %in% names(df)) return(df[[nm]])
# common alternates
alts <- c(est = c("est", "Estimate"),
se = c("se", "S.E."),
p = c("pval","P-Value","PValue"),
z = c("est_se","Z", "Z-value"))
if (nm %in% names(alts)) {
for (cand in alts[[nm]]) if (cand %in% names(df)) return(df[[cand]])
}
rep(default, nrow(df))
}
# Build unstandardized table
if (nrow(par_unstd)) {
U <- data.frame(
Section = if ("paramHeader" %in% names(par_unstd)) par_unstd$paramHeader else "",
Parameter = if ("param" %in% names(par_unstd)) par_unstd$param else "",
Estimate = as.numeric(getcol(par_unstd, "est")),
SE = as.numeric(getcol(par_unstd, "se")),
z = suppressWarnings(as.numeric(getcol(par_unstd, "z"))),
p = as.numeric(getcol(par_unstd, "p")),
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(is.na(U$z)) && any(is.finite(U$Estimate) & is.finite(U$SE) & U$SE > 0)) {
U$z <- with(U, ifelse(is.finite(Estimate) & is.finite(SE) & SE > 0, Estimate/SE, NA_real_))
}
if (all(is.na(U$p)) && any(is.finite(U$z))) {
U$p <- 2 * pnorm(-abs(U$z))
}
} else {
U <- data.frame()
}
# Build STDYX table and merge
if (nrow(par_stdyx)) {
Sx <- data.frame(
Section = if ("paramHeader" %in% names(par_stdyx)) par_stdyx$paramHeader else "",
Parameter = if ("param" %in% names(par_stdyx)) par_stdyx$param else "",
STDYX = as.numeric(getcol(par_stdyx, "est")),
stringsAsFactors = FALSE, check.names = FALSE
)
} else {
Sx <- data.frame(Section=character(), Parameter=character(), STDYX=numeric(), check.names = FALSE)
}
PE <- if (nrow(U)) merge(U, Sx, by = c("Section","Parameter"), all.x = TRUE) else data.frame()
# Keep common, interpretable sections first
section_order <- c("LOADINGS", "BY", "Latent Variable Loadings", "REGRESSIONS", "WITH",
"Latent Variable Covariances", "Residual Variances", "Variances",
"Means", "Intercepts", "Thresholds")
if (nrow(PE)) {
PE$Section <- as.character(PE$Section)
ord <- match(PE$Section, section_order)
PE <- PE[order(ifelse(is.na(ord), length(section_order)+1, ord), PE$Section, PE$Parameter), ]
# Format numbers for APA
PE_out <- transform(PE,
Estimate = fmt(Estimate, 3),
SE = fmt(SE, 3),
z = ifelse(is.na(z), "", fmt(z, 2)),
p = sapply(p, fmt_p),
STDYX = fmt(STDYX, 3))
kable(PE_out[, c("Section","Parameter","Estimate","SE","z","p","STDYX")],
align = "l", caption = "Table 2. Unstandardized and standardized (STDYX) parameter estimates.")
} else {
cat("*(No parameter estimates were found in the output.)*")
}| Section | Parameter | Estimate | SE | z | p | STDYX | |
|---|---|---|---|---|---|---|---|
| 2 | Thresholds | QPLS5MW1$1 | -0.532 | 0.031 | -17.16 | < .001 | -0.470 |
| 1 | QPLS5MW1.ON | RFDENHW1 | 0.149 | 0.009 | 16.56 | < .001 | 0.467 |
# Logit and Probit Regressions in R + Loess regression (Figure 18.5)
``` r
RFDSim <- read.table("plus5RFDSimR.txt", header = TRUE)
RFDSim_subset <- RFDSim[, c("qpls5mw1", "rfdenhw1")]
RFDSim_subset$rfdenhw1 <- pmin(RFDSim_subset$rfdenhw1, 17) # IF rfdenhw1>=17 THEN rfdenhw1=17;
## ---- Libraries ----
suppressPackageStartupMessages({
library(ggplot2)
library(dplyr)
})
## ---- Link helpers ----
logit <- function(p, eps=1e-12){ p <- pmin(pmax(p,eps), 1-eps); log(p/(1-p)) }
inv_logit <- function(eta) plogis(eta)
probit <- function(p, eps=1e-12){ p <- pmin(pmax(p,eps), 1-eps); qnorm(p) }
inv_probit <- function(eta) pnorm(eta)
## ---- Data prep: use qpls5mw1 (binary) and rfdenhw1 from RFDSim_subset ----
stopifnot(exists("RFDSim_subset"))
dat <- RFDSim_subset[, c("qpls5mw1","rfdenhw1")]
# Coerce outcome to 0/1 robustly
y_raw <- dat$qpls5mw1
if (is.factor(y_raw) || is.character(y_raw)) {
levs <- sort(unique(as.character(y_raw)))
dat$y <- as.integer(as.character(y_raw) == levs[length(levs)])
} else {
dat$y <- as.integer(y_raw)
}
if (!all(dat$y %in% c(0L, 1L), na.rm = TRUE)) {
stop("qpls5mw1 must be binary (0/1) or a two-level factor/character.")
}
## ---- Center rfdenhw1
mu_x <- mean(dat$rfdenhw1, na.rm = TRUE)
dat$rfdenhw1_c <- dat$rfdenhw1 - mu_x
sd_x <- sd(dat$rfdenhw1_c, na.rm = TRUE)
## ---- Fit GLMs with centered X ----
library(sandwich); library(lmtest)
fit_logit <- glm(y ~ rfdenhw1_c, data = dat, family = binomial(link = "logit"), na.action = na.omit)
coeftest(fit_logit, vcov = vcovHC(fit_logit, type = "HC1"))##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.882034 0.054039 16.322 < 2.2e-16 ***
## rfdenhw1_c 0.246614 0.014964 16.480 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit_probit <- glm(y ~ rfdenhw1_c, data = dat, family = binomial(link = "probit"), na.action = na.omit)
coeftest(fit_probit, vcov = vcovHC(fit_probit, type = "HC1"))##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5315811 0.0314875 16.882 < 2.2e-16 ***
## rfdenhw1_c 0.1491653 0.0086824 17.180 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ---- Coefficient table with standardizations and thresholds ----
coef_table <- function(fit, link = c("logit","probit"), sd_x){
link <- match.arg(link)
s <- summary(fit)$coefficients
# Pull intercept and slope rows
b0 <- unname(s["(Intercept)", "Estimate"]); se0 <- unname(s["(Intercept)", "Std. Error"])
z0 <- unname(s["(Intercept)", "z value"]); p0 <- unname(s["(Intercept)", "Pr(>|z|)"])
b1 <- unname(s["rfdenhw1_c", "Estimate"]); se1 <- unname(s["rfdenhw1_c", "Std. Error"])
z1 <- unname(s["rfdenhw1_c", "z value"]); p1 <- unname(s["rfdenhw1_c", "Pr(>|z|)"])
# Threshold (Mplus-style): tau = -Intercept because X is centered
tau <- -b0; se_tau <- se0; z_tau <- -z0; p_tau <- p0
# STD.X: multiply by SD(X); SE scales same way (linear rescaling)
beta_stdX <- b1 * sd_x
se_stdX <- se1 * sd_x
z_stdX <- beta_stdX / se_stdX
p_stdX <- 2*pnorm(-abs(z_stdX))
# STDYX: additionally divide by SD(Y*)
sd_y_star <- if (link == "probit") 1 else (pi / sqrt(3))
beta_stdyx <- b1 * sd_x / sd_y_star
se_stdyx <- se1 * sd_x / sd_y_star
z_stdyx <- beta_stdyx / se_stdyx
p_stdyx <- 2*pnorm(-abs(z_stdyx))
tibble::tibble(
Link = link,
Term = c("(Intercept)", "rfdenhw1_c", "Threshold (= -Intercept)", "rfdenhw1_c [STD.X]", "rfdenhw1_c [STDYX]"),
Estimate = c(b0, b1, tau, beta_stdX, beta_stdyx),
SE = c(se0, se1, se_tau, se_stdX, se_stdyx),
z = c(z0, z1, z_tau, z_stdX, z_stdyx),
p = c(p0, p1, p_tau, p_stdX, p_stdyx),
`SD(Y*) used` = c(NA_real_, NA_real_, NA_real_, NA_real_, sd_y_star)
)
}
tab_logit <- coef_table(fit_logit, "logit", sd_x)
tab_probit <- coef_table(fit_probit, "probit", sd_x)
coef_tab <- dplyr::bind_rows(tab_logit, tab_probit) %>%
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 4)))
cat("\nCoefficients (centered X):\n")##
## Coefficients (centered X):print(coef_tab, row.names = FALSE)## # A tibble: 10 × 7
## Link Term Estimate SE z p `SD(Y*) used`
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 logit (Intercept) 0.882 0.0538 16.4 0 NA
## 2 logit rfdenhw1_c 0.247 0.0153 16.1 0 NA
## 3 logit Threshold (= -Intercept) -0.882 0.0538 -16.4 0 NA
## 4 logit rfdenhw1_c [STD.X] 0.874 0.0542 16.1 0 NA
## 5 logit rfdenhw1_c [STDYX] 0.482 0.0299 16.1 0 1.81
## 6 probit (Intercept) 0.532 0.0314 17.0 0 NA
## 7 probit rfdenhw1_c 0.149 0.0089 16.8 0 NA
## 8 probit Threshold (= -Intercept) -0.532 0.0314 -17.0 0 NA
## 9 probit rfdenhw1_c [STD.X] 0.529 0.0315 16.8 0 NA
## 10 probit rfdenhw1_c [STDYX] 0.529 0.0315 16.8 0 1## ---- Predicted probability curves (on raw X axis) + LOESS ----
rng <- range(dat$rfdenhw1, na.rm = TRUE)
grid <- data.frame(rfdenhw1 = seq(rng[1], rng[2], length.out = 300))
grid$rfdenhw1_c <- grid$rfdenhw1 - mu_x
grid$p_logit <- plogis(predict(fit_logit, newdata = grid, type = "link"))
grid$p_probit <- pnorm( predict(fit_probit, newdata = grid, type = "link"))
# LOESS on raw (x, y); clamp to [0,1]
loess_fit <- loess(y ~ rfdenhw1, data = dat, span = 0.75, degree = 1, na.action = na.omit)
grid$p_loess <- pmin(pmax(predict(loess_fit, newdata = grid), 0), 1)
ggplot() +
geom_point(data = dat, aes(rfdenhw1, y),
alpha = 0.25, position = position_jitter(height = 0.02)) +
geom_line(data = grid, aes(rfdenhw1, p_logit, linetype = "logit"), linewidth = 1) +
geom_line(data = grid, aes(rfdenhw1, p_probit, linetype = "probit"), linewidth = 1) +
scale_linetype_manual(values = c(logit = "solid", probit = "dashed", loess = "dotdash"),
name = "Fit") +
coord_cartesian(ylim = c(0, 1)) +
labs(title = "qpls5mw1 ~ rfdenhw1 (centered): logit and probit",
x = "rfdenhw1 (raw scale)", y = "P(qpls5mw1 = 1)") +
theme_minimal(base_size = 13)
ggplot() +
geom_point(data = dat, aes(rfdenhw1, y),
alpha = 0.25, position = position_jitter(height = 0.02)) +
geom_line(data = grid, aes(rfdenhw1, p_loess, linetype = "loess"), linewidth = 1) +
scale_linetype_manual(values = c(logit = "solid", probit = "dashed", loess = "dotdash"),
name = "Fit") +
coord_cartesian(ylim = c(0, 1)) +
labs(title = "qpls5mw1 ~ rfdenhw1 (centered): LOESS by itself (Figure 18.5)",
x = "rfdenhw1 (raw scale)", y = "P(qpls5mw1 = 1)") +
theme_minimal(base_size = 13)
Rasch Model in R Table of percents (Table 18.6)
RFDSim <- read.table("C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/plus5RFDSimR.txt", header = TRUE)
RFDSim_subset <- RFDSim[, c("qpls5mw0", "qpls5mw1","qpls5mw2","qpls5mw3","qpls5mw4","qpls5mw5","qpls5mw6","qpls5mw7","qpls5mw8")]
prop.table(table(RFDSim_subset$qpls5mw4))##
## 0 1
## 0.4833502 0.5166498library(dplyr)
library(tidyr)
library(gt)
# Variables to summarize
vars <- paste0("qpls5mw", 0:8)
# Basic checks
stopifnot(exists("RFDSim_subset"))
stopifnot(all(vars %in% names(RFDSim_subset)))
# Helper: coerce values to 0/1 safely if they are stored as factors/characters/logical
to01 <- function(x) {
if (is.logical(x)) return(as.numeric(x)) # TRUE/FALSE -> 1/0
if (is.numeric(x)) return(x) # assume 0/1/NA
x_chr <- as.character(x)
dplyr::case_when(
x_chr %in% c("1", "TRUE", "T") ~ 1,
x_chr %in% c("0", "FALSE", "F") ~ 0,
TRUE ~ NA_real_
)
}
# Compute counts and percentages among non-missing
apa_tbl <- RFDSim_subset %>%
select(all_of(vars)) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "value") %>%
mutate(value = to01(value)) %>%
group_by(Variable) %>%
summarise(
n = sum(!is.na(value)),
n0 = sum(value == 0, na.rm = TRUE),
n1 = sum(value == 1, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
pct0 = ifelse(n > 0, n0 / n, NA_real_), # proportions for formatting as %
pct1 = ifelse(n > 0, n1 / n, NA_real_)
) %>%
arrange(Variable)
# Build APA-style table
gt(apa_tbl) %>%
cols_label(
Variable = "Variable",
n = html("<em>n</em>"),
n0 = "0s (n)",
n1 = "1s (n)",
pct0 = "0s (%)",
pct1 = "1s (%)"
) %>%
fmt_percent(columns = c(pct0, pct1), decimals = 1) %>% # show as percentages
tab_header(
title = md("**Table 1**"),
subtitle = md("*Counts and Percentages of 0 and 1 (Among Non-Missing)*")
) %>%
tab_source_note(md("Note. Percentages are among non-missing observations; displayed to one decimal place."))| Table 1 | |||||
| Counts and Percentages of 0 and 1 (Among Non-Missing) | |||||
| Variable | n | 0s (n) | 1s (n) | 0s (%) | 1s (%) |
|---|---|---|---|---|---|
| qpls5mw0 | 2561 | 1497 | 1064 | 58.5% | 41.5% |
| qpls5mw1 | 2568 | 1174 | 1394 | 45.7% | 54.3% |
| qpls5mw2 | 2086 | 975 | 1111 | 46.7% | 53.3% |
| qpls5mw3 | 1904 | 890 | 1014 | 46.7% | 53.3% |
| qpls5mw4 | 1982 | 958 | 1024 | 48.3% | 51.7% |
| qpls5mw5 | 1916 | 816 | 1100 | 42.6% | 57.4% |
| qpls5mw6 | 1882 | 825 | 1057 | 43.8% | 56.2% |
| qpls5mw7 | 1896 | 763 | 1133 | 40.2% | 59.8% |
| qpls5mw8 | 1764 | 744 | 1020 | 42.2% | 57.8% |
| Note. Percentages are among non-missing observations; displayed to one decimal place. | |||||
Rasch Model
RFDSim <- read.table("C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/plus5RFDSimR.txt", header = TRUE)
RFDSim_subset <- RFDSim[, c("qpls5mw0", "qpls5mw1","qpls5mw2","qpls5mw3","qpls5mw4","qpls5mw5","qpls5mw6","qpls5mw7","qpls5mw8")]
library(lavaan)
library(dplyr)
library(tidyr)
library(ggplot2)
# --- Items and basic checks ---
vars <- paste0("qpls5mw", 0:8)
stopifnot(exists("RFDSim_subset"), all(vars %in% names(RFDSim_subset)))
# --- Fit Rasch (normal-ogive; equal discrimination) if not already fit ---
if (!exists("fit_rasch")) {
rasch <- paste0("ability =~ ", paste(paste0("a*", vars), collapse = " + "), "\n")
fit_rasch <- cfa(rasch,
data = RFDSim_subset,
ordered = vars,
estimator = "WLSMV",
std.lv = TRUE)
}
summary(fit_rasch)## lavaan 0.6-19 ended normally after 4 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 18
## Number of equality constraints 8
##
## Used Total
## Number of observations 1231 2568
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 343.267 342.910
## Degrees of freedom 35 35
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.027
## Shift parameter 8.776
## simple second-order correction
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## ability =~
## qpls5mw0 (a) 0.891 0.006 139.966 0.000
## qpls5mw1 (a) 0.891 0.006 139.966 0.000
## qpls5mw2 (a) 0.891 0.006 139.966 0.000
## qpls5mw3 (a) 0.891 0.006 139.966 0.000
## qpls5mw4 (a) 0.891 0.006 139.966 0.000
## qpls5mw5 (a) 0.891 0.006 139.966 0.000
## qpls5mw6 (a) 0.891 0.006 139.966 0.000
## qpls5mw7 (a) 0.891 0.006 139.966 0.000
## qpls5mw8 (a) 0.891 0.006 139.966 0.000
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|)
## qpls5mw0|t1 0.427 0.037 11.558 0.000
## qpls5mw1|t1 0.126 0.036 3.503 0.000
## qpls5mw2|t1 0.040 0.036 1.111 0.267
## qpls5mw3|t1 0.007 0.036 0.199 0.842
## qpls5mw4|t1 0.048 0.036 1.339 0.181
## qpls5mw5|t1 -0.093 0.036 -2.592 0.010
## qpls5mw6|t1 -0.081 0.036 -2.250 0.024
## qpls5mw7|t1 -0.185 0.036 -5.154 0.000
## qpls5mw8|t1 -0.148 0.036 -4.130 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .qpls5mw0 0.207
## .qpls5mw1 0.207
## .qpls5mw2 0.207
## .qpls5mw3 0.207
## .qpls5mw4 0.207
## .qpls5mw5 0.207
## .qpls5mw6 0.207
## .qpls5mw7 0.207
## .qpls5mw8 0.207
## ability 1.000# --- Extract thresholds and common discrimination (a) ---
pe <- parameterEstimates(fit_rasch)
thresh <- pe %>%
dplyr::filter(op == "|") %>%
transmute(Item = lhs, threshold = est) %>%
distinct()
# Prefer the labeled 'a' (from a*), else fall back to the (equal) loadings' mean
a_val <- suppressWarnings(unique(pe$est[pe$label == "a"])[1])
if (is.na(a_val)) {
a_val <- mean(pe$est[pe$op == "=~"])
}
# --- Compute model ICCs: P(Y=1 | θ) = Φ(a*θ − τ_i) ---
theta_grid <- seq(-4, 4, length.out = 401)
icc_df <- thresh %>%
tidyr::crossing(theta = theta_grid) %>%
mutate(p = pnorm(a_val * theta - threshold))
# --- Single plot with all items overlaid (no confidence intervals) ---
ggplot(icc_df, aes(x = theta, y = p, color = Item)) +
geom_line(linewidth = 1) +
labs(
title = "Rasch Item Characteristic Curves (normal-ogive)",
x = expression(theta),
y = "Pr(Y = 1 | θ)",
color = "Item"
) +
coord_cartesian(ylim = c(0, 1)) +
theme_minimal(base_size = 13) +
theme(legend.position = "right")
# Optional: save
# ggsave("Rasch_ICCs_all_items.png", width = 7, height = 5, dpi = 300)Two Parameter Logistic
# --- Data and packages (as in your chunk) ---
RFDSim <- read.table("C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/plus5RFDSimR.txt", header = TRUE)
RFDSim_subset <- RFDSim[, c("qpls5mw0","qpls5mw1","qpls5mw2","qpls5mw3",
"qpls5mw4","qpls5mw5","qpls5mw6","qpls5mw7","qpls5mw8")]
library(lavaan)
library(dplyr)
library(tidyr)
library(ggplot2)
library(gt)
vars <- paste0("qpls5mw", 0:8)
stopifnot(all(vars %in% names(RFDSim_subset)))
# --- 2PL model (free discriminations) in lavaan (normal-ogive / probit) ---
twoPL <- paste0("ability =~ ", paste(vars, collapse = " + "), "\n")
fit_2pl <- cfa(twoPL,
data = RFDSim_subset,
ordered = vars, # treat as categorical
estimator= "WLSMV",
std.lv = TRUE) # Var(ability)=1; all loadings free
summary(fit_2pl, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 17 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 18
##
## Used Total
## Number of observations 1231 2568
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 223.026 352.514
## Degrees of freedom 27 27
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 0.642
## Shift parameter 5.147
## simple second-order correction
##
## Model Test Baseline Model:
##
## Test statistic 44712.779 23356.763
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.916
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.996 0.986
## Tucker-Lewis Index (TLI) 0.994 0.981
##
## Robust Comparative Fit Index (CFI) 0.892
## Robust Tucker-Lewis Index (TLI) 0.857
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.077 0.099
## 90 Percent confidence interval - lower 0.068 0.090
## 90 Percent confidence interval - upper 0.086 0.108
## P-value H_0: RMSEA <= 0.050 0.000 0.000
## P-value H_0: RMSEA >= 0.080 0.299 1.000
##
## Robust RMSEA 0.212
## 90 Percent confidence interval - lower 0.186
## 90 Percent confidence interval - upper 0.239
## P-value H_0: Robust RMSEA <= 0.050 0.000
## P-value H_0: Robust RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.056 0.056
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ability =~
## qpls5mw0 0.807 0.019 41.395 0.000 0.807 0.807
## qpls5mw1 0.854 0.016 54.934 0.000 0.854 0.854
## qpls5mw2 0.905 0.012 77.241 0.000 0.905 0.905
## qpls5mw3 0.935 0.009 100.691 0.000 0.935 0.935
## qpls5mw4 0.918 0.011 85.433 0.000 0.918 0.918
## qpls5mw5 0.904 0.011 78.884 0.000 0.904 0.904
## qpls5mw6 0.882 0.013 66.045 0.000 0.882 0.882
## qpls5mw7 0.882 0.014 64.748 0.000 0.882 0.882
## qpls5mw8 0.870 0.014 60.669 0.000 0.870 0.870
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## qpls5mw0|t1 0.427 0.037 11.558 0.000 0.427 0.427
## qpls5mw1|t1 0.126 0.036 3.503 0.000 0.126 0.126
## qpls5mw2|t1 0.040 0.036 1.111 0.267 0.040 0.040
## qpls5mw3|t1 0.007 0.036 0.199 0.842 0.007 0.007
## qpls5mw4|t1 0.048 0.036 1.339 0.181 0.048 0.048
## qpls5mw5|t1 -0.093 0.036 -2.592 0.010 -0.093 -0.093
## qpls5mw6|t1 -0.081 0.036 -2.250 0.024 -0.081 -0.081
## qpls5mw7|t1 -0.185 0.036 -5.154 0.000 -0.185 -0.185
## qpls5mw8|t1 -0.148 0.036 -4.130 0.000 -0.148 -0.148
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .qpls5mw0 0.349 0.349 0.349
## .qpls5mw1 0.270 0.270 0.270
## .qpls5mw2 0.181 0.181 0.181
## .qpls5mw3 0.126 0.126 0.126
## .qpls5mw4 0.158 0.158 0.158
## .qpls5mw5 0.183 0.183 0.183
## .qpls5mw6 0.223 0.223 0.223
## .qpls5mw7 0.223 0.223 0.223
## .qpls5mw8 0.243 0.243 0.243
## ability 1.000 1.000 1.000# --- Extract item parameters: thresholds (tau_i) & discriminations (a_i, probit metric) ---
pe <- parameterEstimates(fit_2pl)
loadings <- pe %>%
filter(op == "=~") %>%
transmute(Item = rhs, a_probit = est)
thresh <- pe %>%
filter(op == "|") %>%
group_by(Item = lhs) %>%
summarise(threshold = first(est), .groups = "drop") # one threshold per binary item
# --- Compute difficulties (b_i = tau_i / a_i) and logistic-scaled a (a_logit = 1.702*a_probit) ---
D <- 1.702
item_parms <- loadings %>%
inner_join(thresh, by = "Item") %>%
mutate(
b_probit = threshold / a_probit, # difficulty on probit/θ scale
a_logit = D * a_probit, # logistic slope
b_logit = b_probit # same location when using D-scaling
) %>%
arrange(Item)
# --- Table of 2PL item parameters ---
gt(item_parms) %>%
cols_label(
Item = "Item",
a_probit = html("<em>a</em> (probit)"),
threshold = html("<em>τ</em> (threshold)"),
b_probit = html("<em>b</em> (difficulty, probit)"),
a_logit = html("<em>a</em> (logit, D=1.702)"),
b_logit = html("<em>b</em> (difficulty, logit)")
) %>%
fmt_number(columns = c(a_probit, threshold, b_probit, a_logit, b_logit), decimals = 3) %>%
tab_header(
title = md("**Two-Parameter Logistic (2PL) Item Parameters**"),
subtitle = md("Estimated with lavaan (normal-ogive); logistic scaling shown via D=1.702")
) %>%
tab_source_note(md("Note. Difficulty computed as τ/α on the probit scale; logistic slopes use a<sub>logit</sub>=1.702·a<sub>probit</sub>."))| Two-Parameter Logistic (2PL) Item Parameters | |||||
| Estimated with lavaan (normal-ogive); logistic scaling shown via D=1.702 | |||||
| Item | a (probit) | τ (threshold) | b (difficulty, probit) | a (logit, D=1.702) | b (difficulty, logit) |
|---|---|---|---|---|---|
| qpls5mw0 | 0.807 | 0.427 | 0.529 | 1.374 | 0.529 |
| qpls5mw1 | 0.854 | 0.126 | 0.147 | 1.454 | 0.147 |
| qpls5mw2 | 0.905 | 0.040 | 0.044 | 1.541 | 0.044 |
| qpls5mw3 | 0.935 | 0.007 | 0.008 | 1.591 | 0.008 |
| qpls5mw4 | 0.918 | 0.048 | 0.052 | 1.562 | 0.052 |
| qpls5mw5 | 0.904 | −0.093 | −0.103 | 1.538 | −0.103 |
| qpls5mw6 | 0.882 | −0.081 | −0.091 | 1.500 | −0.091 |
| qpls5mw7 | 0.882 | −0.185 | −0.210 | 1.500 | −0.210 |
| qpls5mw8 | 0.870 | −0.148 | −0.170 | 1.481 | −0.170 |
| Note. Difficulty computed as τ/α on the probit scale; logistic slopes use alogit=1.702·aprobit. | |||||
# --- ICCs for ALL items on a single plot (logistic display with D=1.702) ---
theta_grid <- seq(-4, 4, length.out = 401)
icc_df <- item_parms %>%
tidyr::crossing(theta = theta_grid) %>%
mutate(
# Logistic ICC: P(Y=1|θ) = logit^{-1}[ D*(a_probit*θ − τ) ] = logit^{-1}[ a_logit*(θ − b) ]
p = plogis(D * (a_probit * theta - threshold))
)
ggplot(icc_df, aes(x = theta, y = p, color = Item)) +
geom_line(linewidth = 1) +
labs(
title = "2PL Item Characteristic Curves (logistic display, D = 1.702)",
x = expression(theta),
y = "Pr(Y = 1 | θ)",
color = "Item"
) +
coord_cartesian(ylim = c(0, 1)) +
theme_minimal(base_size = 13) +
theme(legend.position = "right")
# Optional: save plot
# ggsave("TwoPL_ICCs_all_items.png", width = 7, height = 5, dpi = 300)Orthogonal FCSI (Latent Basis) Growth Curve Model for Ordered Categorical Data
Using R
Path Diagram
Orthogonal FCSI (Latent Basis) for Ordered
Categories Figure
R Code (See Table 18.8)
library(lavaan)
Plus5Sim <- read.table("C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/qpls5_mcR.txt", header = TRUE)
vars <- paste0("qpls5mw", 0:8)
FCSICat <- '
# Thresholds: first fixed at 0; 2–5 constrained equal across items via labels a2..a5
qpls5mw0 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw1 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw2 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw3 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw4 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw5 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw6 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw7 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw8 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
# Slope factor (s): var=1, mean free (label ms)
s =~ NA*qpls5mw0 + L0*qpls5mw0 + L1*qpls5mw1 + L2*qpls5mw2 + L3*qpls5mw3 +
L4*qpls5mw4 + L5*qpls5mw5 + L6*qpls5mw6 + L7*qpls5mw7 + L8*qpls5mw8
s ~~ 1*s
s ~ ms*1
# Intercept factor (i): loadings fixed 1, var free (vi), mean free (mi)
i =~ 1*qpls5mw0 + 1*qpls5mw1 + 1*qpls5mw2 + 1*qpls5mw3 +
1*qpls5mw4 + 1*qpls5mw5 + 1*qpls5mw6 + 1*qpls5mw7 + 1*qpls5mw8
i ~~ vi*i
i ~ mi*1
# Orthogonality
i ~~ 0*s
#Scaling Factors
qpls5mw1 ~*~ NA*qpls5mw1
qpls5mw2 ~*~ NA*qpls5mw2
qpls5mw3 ~*~ NA*qpls5mw3
qpls5mw4 ~*~ NA*qpls5mw4
qpls5mw5 ~*~ NA*qpls5mw5
qpls5mw6 ~*~ NA*qpls5mw6
qpls5mw7 ~*~ NA*qpls5mw7
qpls5mw8 ~*~ NA*qpls5mw8
'
FCSICatfit <- cfa(
FCSICat, data = Plus5Sim,
ordered = vars,
estimator = "WLSMV",
parameterization = "delta",mimic = "Mplus"
)
summary(FCSICatfit, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)## lavaan 0.6-19 ended normally after 83 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 56
## Number of equality constraints 32
##
## Number of observations 1500
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 37.524 57.965
## Degrees of freedom 57 57
## P-value (Chi-square) 0.978 0.439
## Scaling correction factor 0.821
## Shift parameter 12.236
## simple second-order correction (WLSMV)
##
## Model Test Baseline Model:
##
## Test statistic 3687.684 2720.019
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.361
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.003 1.000
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.003
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.000 0.016
## P-value H_0: RMSEA <= 0.050 1.000 1.000
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.022 0.022
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s =~
## qpls5mw0 (L0) 0.500 0.052 9.655 0.000 0.500 0.500
## qpls5mw1 (L1) 0.305 0.047 6.514 0.000 0.305 0.340
## qpls5mw2 (L2) 0.234 0.044 5.307 0.000 0.234 0.273
## qpls5mw3 (L3) 0.202 0.044 4.603 0.000 0.202 0.243
## qpls5mw4 (L4) 0.203 0.045 4.559 0.000 0.203 0.231
## qpls5mw5 (L5) 0.131 0.046 2.810 0.005 0.131 0.150
## qpls5mw6 (L6) 0.069 0.046 1.503 0.133 0.069 0.083
## qpls5mw7 (L7) -0.032 0.052 -0.617 0.537 -0.032 -0.039
## qpls5mw8 (L8) 0.037 0.048 0.769 0.442 0.037 0.042
## i =~
## qpls5mw0 1.000 0.438 0.438
## qpls5mw1 1.000 0.438 0.488
## qpls5mw2 1.000 0.438 0.511
## qpls5mw3 1.000 0.438 0.527
## qpls5mw4 1.000 0.438 0.499
## qpls5mw5 1.000 0.438 0.502
## qpls5mw6 1.000 0.438 0.529
## qpls5mw7 1.000 0.438 0.529
## qpls5mw8 1.000 0.438 0.498
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s ~~
## i 0.000 0.000 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s (ms) -0.594 0.085 -6.991 0.000 -0.594 -0.594
## i (mi) 0.157 0.030 5.161 0.000 0.358 0.358
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## qpls5m0|1 0.000 0.000 0.000
## qpls5m0|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.199
## qpls5m0|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.468
## qpls5m0|4 (a4) 0.902 0.032 27.816 0.000 0.902 0.902
## qpls5m0|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.321
## qpls5m1|1 0.000 0.000 0.000
## qpls5m1|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.222
## qpls5m1|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.521
## qpls5m1|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.005
## qpls5m1|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.472
## qpls5m2|1 0.000 0.000 0.000
## qpls5m2|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.232
## qpls5m2|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.546
## qpls5m2|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.053
## qpls5m2|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.541
## qpls5m3|1 0.000 0.000 0.000
## qpls5m3|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.239
## qpls5m3|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.563
## qpls5m3|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.086
## qpls5m3|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.590
## qpls5m4|1 0.000 0.000 0.000
## qpls5m4|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.227
## qpls5m4|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.533
## qpls5m4|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.027
## qpls5m4|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.504
## qpls5m5|1 0.000 0.000 0.000
## qpls5m5|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.228
## qpls5m5|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.537
## qpls5m5|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.035
## qpls5m5|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.516
## qpls5m6|1 0.000 0.000 0.000
## qpls5m6|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.240
## qpls5m6|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.565
## qpls5m6|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.090
## qpls5m6|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.596
## qpls5m7|1 0.000 0.000 0.000
## qpls5m7|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.240
## qpls5m7|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.565
## qpls5m7|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.090
## qpls5m7|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.596
## qpls5m8|1 0.000 0.000 0.000
## qpls5m8|2 (a2) 0.199 0.008 23.594 0.000 0.199 0.226
## qpls5m8|3 (a3) 0.468 0.017 26.759 0.000 0.468 0.533
## qpls5m8|4 (a4) 0.902 0.032 27.816 0.000 0.902 1.027
## qpls5m8|5 (a5) 1.321 0.047 28.065 0.000 1.321 1.504
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s 1.000 1.000 1.000
## i (vi) 0.192 0.021 9.034 0.000 1.000 1.000
## .qpls5mw0 0.558 0.558 0.558
## .qpls5mw1 0.520 0.520 0.646
## .qpls5mw2 0.488 0.488 0.665
## .qpls5mw3 0.458 0.458 0.663
## .qpls5mw4 0.538 0.538 0.698
## .qpls5mw5 0.551 0.551 0.725
## .qpls5mw6 0.489 0.489 0.713
## .qpls5mw7 0.492 0.492 0.719
## .qpls5mw8 0.579 0.579 0.750
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## qpls5mw1 1.115 0.047 23.911 0.000 1.115 1.000
## qpls5mw2 1.167 0.050 23.404 0.000 1.167 1.000
## qpls5mw3 1.204 0.052 23.126 0.000 1.204 1.000
## qpls5mw4 1.139 0.048 23.530 0.000 1.139 1.000
## qpls5mw5 1.148 0.052 22.254 0.000 1.148 1.000
## qpls5mw6 1.208 0.053 22.737 0.000 1.208 1.000
## qpls5mw7 1.209 0.055 22.036 0.000 1.209 1.000
## qpls5mw8 1.138 0.052 21.923 0.000 1.138 1.000
##
## R-Square:
## Estimate
## qpls5mw0 0.442
## qpls5mw1 0.354
## qpls5mw2 0.335
## qpls5mw3 0.337
## qpls5mw4 0.302
## qpls5mw5 0.275
## qpls5mw6 0.287
## qpls5mw7 0.281
## qpls5mw8 0.250require(lavaangui)
#plot_lavaan(fit)Plot of Free Curve (See Figure 18.8)
library(dplyr)
library(ggplot2)
pe <- parameterEstimates(FCSICatfit)
loadings <- pe %>%
dplyr::filter(lhs == "s", op == "=~", grepl("^qpls5mw\\d+$", rhs)) %>%
dplyr::mutate(wave = as.integer(sub("qpls5mw", "", rhs))) %>%
dplyr::arrange(wave) %>%
dplyr::select(wave, lambda = est) # <-- use est, not .data$est
# Factor scores
lv <- as.data.frame(lavPredict(FCSICatfit, type = "lv"))
if (!all(c("i","s") %in% names(lv))) colnames(lv)[1:2] <- c("i","s")
# Predicted latent response per wave, averaged across people
mu_hat <- outer(lv$s, loadings$lambda, `*`) + lv$i
curve_from_scores <- tibble::tibble(
wave = loadings$wave,
mean_score_based = as.numeric(colMeans(mu_hat, na.rm = TRUE))
)
# Plot
ggplot(curve_from_scores, aes(x = wave, y = mean_score_based)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_x_continuous(breaks = curve_from_scores$wave) +
labs(
x = "Wave",
y = "Mean score-based prediction (i + λ·s)",
title = "Predicted FCSI Average Trajectory (score-based) See Figure 18.8"
) +
theme_minimal(base_size = 13)
Using MplusAutomation
<a https://drive.google.com/file/d/1uc-D22iyKqYONW61EqsGF_GNlsXqeLn0/view?usp=drive_link
<a https://drive.google.com/file/d/1VcuMP2jR8Lv5O_q4dQCRSSuh9FpQC81a/view?usp=drive_link
library(MplusAutomation)
library(knitr)
# ---- Make knitting behave like your Rmd folder ----
if (!is.null(knitr::current_input())) {
knitr::opts_knit$set(root.dir = dirname(knitr::current_input()))
}
# ---- Tell MplusAutomation exactly where Mplus is (edit if needed) ----
if (.Platform$OS.type == "windows") {
options(Mplus_command = "C:/Program Files/Mplus/Mplus.exe")
} else if (Sys.info()[["sysname"]] == "Darwin") {
# Common macOS install path (adjust if different)
options(Mplus_command = "/Applications/Mplus/mplus")
} else {
# Linux/other (adjust if different)
options(Mplus_command = "/usr/local/bin/mplus")
}
stopifnot(!is.null(getOption("Mplus_command")),
nzchar(getOption("Mplus_command")),
file.exists(getOption("Mplus_command")))
# ---- File names ----
inp_file <- "C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/montecarlorunFCSI.inp"
out_file <- sub("\\.inp$", ".out", inp_file)
stopifnot(file.exists(inp_file))
# ---- Run & read model ----
runModels(target = inp_file, showOutput = FALSE, replaceOutfile = "always")
mod <- readModels(target = out_file)
# ---- Helpers ----
fmt <- function(x, digits = 3) ifelse(is.na(x), "", formatC(x, format = "f", digits = digits))
fmt_p <- function(p) {
ifelse(is.na(p), "",
ifelse(p < .001, "< .001", paste0("= ", sub("^0", ".", formatC(p, format="f", digits=3)))))
}
# ---- Table 1: Model fit ----
S <- if (is.data.frame(mod$summaries)) mod$summaries[1, ] else as.data.frame(mod$summaries)
getv <- function(nm) if (nm %in% names(S)) as.numeric(S[[nm]]) else NA_real_
gets <- function(nm) if (nm %in% names(S)) as.character(S[[nm]]) else ""
chisq <- getv("ChiSqM_Value")
df <- getv("ChiSqM_DF")
pp <- getv("ChiSqM_PValue")
cfi <- getv("CFI")
tli <- getv("TLI")
rmsea <- getv("RMSEA_Estimate")
rm_lb <- getv("RMSEA_90CI_LB")
rm_ub <- getv("RMSEA_90CI_UB")
srmr <- getv("SRMR")
aic <- getv("AIC")
bic <- getv("BIC")
estm <- gets("Estimator")
atype <- gets("AnalysisType")
fit_row <- data.frame(
`χ² (df)` = if (is.na(chisq) | is.na(df)) "" else paste0(fmt(chisq, 2), " (", df, ")"),
`p` = fmt_p(pp),
`CFI` = fmt(cfi, 3),
`TLI` = fmt(tli, 3),
`RMSEA [90% CI]` = if (is.na(rmsea)) "" else paste0(fmt(rmsea, 3), " [", fmt(rm_lb, 3), ", ", fmt(rm_ub, 3), "]"),
`SRMR` = fmt(srmr, 3),
`AIC` = fmt(aic, 2),
`BIC` = fmt(bic, 2),
`Estimator` = estm,
`Analysis` = atype,
check.names = FALSE
)
kable(fit_row, align = "l", caption = "Table 1. Overall model fit.")| χ² (df) | p | CFI | TLI | RMSEA [90% CI] | SRMR | AIC | BIC | Estimator | Analysis |
|---|---|---|---|---|---|---|---|---|---|
| 57.97 (57) | = ..440 | 1.000 | 1.000 | 0.003 [0.000, 0.016] | 0.016 | WLSMV | GENERAL |
# ---- Table 2: Unstandardized & Standardized (STDYX) estimates ----
# Safely grab parameter frames if present
par_unstd <- try(mod$parameters$unstandardized, silent = TRUE)
if (inherits(par_unstd, "try-error") || is.null(par_unstd)) par_unstd <- data.frame()
par_stdyx <- try(mod$parameters$stdyx.standardized, silent = TRUE)
if (inherits(par_stdyx, "try-error") || is.null(par_stdyx)) par_stdyx <- data.frame()
# Helper accessors for column names that vary across Mplus versions
getcol <- function(df, nm, default = NA_real_) {
if (length(df) == 0) return(rep(default, 0))
if (nm %in% names(df)) return(df[[nm]])
# common alternates
alts <- c(est = c("est", "Estimate"),
se = c("se", "S.E."),
p = c("pval","P-Value","PValue"),
z = c("est_se","Z", "Z-value"))
if (nm %in% names(alts)) {
for (cand in alts[[nm]]) if (cand %in% names(df)) return(df[[cand]])
}
rep(default, nrow(df))
}
# Build unstandardized table
if (nrow(par_unstd)) {
U <- data.frame(
Section = if ("paramHeader" %in% names(par_unstd)) par_unstd$paramHeader else "",
Parameter = if ("param" %in% names(par_unstd)) par_unstd$param else "",
Estimate = as.numeric(getcol(par_unstd, "est")),
SE = as.numeric(getcol(par_unstd, "se")),
z = suppressWarnings(as.numeric(getcol(par_unstd, "z"))),
p = as.numeric(getcol(par_unstd, "p")),
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(is.na(U$z)) && any(is.finite(U$Estimate) & is.finite(U$SE) & U$SE > 0)) {
U$z <- with(U, ifelse(is.finite(Estimate) & is.finite(SE) & SE > 0, Estimate/SE, NA_real_))
}
if (all(is.na(U$p)) && any(is.finite(U$z))) {
U$p <- 2 * pnorm(-abs(U$z))
}
} else {
U <- data.frame()
}
# Build STDYX table and merge
if (nrow(par_stdyx)) {
Sx <- data.frame(
Section = if ("paramHeader" %in% names(par_stdyx)) par_stdyx$paramHeader else "",
Parameter = if ("param" %in% names(par_stdyx)) par_stdyx$param else "",
STDYX = as.numeric(getcol(par_stdyx, "est")),
stringsAsFactors = FALSE, check.names = FALSE
)
} else {
Sx <- data.frame(Section=character(), Parameter=character(), STDYX=numeric(), check.names = FALSE)
}
PE <- if (nrow(U)) merge(U, Sx, by = c("Section","Parameter"), all.x = TRUE) else data.frame()
# Keep common, interpretable sections first
section_order <- c("LOADINGS", "BY", "Latent Variable Loadings", "REGRESSIONS", "WITH",
"Latent Variable Covariances", "Residual Variances", "Variances",
"Means", "Intercepts", "Thresholds")
if (nrow(PE)) {
PE$Section <- as.character(PE$Section)
ord <- match(PE$Section, section_order)
PE <- PE[order(ifelse(is.na(ord), length(section_order)+1, ord), PE$Section, PE$Parameter), ]
# Format numbers for APA
PE_out <- transform(PE,
Estimate = fmt(Estimate, 3),
SE = fmt(SE, 3),
z = ifelse(is.na(z), "", fmt(z, 2)),
p = sapply(p, fmt_p),
STDYX = fmt(STDYX, 3))
kable(PE_out[, c("Section","Parameter","Estimate","SE","z","p","STDYX")],
align = "l", caption = "Table 2. Unstandardized and standardized (STDYX) parameter estimates.")
} else {
cat("*(No parameter estimates were found in the output.)*")
}| Section | Parameter | Estimate | SE | z | p | STDYX | |
|---|---|---|---|---|---|---|---|
| 76 | Variances | I | 0.192 | 0.021 | 9.14 | < .001 | 1.000 |
| 77 | Variances | S | 1.000 | 0.000 | 1.000 | ||
| 11 | Means | I | 0.157 | 0.030 | 5.23 | < .001 | 0.358 |
| 12 | Means | S | 0.594 | 0.085 | 6.99 | < .001 | 0.594 |
| 31 | Thresholds | QPLS5MW0$1 | 0.000 | 0.000 | 0.000 | ||
| 32 | Thresholds | QPLS5MW0$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.199 |
| 33 | Thresholds | QPLS5MW0$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.468 |
| 34 | Thresholds | QPLS5MW0$4 | 0.902 | 0.032 | 28.19 | < .001 | 0.902 |
| 35 | Thresholds | QPLS5MW0$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.321 |
| 36 | Thresholds | QPLS5MW1$1 | 0.000 | 0.000 | 0.000 | ||
| 37 | Thresholds | QPLS5MW1$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.222 |
| 38 | Thresholds | QPLS5MW1$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.521 |
| 39 | Thresholds | QPLS5MW1$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.005 |
| 40 | Thresholds | QPLS5MW1$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.472 |
| 41 | Thresholds | QPLS5MW2$1 | 0.000 | 0.000 | 0.000 | ||
| 42 | Thresholds | QPLS5MW2$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.232 |
| 43 | Thresholds | QPLS5MW2$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.546 |
| 44 | Thresholds | QPLS5MW2$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.053 |
| 45 | Thresholds | QPLS5MW2$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.541 |
| 46 | Thresholds | QPLS5MW3$1 | 0.000 | 0.000 | 0.000 | ||
| 47 | Thresholds | QPLS5MW3$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.239 |
| 48 | Thresholds | QPLS5MW3$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.563 |
| 49 | Thresholds | QPLS5MW3$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.086 |
| 50 | Thresholds | QPLS5MW3$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.590 |
| 51 | Thresholds | QPLS5MW4$1 | 0.000 | 0.000 | 0.000 | ||
| 52 | Thresholds | QPLS5MW4$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.227 |
| 53 | Thresholds | QPLS5MW4$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.533 |
| 54 | Thresholds | QPLS5MW4$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.027 |
| 55 | Thresholds | QPLS5MW4$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.504 |
| 56 | Thresholds | QPLS5MW5$1 | 0.000 | 0.000 | 0.000 | ||
| 57 | Thresholds | QPLS5MW5$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.228 |
| 58 | Thresholds | QPLS5MW5$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.537 |
| 59 | Thresholds | QPLS5MW5$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.035 |
| 60 | Thresholds | QPLS5MW5$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.516 |
| 61 | Thresholds | QPLS5MW6$1 | 0.000 | 0.000 | 0.000 | ||
| 62 | Thresholds | QPLS5MW6$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.240 |
| 63 | Thresholds | QPLS5MW6$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.565 |
| 64 | Thresholds | QPLS5MW6$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.090 |
| 65 | Thresholds | QPLS5MW6$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.596 |
| 66 | Thresholds | QPLS5MW7$1 | 0.000 | 0.000 | 0.000 | ||
| 67 | Thresholds | QPLS5MW7$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.240 |
| 68 | Thresholds | QPLS5MW7$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.565 |
| 69 | Thresholds | QPLS5MW7$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.090 |
| 70 | Thresholds | QPLS5MW7$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.596 |
| 71 | Thresholds | QPLS5MW8$1 | 0.000 | 0.000 | 0.000 | ||
| 72 | Thresholds | QPLS5MW8$2 | 0.199 | 0.008 | 24.88 | < .001 | 0.226 |
| 73 | Thresholds | QPLS5MW8$3 | 0.468 | 0.017 | 27.53 | < .001 | 0.533 |
| 74 | Thresholds | QPLS5MW8$4 | 0.902 | 0.032 | 28.19 | < .001 | 1.027 |
| 75 | Thresholds | QPLS5MW8$5 | 1.321 | 0.047 | 28.11 | < .001 | 1.504 |
| 1 | I.BY | QPLS5MW0 | 1.000 | 0.000 | 0.438 | ||
| 2 | I.BY | QPLS5MW1 | 1.000 | 0.000 | 0.488 | ||
| 3 | I.BY | QPLS5MW2 | 1.000 | 0.000 | 0.511 | ||
| 4 | I.BY | QPLS5MW3 | 1.000 | 0.000 | 0.527 | ||
| 5 | I.BY | QPLS5MW4 | 1.000 | 0.000 | 0.499 | ||
| 6 | I.BY | QPLS5MW5 | 1.000 | 0.000 | 0.502 | ||
| 7 | I.BY | QPLS5MW6 | 1.000 | 0.000 | 0.529 | ||
| 8 | I.BY | QPLS5MW7 | 1.000 | 0.000 | 0.529 | ||
| 9 | I.BY | QPLS5MW8 | 1.000 | 0.000 | 0.498 | ||
| 10 | I.WITH | S | 0.000 | 0.000 | 0.000 | ||
| 13 | S.BY | QPLS5MW0 | -0.500 | 0.052 | -9.62 | < .001 | -0.500 |
| 14 | S.BY | QPLS5MW1 | -0.305 | 0.047 | -6.49 | < .001 | -0.340 |
| 15 | S.BY | QPLS5MW2 | -0.234 | 0.044 | -5.32 | < .001 | -0.273 |
| 16 | S.BY | QPLS5MW3 | -0.202 | 0.044 | -4.59 | < .001 | -0.243 |
| 17 | S.BY | QPLS5MW4 | -0.203 | 0.044 | -4.61 | < .001 | -0.231 |
| 18 | S.BY | QPLS5MW5 | -0.131 | 0.046 | -2.85 | = ..004 | -0.150 |
| 19 | S.BY | QPLS5MW6 | -0.069 | 0.046 | -1.50 | = ..134 | -0.083 |
| 20 | S.BY | QPLS5MW7 | 0.032 | 0.052 | 0.62 | = ..538 | 0.039 |
| 21 | S.BY | QPLS5MW8 | -0.037 | 0.048 | -0.77 | = ..441 | -0.042 |
| 22 | Scales | QPLS5MW0 | 1.000 | 0.000 | 1.000 | ||
| 23 | Scales | QPLS5MW1 | 1.115 | 0.047 | 23.72 | < .001 | 1.000 |
| 24 | Scales | QPLS5MW2 | 1.167 | 0.050 | 23.34 | < .001 | 1.000 |
| 25 | Scales | QPLS5MW3 | 1.204 | 0.052 | 23.15 | < .001 | 1.000 |
| 26 | Scales | QPLS5MW4 | 1.139 | 0.048 | 23.73 | < .001 | 1.000 |
| 27 | Scales | QPLS5MW5 | 1.148 | 0.052 | 22.08 | < .001 | 1.000 |
| 28 | Scales | QPLS5MW6 | 1.208 | 0.053 | 22.79 | < .001 | 1.000 |
| 29 | Scales | QPLS5MW7 | 1.209 | 0.055 | 21.98 | < .001 | 1.000 |
| 30 | Scales | QPLS5MW8 | 1.139 | 0.052 | 21.90 | < .001 | 1.000 |
Orthogonal Linear Growth Curve Model for Ordered Categorical Data See Figure 18.8
Using R
Path Diagram
### Fit
Orthogonal Linear Model (See Table 18.9)
library(lavaan)
Plus5Sim <- read.table("C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/qpls5_mcR.txt", header = TRUE)
vars <- paste0("qpls5mw", 0:8)
LinearCat <- '
# Thresholds: first fixed at 0; 2–5 constrained equal across items via labels a2..a5
qpls5mw0 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw1 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw2 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw3 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw4 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw5 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw6 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw7 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
qpls5mw8 | 0*t1 + a2*t2 + a3*t3 + a4*t4 + a5*t5
# Slope factor (s): var=1, mean free (label ms)
s =~ NA*qpls5mw0 + L0*qpls5mw0 + L1*qpls5mw1 + L2*qpls5mw2 + L3*qpls5mw3 +
L4*qpls5mw4 + L5*qpls5mw5 + L6*qpls5mw6 + L7*qpls5mw7 + L8*qpls5mw8
s ~~ 1*s
s ~ ms*1
# Intercept factor (i): loadings fixed 1, var free (vi), mean free (mi)
i =~ 1*qpls5mw0 + 1*qpls5mw1 + 1*qpls5mw2 + 1*qpls5mw3 +
1*qpls5mw4 + 1*qpls5mw5 + 1*qpls5mw6 + 1*qpls5mw7 + 1*qpls5mw8
i ~~ vi*i
i ~ mi*1
# Orthogonality
i ~~ 0*s
#Scaling Factors
qpls5mw1 ~*~ NA*qpls5mw1
qpls5mw2 ~*~ NA*qpls5mw2
qpls5mw3 ~*~ NA*qpls5mw3
qpls5mw4 ~*~ NA*qpls5mw4
qpls5mw5 ~*~ NA*qpls5mw5
qpls5mw6 ~*~ NA*qpls5mw6
qpls5mw7 ~*~ NA*qpls5mw7
qpls5mw8 ~*~ NA*qpls5mw8
#Constraints
# constraints (lavaan uses "==")
L1 == 0.10*L8 + 0.90*L0
L2 == 0.21*L8 + 0.79*L0
L3 == 0.36*L8 + 0.64*L0
L4 == 0.47*L8 + 0.53*L0
L5 == 0.63*L8 + 0.37*L0
L6 == 0.73*L8 + 0.27*L0
L7 == 0.89*L8 + 0.11*L0
'
LinearFit <- cfa(
LinearCat, data = Plus5Sim,
ordered = vars,
estimator = "WLSMV",
parameterization = "delta",mimic = "Mplus"
)
summary(LinearFit, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)## lavaan 0.6-19 ended normally after 65 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 56
## Number of equality constraints 39
##
## Number of observations 1500
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 62.471 80.414
## Degrees of freedom 64 64
## P-value (Chi-square) 0.531 0.081
## Scaling correction factor 0.943
## Shift parameter 14.189
## simple second-order correction (WLSMV)
##
## Model Test Baseline Model:
##
## Test statistic 3687.684 2720.019
## Degrees of freedom 36 36
## P-value 0.000 0.000
## Scaling correction factor 1.361
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 0.994
## Tucker-Lewis Index (TLI) 1.000 0.997
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.013
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.015 0.021
## P-value H_0: RMSEA <= 0.050 1.000 1.000
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.025 0.025
##
## Parameter Estimates:
##
## Parameterization Delta
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s =~
## qpls5mw0 (L0) 0.406 0.039 10.352 0.000 0.406 0.406
## qpls5mw1 (L1) 0.365 0.039 9.454 0.000 0.365 0.369
## qpls5mw2 (L2) 0.320 0.039 8.282 0.000 0.320 0.332
## qpls5mw3 (L3) 0.259 0.040 6.495 0.000 0.259 0.283
## qpls5mw4 (L4) 0.214 0.042 5.153 0.000 0.214 0.229
## qpls5mw5 (L5) 0.149 0.045 3.306 0.001 0.149 0.159
## qpls5mw6 (L6) 0.108 0.048 2.265 0.024 0.108 0.120
## qpls5mw7 (L7) 0.043 0.053 0.812 0.417 0.043 0.046
## qpls5mw8 (L8) -0.002 0.056 -0.039 0.969 -0.002 -0.002
## i =~
## qpls5mw0 1.000 0.466 0.466
## qpls5mw1 1.000 0.466 0.472
## qpls5mw2 1.000 0.466 0.484
## qpls5mw3 1.000 0.466 0.510
## qpls5mw4 1.000 0.466 0.498
## qpls5mw5 1.000 0.466 0.499
## qpls5mw6 1.000 0.466 0.517
## qpls5mw7 1.000 0.466 0.502
## qpls5mw8 1.000 0.466 0.516
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s ~~
## i 0.000 0.000 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s (ms) -0.632 0.102 -6.167 0.000 -0.632 -0.632
## i (mi) 0.184 0.042 4.423 0.000 0.395 0.395
##
## Thresholds:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## qpls5m0|1 0.000 0.000 0.000
## qpls5m0|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.215
## qpls5m0|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.505
## qpls5m0|4 (a4) 0.974 0.028 34.505 0.000 0.974 0.974
## qpls5m0|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.425
## qpls5m1|1 0.000 0.000 0.000
## qpls5m1|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.217
## qpls5m1|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.511
## qpls5m1|4 (a4) 0.974 0.028 34.505 0.000 0.974 0.985
## qpls5m1|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.442
## qpls5m2|1 0.000 0.000 0.000
## qpls5m2|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.223
## qpls5m2|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.524
## qpls5m2|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.011
## qpls5m2|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.480
## qpls5m3|1 0.000 0.000 0.000
## qpls5m3|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.235
## qpls5m3|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.552
## qpls5m3|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.064
## qpls5m3|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.557
## qpls5m4|1 0.000 0.000 0.000
## qpls5m4|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.229
## qpls5m4|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.539
## qpls5m4|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.040
## qpls5m4|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.522
## qpls5m5|1 0.000 0.000 0.000
## qpls5m5|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.230
## qpls5m5|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.540
## qpls5m5|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.042
## qpls5m5|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.525
## qpls5m6|1 0.000 0.000 0.000
## qpls5m6|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.238
## qpls5m6|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.559
## qpls5m6|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.079
## qpls5m6|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.579
## qpls5m7|1 0.000 0.000 0.000
## qpls5m7|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.231
## qpls5m7|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.543
## qpls5m7|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.048
## qpls5m7|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.533
## qpls5m8|1 0.000 0.000 0.000
## qpls5m8|2 (a2) 0.215 0.008 27.003 0.000 0.215 0.238
## qpls5m8|3 (a3) 0.505 0.016 32.212 0.000 0.505 0.559
## qpls5m8|4 (a4) 0.974 0.028 34.505 0.000 0.974 1.078
## qpls5m8|5 (a5) 1.425 0.041 34.686 0.000 1.425 1.577
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## s 1.000 1.000 1.000
## i (vi) 0.218 0.024 8.995 0.000 1.000 1.000
## .qpls5mw0 0.618 0.618 0.618
## .qpls5mw1 0.626 0.626 0.641
## .qpls5mw2 0.608 0.608 0.655
## .qpls5mw3 0.553 0.553 0.660
## .qpls5mw4 0.613 0.613 0.699
## .qpls5mw5 0.633 0.633 0.725
## .qpls5mw6 0.585 0.585 0.719
## .qpls5mw7 0.644 0.644 0.746
## .qpls5mw8 0.599 0.599 0.734
##
## Scales y*:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## qpls5mw1 1.012 0.033 30.592 0.000 1.012 1.000
## qpls5mw2 1.038 0.035 29.607 0.000 1.038 1.000
## qpls5mw3 1.092 0.037 29.556 0.000 1.092 1.000
## qpls5mw4 1.068 0.037 28.674 0.000 1.068 1.000
## qpls5mw5 1.070 0.041 26.272 0.000 1.070 1.000
## qpls5mw6 1.108 0.041 26.831 0.000 1.108 1.000
## qpls5mw7 1.076 0.043 25.141 0.000 1.076 1.000
## qpls5mw8 1.107 0.044 24.989 0.000 1.107 1.000
##
## R-Square:
## Estimate
## qpls5mw0 0.382
## qpls5mw1 0.359
## qpls5mw2 0.345
## qpls5mw3 0.340
## qpls5mw4 0.301
## qpls5mw5 0.275
## qpls5mw6 0.281
## qpls5mw7 0.254
## qpls5mw8 0.266
##
## Constraints:
## |Slack|
## L1 - (0.10*L8+0.90*L0) 0.000
## L2 - (0.21*L8+0.79*L0) 0.000
## L3 - (0.36*L8+0.64*L0) 0.000
## L4 - (0.47*L8+0.53*L0) 0.000
## L5 - (0.63*L8+0.37*L0) 0.000
## L6 - (0.73*L8+0.27*L0) 0.000
## L7 - (0.89*L8+0.11*L0) 0.000require(lavaangui)
#plot_lavaan(LinearFit)Plot Orthogonal Linear Model (See Figure 18.8)
library(dplyr)
library(ggplot2)
pe <- parameterEstimates(LinearFit)
loadings <- pe %>%
dplyr::filter(lhs == "s", op == "=~", grepl("^qpls5mw\\d+$", rhs)) %>%
dplyr::mutate(wave = as.integer(sub("qpls5mw", "", rhs))) %>%
dplyr::arrange(wave) %>%
dplyr::select(wave, lambda = est) # <-- use est, not .data$est
# Factor scores
lv <- as.data.frame(lavPredict(LinearFit, type = "lv"))
if (!all(c("i","s") %in% names(lv))) colnames(lv)[1:2] <- c("i","s")
# Predicted latent response per wave, averaged across people
mu_hat <- outer(lv$s, loadings$lambda, `*`) + lv$i
curve_from_scores <- tibble::tibble(
wave = loadings$wave,
mean_score_based = as.numeric(colMeans(mu_hat, na.rm = TRUE))
)
# Plot
ggplot(curve_from_scores, aes(x = wave, y = mean_score_based)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_x_continuous(breaks = curve_from_scores$wave) +
labs(
x = "Wave",
y = "Mean score-based prediction (i + λ·s)",
title = "Linear Predicted Average Trajectory (score-based)"
) +
theme_minimal(base_size = 13)
Using MplusAutomation
<a https://drive.google.com/file/d/1dnzD7WmvYD9FYV-saq9HHGGtax8Q4bwW/view?usp=drive_link
<a https://drive.google.com/file/d/14yGVcfndep90cTY9yqEWPGuQWI44njz2/view?usp=drive_link
Executing Program and Fit Statistics
library(MplusAutomation)
library(knitr)
# ---- Make knitting behave like your Rmd folder ----
if (!is.null(knitr::current_input())) {
knitr::opts_knit$set(root.dir = dirname(knitr::current_input()))
}
# ---- Tell MplusAutomation exactly where Mplus is (edit if needed) ----
if (.Platform$OS.type == "windows") {
options(Mplus_command = "C:/Program Files/Mplus/Mplus.exe")
} else if (Sys.info()[["sysname"]] == "Darwin") {
# Common macOS install path (adjust if different)
options(Mplus_command = "/Applications/Mplus/mplus")
} else {
# Linux/other (adjust if different)
options(Mplus_command = "/usr/local/bin/mplus")
}
stopifnot(!is.null(getOption("Mplus_command")),
nzchar(getOption("Mplus_command")),
file.exists(getOption("Mplus_command")))
# ---- File names ----
inp_file <- "C:/Users/woodph/OneDrive - University of Missouri/Documents/SEMR/Ch18/montecarlorunlinear.inp"
out_file <- sub("\\.inp$", ".out", inp_file)
stopifnot(file.exists(inp_file))
# ---- Run & read model ----
runModels(target = inp_file, showOutput = FALSE, replaceOutfile = "always")
mod <- readModels(target = out_file)
# ---- Helpers ----
fmt <- function(x, digits = 3) ifelse(is.na(x), "", formatC(x, format = "f", digits = digits))
fmt_p <- function(p) {
ifelse(is.na(p), "",
ifelse(p < .001, "< .001", paste0("= ", sub("^0", ".", formatC(p, format="f", digits=3)))))
}
# ---- Table 1: Model fit ----
S <- if (is.data.frame(mod$summaries)) mod$summaries[1, ] else as.data.frame(mod$summaries)
getv <- function(nm) if (nm %in% names(S)) as.numeric(S[[nm]]) else NA_real_
gets <- function(nm) if (nm %in% names(S)) as.character(S[[nm]]) else ""
chisq <- getv("ChiSqM_Value")
df <- getv("ChiSqM_DF")
pp <- getv("ChiSqM_PValue")
cfi <- getv("CFI")
tli <- getv("TLI")
rmsea <- getv("RMSEA_Estimate")
rm_lb <- getv("RMSEA_90CI_LB")
rm_ub <- getv("RMSEA_90CI_UB")
srmr <- getv("SRMR")
aic <- getv("AIC")
bic <- getv("BIC")
estm <- gets("Estimator")
atype <- gets("AnalysisType")
fit_row <- data.frame(
`χ² (df)` = if (is.na(chisq) | is.na(df)) "" else paste0(fmt(chisq, 2), " (", df, ")"),
`p` = fmt_p(pp),
`CFI` = fmt(cfi, 3),
`TLI` = fmt(tli, 3),
`RMSEA [90% CI]` = if (is.na(rmsea)) "" else paste0(fmt(rmsea, 3), " [", fmt(rm_lb, 3), ", ", fmt(rm_ub, 3), "]"),
`SRMR` = fmt(srmr, 3),
`AIC` = fmt(aic, 2),
`BIC` = fmt(bic, 2),
`Estimator` = estm,
`Analysis` = atype,
check.names = FALSE
)
kable(fit_row, align = "l", caption = "Table 1. Overall model fit.")| χ² (df) | p | CFI | TLI | RMSEA [90% CI] | SRMR | AIC | BIC | Estimator | Analysis |
|---|---|---|---|---|---|---|---|---|---|
| 80.41 (64) | = ..081 | 0.994 | 0.997 | 0.013 [0.000, 0.021] | 0.019 | WLSMV | GENERAL |
Parameter Estimates
# ---- Table 2: Unstandardized & Standardized (STDYX) estimates ----
# Safely grab parameter frames if present
par_unstd <- try(mod$parameters$unstandardized, silent = TRUE)
if (inherits(par_unstd, "try-error") || is.null(par_unstd)) par_unstd <- data.frame()
par_stdyx <- try(mod$parameters$stdyx.standardized, silent = TRUE)
if (inherits(par_stdyx, "try-error") || is.null(par_stdyx)) par_stdyx <- data.frame()
# Helper accessors for column names that vary across Mplus versions
getcol <- function(df, nm, default = NA_real_) {
if (length(df) == 0) return(rep(default, 0))
if (nm %in% names(df)) return(df[[nm]])
# common alternates
alts <- c(est = c("est", "Estimate"),
se = c("se", "S.E."),
p = c("pval","P-Value","PValue"),
z = c("est_se","Z", "Z-value"))
if (nm %in% names(alts)) {
for (cand in alts[[nm]]) if (cand %in% names(df)) return(df[[cand]])
}
rep(default, nrow(df))
}
# Build unstandardized table
if (nrow(par_unstd)) {
U <- data.frame(
Section = if ("paramHeader" %in% names(par_unstd)) par_unstd$paramHeader else "",
Parameter = if ("param" %in% names(par_unstd)) par_unstd$param else "",
Estimate = as.numeric(getcol(par_unstd, "est")),
SE = as.numeric(getcol(par_unstd, "se")),
z = suppressWarnings(as.numeric(getcol(par_unstd, "z"))),
p = as.numeric(getcol(par_unstd, "p")),
stringsAsFactors = FALSE, check.names = FALSE
)
if (all(is.na(U$z)) && any(is.finite(U$Estimate) & is.finite(U$SE) & U$SE > 0)) {
U$z <- with(U, ifelse(is.finite(Estimate) & is.finite(SE) & SE > 0, Estimate/SE, NA_real_))
}
if (all(is.na(U$p)) && any(is.finite(U$z))) {
U$p <- 2 * pnorm(-abs(U$z))
}
} else {
U <- data.frame()
}
# Build STDYX table and merge
if (nrow(par_stdyx)) {
Sx <- data.frame(
Section = if ("paramHeader" %in% names(par_stdyx)) par_stdyx$paramHeader else "",
Parameter = if ("param" %in% names(par_stdyx)) par_stdyx$param else "",
STDYX = as.numeric(getcol(par_stdyx, "est")),
stringsAsFactors = FALSE, check.names = FALSE
)
} else {
Sx <- data.frame(Section=character(), Parameter=character(), STDYX=numeric(), check.names = FALSE)
}
PE <- if (nrow(U)) merge(U, Sx, by = c("Section","Parameter"), all.x = TRUE) else data.frame()
# Keep common, interpretable sections first
section_order <- c("LOADINGS", "BY", "Latent Variable Loadings", "REGRESSIONS", "WITH",
"Latent Variable Covariances", "Residual Variances", "Variances",
"Means", "Intercepts", "Thresholds")
if (nrow(PE)) {
PE$Section <- as.character(PE$Section)
ord <- match(PE$Section, section_order)
PE <- PE[order(ifelse(is.na(ord), length(section_order)+1, ord), PE$Section, PE$Parameter), ]
# Format numbers for APA
PE_out <- transform(PE,
Estimate = fmt(Estimate, 3),
SE = fmt(SE, 3),
z = ifelse(is.na(z), "", fmt(z, 2)),
p = sapply(p, fmt_p),
STDYX = fmt(STDYX, 3))
kable(PE_out[, c("Section","Parameter","Estimate","SE","z","p","STDYX")],
align = "l", caption = "Unstandardized and standardized (STDYX) parameter estimates Orthogonal Linear.")
} else {
cat("*(No parameter estimates were found in the output.)*")
}| Section | Parameter | Estimate | SE | z | p | STDYX | |
|---|---|---|---|---|---|---|---|
| 76 | Variances | I | 0.218 | 0.024 | 9.08 | < .001 | 1.000 |
| 77 | Variances | S | 1.000 | 0.000 | 1.000 | ||
| 11 | Means | I | 0.184 | 0.042 | 4.38 | < .001 | 0.395 |
| 12 | Means | S | 0.632 | 0.102 | 6.20 | < .001 | 0.632 |
| 31 | Thresholds | QPLS5MW0$1 | 0.000 | 0.000 | 0.000 | ||
| 32 | Thresholds | QPLS5MW0$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.215 |
| 33 | Thresholds | QPLS5MW0$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.505 |
| 34 | Thresholds | QPLS5MW0$4 | 0.974 | 0.028 | 34.79 | < .001 | 0.974 |
| 35 | Thresholds | QPLS5MW0$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.425 |
| 36 | Thresholds | QPLS5MW1$1 | 0.000 | 0.000 | 0.000 | ||
| 37 | Thresholds | QPLS5MW1$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.217 |
| 38 | Thresholds | QPLS5MW1$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.511 |
| 39 | Thresholds | QPLS5MW1$4 | 0.974 | 0.028 | 34.79 | < .001 | 0.985 |
| 40 | Thresholds | QPLS5MW1$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.442 |
| 41 | Thresholds | QPLS5MW2$1 | 0.000 | 0.000 | 0.000 | ||
| 42 | Thresholds | QPLS5MW2$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.223 |
| 43 | Thresholds | QPLS5MW2$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.524 |
| 44 | Thresholds | QPLS5MW2$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.011 |
| 45 | Thresholds | QPLS5MW2$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.480 |
| 46 | Thresholds | QPLS5MW3$1 | 0.000 | 0.000 | 0.000 | ||
| 47 | Thresholds | QPLS5MW3$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.235 |
| 48 | Thresholds | QPLS5MW3$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.552 |
| 49 | Thresholds | QPLS5MW3$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.064 |
| 50 | Thresholds | QPLS5MW3$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.557 |
| 51 | Thresholds | QPLS5MW4$1 | 0.000 | 0.000 | 0.000 | ||
| 52 | Thresholds | QPLS5MW4$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.229 |
| 53 | Thresholds | QPLS5MW4$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.539 |
| 54 | Thresholds | QPLS5MW4$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.040 |
| 55 | Thresholds | QPLS5MW4$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.522 |
| 56 | Thresholds | QPLS5MW5$1 | 0.000 | 0.000 | 0.000 | ||
| 57 | Thresholds | QPLS5MW5$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.230 |
| 58 | Thresholds | QPLS5MW5$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.540 |
| 59 | Thresholds | QPLS5MW5$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.042 |
| 60 | Thresholds | QPLS5MW5$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.525 |
| 61 | Thresholds | QPLS5MW6$1 | 0.000 | 0.000 | 0.000 | ||
| 62 | Thresholds | QPLS5MW6$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.238 |
| 63 | Thresholds | QPLS5MW6$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.559 |
| 64 | Thresholds | QPLS5MW6$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.079 |
| 65 | Thresholds | QPLS5MW6$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.579 |
| 66 | Thresholds | QPLS5MW7$1 | 0.000 | 0.000 | 0.000 | ||
| 67 | Thresholds | QPLS5MW7$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.231 |
| 68 | Thresholds | QPLS5MW7$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.543 |
| 69 | Thresholds | QPLS5MW7$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.048 |
| 70 | Thresholds | QPLS5MW7$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.533 |
| 71 | Thresholds | QPLS5MW8$1 | 0.000 | 0.000 | 0.000 | ||
| 72 | Thresholds | QPLS5MW8$2 | 0.215 | 0.008 | 26.88 | < .001 | 0.238 |
| 73 | Thresholds | QPLS5MW8$3 | 0.505 | 0.016 | 31.56 | < .001 | 0.559 |
| 74 | Thresholds | QPLS5MW8$4 | 0.974 | 0.028 | 34.79 | < .001 | 1.078 |
| 75 | Thresholds | QPLS5MW8$5 | 1.425 | 0.041 | 34.76 | < .001 | 1.577 |
| 1 | I.BY | QPLS5MW0 | 1.000 | 0.000 | 0.466 | ||
| 2 | I.BY | QPLS5MW1 | 1.000 | 0.000 | 0.472 | ||
| 3 | I.BY | QPLS5MW2 | 1.000 | 0.000 | 0.484 | ||
| 4 | I.BY | QPLS5MW3 | 1.000 | 0.000 | 0.510 | ||
| 5 | I.BY | QPLS5MW4 | 1.000 | 0.000 | 0.498 | ||
| 6 | I.BY | QPLS5MW5 | 1.000 | 0.000 | 0.499 | ||
| 7 | I.BY | QPLS5MW6 | 1.000 | 0.000 | 0.517 | ||
| 8 | I.BY | QPLS5MW7 | 1.000 | 0.000 | 0.502 | ||
| 9 | I.BY | QPLS5MW8 | 1.000 | 0.000 | 0.516 | ||
| 10 | I.WITH | S | 0.000 | 0.000 | 0.000 | ||
| 13 | S.BY | QPLS5MW0 | -0.406 | 0.039 | -10.41 | < .001 | -0.406 |
| 14 | S.BY | QPLS5MW1 | -0.365 | 0.039 | -9.36 | < .001 | -0.369 |
| 15 | S.BY | QPLS5MW2 | -0.320 | 0.039 | -8.21 | < .001 | -0.332 |
| 16 | S.BY | QPLS5MW3 | -0.259 | 0.040 | -6.47 | < .001 | -0.283 |
| 17 | S.BY | QPLS5MW4 | -0.214 | 0.042 | -5.10 | < .001 | -0.229 |
| 18 | S.BY | QPLS5MW5 | -0.149 | 0.045 | -3.31 | < .001 | -0.159 |
| 19 | S.BY | QPLS5MW6 | -0.108 | 0.048 | -2.25 | = ..024 | -0.120 |
| 20 | S.BY | QPLS5MW7 | -0.043 | 0.053 | -0.81 | = ..417 | -0.046 |
| 21 | S.BY | QPLS5MW8 | 0.002 | 0.056 | 0.04 | = ..972 | 0.002 |
| 22 | Scales | QPLS5MW0 | 1.000 | 0.000 | 1.000 | ||
| 23 | Scales | QPLS5MW1 | 1.012 | 0.033 | 30.67 | < .001 | 1.000 |
| 24 | Scales | QPLS5MW2 | 1.038 | 0.035 | 29.66 | < .001 | 1.000 |
| 25 | Scales | QPLS5MW3 | 1.092 | 0.037 | 29.51 | < .001 | 1.000 |
| 26 | Scales | QPLS5MW4 | 1.068 | 0.037 | 28.86 | < .001 | 1.000 |
| 27 | Scales | QPLS5MW5 | 1.070 | 0.041 | 26.10 | < .001 | 1.000 |
| 28 | Scales | QPLS5MW6 | 1.108 | 0.041 | 27.02 | < .001 | 1.000 |
| 29 | Scales | QPLS5MW7 | 1.076 | 0.043 | 25.02 | < .001 | 1.000 |
| 30 | Scales | QPLS5MW8 | 1.107 | 0.044 | 25.16 | < .001 | 1.000 |