add sex / age info to nlsy97; restrict tests to 1983 cohort like esta…

…brook 2015
khusmann · Mar 16, 2020 · 9b004db · 9b004db
1 parent 2f51b18
commit 9b004db
Show file tree

Hide file tree

Showing 6 changed files with 45 additions and 27 deletions.
diff --git a/R/data.R b/R/data.R
@@ -8,6 +8,9 @@
 #' @format A data frame with 26952 rows and 7 variables:
 #' \describe{
 #'   \item{pid}{Unique ID of participant}
+#'   \item{sex}{Sex of participant}
+#'   \item{birth_m}{Birth Month}
+#'   \item{birth_y}{Birth Year}
 #'   \item{occasion}{Measurement occasion}
 #'   \item{nervous}{How often participant felt 'like a nervous person'}
 #'   \item{calm}{How often participant felt 'calm and peaceful'}

diff --git a/data-raw/nlsy97depression.R b/data-raw/nlsy97depression.R
@@ -5,9 +5,9 @@
 
 var_list <- list(
   pid = 'R0000100',
-  # sex = 'R0536300',
-  # birth_m = 'R0536401',
-  # birth_y = 'R0536402',
+  sex = 'R0536300',
+  birth_m = 'R0536401',
+  birth_y = 'R0536402',
   # sample_type = 'R1235800',
   # eth = 'R1482600',
   nervous0 = 'R4893600',
@@ -30,12 +30,19 @@ var_list <- list(
 nlsy97depression <- read.csv('data-raw/nlsy97_subset.csv')
 nlsy97depression <- nlsy97depression[unlist(var_list)]
 colnames(nlsy97depression) <- names(var_list)
+
+# Add NAs
 nlsy97depression[nlsy97depression < 0] = NA
+
+# Wide -> Long
 nlsy97depression <- reshape(nlsy97depression,
                             idvar='pid', timevar='occasion',
-                            varying=names(var_list[-1]),
+                            varying=5:19,
                             direction='long', sep='')
+
+# Fix df attributes
 nlsy97depression['pid'] <- factor(nlsy97depression$pid)
+nlsy97depression['sex'] <- factor(nlsy97depression$sex, labels=c('M', 'F'))
 attr(nlsy97depression, "reshapeLong") <- NULL
 rownames(nlsy97depression) <- NULL
 

diff --git a/data/nlsy97depression.rda b/data/nlsy97depression.rda
diff --git a/man/nlsy97depression.Rd b/man/nlsy97depression.Rd
diff --git a/tests/testthat/test-refimpl.R b/tests/testthat/test-refimpl.R
@@ -1,5 +1,10 @@
 context('Reference implementation comparison')
 
+nlsy97_1983_cohort <- function() {
+  data(nlsy97depression)
+  nlsy97depression[nlsy97depression$birth_y==1983,c('pid', 'occasion', 'nervous', 'calm', 'down', 'happy', 'depressed')]
+}
+
 mxmmod_ref <- function(df, do_fiml=F) {
   require(OpenMx)
   require(mxmmod)
@@ -47,7 +52,7 @@ estabrook_ref <- function(df) {
   longWeight <- as.vector(t(solve(weight)))
 
   mmod2 <- mxModel("One Factor MMOD 2",
-                    mxData(covData, type="cov", numObs=nrow(na.omit(rawData))),
+                    mxData(covData, type="cov", numObs=1397),
                     type="RAM",
                     manifestVars=colnames(covData),
                     latentVars=c("F0", "F1", "F2", latDeriv),
@@ -95,17 +100,18 @@ estabrook_ref <- function(df) {
 }
 
 test_that('Same results as Estabrook 2015', {
-  data(nlsy97depression)
-  a <- mxmmod_ref(nlsy97depression)
-  b <- estabrook_ref(nlsy97depression)
+  df <- nlsy97_1983_cohort()
+  a <- mxmmod_ref(df)
+  b <- estabrook_ref(df)
   expect_equal(a$parameters$Estimate, b$parameters$Estimate, tolerance = .001)
   expect_equal(a$parameters$Std.Error, b$parameters$Std.Error, tolerance = .001)
 })
 
 test_that('Floating point occasions', {
-  data(nlsy97depression)
+  df <- nlsy97_1983_cohort()
+
   scale_val <- 2
-  scaled_nlsy <- nlsy97depression
+  scaled_nlsy <- df
   scaled_nlsy$occasion <- scaled_nlsy$occasion / scale_val
 
   trans <- rep(1, 33)
@@ -115,7 +121,7 @@ test_that('Floating point occasions', {
   trans[29:33] <- scale_val^4
 
   a <- mxmmod_ref(scaled_nlsy)
-  b <- estabrook_ref(nlsy97depression)
+  b <- estabrook_ref(df)
 
   a$parameters$Estimate <- a$parameters$Estimate / trans
   a$parameters$Std.Error <- a$parameters$Std.Error / trans
@@ -125,28 +131,27 @@ test_that('Floating point occasions', {
 })
 
 test_that('FIML', {
-  data(nlsy97depression)
-  nlsy97depression <- na.omit(nlsy97depression)
-  a <- mxmmod_ref(nlsy97depression, do_fiml=T)
-  b <- estabrook_ref(nlsy97depression)
-  expect_equal(a$parameters$Estimate[1:33], b$parameters$Estimate, tolerance = .02)
-  expect_equal(a$parameters$Std.Error[1:33], b$parameters$Std.Error, tolerance = .001)
+  df <- na.omit(nlsy97_1983_cohort())
+  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('Argument check timevar numeric', {
-  data(nlsy97depression)
-  nlsy97depression$occasion <- as.factor(nlsy97depression$occasion)
-  expect_error(mxmmod_ref(nlsy97depression), 'timevar must be a numeric type')
+  df <- nlsy97_1983_cohort()
+  df$occasion <- as.factor(df$occasion)
+  expect_error(mxmmod_ref(df), 'timevar must be a numeric type')
 })
 
 test_that('Argument check timevar', {
-  data(nlsy97depression)
-  nlsy97depression$occasion <- NULL
-  expect_error(mxmmod_ref(nlsy97depression), "Column 'occasion' not found")
+  df <- nlsy97_1983_cohort()
+  df$occasion <- NULL
+  expect_error(mxmmod_ref(df), "Column 'occasion' not found")
 })
 
 test_that('Argument check idvar', {
-  data(nlsy97depression)
-  nlsy97depression$pid <- NULL
-  expect_error(mxmmod_ref(nlsy97depression), "Column 'pid' not found")
+  df <- nlsy97_1983_cohort()
+  df$pid <- NULL
+  expect_error(mxmmod_ref(df), "Column 'pid' not found")
 })
diff --git a/vignettes/mmod_tutorial.Rmd b/vignettes/mmod_tutorial.Rmd
@@ -60,7 +60,7 @@ 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:
 
 ```{r}
-data("nlsy97depression")
+data(nlsy97depression)
 summary(nlsy97depression)
 ```