# principles and practice of sem, 5th ed.
# rex b. kline, guilford press, 2023
# chapter 8, table 8.1, analysis 4
# bias-corrected bootstapped 95% confidence intervals for
# indirect and other effects in the roth et al. (1989)
# parametric path model of illness
# to avoid the problem that some R packages share the same name
# for different functions, all functions are specified next as
# package::function
# which prevents masking, or the default hiding of
# a function with a redundant name from a package used next
date()
v <- R.Version()
print(paste0(v$language, " version ", v$major, ".",
v$minor, " (", v$year, "-", v$month, "-", v$day, ")"))
library(semTools)
library(lavaan)
library(bmem)
library(sem)
# packages semTools and lavaan are used to generate raw
# scores for analyses in packages bmem and sem
# get citation information
citation("semTools", auto = TRUE)
citation("lavaan", auto = TRUE)
citation("bmem", auto = TRUE)
citation("sem", auto = TRUE)
# generate raw scores with descriptive statistics that
# exactly match those from the summarey data in table 4.3
# reuses syntax from analysis 1 ("roth-generate-scores.r") to
# avoid reading in an external raw data file
# functions from semTools and lavaan
roth.cor <- matrix(c(1.00, -.03, .39, -.05, -.08,
-.03, 1.00, .07, -.23, -.16,
.39, .07, 1.00, -.13, -.29,
-.05, -.23, -.13, 1.00, .34,
-.08, -.16, -.29, .34, 1.00), ncol = 5, nrow = 5)
roth.sd <- c(66.5, 38.0, 18.4, 33.5, 62.48)
roth.cov <- diag(roth.sd) %*% roth.cor %*% diag(roth.sd)
c = 372/373
roth.cov = c * roth.cov
roth.data <- semTools::kd(roth.cov, 373, type = "exact")
names(roth.data) <- c("exercise", "hardy", "fitness", "stress", "illness")
roth.data <- transform(roth.data, exercise = exercise + 40.9)
roth.data <- transform(roth.data, fitness = fitness + 67.1)
roth.data <- transform(roth.data, stress = stress + 24.0)
roth.data <- transform(roth.data, illness = illness + 71.67)
# specify path model in package sem
# exercise and hardy covary
# direct effects of exercise, hardy on illness are
# constrained to zero
# disturbance variances are free parameters
# path coefficients are labeled
roth.model <- sem::specifyEquations(covs = "exercise, hardy",
endog.variances = TRUE,
text = "
fitness = a*exercise
stress = c*hardy
illness = b*fitness + d*stress ")
# define indirect effects and copy to vector
indirect <- c("a*b", "c*d")
# 95% bias-corrected confidence intervals for
# indirect effects and all other model parameters
# estimated in bmem package
# defaults are the .95 confidence level and 1,000 generated
# samples
# intercepts are not estimated (the model has a covariance
# structure only)
# automatically generates a message about the bollen-stine
# boostrap method, which is described in chapter 9
# there is no option in bmem to specify a seed,
# so results will probably differ from those reported
# in table 8.5 if this syntax is rerun
bmem::bmem(roth.data, roth.model, indirect, ci = "bc",
moment = FALSE)