Getting Started with mxmmod

Kyle Husmann

2020-03-26

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:

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      
#>  Min.   :1.000   Min.   :1.000   Min.   :1.00   Min.   :1.000  
#>  1st Qu.:3.000   1st Qu.:2.000   1st Qu.:3.00   1st Qu.:2.000  
#>  Median :3.000   Median :2.000   Median :3.00   Median :2.000  
#>  Mean   :3.187   Mean   :2.386   Mean   :3.15   Mean   :2.205  
#>  3rd Qu.:4.000   3rd Qu.:3.000   3rd Qu.:4.00   3rd Qu.:3.000  
#>  Max.   :4.000   Max.   :4.000   Max.   :4.00   Max.   :4.000  
#>  NA's   :3678    NA's   :3646    NA's   :3716   NA's   :3631   
#>    depressed    
#>  Min.   :1.000  
#>  1st Qu.:3.000  
#>  Median :4.000  
#>  Mean   :3.572  
#>  3rd Qu.:4.000  
#>  Max.   :4.000  
#>  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

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   Std.Error A
#> 1  1 Factor MMOD.A[19,16]      A  19  16  0.34410298 0.006321527  
#> 2  1 Factor MMOD.A[20,16]      A  20  16  0.41640197 0.005922426  
#> 3  1 Factor MMOD.A[21,16]      A  21  16  0.32588076 0.005705342  
#> 4  1 Factor MMOD.A[22,16]      A  22  16 -0.37367085 0.006246139  
#> 5  1 Factor MMOD.A[23,16]      A  23  16 -0.38437702 0.005947542  
#> 6  1 Factor MMOD.A[24,17]      A  24  17  0.15339097 0.005712813  
#> 7  1 Factor MMOD.A[25,17]      A  25  17  0.26887738 0.005963080  
#> 8  1 Factor MMOD.A[26,17]      A  26  17  0.21208038 0.005478890  
#> 9  1 Factor MMOD.A[27,17]      A  27  17 -0.21382262 0.006227886  
#> 10 1 Factor MMOD.A[28,17]      A  28  17 -0.25476963 0.005893632  
#> 11 1 Factor MMOD.A[29,18]      A  29  18  0.27951265 0.010202608  
#> 12 1 Factor MMOD.A[30,18]      A  30  18  0.40806264 0.010590678  
#> 13 1 Factor MMOD.A[31,18]      A  31  18  0.33052611 0.009869613  
#> 14 1 Factor MMOD.A[32,18]      A  32  18 -0.34191912 0.010786454  
#> 15 1 Factor MMOD.A[33,18]      A  33  18 -0.38504243 0.010290284  
#> 16 1 Factor MMOD.S[16,17]      S  16  17 -0.03903612 0.016001941  
#> 17 1 Factor MMOD.S[16,18]      S  16  18 -0.05198218 0.016533674  
#> 18 1 Factor MMOD.S[17,18]      S  17  18 -0.03198219 0.018222327  
#> 19 1 Factor MMOD.S[19,19]      S  19  19  0.16422156 0.003320386  
#> 20 1 Factor MMOD.S[20,20]      S  20  20  0.10017308 0.002764829  
#> 21 1 Factor MMOD.S[21,21]      S  21  21  0.12320016 0.002654192  
#> 22 1 Factor MMOD.S[22,22]      S  22  22  0.13467958 0.003149646  
#> 23 1 Factor MMOD.S[23,23]      S  23  23  0.11556750 0.002814601  
#> 24 1 Factor MMOD.S[24,24]      S  24  24  0.13501631 0.002553492  
#> 25 1 Factor MMOD.S[25,25]      S  25  25  0.09921611 0.002705456  
#> 26 1 Factor MMOD.S[26,26]      S  26  26  0.10089526 0.002262878  
#> 27 1 Factor MMOD.S[27,27]      S  27  27  0.13235628 0.002847126  
#> 28 1 Factor MMOD.S[28,28]      S  28  28  0.10017537 0.002613517  
#> 29 1 Factor MMOD.S[29,29]      S  29  29  0.38803367 0.007620135  
#> 30 1 Factor MMOD.S[30,30]      S  30  30  0.31074598 0.008009581  
#> 31 1 Factor MMOD.S[31,31]      S  31  31  0.31667128 0.006957156  
#> 32 1 Factor MMOD.S[32,32]      S  32  32  0.36848071 0.008084456  
#> 33 1 Factor MMOD.S[33,33]      S  33  33  0.30162331 0.007521336  
#> 
#> 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: 2020-03-26 13:56:45 
#> Wall clock time: 0.1772234 secs 
#> optimizer:  CSOLNP 
#> OpenMx version number: 2.15.4.10 
#> 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   Std.Error A
#> 1  2 Factor MMOD.A[22,16]      A  22  16  0.34208388 0.006429899  
#> 2  2 Factor MMOD.A[23,16]      A  23  16  0.44449741 0.005817207  
#> 3  2 Factor MMOD.A[24,16]      A  24  16  0.34568953 0.005541246  
#> 4  2 Factor MMOD.A[25,17]      A  25  17  0.43364577 0.005941030  
#> 5  2 Factor MMOD.A[26,17]      A  26  17  0.40726028 0.006126788  
#> 6  2 Factor MMOD.A[27,18]      A  27  18  0.15141701 0.005821460  
#> 7  2 Factor MMOD.A[28,18]      A  28  18  0.29720491 0.006181797  
#> 8  2 Factor MMOD.A[29,18]      A  29  18  0.22737258 0.005478688  
#> 9  2 Factor MMOD.A[30,19]      A  30  19  0.31078553 0.006939844  
#> 10 2 Factor MMOD.A[31,19]      A  31  19  0.23510537 0.006347763  
#> 11 2 Factor MMOD.A[32,20]      A  32  20  0.27996023 0.010462904  
#> 12 2 Factor MMOD.A[33,20]      A  33  20  0.45297787 0.011082778  
#> 13 2 Factor MMOD.A[34,20]      A  34  20  0.35787913 0.009959036  
#> 14 2 Factor MMOD.A[35,21]      A  35  21  0.46634674 0.012210854  
#> 15 2 Factor MMOD.A[36,21]      A  36  21  0.38075360 0.011297882  
#> 16 2 Factor MMOD.S[16,17]      S  16  17 -0.78691973 0.008525433  
#> 17 2 Factor MMOD.S[16,18]      S  16  18 -0.06014055 0.017059958  
#> 18 2 Factor MMOD.S[17,18]      S  17  18  0.03925656 0.017223501  
#> 19 2 Factor MMOD.S[16,19]      S  16  19  0.01587567 0.016975544  
#> 20 2 Factor MMOD.S[17,19]      S  17  19 -0.01766174 0.017100865  
#> 21 2 Factor MMOD.S[18,19]      S  18  19 -0.69546523 0.016312123  
#> 22 2 Factor MMOD.S[16,20]      S  16  20 -0.06789090 0.017903367  
#> 23 2 Factor MMOD.S[17,20]      S  17  20  0.01426604 0.018092996  
#> 24 2 Factor MMOD.S[18,20]      S  18  20 -0.02309915 0.020277245  
#> 25 2 Factor MMOD.S[19,20]      S  19  20  0.03211214 0.020082629  
#> 26 2 Factor MMOD.S[16,21]      S  16  21  0.04478213 0.017994100  
#> 27 2 Factor MMOD.S[17,21]      S  17  21 -0.03333397 0.018147184  
#> 28 2 Factor MMOD.S[18,21]      S  18  21  0.01278158 0.020314127  
#> 29 2 Factor MMOD.S[19,21]      S  19  21 -0.03144049 0.020155622  
#> 30 2 Factor MMOD.S[20,21]      S  20  21 -0.68798156 0.018798232  
#> 31 2 Factor MMOD.S[22,22]      S  22  22  0.16560698 0.003424758  
#> 32 2 Factor MMOD.S[23,23]      S  23  23  0.07598569 0.002732376  
#> 33 2 Factor MMOD.S[24,24]      S  24  24  0.10989722 0.002450649  
#> 34 2 Factor MMOD.S[25,25]      S  25  25  0.07526435 0.002980122  
#> 35 2 Factor MMOD.S[26,26]      S  26  26  0.10844852 0.003024748  
#> 36 2 Factor MMOD.S[27,27]      S  27  27  0.13561787 0.002579903  
#> 37 2 Factor MMOD.S[28,28]      S  28  28  0.08318042 0.002959210  
#> 38 2 Factor MMOD.S[29,29]      S  29  29  0.09417507 0.002264613  
#> 39 2 Factor MMOD.S[30,30]      S  30  30  0.06849524 0.003628240  
#> 40 2 Factor MMOD.S[31,31]      S  31  31  0.12280186 0.002904601  
#> 41 2 Factor MMOD.S[32,32]      S  32  32  0.38778324 0.007729036  
#> 42 2 Factor MMOD.S[33,33]      S  33  33  0.27207204 0.008747856  
#> 43 2 Factor MMOD.S[34,34]      S  34  34  0.29784116 0.007042508  
#> 44 2 Factor MMOD.S[35,35]      S  35  35  0.23240173 0.010047949  
#> 45 2 Factor MMOD.S[36,36]      S  36  36  0.34041591 0.008535253  
#> 
#> 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: 2020-03-26 13:56:45 
#> Wall clock time: 0.2176533 secs 
#> optimizer:  CSOLNP 
#> OpenMx version number: 2.15.4.10 
#> 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.