This tutorial introduces the mxmmod
package, for building Measurement Model of Derivatives (MMOD; Estabrook, 2015) with OpenMx.
A. Introduction to the Measurement Model of Derivatives (MMOD)
B. Data set
C. Example 1: One factor model
D. Example 2: Two factor model
E. Discussion
First, let’s load the required libraries for this tutorial:
library(tidyverse)
library(OpenMx)
library(mxmmod)
The Measurement Model of Derivatives (MMOD; Estabrook, 2015) is a method for evaluating test item structures that includes the temporal dynamics of item responses. Unlike traditional confirmatory factor analysis which only evaluates factor structures cross-sectionally at a single time point, an MMOD operates longitudinally, taking into account how latent factors and their associated items change over time. The MMOD makes the assumption that items from the same construct will exhibit similar temporal dynamcs (as defined by their deratives). In doing so, the MMOD can uniquely identify factor structures that would otherwise be indistinguishable cross-sectionally. By reducing the ambiguity in factor structure, the MMOD is a powerful tool to validate and sharpen theoretical distinctions among constructs in longitudinal data.
This tutorial follows the example in Estabrook (2015) and makes use of data from the National Longitudinal Survey of Youth. The NLSY 1997 sample (NLSY97) has a 5-item depression scale that was administered at three occasions. The five items are all on a 4-point Likert scale. Participants were asked how often they felt “like a nervous person”, “calm and peaceful”, “down or blue”, “like a happy person”, and “depressed” in the last month. These example data are included in the mxmmod
package:
data(nlsy97depression)
summary(nlsy97depression)
#> pid sex birth_m birth_y occasion
#> 1 : 3 M:13797 Min. : 1.000 Min. :1980 Min. :0
#> 2 : 3 F:13155 1st Qu.: 3.000 1st Qu.:1981 1st Qu.:0
#> 3 : 3 Median : 7.000 Median :1982 Median :1
#> 4 : 3 Mean : 6.556 Mean :1982 Mean :1
#> 5 : 3 3rd Qu.:10.000 3rd Qu.:1983 3rd Qu.:2
#> 6 : 3 Max. :12.000 Max. :1984 Max. :2
#> (Other):26934
#> nervous calm down happy depressed
#> Min. :1.000 Min. :1.000 Min. :1.00 Min. :1.000 Min. :1.000
#> 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:3.00 1st Qu.:2.000 1st Qu.:3.000
#> Median :3.000 Median :2.000 Median :3.00 Median :2.000 Median :4.000
#> Mean :3.187 Mean :2.386 Mean :3.15 Mean :2.205 Mean :3.572
#> 3rd Qu.:4.000 3rd Qu.:3.000 3rd Qu.:4.00 3rd Qu.:3.000 3rd Qu.:4.000
#> Max. :4.000 Max. :4.000 Max. :4.00 Max. :4.000 Max. :4.000
#> NA's :3678 NA's :3646 NA's :3716 NA's :3631 NA's :3730
Before building any models, we first plot a few example trajectories and mean trajectories for the five items assessed:
set.seed(1000)
<- sample(unique(nlsy97depression$pid), 9)
subset
%>%
nlsy97depression filter(pid %in% subset) %>%
gather(measure, val, -pid, -occasion) %>%
ggplot(aes(x=occasion, group=measure, color=measure, y=val)) +
geom_line(position=position_jitter(w=0.1, h=0.1)) +
facet_wrap(~pid)
#> Warning: attributes are not identical across measure variables;
#> they will be dropped
%>%
nlsy97depression gather(measure, val, -occasion, -pid) %>%
na.omit() %>%
ggplot(aes(x=occasion, color=measure, y=val)) +
stat_summary(fun.y = mean, geom='line') +
stat_summary(fun.y = mean, geom='point') +
stat_summary(fun.data = mean_se, geom='errorbar', width=0.2)
#> Warning: attributes are not identical across measure variables;
#> they will be dropped
#> Warning: `fun.y` is deprecated. Use `fun` instead.
#> Warning: `fun.y` is deprecated. Use `fun` instead.
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
We’ll start by building a 1-factor MMOD, with all items loading onto a single latent factor.
<- list(
structure F1 = c('nervous', 'down', 'depressed', 'calm', 'happy')
)<- mxMmodModel(data=nlsy97depression,
mmod_model modelName='1 Factor MMOD',
idvar='pid', timevar='occasion', structure=structure, fiml=F)
#> Warning in mxMmodModel(data = nlsy97depression, modelName = "1 Factor MMOD", :
#> Missing values detected; omitting them.
<- mxRun(mmod_model)
mmod_fit #> Running 1 Factor MMOD with 33 parameters
<- summary(mmod_fit))
(mmod_summary #> Summary of 1 Factor MMOD
#>
#> free parameters:
#> name matrix row col Estimate
#> 1 1 Factor MMOD.A[19,16] A dnervous_1 F1_1 0.34410308
#> 2 1 Factor MMOD.A[20,16] A ddown_1 F1_1 0.41640202
#> 3 1 Factor MMOD.A[21,16] A ddepressed_1 F1_1 0.32588085
#> 4 1 Factor MMOD.A[22,16] A dcalm_1 F1_1 -0.37367083
#> 5 1 Factor MMOD.A[23,16] A dhappy_1 F1_1 -0.38437697
#> 6 1 Factor MMOD.A[24,17] A dnervous_2 F1_2 0.15339125
#> 7 1 Factor MMOD.A[25,17] A ddown_2 F1_2 0.26887741
#> 8 1 Factor MMOD.A[26,17] A ddepressed_2 F1_2 0.21208041
#> 9 1 Factor MMOD.A[27,17] A dcalm_2 F1_2 -0.21382244
#> 10 1 Factor MMOD.A[28,17] A dhappy_2 F1_2 -0.25476948
#> 11 1 Factor MMOD.A[29,18] A dnervous_3 F1_3 0.27951244
#> 12 1 Factor MMOD.A[30,18] A ddown_3 F1_3 0.40806283
#> 13 1 Factor MMOD.A[31,18] A ddepressed_3 F1_3 0.33052618
#> 14 1 Factor MMOD.A[32,18] A dcalm_3 F1_3 -0.34191883
#> 15 1 Factor MMOD.A[33,18] A dhappy_3 F1_3 -0.38504213
#> 16 1 Factor MMOD.S[16,17] S F1_1 F1_2 -0.03903613
#> 17 1 Factor MMOD.S[16,18] S F1_1 F1_3 -0.05198215
#> 18 1 Factor MMOD.S[17,18] S F1_2 F1_3 -0.03198260
#> 19 1 Factor MMOD.S[19,19] S dnervous_1 dnervous_1 0.16422156
#> 20 1 Factor MMOD.S[20,20] S ddown_1 ddown_1 0.10017311
#> 21 1 Factor MMOD.S[21,21] S ddepressed_1 ddepressed_1 0.12320024
#> 22 1 Factor MMOD.S[22,22] S dcalm_1 dcalm_1 0.13467963
#> 23 1 Factor MMOD.S[23,23] S dhappy_1 dhappy_1 0.11556757
#> 24 1 Factor MMOD.S[24,24] S dnervous_2 dnervous_2 0.13501621
#> 25 1 Factor MMOD.S[25,25] S ddown_2 ddown_2 0.09921615
#> 26 1 Factor MMOD.S[26,26] S ddepressed_2 ddepressed_2 0.10089529
#> 27 1 Factor MMOD.S[27,27] S dcalm_2 dcalm_2 0.13235635
#> 28 1 Factor MMOD.S[28,28] S dhappy_2 dhappy_2 0.10017545
#> 29 1 Factor MMOD.S[29,29] S dnervous_3 dnervous_3 0.38803391
#> 30 1 Factor MMOD.S[30,30] S ddown_3 ddown_3 0.31074588
#> 31 1 Factor MMOD.S[31,31] S ddepressed_3 ddepressed_3 0.31667122
#> 32 1 Factor MMOD.S[32,32] S dcalm_3 dcalm_3 0.36848075
#> 33 1 Factor MMOD.S[33,33] S dhappy_3 dhappy_3 0.30162336
#> Std.Error A
#> 1 0.006321509
#> 2 0.005922424
#> 3 0.005705315
#> 4 0.006246117
#> 5 0.005947534
#> 6 0.005712803
#> 7 0.005963130
#> 8 0.005478888
#> 9 0.006227906
#> 10 0.005893614
#> 11 0.010202523
#> 12 0.010590604
#> 13 0.009869528
#> 14 0.010786293
#> 15 0.010290236
#> 16 0.016002546
#> 17 0.016533647
#> 18 0.018222365
#> 19 0.003320392
#> 20 0.002764820
#> 21 0.002654193
#> 22 0.003149643
#> 23 0.002814603
#> 24 0.002553491
#> 25 0.002705463
#> 26 0.002262882
#> 27 0.002847131
#> 28 0.002613522
#> 29 0.007620039
#> 30 0.008009461
#> 31 0.006957048
#> 32 0.008084398
#> 33 0.007521245
#>
#> Model Statistics:
#> | Parameters | Degrees of Freedom | Fit (-2lnL units)
#> Model: 33 87 -4384.387
#> Saturated: 120 0 -6755.909
#> Independence: 15 105 25163.633
#> Number of observations/statistics: 6566/120
#>
#> chi-square: χ² ( df=87 ) = 2371.522, p = 0
#> Information Criteria:
#> | df Penalty | Parameters Penalty | Sample-Size Adjusted
#> AIC: 2197.522 2437.522 2437.866
#> BIC: 1606.822 2661.581 2556.715
#> CFI: 0.9281925
#> TLI: 0.9133358 (also known as NNFI)
#> RMSEA: 0.06323938 [95% CI (0.06063494, 0.06586998)]
#> Prob(RMSEA <= 0.05): 0
#> timestamp: 2021-05-18 11:24:21
#> Wall clock time: 0.1550233 secs
#> optimizer: SLSQP
#> OpenMx version number: 2.18.1
#> Need help? See help(mxSummary)
The path diagram of this MMOD can be rendered by semPlot::semPaths
# Note: This can take a while to draw...
::semPaths(mmod_fit, 'est') semPlot
Next, let’s build a two-factor MMOD with one latent factor for negative items (nervous, down, depressed), and the other for positive items (happy, calm):
<- list(
structure2 F1 = c('nervous', 'down', 'depressed'),
F2 = c('happy', 'calm')
)<- mxMmodModel(data=nlsy97depression,
mmod_model2 modelName='2 Factor MMOD',
idvar='pid', timevar='occasion', structure=structure2)
#> Warning in mxMmodModel(data = nlsy97depression, modelName = "2 Factor MMOD", :
#> Missing values detected; omitting them.
<- mxRun(mmod_model2)
mmod_fit2 #> Running 2 Factor MMOD with 45 parameters
<- summary(mmod_fit2))
(mmod_summary2 #> Summary of 2 Factor MMOD
#>
#> free parameters:
#> name matrix row col Estimate
#> 1 2 Factor MMOD.A[22,16] A dnervous_1 F1_1 0.34208391
#> 2 2 Factor MMOD.A[23,16] A ddown_1 F1_1 0.44449759
#> 3 2 Factor MMOD.A[24,16] A ddepressed_1 F1_1 0.34568966
#> 4 2 Factor MMOD.A[25,17] A dhappy_1 F2_1 -0.43364599
#> 5 2 Factor MMOD.A[26,17] A dcalm_1 F2_1 -0.40726033
#> 6 2 Factor MMOD.A[27,18] A dnervous_2 F1_2 -0.15141710
#> 7 2 Factor MMOD.A[28,18] A ddown_2 F1_2 -0.29720493
#> 8 2 Factor MMOD.A[29,18] A ddepressed_2 F1_2 -0.22737277
#> 9 2 Factor MMOD.A[30,19] A dhappy_2 F2_2 -0.31078549
#> 10 2 Factor MMOD.A[31,19] A dcalm_2 F2_2 -0.23510537
#> 11 2 Factor MMOD.A[32,20] A dnervous_3 F1_3 0.27996036
#> 12 2 Factor MMOD.A[33,20] A ddown_3 F1_3 0.45297788
#> 13 2 Factor MMOD.A[34,20] A ddepressed_3 F1_3 0.35787907
#> 14 2 Factor MMOD.A[35,21] A dhappy_3 F2_3 0.46634671
#> 15 2 Factor MMOD.A[36,21] A dcalm_3 F2_3 0.38075384
#> 16 2 Factor MMOD.S[16,17] S F1_1 F2_1 0.78691979
#> 17 2 Factor MMOD.S[16,18] S F1_1 F1_2 0.06014108
#> 18 2 Factor MMOD.S[17,18] S F2_1 F1_2 0.03925727
#> 19 2 Factor MMOD.S[16,19] S F1_1 F2_2 -0.01587564
#> 20 2 Factor MMOD.S[17,19] S F2_1 F2_2 -0.01766206
#> 21 2 Factor MMOD.S[18,19] S F1_2 F2_2 -0.69546524
#> 22 2 Factor MMOD.S[16,20] S F1_1 F1_3 -0.06789042
#> 23 2 Factor MMOD.S[17,20] S F2_1 F1_3 -0.01426543
#> 24 2 Factor MMOD.S[18,20] S F1_2 F1_3 0.02309934
#> 25 2 Factor MMOD.S[19,20] S F2_2 F1_3 -0.03211301
#> 26 2 Factor MMOD.S[16,21] S F1_1 F2_3 0.04478221
#> 27 2 Factor MMOD.S[17,21] S F2_1 F2_3 0.03333399
#> 28 2 Factor MMOD.S[18,21] S F1_2 F2_3 -0.01278190
#> 29 2 Factor MMOD.S[19,21] S F2_2 F2_3 0.03144088
#> 30 2 Factor MMOD.S[20,21] S F1_3 F2_3 -0.68798195
#> 31 2 Factor MMOD.S[22,22] S dnervous_1 dnervous_1 0.16560711
#> 32 2 Factor MMOD.S[23,23] S ddown_1 ddown_1 0.07598566
#> 33 2 Factor MMOD.S[24,24] S ddepressed_1 ddepressed_1 0.10989724
#> 34 2 Factor MMOD.S[25,25] S dhappy_1 dhappy_1 0.07526441
#> 35 2 Factor MMOD.S[26,26] S dcalm_1 dcalm_1 0.10844859
#> 36 2 Factor MMOD.S[27,27] S dnervous_2 dnervous_2 0.13561792
#> 37 2 Factor MMOD.S[28,28] S ddown_2 ddown_2 0.08318044
#> 38 2 Factor MMOD.S[29,29] S ddepressed_2 ddepressed_2 0.09417501
#> 39 2 Factor MMOD.S[30,30] S dhappy_2 dhappy_2 0.06849531
#> 40 2 Factor MMOD.S[31,31] S dcalm_2 dcalm_2 0.12280189
#> 41 2 Factor MMOD.S[32,32] S dnervous_3 dnervous_3 0.38778330
#> 42 2 Factor MMOD.S[33,33] S ddown_3 ddown_3 0.27207224
#> 43 2 Factor MMOD.S[34,34] S ddepressed_3 ddepressed_3 0.29784132
#> 44 2 Factor MMOD.S[35,35] S dhappy_3 dhappy_3 0.23240185
#> 45 2 Factor MMOD.S[36,36] S dcalm_3 dcalm_3 0.34041584
#> Std.Error A
#> 1 0.006429926
#> 2 0.005817232
#> 3 0.005541267
#> 4 0.005941098
#> 5 0.006126876
#> 6 0.005821508
#> 7 0.006181826
#> 8 0.005478696
#> 9 0.006940047
#> 10 0.006347826
#> 11 0.010462995
#> 12 0.011082682
#> 13 0.009959067
#> 14 0.012211544
#> 15 0.011298501
#> 16 0.008525574
#> 17 0.017062146
#> 18 0.017225826
#> 19 0.016978516
#> 20 0.017103537
#> 21 0.016313765
#> 22 0.017904571
#> 23 0.018094869
#> 24 0.020276675
#> 25 0.020086240
#> 26 0.017996622
#> 27 0.018150032
#> 28 0.020313810
#> 29 0.020158937
#> 30 0.018799567
#> 31 0.003424771
#> 32 0.002732386
#> 33 0.002450655
#> 34 0.002980145
#> 35 0.003024774
#> 36 0.002579910
#> 37 0.002959253
#> 38 0.002264617
#> 39 0.003628395
#> 40 0.002904624
#> 41 0.007729192
#> 42 0.008747754
#> 43 0.007042538
#> 44 0.010048715
#> 45 0.008535637
#>
#> Model Statistics:
#> | Parameters | Degrees of Freedom | Fit (-2lnL units)
#> Model: 45 75 -5917.303
#> Saturated: 120 0 -6755.909
#> Independence: 15 105 25163.633
#> Number of observations/statistics: 6566/120
#>
#> chi-square: χ² ( df=75 ) = 838.6063, p = 2.030077e-129
#> Information Criteria:
#> | df Penalty | Parameters Penalty | Sample-Size Adjusted
#> AIC: 688.6063 928.6063 929.2412
#> BIC: 179.3818 1234.1410 1091.1423
#> CFI: 0.9759982
#> TLI: 0.9663975 (also known as NNFI)
#> RMSEA: 0.039378 [95% CI (0.03654012, 0.04226008)]
#> Prob(RMSEA <= 0.05): 1
#> timestamp: 2021-05-18 11:24:22
#> Wall clock time: 0.2585125 secs
#> optimizer: SLSQP
#> OpenMx version number: 2.18.1
#> Need help? See help(mxSummary)
The path diagram of this MMOD can be rendered by semPlot::semPaths
# Note: This can take a while to draw...
::semPaths(mmod_fit2, 'est') semPlot
Finally, let’s create a summary table of the fits from the two models so we can compare:
<- list(mmod_summary, mmod_summary2)
fits
<- tibble(
(compare_models name=map_chr(fits, 'modelName'),
chisq=map_dbl(fits, 'Chi'),
dof=map_dbl(fits, 'ChiDoF'),
`-2ll`=map_dbl(fits, 'Minus2LogLikelihood'),
aic=map_dbl(fits, 'AIC.Mx'),
bic=map_dbl(fits, 'BIC.Mx'),
rmsea=map_dbl(fits, 'RMSEA'),
cfi=map_dbl(fits, 'CFI'),
tli=map_dbl(fits, 'TLI')
))#> # A tibble: 2 x 9
#> name chisq dof `-2ll` aic bic rmsea cfi tli
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 Factor MMOD 2372. 87 -4384. 2198. 1607. 0.0632 0.928 0.913
#> 2 2 Factor MMOD 839. 75 -5917. 689. 179. 0.0394 0.976 0.966
The two-factor model is superior, across every fit metric.