# Getting Started with mxmmod

## Overview

This tutorial introduces the mxmmod package, for building Measurement Model of Derivatives (MMOD; Estabrook, 2015) with OpenMx.

## Outline

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

### Prelim: Prepare environment

First, let’s load the required libraries for this tutorial:

library(tidyverse)
library(OpenMx)
library(mxmmod)

## A. Introduction to the Measurement Model of Derivatives (MMOD)

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.

## B. Data set

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)
subset <- sample(unique(nlsy97depression\$pid), 9)

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?

## C. Example 1: One factor model

We’ll start by building a 1-factor MMOD, with all items loading onto a single latent factor.

structure <- list(
F1 = c('nervous', 'down', 'depressed', 'calm', 'happy')
)
mmod_model <- mxMmodModel(data=nlsy97depression,
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.
mmod_fit <- mxRun(mmod_model)
#> Running 1 Factor MMOD with 33 parameters
(mmod_summary <- summary(mmod_fit))
#> 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...
semPlot::semPaths(mmod_fit, 'est')

## D. Example 2: Two factor model

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):

structure2 <- list(
F1 = c('nervous', 'down', 'depressed'),
F2 = c('happy', 'calm')
)
mmod_model2 <- mxMmodModel(data=nlsy97depression,
modelName='2 Factor MMOD',
idvar='pid', timevar='occasion', structure=structure2)
#> Warning in mxMmodModel(data = nlsy97depression, modelName = "2 Factor MMOD", :
#> Missing values detected; omitting them.
mmod_fit2 <- mxRun(mmod_model2)
#> Running 2 Factor MMOD with 45 parameters
(mmod_summary2 <- summary(mmod_fit2))
#> 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...
semPlot::semPaths(mmod_fit2, 'est')

## E. Discussion

Finally, let’s create a summary table of the fits from the two models so we can compare:

fits <- list(mmod_summary, mmod_summary2)

(compare_models <- tibble(
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.