Load necessary packages
library(lavaan)
library(stargazer)
Generate data - four factors, each with 3 indicators
# factor X
X <- rnorm(300, 0, 1)
x1 <- X + rnorm(300, 0, 1)
x2 <- X + rnorm(300, 0, 1)
x3 <- X + rnorm(300, 0, 1)
# factor W
W <- rnorm(300, 0, 1)
w1 <- W + rnorm(300, 0, 1)
w2 <- W + rnorm(300, 0, 1)
w3 <- W + rnorm(300, 0, 1)
# factor Z
Z <- rnorm(300, 0, 1)
z1 <- Z + rnorm(300, 0, 1)
z2 <- Z + rnorm(300, 0, 1)
z3 <- Z + rnorm(300, 0, 1)
# factor Y
Y <- rnorm(300, 0, 1)
y1 <- Y + X + W + Z + rnorm(300, 0, 1)
y2 <- Y + X + W + Z + rnorm(300, 0, 1)
y3 <- Y + X + W + Z + rnorm(300, 0, 1)
# compile into data frame
data <- data.frame(x1, x2, x3, w1, w2, w3, z1, z2, z3, y1, y2, y3)
Hide upper portion of correlation matrix. Specifying ‘as.numeric(“”)’ prevents an issue with the ‘initial.zero’ argument below.
cor <- cor(data)
cor[upper.tri(cor, diag=TRUE)] <- as.numeric("")
stargazer(cor[,1:11],
summary = FALSE,
type = 'html',
initial.zero = FALSE,
digits=2)
| x1 | x2 | x3 | w1 | w2 | w3 | z1 | z2 | z3 | y1 | y2 | |
| x1 | |||||||||||
| x2 | .50 | ||||||||||
| x3 | .48 | .49 | |||||||||
| w1 | .01 | -.02 | .01 | ||||||||
| w2 | -.03 | -.04 | -.07 | .44 | |||||||
| w3 | .02 | -.03 | -.001 | .45 | .50 | ||||||
| z1 | -.08 | .02 | -.01 | .03 | .05 | .02 | |||||
| z2 | -.03 | -.10 | .08 | .05 | -.04 | -.01 | .45 | ||||
| z3 | .02 | -.03 | .08 | .12 | .09 | .10 | .49 | .49 | |||
| y1 | .29 | .38 | .34 | .35 | .35 | .34 | .22 | .25 | .31 | ||
| y2 | .30 | .34 | .36 | .35 | .32 | .32 | .23 | .24 | .33 | .80 | |
| y3 | .31 | .35 | .32 | .35 | .33 | .33 | .28 | .25 | .33 | .82 | .79 |
Specify and estimate model
m1.syntax <-'
# exogenous factors
x =~ x1 + x2 + x3
w =~ w1 + w2 + w3
z =~ z1 + z2 + z3
# endogenous factor
y =~ y1 + y2 + y3
# structural model
y ~ x + w + z
'
m1.fit <- sem(m1.syntax, std.lv=TRUE, data=data)
Extract parameter estimates
pars.factors <- standardizedSolution(m1.fit)[ standardizedSolution(m1.fit)[,'op']=='=~', c(1:5)]
pars.regressions <- standardizedSolution(m1.fit)[ standardizedSolution(m1.fit)[,'op']=='~', c(3,4,5,7)]
Use ‘stargazer()’ to save the table as an html file. Once saved, you can open the file using Word and make some minor edits (e.g., labels, column numbers, etc.).
stargazer(pars.factors, summary=FALSE, type='html', rownames=FALSE, initial.zero=FALSE, digits=3, title='Factor loadings')
| lhs | op | rhs | est.std | se |
| x | =~ | x1 | .669 | .043 |
| x | =~ | x2 | .755 | .039 |
| x | =~ | x3 | .669 | .043 |
| w | =~ | w1 | .645 | .045 |
| w | =~ | w2 | .705 | .043 |
| w | =~ | w3 | .686 | .043 |
| z | =~ | z1 | .659 | .046 |
| z | =~ | z2 | .686 | .045 |
| z | =~ | z3 | .728 | .043 |
| y | =~ | y1 | .911 | .014 |
| y | =~ | y2 | .882 | .016 |
| y | =~ | y3 | .900 | .014 |
stargazer(pars.regressions, summary=FALSE, type='html', rownames=FALSE, initial.zero=FALSE, digits=3, title='Predicting Y')
| rhs | est.std | se | pvalue |
| x | .555 | .049 | 0 |
| w | .536 | .050 | 0 |
| z | .397 | .051 | 0 |
And that’s it!