Bring branch up to date with master

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: chartplotter
 Type: Package
 Title: Analysing and Plotting Growth Curves
-Version: 0.31.2
+Version: 0.33.0
 Authors@R: person("Stef", "van Buuren", email = "stef.vanbuuren@tno.nl", role = c("aut", "cre"))
 Maintainer: Stef van Buuren <stef.vanbuuren@tno.nl>
 Description: This package plots the target child's growth curves on 

NEWS.md

@@ -1,3 +1,15 @@
+# chartplotter 0.33.0
+
+- Resolves the incorrect time sequence of the WFH curve when height at a later age is lower. See <https://github.com/growthcharts/james/issues/24> It works for `curve_interpolation` is `FALSE`. Some approximation error remain for `curve_interpolation` is `TRUE`.
+
+# chartplotter 0.32.1
+
+- Set color blue for preterm D-score chart
+
+# chartplotter 0.32.0
+
+- Solves a problem resulted in `Error in eval(predvars, data, env) : object 'hgt_z_0' not found`. The problem occured in `find_matches()` when `hat` contains no brokenstick estimates. This may happen when child observations do not start at birth. The function now removes names of earlier predictors and matches only on covariates.
+
 # chartplotter 0.31.2
 
 - Sets colors to green for new D-score references (despite their WHOpink/blue defaults)
@@ -119,7 +131,7 @@ This version major extension of the chartplotter package.
 
 # chartplotter 0.13.1
 
-# Repairs prediction line plotting
+* Repairs prediction line plotting
 
 # chartplotter 0.13.0
 

R/find_matches.R

@@ -102,7 +102,8 @@ find_matches <- function(target,
                      make_xname(x,
                                 xnames_complete,
                                 user_model = user_model,
-                                current_age = period[1L]
+                                current_age = period[1L],
+                                hat
                      )
                    },
                    simplify = FALSE

R/internal.R

