.Rbuildignore

@@ -1,3 +1,5 @@
 ^data-raw$
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^cran-comments\.md$
+^CRAN-RELEASE$

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mxmmod
 Type: Package
 Title: Measurement Model of Derivatives in 'OpenMx'
-Version: 1.0.1
+Version: 1.1.0
 Authors@R: c(
     person(
       given = 'Kyle D.',
@@ -34,5 +34,5 @@ Suggests:
     testthat,
     tidyverse
 VignetteBuilder: knitr
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.1
 Depends: R (>= 2.10)

NEWS.md

@@ -0,0 +1,13 @@
+# mxmmod 1.1.0
+
+- New feature: orthogonal structures
+- New feature: time delay embedding
+- Add testthat for two factor models
+
+# mxmmod 1.0.1
+
+Minor formatting changes for CRAN
+
+# mxmmod 1.0.0
+
+First release

R/mxMmodModel.R

@@ -8,6 +8,9 @@
 #' @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
+#'                   (A factor and its derivatives will still intercorrelate)
+#' @param embed_dim time delay embedding dimension
 #' @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 +39,13 @@
 #' summary(mmod_fit)
 #' @export
 
-mxMmodModel <- function(data, modelName, idvar, timevar, structure, fiml=F) {
+mxMmodModel <- function(data, modelName, idvar, timevar, structure, orthogonal=F, embed_dim=NULL, fiml=F) {
+  if (!is.null(embed_dim)) {
+    data <- time_delay_embed(data, idvar, timevar, embed_dim)
+    idvar <- paste0(idvar, "_embed")
+    timevar <- paste0(timevar, "_embed")
+  }
+
   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
@@ -95,7 +104,7 @@ mxMmodModel <- function(data, modelName, idvar, timevar, structure, fiml=F) {
   #     d3anxious = c('anxious_1', 'anxious_2', 'anxious_3')
   #   )
   derivStruct <- lapply(occasions, function(o) {
-    measures_flat <- unlist(structure, use.names=F)
+    measures_flat <- unique(unlist(structure, use.names=F))
     tmp <- lapply(measures_flat, function(m) {
       sapply(occasions, function(oo) {itemName(oo, m)})
     })
@@ -104,10 +113,10 @@ mxMmodModel <- function(data, modelName, idvar, timevar, structure, fiml=F) {
   })
 
   factors <- names(factorStruct)
-  derivatives <- unlist(factorStruct, use.names=F)
+  derivatives <- unique(unlist(factorStruct, use.names=F))
   manifests <- unique(unlist(derivStruct))
 
-  data <- data[c(idvar, timevar, unlist(structure))]
+  data <- data[c(idvar, timevar, unique(unlist(structure)))]
   data <- stats::reshape(as.data.frame(data), timevar=timevar, idvar=idvar, direction='wide', sep='_')[-1]
   stopifnot(setequal(manifests, names(data))) # Sanity check
 
@@ -140,7 +149,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
@@ -174,3 +189,31 @@ ContrastsGOLD <- function(T, max) {
   }
   return(Xi[1:(max+1),])
 }
+
+#' Generate Time Delay Embeddings
+#'
+#' @param data a data frame with measurements in long format
+#' @param idvar name of column for subject IDs
+#' @param timevar name of column for measurement occasion
+#' @param n_embed embedding dimension
+#'
+#' @keywords internal
+time_delay_embed <- function(data, idvar, timevar, n_embed) {
+  unique_occasions <- sort(unique(data[[timevar]]))
+  unique_ids <- unique(data[[idvar]])
+
+  n_copies <- length(unique_occasions) - n_embed + 1
+
+  embedings <- unlist(lapply(1:n_copies, function (i) unique_occasions[i:(i+n_embed-1)]))
+
+  embed_map <- data.frame(
+    idvar = rep(unique_ids, each=length(embedings)),
+    timevar = rep(embedings, length(unique_ids)),
+    embedvar = rep(rep(1:n_embed, n_copies), length(unique_ids))
+  )
+  names(embed_map) <- c(idvar, timevar, paste0(timevar, "_embed"))
+  embed_map[paste0(idvar, "_embed")] <-
+    paste0(embed_map[[idvar]], "_", rep(rep(1:n_copies, each=n_embed), length(unique_ids)))
+
+  merge(data, embed_map)
+}

cran-comments.md

@@ -0,0 +1,14 @@
+## Notes
+
+CRAN package check on solaris fails: (https://cloud.r-project.org/web/checks/check\_results\_mxmmod.html)
+
+But this is due to a segfault in the OpenMx dependency: (https://cran.r-project.org/web/checks/check\_results\_OpenMx.html)
+
+## Test environments
+* local R installation, R 3.6.3
+* ubuntu 16.04 (on travis-ci), R 3.6.3
+* win-builder (devel)
+
+## R CMD check results
+
+0 errors | 0 warnings | 0 notes

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)