-
Notifications
You must be signed in to change notification settings - Fork 1
/
semlrtp.Rmd
266 lines (215 loc) · 5.27 KB
/
semlrtp.Rmd
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
---
title: "Get Started"
author: "Shu Fai Cheung & Mark Hok Chio Lai"
date: "`r Sys.Date()`"
output:
rmarkdown::html_vignette:
number_sections: true
vignette: >
%\VignetteIndexEntry{Get Started}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: references.bib
csl: apa.csl
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 6,
fig.height = 6,
fig.align = "center",
fig.path = ""
)
```
# Introduction
This article illustrates how to use
[`semlrtp`](https://sfcheung.github.io/semlrtp/) to
compute LRT *p*-values for selected free parameters in
a model fitted by `lavaan`.
These are the packages needed for this illustration:
```{r}
library(lavaan)
library(semlrtp)
```
# LRT *p*-Values
What we called the LRT *p*-value of a parameter is
simply the *p*-value of the likelihood
ratio test (LRT, a.k.a. $\chi^2$ difference test)
comparing the fitted model with the model with
this parameter fixed to zero. It is not new but
we are not aware of packages using it as the default
*p*-value in testing model parameters.
The package
[`semlrtp`](https://sfcheung.github.io/semlrtp/)
is developed to facilitate the use of
this *p*-value.
# Basic Workflow
## Data
This is the sample dataset from the package,
with 16 variables and a group variable:
```{r}
data(data_sem16)
print(head(data_sem16), digits = 2)
```
## Model
This is a model to be fitted, with
four latent factors, and a structural
model for the factors:
```{r}
mod <-
"
f1 =~ x1 + x2 + x3
f2 =~ x5 + x6 + x7
f3 =~ x9 + x10 + x11
f4 =~ x13 + x14 + x15
f3 ~ f1 + f2
f4 ~ f3
"
```
We first fit a model to the whole sample:
```{r}
fit <- sem(mod,
data_sem16)
```
## LRT *p*-Values For Selected Parameters
If we use the default settings,
we can compute the LRT *p*-values just
by calling `lrtp()` on the `lavann`
output.
By default, LRT *p*-values will be computed
only for regression paths
(`"~"`) and covariances (`"~~"`),
excluding variances and error covariances.
```{r}
fit_lrtp <- lrtp(fit)
```
By default, the output will be printed
in `lavaan` style, and only parameters
with LRT *p*-values are printed:
```{r}
fit_lrtp
```
The column `Chisq` shows the $\chi^2$
difference of the likelihood ratio test
when a parameter is fixed to zero.
The column `LRTp` shows the LRT *p*-value.
## How LRT *p*-Values Are Computed
It can be verified that the LRT *p*-value
of a parameter is the likelihood ratio
test (LRT) *p*-value (a.k.a. the $\chi^2$
difference test) when this parameter is
fitted to zero.
For example, we fix the path from `f2`
to `f3` to zero and then do an LR test:
```{r}
mod_f3_f2 <-
"
f1 =~ x1 + x2 + x3
f2 =~ x5 + x6 + x7
f3 =~ x9 + x10 + x11
f4 =~ x13 + x14 + x15
f3 ~ f1 + 0*f2
f4 ~ f3
"
fit_f3_f2 <- sem(mod_f3_f2,
data_sem16)
lavTestLRT(fit_f3_f2,
fit)
```
```{r echo = FALSE}
est <- parameterEstimates(fit)
f3_f2_wald_p <- est[14, "pvalue"]
f3_f2_lrt_p <- fit_lrtp[14, "LRTp"]
```
Unlike the original *p*-value of
this path,
`r formatC(f3_f2_wald_p, digits = 3, format = "f")`,
the LRT *p*-value is
`r formatC(f3_f2_lrt_p, digits = 3, format = "f")`,
suggesting that the path from `f2` to `f2`
is significant.
# Why LRT *p*-Value
The example above illustrates the importance
of the LRT *p*-value. The usual *p*-values
in `lavaan` (and many other SEM
programs) are Wald-based *p*-values.
The Wald-based p-value is an approximation
of the
LRT *p*-value when a parameter is fixed
to zero. It is an approximation and
so can be different from the
LRT *p*-value. Moreover, they may also
depend on
the parameterization [@gonzalez_testing_2001].
For example, we can fit the same model
by changing the indicators being fixed to 1.
```{r}
mod2 <-
"
f1 =~ x2 + x1 + x3
f2 =~ x7 + x5 + x6
f3 =~ x11 + x9 + x10
f4 =~ x14 + x13 + x15
f3 ~ f1 + f2
f4 ~ f3
"
fit2 <- sem(mod2,
data_sem16)
```
This model and the previous one have exactly
identical model fit, as expected:
```{r}
fitMeasures(fit, c("chisq", "df"))
fitMeasures(fit2, c("chisq", "df"))
```
This is the parameter estimates with
Wald *p*-values (only those for
the regression paths are displayed)
```{r echo = FALSE}
est2 <- parameterEstimates(fit2, output = "text")
tmp <- capture.output(est2)
f3_f2_wald_p2 <- est2[14, "pvalue"]
```
```r
parameterEstimates(fit2,
output = "text")
```
```{r echo = FALSE}
cat(tmp[21:27], sep = "\n")
```
The Wald *p*-value is
`r formatC(f3_f2_wald_p2, digits = 3, format = "f")`,
even larger than
`r formatC(f3_f2_wald_p, digits = 3, format = "f")`
in the original model.
However, the LRT *p*-values are the same:
```{r}
fit2_lrtp <- lrtp(fit2)
fit2_lrtp
```
It is because LRT *p*-value is invariant to
parameterization.
# Limitations
Because LRT *p*-value are computed by
fixing a parameter to zero, there are
parameter for which the LRT *p*-value
cannot be computed. For example,
suppose we request LRT *p*-values for
factor loadings using the argument
`op`:
```{r}
fit_lrtp_loadings <- lrtp(fit,
op = "=~")
fit_lrtp_loadings
```
As shown above, LRT *p*-values are not
computed for indicators with loadings
fixed to zero.
# Further Information
Please refer to the help page of
`lrtp()` for other arguments, and
the print method of `lrtp()` output
(`print.lrtp()`) for options
in printing.
# References