1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/ordMoments.R
\name{lrv2ord}
\alias{lrv2ord}
\title{Calculate Population Moments for Ordinal Data Treated as Numeric}
\usage{
lrv2ord(Sigma, Mu, thresholds, cWts)
}
\arguments{
\item{Sigma}{Population covariance \code{\link[=matrix]{matrix()}}, with variable names
saved in the \code{\link[=dimnames]{dimnames()}} attribute.}
\item{Mu}{Optional \code{numeric} vector of population means. If missing,
all means will be set to zero.}
\item{thresholds}{Either a single \code{numeric} vector of population
thresholds used to discretize each normally distributed variable, or a
named \code{list} of each discretized variable's vector of thresholds.
The discretized variables may be a subset of all variables in \code{Sigma}
if the remaining variables are intended to be observed rather than latent
normally distributed variables.}
\item{cWts}{Optional (default when missing is to use 0 for the lowest
category, followed by successive integers for each higher category).
Either a single \code{numeric} vector of category weights (if they are
identical across all variables) or a named \code{list} of each
discretized variable's vector of category weights.}
}
\value{
A \code{list} including the LRV-scale population moments (means,
covariance matrix, correlation matrix, and thresholds), the category
weights, a \code{data.frame} of implied univariate moments (means,
\emph{SD}s, skewness, and excess kurtosis (i.e., in excess of 3, which is
the kurtosis of the normal distribution) for discretized data treated as
\code{numeric}, and the implied covariance and correlation matrix of
discretized data treated as \code{numeric}.
}
\description{
This function calculates ordinal-scale moments implied by LRV-scale moments
}
\details{
Binary and ordinal data are frequently accommodated in SEM by incorporating
a threshold model that links each observed categorical response variable to
a corresponding latent response variable that is typically assumed to be
normally distributed (Kamata & Bauer, 2008; Wirth & Edwards, 2007).
This function can be useful for real-data analysis or for designing
Monte Carlo simulations, as described by Jorgensen and Johnson (2022).
}
\examples{
## SCENARIO 1: DIRECTLY SPECIFY POPULATION PARAMETERS
## specify population model in LISREL matrices
Nu <- rep(0, 4)
Alpha <- c(1, -0.5)
Lambda <- matrix(c(1, 1, 0, 0, 0, 0, 1, 1), nrow = 4, ncol = 2,
dimnames = list(paste0("y", 1:4), paste0("eta", 1:2)))
Psi <- diag(c(1, .75))
Theta <- diag(4)
Beta <- matrix(c(0, .5, 0, 0), nrow = 2, ncol = 2)
## calculate model-implied population means and covariance matrix
## of latent response variables (LRVs)
IB <- solve(diag(2) - Beta) # to save time and space
Mu_LRV <- Nu + Lambda \%*\% IB \%*\% Alpha
Sigma_LRV <- Lambda \%*\% IB \%*\% Psi \%*\% t(IB) \%*\% t(Lambda) + Theta
## Specify (unstandardized) thresholds to discretize normally distributed data
## generated from Mu_LRV and Sigma_LRV, based on marginal probabilities
PiList <- list(y1 = c(.25, .5, .25),
y2 = c(.17, .33, .33, .17),
y3 = c(.1, .2, .4, .2, .1),
## make final variable highly asymmetric
y4 = c(.33, .25, .17, .12, .08, .05))
sapply(PiList, sum) # all sum to 100\%
CumProbs <- sapply(PiList, cumsum)
## unstandardized thresholds
TauList <- mapply(qnorm, p = lapply(CumProbs, function(x) x[-length(x)]),
m = Mu_LRV, sd = sqrt(diag(Sigma_LRV)))
for (i in 1:4) names(TauList[[i]]) <- paste0(names(TauList)[i], "|t",
1:length(TauList[[i]]))
## assign numeric weights to each category (optional, see default)
NumCodes <- list(y1 = c(-0.5, 0, 0.5), y2 = 0:3, y3 = 1:5, y4 = 1:6)
## Calculate Population Moments for Numerically Coded Ordinal Variables
lrv2ord(Sigma = Sigma_LRV, Mu = Mu_LRV, thresholds = TauList, cWts = NumCodes)
## SCENARIO 2: USE ESTIMATED PARAMETERS AS POPULATION
data(datCat) # already stored as c("ordered","factor")
fit <- cfa(' f =~ 1*u1 + 1*u2 + 1*u3 + 1*u4 ', data = datCat)
lrv2ord(Sigma = fit, thresholds = fit) # use same fit for both
## or use estimated thresholds with specified parameters, but note that
## lrv2ord() will only extract standardized thresholds
dimnames(Sigma_LRV) <- list(paste0("u", 1:4), paste0("u", 1:4))
lrv2ord(Sigma = cov2cor(Sigma_LRV), thresholds = fit)
}
\references{
Jorgensen, T. D., & Johnson, A. R. (2022). How to derive expected values of
structural equation model parameters when treating discrete data as
continuous. \emph{Structural Equation Modeling, 29}(4), 639--650.
\doi{10.1080/10705511.2021.1988609}
Kamata, A., & Bauer, D. J. (2008). A note on the relation between factor
analytic and item response theory models.
\emph{Structural Equation Modeling, 15}(1), 136--153.
\doi{10.1080/10705510701758406}
Wirth, R. J., & Edwards, M. C. (2007). Item factor analysis: Current
approaches and future directions. \emph{Psychological Methods, 12}(1),
58--79. \doi{10.1037/1082-989X.12.1.58}
}
\author{
Terrence D. Jorgensen (University of Amsterdam; \email{TJorgensen314@gmail.com})
Andrew R. Johnson (Curtin University; \email{andrew.johnson@curtin.edu.au})
}
|