@@ -25,18 +25,12 @@ set_xout <- function(chartcode, yname) {
   if (design == "A") {
     return(round(seq(0.5, 15, 0.5) / 12, 4L))
   }
-  if (design == "B" & yname == "wfh") {
-    return(round(seq(50, 120, by = 2), 4L))
-  }
   if (design == "B" & yname %in% c("hgt", "dsc")) {
     return(round(c(0.5, 0.75, 1:48) / 12, 4L))
   }
-  if (design == "B" & yname == "hdc") {
+  if (design == "B" & yname %in% c("hdc", "wfh")) {
     return(round(seq(0.1, 4, by = 0.1), 4L))
   }
-  if (design == "C" & yname == "wfh") {
-    return(round(seq(60, 184, by = 4), 4L))
-  }
   if (design == "C") {
     return(round(seq(1, 21, by = 0.5), 4L))
   }

R/make_xname.R

@@ -7,6 +7,11 @@
 #'   `paste(yname, d, sep = "_")`, where `d` is a vector
 #'   contuining decimal ages.
 #' @param current_age Current age (in decimal years)
+#' @param hat A list with brokenstick prediction for the target that
+#' contains the names of the brokenstick proedictions. The entry is
+#' `NULL` if these estimates cannot be made. In that case, this
+#' function does not include the previous yname scores into the
+#' matching.
 #' @inheritParams process_chart
 #' @return A character vector with the `xname` specification
 #' @examples
@@ -23,15 +28,22 @@
 make_xname <- function(yname,
                        xnames,
                        user_model,
-                       current_age) {
+                       current_age,
+                       hat = NULL) {
   covariates <- c(
     "sex", "etn", "ga", "bw", "twin", "smo",
     "edu", "agem", "hgtf", "wgtf", "hgtm", "wgtm",
     "durbrst"
   )
   ynames <- c("hgt", "wgt", "hdc", "bmi", "wfh", "dsc")
   yname <- match.arg(yname, ynames)
-  xn <- get_xname(yname, xnames)
+  # read names from hat if brokenstick estimates are present
+  # else create names
+  if (missing(hat)) {
+    xn <- get_xname(yname, xnames)
+  } else {
+    xn <- names(hat[[yname]])
+  }
   xa <- get_age(xn)
 
   idx <- xa <= current_age

R/process_chart.R

@@ -93,11 +93,12 @@ process_chart <- function(target,
   }
 
   # set the palette
-  # set green for D-score plots despite their WHOblue and WHOpink defaults
+  # set green/blue for D-score plots despite their WHOblue and WHOpink defaults
   palettes <- chartbox::palettes
   parsed <- parse_chartcode(chartcode)
   if (parsed$side == "dsc") {
-    palette(palettes["NL",])
+    if (parsed$week == "40") palette(palettes["NL", ])
+    else palette(palettes["PT", ])
   } else {
     palette(palettes[parsed$population,])
   }

R/set_curves.R

@@ -37,9 +37,13 @@ set_curves <- function(g,
     mutate(
       id = -1,
       sex = child$sex,
-      ga = child$ga
+      ga = child$ga,
+      x2 = .data$x
     ) %>%
-    select(all_of(c("id", "age", "xname", "yname", "x", "y", "sex", "ga")))
+    select(all_of(c("id", "age", "xname", "yname", "x", "y", "sex", "ga", "x2")))
+  # For WFH, temporary sort on age to get correct time sequence, use x2 to store x
+  idx <- data$yname == "wfh"
+  data$x[idx] <- data$age[idx]
 
   # get data of matches
   time <- vector("list", length(ynames))
@@ -107,7 +111,17 @@ set_curves <- function(g,
   # append synthetic data
   data <- data %>%
     bind_rows(synt) %>%
-    arrange(.data$id, .data$yname, .data$x, .data$age) %>%
+    arrange(.data$id, .data$yname, .data$x)
+
+  # For wfh, interpolate hgt from age, and overwrite data$x with hgt
+  idx <- data$yname == "wfh"
+  if (any(idx)) {
+    data$x[idx] <- safe_approx(x = data$x[idx], y = data$x2[idx],
+                               xout = data$x[idx], ties = list("ordered", mean))$y
+  }
+
+  # add Z-scores
+  data <- data %>%
     select(-"age") %>%  # fool set_refcodes()
     mutate(refcode_z = nlreferences::set_refcodes(.)) %>%
     mutate(
@@ -133,8 +147,7 @@ set_curves <- function(g,
     group_by(.data$id, .data$yname, .data$pred) %>%
     mutate(z = safe_approx(
       x = .data$x, y = .data$z, xout = .data$x,
-      ties = list("ordered", mean)
-    )$y) %>%
+      ties = mean)$y) %>%
     ungroup()
 
   # set refcode as target's sex and ga

tests/testthat/test-process_chart.R

@@ -4,7 +4,7 @@ library(brokenstick)
 library(nlreferences)
 
 test_that("returns empty chart if not a list", {
-          expect_silent(process_chart(NULL, chartcode = "NJAA"))
+  expect_silent(process_chart(NULL, chartcode = "NJAA"))
 })
 
 target <- list(
@@ -20,21 +20,43 @@ fn <- system.file("extdata", "bds_v1.0", "smocc", "Laura_S.json", package = "jam
 ind <- bdsreader::read_bds(fn)
 
 g <- process_chart(ind, chartcode = "NJAA",
-                  dnr = "smocc", period = c(0.5, 1.1667), nmatch = 10, break_ties = TRUE)
+                   dnr = "smocc", period = c(0.5, 1.1667), nmatch = 10, break_ties = TRUE)
+
+# incorrect order of observations for WFH when hgt is not monotone
+day <- c(0, 13, 42, 91, 152, 287, 336, 434, 541, 632, 744, 905)
+hgt <- c(NA, NA, 56, 61.5, 67, 72.5, 74, 78, 83, 84, 89, 88)
+wgt <- c(2.879, 3.14, 4.4, 6.055, 7.15, 7.915, 8.25, 9.45, 10.8, 10.45, 10.8, 11.9)
+df <- data.frame(age = round(day / 365.25, 4), hgt, wgt)
+xyz <- ind$xyz[ind$xyz$yname == "wfh", ]
+xyz$age <- df$age[-(1:2)]
+xyz$x <- df$hgt[-(1:2)]
+xyz$y <- df$wgt[-(1:2)]
+xyz$z <- centile::y2z(y = xyz$y, x = xyz$x, refcode = "nl_1997_wfh_female_nla",
+                      pkg = "nlreferences", rule = 2L)
+data <- ind
+data$xyz <- xyz
+
+# Note WFH curve age sequence: correct for curve_interpolation is FALSE
+# But has approximation errors when curve_interpolation is TRUE
+test_that("Weight for height curve has correct time sequence", {
+  expect_silent(process_chart(data, chartcode = "NJBR",
+                              dnr = "smocc", period = c(0.5, 1.1667), nmatch = 10,
+                              break_ties = TRUE, curve_interpolation = FALSE))
+})
 
 test_that("terneuzen donordata yields matches", {
   expect_silent(process_chart(ind, chartcode = "NJCH",
-                   dnr = "terneuzen", period = c(2, 18),
-                   nmatch = 10, break_ties = TRUE))
+                              dnr = "terneuzen", period = c(2, 18),
+                              nmatch = 10, break_ties = TRUE))
 })
 
 test_that("prediction line connects last observation to prediction", {
   # does not warn anymore for mutate_ in curvematching::calculate_matches()
   expect_silent(process_chart(ind,
-    chartcode = "NJCH",
-    dnr = "terneuzen", period = c(3, 10),
-    nmatch = 25, break_ties = TRUE,
-    show_realized = TRUE, show_future = TRUE
+                              chartcode = "NJCH",
+                              dnr = "terneuzen", period = c(3, 10),
+                              nmatch = 25, break_ties = TRUE,
+                              show_realized = TRUE, show_future = TRUE
   ))
 })
 
@@ -43,17 +65,17 @@ ind <- bdsreader::read_bds(fn)
 test_that("Kevin S is drawn silently", {
   # warns for mutate_ in curvematching::calculate_matches()
   expect_silent(process_chart(ind,
-    chartcode = "PJAAN34", dnr = "smocc", period = c(0.6, 1.1667),
-    nmatch = 10, exact_ga = FALSE, break_ties = TRUE,
-    show_future = TRUE, show_realized = TRUE, curve_interpolation = TRUE
+                              chartcode = "PJAAN34", dnr = "smocc", period = c(0.6, 1.1667),
+                              nmatch = 10, exact_ga = FALSE, break_ties = TRUE,
+                              show_future = TRUE, show_realized = TRUE, curve_interpolation = TRUE
   ))
 })
 
 test_that("Kevin S predict hdc using lollypop", {
   # warns for mutate_ in curvematching::calculate_matches()
   expect_silent(process_chart(ind,
-    chartcode = "PJAAN34", dnr = "lollypop", period = c(0.6, 1.1667),
-    nmatch = 10, show_future = TRUE, show_realized = TRUE
+                              chartcode = "PJAAN34", dnr = "lollypop", period = c(0.6, 1.1667),
+                              nmatch = 10, show_future = TRUE, show_realized = TRUE
   ))
 })
 
@@ -84,8 +106,8 @@ fn <- system.file("extdata", "bds_v1.0", "smocc", "Laura_S.json", package = "jam
 ind <- bdsreader::read_bds(fn)
 test_that("D-score prediction does not go beyond 24 months", {
   expect_silent(g <- process_chart(ind,
-    chartcode = "NMBD", period = c(1, 3),
-    nmatch = 10, show_realized = TRUE, show_future = TRUE
+                                   chartcode = "NMBD", period = c(1, 3),
+                                   nmatch = 10, show_realized = TRUE, show_future = TRUE
   ))
 })
 
@@ -95,7 +117,7 @@ ind <- bdsreader::read_bds(jtf[5])
 
 test_that("test5.json passes individual_to_donordata()", {
   expect_silent(process_chart(ind, chartcode = "NJBA", nmatch = 1, period = c(0, 1)))
-  })
+})
 
 
 # g <- process_chart(ind, chartcode = "NMBD", period = c(1, 2),
@@ -147,3 +169,17 @@ test_that("Height plots on NMCH", {
 # ind <- bdsreader::read_bds(fn)
 # g <- process_chart(ind, chartcode = "NMAH")
 # grid::grid.draw(g)
+
+# Mar 2024
+# Error in eval(predvars, data, env) : object 'hgt_z_0' not found
+# Occurs when period[1] < min(age)
+# Solution: predict from the population mean/no prediction/random sample?
+fn <- system.file("examples/maria.json", package = "bdsreader")
+tgt <- bdsreader::read_bds(fn)
+test_that("Matches do not condition on yname when there are no brokenstick estimates", {
+  expect_silent(process_chart(
+    target = tgt,
+    chartcode = "PMAHN27",
+    nmatch = 10,
+    period = c(0.01, 1.1667)))
+})