add option for orthogonal factors

khusmann · Apr 3, 2021 · d3bca3f · d3bca3f
1 parent 55eb1ee
commit d3bca3f
Show file tree

Hide file tree

Showing 3 changed files with 122 additions and 7 deletions.
diff --git a/R/mxMmodModel.R b/R/mxMmodModel.R
@@ -8,6 +8,8 @@
 #' @param idvar name of column for subject IDs
 #' @param timevar name of column for measurement occasion
 #' @param structure factor structure, see 'Details'
+#' @param orthogonal if true, fix correlations between factors to 0
+#'                   (correlations within factor derivatives are still estimated)
 #' @param fiml if true, use raw data to fit model with FIML. Otherwise, fit using cov matrix
 #'             (dropping missing values if necessary).
 #' @return an MMOD as an mxModel object
@@ -36,7 +38,7 @@
 #' summary(mmod_fit)
 #' @export
 
-mxMmodModel <- function(data, modelName, idvar, timevar, structure, fiml=F) {
+mxMmodModel <- function(data, modelName, idvar, timevar, structure, orthogonal=F, fiml=F) {
   derivName <- function(o, m) {paste0('d', m, '_', o)} # derivName(1, 'nervous') -> dnervous_1
   itemName <- function(o, m) {paste0(m, '_', o)} # itemName(1, 'nervous') -> nervous_1
   factorName <- function(o, f) {paste0(f, '_', o)} # factorName(1, 'F') -> F_1
@@ -140,7 +142,13 @@ mxMmodModel <- function(data, modelName, idvar, timevar, structure, fiml=F) {
     # factor variances
     OpenMx::mxPath(from=factors, arrows=2, values=1, free=F),
     # factor correlations
-    OpenMx::mxPath(from=factors, arrows=2, connect="unique.bivariate", free=T),
+    if (orthogonal) {
+      lapply(names(structure), function(f) {
+        OpenMx::mxPath(from=factorName(occasions, f), arrows=2, connect='unique.bivariate', free=T)
+      })
+    } else {
+      OpenMx::mxPath(from=factors, arrows=2, connect="unique.bivariate", free=T)
+    },
     # residual variances(only for latent derivatives !)
     OpenMx::mxPath(from=derivatives, arrows=2, values=1)),
     # transformation

diff --git a/tests/testthat/test-refimpl.R b/tests/testthat/test-refimpl.R
@@ -1,8 +1,16 @@
 context('Reference implementation comparison')
 
-nlsy97_1983_cohort <- function() {
+nlsy97_1983_cohort <- function(omit.na=T) {
   data(nlsy97depression)
-  nlsy97depression[nlsy97depression$birth_y==1983,c('pid', 'occasion', 'nervous', 'calm', 'down', 'happy', 'depressed')]
+  df <- nlsy97depression[nlsy97depression$birth_y==1983,c('pid', 'occasion', 'nervous', 'calm', 'down', 'happy', 'depressed')]
+  if (omit.na) {
+    df_wide <- na.omit(reshape(df, timevar='occasion', idvar='pid', direction='wide',
+                               v.names=c('nervous', 'calm', 'down', 'happy', 'depressed')))
+    expect_equal(nrow(df_wide), 1397)
+    reshape(df_wide)
+  } else {
+    df
+  }
 }
 
 mxmmod_ref <- function(df, do_fiml=F) {
@@ -131,13 +139,18 @@ test_that('Floating point occasions', {
 })
 
 test_that('FIML', {
-  df <- na.omit(nlsy97_1983_cohort())
+  df <- na.omit(nlsy97_1983_cohort(omit.na=F))
   a <- mxmmod_ref(df, do_fiml=T)
   b <- estabrook_ref(df)
   expect_equal(a$parameters$Estimate[1:33], b$parameters$Estimate, tolerance = .05)
   expect_equal(a$parameters$Std.Error[1:33], b$parameters$Std.Error, tolerance = .01)
 })
 
+test_that('Omit missing values warning', {
+  df <- na.omit(nlsy97_1983_cohort(omit.na=F))
+  expect_warning(mxmmod_ref(df, do_fiml=F), regexp="Missing values")
+})
+
 test_that('Argument check timevar numeric', {
   df <- nlsy97_1983_cohort()
   df$occasion <- as.factor(df$occasion)

diff --git a/tests/testthat/test-twofactor.R b/tests/testthat/test-twofactor.R
@@ -1,8 +1,16 @@
 context('Two factor tests')
 
-nlsy97_1983_cohort <- function() {
+nlsy97_1983_cohort <- function(omit.na=T) {
   data(nlsy97depression)
-  nlsy97depression[nlsy97depression$birth_y==1983,c('pid', 'occasion', 'nervous', 'calm', 'down', 'happy', 'depressed')]
+  df <- nlsy97depression[nlsy97depression$birth_y==1983,c('pid', 'occasion', 'nervous', 'calm', 'down', 'happy', 'depressed')]
+  if (omit.na) {
+    df_wide <- na.omit(reshape(df, timevar='occasion', idvar='pid', direction='wide',
+                               v.names=c('nervous', 'calm', 'down', 'happy', 'depressed')))
+    expect_equal(nrow(df_wide), 1397)
+    reshape(df_wide)
+  } else {
+    df
+  }
 }
 
 test_that("basic two factor model works", {
@@ -170,3 +178,89 @@ test_that("factors can have crossloadings", {
   expect_equal(mmod_model$S, mmod_model_ref$S)
   expect_equal(mmod_model$F, mmod_model_ref$F)
 })
+
+test_that("orthogonal factors", {
+  require(OpenMx)
+  require(mxmmod)
+
+  df <- nlsy97_1983_cohort()
+
+  structure <- list(
+    F1 = c('nervous', 'down', 'depressed'),
+    F2 = c('calm', 'happy', 'nervous')
+  )
+  mmod_model <- mxMmodModel(data=df,
+                            modelName='2 Factor orthogonal MMOD',
+                            idvar='pid', timevar='occasion', structure=structure, fiml=F,
+                            orthogonal=T)
+  summary(mxRun(mmod_model))
+
+  dval1 <- c(1, 1, 1)
+  dval2 <- c(-1, 0, 1)
+  dval3 <- c(1/3, -2/3, 1/3)
+
+  mmod_model_ref <- mxModel(
+    type='RAM',
+    manifestVars=c('nervous_1', 'nervous_2', 'nervous_3',
+                   'down_1', 'down_2', 'down_3',
+                   'depressed_1', 'depressed_2', 'depressed_3',
+                   'calm_1', 'calm_2', 'calm_3',
+                   'happy_1', 'happy_2', 'happy_3'
+    ),
+    latentVars=  c('F1_1', 'F2_1', 'F1_2', 'F2_2', 'F1_3', 'F2_3',
+                   'dnervous_1', 'ddown_1', 'ddepressed_1', 'dcalm_1', 'dhappy_1',
+                   'dnervous_2', 'ddown_2', 'ddepressed_2', 'dcalm_2', 'dhappy_2',
+                   'dnervous_3', 'ddown_3', 'ddepressed_3', 'dcalm_3', 'dhappy_3'
+    ),
+
+    # Factor structures
+    mxPath(from='F1_1', to=c('dnervous_1', 'ddown_1', 'ddepressed_1'), values=0.5),
+    mxPath(from='F1_2', to=c('dnervous_2', 'ddown_2', 'ddepressed_2'), values=0.5),
+    mxPath(from='F1_3', to=c('dnervous_3', 'ddown_3', 'ddepressed_3'), values=0.5),
+
+    mxPath(from='F2_1', to=c('dcalm_1', 'dhappy_1', 'dnervous_1'), values=0.5), # (crossloading!)
+    mxPath(from='F2_2', to=c('dcalm_2', 'dhappy_2', 'dnervous_2'), values=0.5), # (crossloading!)
+    mxPath(from='F2_3', to=c('dcalm_3', 'dhappy_3', 'dnervous_3'), values=0.5), # (crossloading!)
+
+    # Factor Covariances (orthogonal => within cov, NO BETWEEN COV)
+    mxPath(from=c('F1_1', 'F1_2', 'F1_3'),
+           to=c('F1_1', 'F1_2', 'F1_3'), 'unique.bivariate', arrows=2),
+    mxPath(from=c('F2_1', 'F2_2', 'F2_3'),
+           to=c('F2_1', 'F2_2', 'F2_3'), 'unique.bivariate', arrows=2),
+
+    # Factor variances fixed @ 1
+    mxPath(from=c('F1_1', 'F1_2', 'F1_3', 'F2_1', 'F2_2', 'F2_3'),
+           to=c('F1_1', 'F1_2', 'F1_3', 'F2_1', 'F2_2', 'F2_3'), 'single', arrows=2, values=1, free=F),
+
+    # residual variances (for latent derivs)
+    mxPath(from=c('dnervous_1', 'ddown_1', 'ddepressed_1', 'dcalm_1', 'dhappy_1',
+                  'dnervous_2', 'ddown_2', 'ddepressed_2', 'dcalm_2', 'dhappy_2',
+                  'dnervous_3', 'ddown_3', 'ddepressed_3', 'dcalm_3', 'dhappy_3'
+    ), arrows=2, values=1),
+
+    # Build latent derivs
+    mxPath(from='dnervous_1', to=c('nervous_1', 'nervous_2', 'nervous_3'), free=F, values=dval1),
+    mxPath(from='dnervous_2', to=c('nervous_1', 'nervous_2', 'nervous_3'), free=F, values=dval2),
+    mxPath(from='dnervous_3', to=c('nervous_1', 'nervous_2', 'nervous_3'), free=F, values=dval3),
+
+    mxPath(from='ddown_1', to=c('down_1', 'down_2', 'down_3'), free=F, values=dval1),
+    mxPath(from='ddown_2', to=c('down_1', 'down_2', 'down_3'), free=F, values=dval2),
+    mxPath(from='ddown_3', to=c('down_1', 'down_2', 'down_3'), free=F, values=dval3),
+
+    mxPath(from='ddepressed_1', to=c('depressed_1', 'depressed_2', 'depressed_3'), free=F, values=dval1),
+    mxPath(from='ddepressed_2', to=c('depressed_1', 'depressed_2', 'depressed_3'), free=F, values=dval2),
+    mxPath(from='ddepressed_3', to=c('depressed_1', 'depressed_2', 'depressed_3'), free=F, values=dval3),
+
+    mxPath(from='dcalm_1', to=c('calm_1', 'calm_2', 'calm_3'), free=F, values=dval1),
+    mxPath(from='dcalm_2', to=c('calm_1', 'calm_2', 'calm_3'), free=F, values=dval2),
+    mxPath(from='dcalm_3', to=c('calm_1', 'calm_2', 'calm_3'), free=F, values=dval3),
+
+    mxPath(from='dhappy_1', to=c('happy_1', 'happy_2', 'happy_3'), free=F, values=dval1),
+    mxPath(from='dhappy_2', to=c('happy_1', 'happy_2', 'happy_3'), free=F, values=dval2),
+    mxPath(from='dhappy_3', to=c('happy_1', 'happy_2', 'happy_3'), free=F, values=dval3)
+  )
+
+  expect_equal(mmod_model$A, mmod_model_ref$A)
+  expect_equal(mmod_model$S, mmod_model_ref$S)
+  expect_equal(mmod_model$F, mmod_model_ref$F)
+})