diff --git a/.Rbuildignore b/.Rbuildignore
new file mode 100644
index 0000000..45a7797
--- /dev/null
+++ b/.Rbuildignore
@@ -0,0 +1,4 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^README\.Rmd$
+^data-raw$
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5b6a065
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..85399f7
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,29 @@
+Package: growthscreener
+Title: Finding Children with Unusual Growth Patterns
+Version: 1.0.0
+Authors@R: 
+    c(person(given = "Paula",
+           family = "van Dommelen",
+           role = c("aut", "cre"),
+           email = "paula.vandommelen@tno.nl"),
+       person(given = "Stef",
+           family = "van Buuren",
+           role = c("aut"),
+           email = "stef.vanbuuren@tno.nl"))
+Description: Unusual child growth may be a sign of an underlying 
+    disease or condition. This package implements functions that 
+    evaluate the growth pattern of a child relative to Dutch 
+    guidelines.
+Depends: 
+    R (>= 2.10)
+Imports: 
+    AGD (>= 0.41.0)
+Remotes: 
+    stefvanbuuren/AGD
+License: AGPL-3 + file LICENSE
+Encoding: UTF-8
+LazyData: true
+Roxygen: list(markdown = TRUE)
+RoxygenNote: 6.1.1
+Suggests: 
+    testthat (>= 2.1.0)
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..d82115c
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,31 @@
+    growthscreener: Finding Children with Unusual Growth Patterns
+    Copyright (C) 2019 Paula van Dommelen, Stef van Buuren
+
+    Source: <https://github.com/stefvanbuuren/growthscreener>
+    Contact: paula.vandommelen at tno.nl
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Affero General Public License as
+    published by the Free Software Foundation, either version 3 of the
+    License, or (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Affero General Public License for more details.
+
+    You should have received a copy of the GNU Affero General Public License
+    along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+    If your software can interact with users remotely through a computer
+    network, you should also make sure that it provides a way for users
+    to get its source. For example, if your program is a web application,
+    its interface could display a "Source" link that leads users to an
+    archive of the code. There are many ways you could offer source,
+    and different solutions will be better for different programs; see
+    section 13 for the specific requirements.
+
+    You should also get your employer (if you work as a programmer) or
+    school, if any, to sign a "copyright disclaimer" for the program,
+    if necessary. For more information on this, and how to apply and
+    follow the GNU AGPL, see <https://www.gnu.org/licenses/>.
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..94d4cb2
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,9 @@
+# Generated by roxygen2: do not edit by hand
+
+export(calculate_advice_hgt)
+export(calculate_birth_z)
+export(calculate_helpers_hgt)
+export(calculate_th)
+export(fold)
+export(msg)
+importFrom(AGD,y2z)
diff --git a/R/calculate_advice_hgt.R b/R/calculate_advice_hgt.R
new file mode 100644
index 0000000..4bfd6a8
--- /dev/null
+++ b/R/calculate_advice_hgt.R
@@ -0,0 +1,148 @@
+#' Referral advice for body height
+#'
+#' This function traverses the decision tree of the
+#' "JGZ-Richtlijn Lengtegroei 2019".
+#'
+#' The decision tree assesses both single and paired measurements.
+#' The last observations (\code{y1}) is generally taken as the
+#' last measurement, whereas \code{y0} can be one of the previous
+#' measurements. For more than two measurements, there are many
+#' pairs possible, and these pairs need not be consecutive.
+#' The \code{y0} measurement needs to be defined by the user,
+#' and is informally taken as an earlier measurement that maximumizes
+#' the referal probability. On the other hand, defining pairs that are
+#' remote in ages (e.g. period between 1 month and 14 years) is probably
+#' not that useful. In practice, we may be interested in setting the
+#' maximum period to, say, five years.
+#' @param sex   Character, either \code{"M"} (male) or \code{"F"} (female)
+#' @param dob   Date of birth (class Date)
+#' @param bw    Birth weight (grammes)
+#' @param bl    Birth length (cm)
+#' @param ga    Gestational age, completed weeks (Integer or character)
+#' @param etn   Etnicity, one of \code{"N"} (dutch), \code{"T"} (turkish),
+#'              \code{"M"} (moroccan) or \code{"H"} (hindustani).
+#' @param hgtf  Height of father (cm)
+#' @param hgtm  Height of mother (cm)
+#' @param dom1  Date of last measurement (Date)
+#' @param y1    Height at last measurement (cm)
+#' @param dom0  Date of previous measurement (Date)
+#' @param y0    Height at previous measurement (cm)
+#' @param d     Optional, list of derived variables, obtained by
+#'              \code{calculate_derived_variables()}
+#' @return \code{calculate_advice_hgt} returns an integer, the \code{msgcode}
+#' @author Paula van Dommelen, Stef van Buuren, 2019
+#' @rdname advice_hgt
+#' @examples
+#' msg(calculate_advice_hgt())
+#' msgcode <- calculate_advice_hgt(sex = "M", etn = "N",
+#'                                 dob = as.Date("2018-07-31"),
+#'                                 dom1 = as.Date("2018-12-12"), y1 = 64)
+#' msg(msgcode)
+#' @export
+calculate_advice_hgt <- function(sex = NA_character_, dob = as.Date(NA),
+                                 bw = NA, bl = NA, ga = NA,
+                                 etn = NA_character_, hgtf = NA, hgtm = NA,
+                                 dom1 = as.Date(NA), y1 = NA,
+                                 dom0 = as.Date(NA), y0 = NA,
+                                 d = NULL) {
+
+  if (is.null(d))
+    d <- calculate_helpers_hgt(sex, dob, bw, bl, ga, etn,
+                               hgtf, hgtm, dom1, y1, dom0, y0)
+  bw_z <- d$bw_z
+  bl_z <- d$bl_z
+  th_z <- d$th_z
+  age1 <- d$age1
+  age0 <- d$age0
+  z1   <- d$z1
+  z0   <- d$z0
+
+  # start the sieve
+
+  # return early if data are insufficient
+  # if (is.na(sex)) return(19)
+  if (!sex %in% c("M", "F")) return(19)
+  if (is.na(dob)) return(16)
+  if (is.na(dom1)) return(15)
+  if (is.na(y1)) return(ifelse(age1 < 18.0, 18, 21))
+  if (!etn %in% c("N", "T", "M", "H")) return(20)
+
+  # outside age range
+  if (age1 >= 18.0) return(21)
+  if (age1 < 0.0833) return(22)
+
+  # check single measurement
+  if (age1 < 3.0) {
+    # short
+    if (z1 < -2.5 & is.na(bw)) return(13)
+    if (z1 < -3.0 & bw >= 2500) return(45)
+
+    # tall
+    if (z1 > 3.0) return(48)
+    if (z1 > 2.5 & is.na(z0)) return(11)
+    if (z1 > 1.0) return(77)
+  }
+
+  if (age1 >= 3.0 & age1 < 10.0) {
+
+    # short
+    if (z1 < -2.5) return(44)
+    if (z1 < -2.0 & !is.na(bw_z)) if (bw_z < -2.0) return(42)
+    if (z1 < -2.0 & !is.na(bl_z)) if (bl_z < -2.0) return(41)
+    if (!is.na(th_z)) {
+      if (z1 <  -2.0 & (z1 - th_z) < -1.6) return(43)
+      if (z1 >= -2.0 & z1 < -1.0 & (z1 - th_z) < -2.0) return(53)
+    }
+    if (z1 < -2.0 & is.na(bw)) return(13)
+    if (z1 < -2.0 & is.na(bl)) return(12)
+    if (z1 < -1.0 & is.na(th_z)) return(14)
+
+    # tall
+    if (z1 > 2.5) return(47)
+    if (!is.na(th_z)) if (z1 > 2.0 & (z1 - th_z) > 2.0) return(46)
+    if (z1 > 2.0 & is.na(th_z)) return(82)
+    if (z1 > 1.0 &
+        ((age1 < 8.0 & sex == "F") | (age1 < 9.0 & sex == "M"))) return(79)
+    if (z1 > 1.0) return(81)
+
+    if (is.na(z0)) return(11)
+  }
+
+  if (age1 >= 10.0 & age1 < 18.0) {
+    # short
+    if (z1 < -2.5) return(44)
+
+    # tall
+    if (z1 > 2.5) return(47)
+    if (z1 > 2.0 & sex == "F" & y1 >= 170) return(71)
+    if (z1 > 2.0 & sex == "M" & y1 >= 185) return(72)
+    if (z1 > 2.0) return(73)
+    if (sex == "F" & y1 >= 170.0) return(74)
+    if (sex == "M" & y1 >= 185.0) return(75)
+  }
+
+  # check for gain z1 - z0
+  if (!is.na(z0)) {
+    if (age1 < 3.0) {
+      # short
+      if (z1 < -2.5 & z0 < -2.5 & is.na(bw)) return(13)
+      if (z1 < -2.5 & z0 < -2.5 & bw >= 2500) return(49)
+
+      # tall
+      if (z1 > 2.5 & z0 > 2.5) return(50)
+    }
+
+    if (age1 >= 3.0 & age1 < 10.0) {
+      # short
+      if ((z1 - z0) < -2.0) return(55)
+      if (!is.na(th_z))
+        if (z1 >= -2.0 & (z1 - z0) < -1.0 & (z1 - th_z) < -1.0) return(76)
+
+      # tall
+      if ((z1 - z0) > 2.0) return(54)
+    }
+  }
+
+  # signal everthing is OK
+  return(31)
+}
diff --git a/R/calculate_birth_z.R b/R/calculate_birth_z.R
new file mode 100644
index 0000000..9a434aa
--- /dev/null
+++ b/R/calculate_birth_z.R
@@ -0,0 +1,43 @@
+#' Calculate birth weight SDS relative to Dutch references
+#'
+#' @param y     Birth weight (grammes) or birth length (cm).
+#'              May be a vector. Converted to numeric.
+#' @param sex   Character, either \code{"M"} (male) or \code{"F"} (female)
+#' @param ga    Gestational age, completed week (Integer or character)
+#' @param yname Either \code{"wgt"} (for birth weight) or \code{"hgt"}
+#'              (for birth length)
+#' @param dec   Number of decimals for rounding
+#' @return      Numeric vector of \code{length(bw)} elements with
+#'              standard deviation scores relative to Dutch birth
+#'              weight references
+#' @author      Stef van Buuren, 2019
+#' @examples
+#' calculate_birth_z(c(2500, 3000), sex = "M", ga = 36)
+#' @export
+calculate_birth_z <- function(y, sex, ga, yname = "wgt",
+                              dec = 3) {
+  # convert inputs
+  y <- suppressWarnings(as.numeric(y))
+  sex <- as.character(sex[1L])
+  ga <- suppressWarnings(as.integer(ga[1L]))
+  if (is.na(ga)) return(rep(NA, length(y)))
+
+  # find reference
+  ref <- NULL
+  if (yname == "wgt") {
+    ref <- growthscreener::ref.nl5defSGAgewicht
+    y <- y / 1000
+  }
+  if (yname == "hgt") ref <- growthscreener::ref.nl5defSGAlengte
+  if (is.null(ref)) return(rep(NA, length(y)))
+
+  # find proper row
+  idx <- ref$sex == sex & ref$ga == ga
+  if (any(is.na(idx))) return(rep(NA, length(y)))
+  if (!any(idx)) return(rep(NA, length(y)))
+
+  # calculate
+  z <- round(((y / ref$M[idx])^ref$L[idx] - 1) / (ref$L[idx] * ref$S[idx]), dec)
+  z[is.nan(z)] <- NA
+  z
+}
diff --git a/R/calculate_helpers_hgt.R b/R/calculate_helpers_hgt.R
new file mode 100644
index 0000000..1a04993
--- /dev/null
+++ b/R/calculate_helpers_hgt.R
@@ -0,0 +1,39 @@
+#' @details \code{calculate_helpers_hgt()} provides an optional
+#' pre-calculation for \code{calculate_advice_hgt()}. The user may
+#' wish to divide up calculations into two steps if intermediate
+#' results are needed.
+#' @rdname advice_hgt
+#' @inheritParams calculate_advice_hgt
+#' @return \code{calculate_helpers_hgt()} returns a \code{list} with
+#' the following elements:
+#' \describe{
+#' \item{\code{bw_z}}{Birth weight SDS}
+#' \item{\code{bl_z}}{Birth length SDS}
+#' \item{\code{th}}{Target height (cm)}
+#' \item{\code{th_z}}{Target height SDS}
+#' \item{\code{age1}}{Age at last measurement}
+#' \item{\code{age0}}{Age at previous measurement}
+#' \item{\code{z1}}{Height SDS at \code{age1}}
+#' \item{\code{z0}}{Height SDS at \code{age0}}
+#' }
+#' @export
+calculate_helpers_hgt <- function(sex = NA_character_, dob = as.Date(NA),
+                                  bw = NA, bl = NA, ga = NA,
+                                  etn = NA_character_, hgtf = NA, hgtm = NA,
+                                  dom1 = as.Date(NA), y1 = NA,
+                                  dom0 = as.Date(NA), y0 = NA) {
+  bw_z <- calculate_birth_z(bw, sex, ga, yname = "wgt")
+  bl_z <- calculate_birth_z(bl, sex, ga, yname = "hgt")
+  thl  <- calculate_th(hgtf, hgtm, sex = sex, etn = etn)
+  th   <- thl[1L]
+  th_z <- thl[2L]
+  age1 <- as.integer(dom1 - dob)/365.25
+  age0 <- as.integer(dom0 - dob)/365.25
+  z1 <- y2z(y = y1, x = age1, sex = sex, sub = etn,
+            ref = growthscreener::nl5.hgt)
+  z0 <- y2z(y = y0, x = age0, sex = sex, sub = etn,
+            ref = growthscreener::nl5.hgt)
+
+  list(bw_z = bw_z, bl_z = bl_z, th = th, th_z = th_z,
+       age1 = age1, age0 = age0, z1 = z1, z0 = z0)
+}
diff --git a/R/calculate_th.R b/R/calculate_th.R
new file mode 100644
index 0000000..ed8be2b
--- /dev/null
+++ b/R/calculate_th.R
@@ -0,0 +1,60 @@
+#' Calculate target height SDS
+#'
+#' @param hgtf  Length of biological father (cm)
+#' @param hgtm  Length of biological mother (cm)
+#' @param sex   Character, either \code{"M"} (male) or \code{"F"} (female)
+#' @param etn   Etnicity, one of \code{"N"} (dutch), \code{"T"} (turkish),
+#'              \code{"M"} (moroccan) or \code{"H"} (hindustani).
+#' @param dec   Integer vector, length 2, indicating rounding for
+#'              th and th_z, respectively
+#' @return      Numeric, length 2: target height (cm) and target height
+#'              standard deviation score (z), relative to populations
+#'              living in The Netherlands
+#' @author      Stef van Buuren, 2019
+#' @examples
+#' calculate_th(180, 170, "M", "N")
+#' @export
+calculate_th <- function(hgtf, hgtm,
+                         sex = NULL, etn = NULL,
+                         dec = c(1L, 3L)) {
+
+  # don't calculate if we're missing proper sex or etn
+  if (is.null(sex) | is.null(etn)) return(c(NA, NA))
+  if (!sex %in% c("M", "F")) return(c(NA, NA))
+  if (!etn %in% c("N", "T", "M", "H")) return(c(NA, NA))
+
+  if (sex == "M") {
+    th <- switch(EXPR = etn,
+                 "N" = 44.5 + 0.376 * hgtf + 0.411 * hgtm,
+                 "T" = 29.6 + 0.441 * hgtf + 0.465 * hgtm,
+                 "M" = 22.4 + 0.439 * hgtf + 0.508 * hgtm,
+                 "H" = 43.6 + 0.366 * hgtf + 0.431 * hgtm,
+                 NA)
+    th_z <- switch(EXPR = etn,
+                   "N" = (th - 183.8) / 7.1,
+                   "T" = (th - 176.8) / 6.8,
+                   "M" = (th - 177.2) / 7.7,
+                   "H" = (th - 174.3) / 7.0,
+                   NA)
+  }
+
+  if (sex == "F") {
+    th <- switch(EXPR = etn,
+                 "N" = 47.1 + 0.334 * hgtf + 0.364 * hgtm,
+                 "T" = 32.8 + 0.389 * hgtf + 0.410 * hgtm,
+                 "M" = 32.1 + 0.370 * hgtf + 0.429 * hgtm,
+                 "H" = 49.4 + 0.308 * hgtf + 0.364 * hgtm,
+                 NA)
+    th_z <- switch(EXPR = etn,
+                   "N" = (th - 170.7) / 6.3,
+                   "T" = (th - 162.6) / 6.0,
+                   "M" = (th - 162.8) / 6.5,
+                   "H" = (th - 159.6) / 5.9,
+                   NA)
+  }
+
+  th <- round(th, dec[1])
+  th_z <- round(th_z, dec[2])
+
+  c(th, th_z)
+}
diff --git a/R/datasets.R b/R/datasets.R
new file mode 100644
index 0000000..b2d070c
--- /dev/null
+++ b/R/datasets.R
@@ -0,0 +1,15 @@
+#' Table of Dutch birth weight references
+#' @name datasets
+"ref.nl5defSGAgewicht"
+
+#' Table of Dutch birth length references
+#' @rdname datasets
+"ref.nl5defSGAlengte"
+
+#' Table of Dutch nl5 references
+#' @rdname datasets
+"nl5.hgt"
+
+#' Table of return message and return code
+#' @rdname datasets
+"messages_hgt"
diff --git a/R/fold.R b/R/fold.R
new file mode 100644
index 0000000..016dbc4
--- /dev/null
+++ b/R/fold.R
@@ -0,0 +1,8 @@
+#' Fold long string across multiple lines
+#'
+#' @param s string
+#' @seealso https://stackoverflow.com/questions/2351744/insert-line-breaks-in-long-string-word-wrap
+#' @export
+fold <- function(s) {
+  gsub('(.{1,80})(\\s|$)', '\\1\n', s)
+}
\ No newline at end of file
diff --git a/R/import.R b/R/import.R
new file mode 100644
index 0000000..dec48cf
--- /dev/null
+++ b/R/import.R
@@ -0,0 +1,2 @@
+#' @importFrom AGD            y2z
+NULL
diff --git a/R/msg.R b/R/msg.R
new file mode 100644
index 0000000..adf61ea
--- /dev/null
+++ b/R/msg.R
@@ -0,0 +1,11 @@
+#' Find message string
+#'
+#' @param msgcode Integer with message code
+#' @return A string with the message code
+#' @examples
+#' msg(31)
+#' @export
+msg <- function(msgcode) {
+  mess <- growthscreener::messages_hgt
+  mess[mess[, "msgcode"] == msgcode, "msg"]
+}
diff --git a/README.Rmd b/README.Rmd
new file mode 100644
index 0000000..f3bffe9
--- /dev/null
+++ b/README.Rmd
@@ -0,0 +1,70 @@
+---
+output: github_document
+---
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.path = "man/figures/README-",
+  out.width = "100%"
+)
+options(width = 80)
+```
+# growthscreener
+
+<!-- badges: start -->
+<!-- badges: end -->
+
+The `growthscreener` package implements tools to evaluate child 
+growth with respect to Dutch criteria for unusual growth. 
+Application of these tools helps to identify children that meet 
+criteria for criteria for referral from youth health care (JGZ) for
+follow-up with a general physician or pediatrician.
+
+Dutch guidelines for unusual height growth are currently implemented.
+
+## Installation
+
+The following statements will install the `growthscreener` package
+
+```{r eval = FALSE}
+install.packages("remotes")
+remotes::install_github("stefvanbuuren/growthscreener")
+```
+
+## Example
+
+Find the advice for a very short girl:
+
+```{r example}
+library(growthscreener)
+
+# a very short girl, 4 months old
+msgcode <- calculate_advice_hgt(sex = "F", etn = "N", 
+                                bw = 3250, ga = 40,
+                                dob = as.Date("2018-07-31"),
+                                dom1 = as.Date("2018-12-12"), 
+                                y1 = 55)
+msgcode
+cat(fold(msg(msgcode)))
+
+# some more details
+d <- calculate_helpers_hgt(sex = "F", etn = "N", 
+                           bw = 3250, ga = 40,
+                           dob = as.Date("2018-07-31"),
+                           dom1 = as.Date("2018-12-12"), 
+                           y1 = 55)
+d
+```
+
+The height SDS at the age of 4 months is equal to -3.26, which is 
+the reason for referral. There are 40 different messages for height.
+
+## Background
+
+The calculations follow to the "JGZ-Richtlijn Lengtegroei 2019". See
+<https://www.ncj.nl/richtlijnen/alle-richtlijnen/richtlijn/lengtegroei-2019>
+for more details.
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..ba2a142
--- /dev/null
+++ b/README.md
@@ -0,0 +1,86 @@
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+# growthscreener
+
+<!-- badges: start -->
+
+<!-- badges: end -->
+
+The `growthscreener` package implements tools to evaluate child growth
+with respect to Dutch criteria for unusual growth. Application of these
+tools helps to identify children that meet criteria for criteria for
+referral from youth health care (JGZ) for follow-up with a general
+physician or pediatrician.
+
+Dutch guidelines for unusual height growth are currently implemented.
+
+## Installation
+
+The following statements will install the `growthscreener` package
+
+``` r
+install.packages("remotes")
+remotes::install_github("stefvanbuuren/growthscreener")
+```
+
+## Example
+
+Find the advice for a very short girl:
+
+``` r
+library(growthscreener)
+
+# a very short girl, 4 months old
+msgcode <- calculate_advice_hgt(sex = "F", etn = "N", 
+                                bw = 3250, ga = 40,
+                                dob = as.Date("2018-07-31"),
+                                dom1 = as.Date("2018-12-12"), 
+                                y1 = 55)
+msgcode
+#> [1] 45
+cat(fold(msg(msgcode)))
+#> Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar
+#> huisarts/kinderarts, omdat de lengte < -3 SDS is en het geboortegewicht >= 2500
+#> gram is.
+
+# some more details
+d <- calculate_helpers_hgt(sex = "F", etn = "N", 
+                           bw = 3250, ga = 40,
+                           dob = as.Date("2018-07-31"),
+                           dom1 = as.Date("2018-12-12"), 
+                           y1 = 55)
+d
+#> $bw_z
+#> [1] -0.474
+#> 
+#> $bl_z
+#> [1] NA
+#> 
+#> $th
+#> [1] NA
+#> 
+#> $th_z
+#> [1] NA
+#> 
+#> $age1
+#> [1] 0.367
+#> 
+#> $age0
+#> [1] NA
+#> 
+#> $z1
+#> [1] -3.26
+#> 
+#> $z0
+#> [1] NA
+```
+
+The height SDS at the age of 4 months is equal to -3.26, which is the
+reason for referral. There are 40 different messages for height.
+
+## Background
+
+The calculations follow to the “JGZ-Richtlijn Lengtegroei 2019”. See
+<https://www.ncj.nl/richtlijnen/alle-richtlijnen/richtlijn/lengtegroei-2019>
+for more details.
diff --git a/data-raw/R/create_ref.nl5def.R b/data-raw/R/create_ref.nl5def.R
new file mode 100644
index 0000000..28832ac
--- /dev/null
+++ b/data-raw/R/create_ref.nl5def.R
@@ -0,0 +1,175 @@
+#' Create table of Dutch height reference 2010
+#' @export
+create_ref.nl5def <- function() {
+  ref.nl5def<-data.frame(matrix(NA,570,7))
+  ref.nl5def[,1]<-rep("nl5",570)
+  ref.nl5def[,2]<-c(rep("N",190),rep("T",190),rep("M",190))
+  ref.nl5def[,3]<-c(rep("M",95),rep("F",95),rep("M",95),rep("F",95),rep("M",95),rep("F",95))
+  ref.nl5def[,4]<-c(0.0000,0.0027,0.0055,0.0082,0.0110,0.0137,0.0164,0.0192,0.0219,0.0246,
+                    0.0274,0.0301,0.0329,0.0356,0.0383,0.0575,0.0767,0.0958,0.1150,0.1342,
+                    0.1533,0.1725,0.1916,0.2108,0.2300,0.2491,0.2500,0.2917,0.3333,0.3750,
+                    0.4167,0.4583,0.5000,0.5417,0.5833,0.6250,0.6667,0.7083,0.7500,0.7917,
+                    0.8333,0.8750,0.9167,0.9583,1.0000,1.0833,1.1667,1.2500,1.3333,1.4167,
+                    1.5000,1.5833,1.6667,1.7500,1.8333,1.9167,2.0000,2.5000,3.0000,3.5000,
+                    4.0000,4.5000,5.0000,5.5000,6.0000,6.5000,7.0000,7.5000,8.0000,8.5000,
+                    9.0000,9.5000,10.0000,10.5000,11.0000,11.5000,12.0000,12.5000,13.0000,13.5000,
+                    14.0000,14.5000,15.0000,15.5000,16.0000,16.5000,17.0000,17.5000,18.0000,18.5000,
+                    19.0000,19.5000,20.0000,20.5000,21.0000,0.0000,0.0027,0.0055,0.0082,0.0110,
+                    0.0137,0.0164,0.0192,0.0219,0.0246,0.0274,0.0301,0.0329,0.0356,0.0383,
+                    0.0575,0.0767,0.0958,0.1150,0.1342,0.1533,0.1725,0.1916,0.2108,0.2300,
+                    0.2491,0.2500,0.2917,0.3333,0.3750,0.4167,0.4583,0.5000,0.5417,0.5833,
+                    0.6250,0.6667,0.7083,0.7500,0.7917,0.8333,0.8750,0.9167,0.9583,1.0000,
+                    1.0833,1.1667,1.2500,1.3333,1.4167,1.5000,1.5833,1.6667,1.7500,1.8333,
+                    1.9167,2.0000,2.5000,3.0000,3.5000,4.0000,4.5000,5.0000,5.5000,6.0000,
+                    6.5000,7.0000,7.5000,8.0000,8.5000,9.0000,9.5000,10.0000,10.5000,11.0000,
+                    11.5000,12.0000,12.5000,13.0000,13.5000,14.0000,14.5000,15.0000,15.5000,16.0000,
+                    16.5000,17.0000,17.5000,18.0000,18.5000,19.0000,19.5000,20.0000,20.5000,21.0000,
+                    0.0000,0.0027,0.0055,0.0082,0.0110,0.0137,0.0164,0.0192,0.0219,0.0246,
+                    0.0274,0.0301,0.0329,0.0356,0.0383,0.0575,0.0767,0.0958,0.1150,0.1342,
+                    0.1533,0.1725,0.1916,0.2108,0.2300,0.2491,0.2500,0.2917,0.3333,0.3750,
+                    0.4167,0.4583,0.5000,0.5417,0.5833,0.6250,0.6667,0.7083,0.7500,0.7917,
+                    0.8333,0.8750,0.9167,0.9583,1.0000,1.0833,1.1667,1.2500,1.3333,1.4167,
+                    1.5000,1.5833,1.6667,1.7500,1.8333,1.9167,2.0000,2.5000,3.0000,3.5000,
+                    4.0000,4.5000,5.0000,5.5000,6.0000,6.5000,7.0000,7.5000,8.0000,8.5000,
+                    9.0000,9.5000,10.0000,10.5000,11.0000,11.5000,12.0000,12.5000,13.0000,13.5000,
+                    14.0000,14.5000,15.0000,15.5000,16.0000,16.5000,17.0000,17.5000,18.0000,18.5000,
+                    19.0000,19.5000,20.0000,20.5000,21.0000,0.0000,0.0027,0.0055,0.0082,0.0110,
+                    0.0137,0.0164,0.0192,0.0219,0.0246,0.0274,0.0301,0.0329,0.0356,0.0383,
+                    0.0575,0.0767,0.0958,0.1150,0.1342,0.1533,0.1725,0.1916,0.2108,0.2300,
+                    0.2491,0.2500,0.2917,0.3333,0.3750,0.4167,0.4583,0.5000,0.5417,0.5833,
+                    0.6250,0.6667,0.7083,0.7500,0.7917,0.8333,0.8750,0.9167,0.9583,1.0000,
+                    1.0833,1.1667,1.2500,1.3333,1.4167,1.5000,1.5833,1.6667,1.7500,1.8333,
+                    1.9167,2.0000,2.5000,3.0000,3.5000,4.0000,4.5000,5.0000,5.5000,6.0000,
+                    6.5000,7.0000,7.5000,8.0000,8.5000,9.0000,9.5000,10.0000,10.5000,11.0000,
+                    11.5000,12.0000,12.5000,13.0000,13.5000,14.0000,14.5000,15.0000,15.5000,16.0000,
+                    16.5000,17.0000,17.5000,18.0000,18.5000,19.0000,19.5000,20.0000,20.5000,21.0000,
+                    0.0000,0.0027,0.0055,0.0082,0.0110,0.0137,0.0164,0.0192,0.0219,0.0246,
+                    0.0274,0.0301,0.0329,0.0356,0.0383,0.0575,0.0767,0.0958,0.1150,0.1342,
+                    0.1533,0.1725,0.1916,0.2108,0.2300,0.2491,0.2500,0.2917,0.3333,0.3750,
+                    0.4167,0.4583,0.5000,0.5417,0.5833,0.6250,0.6667,0.7083,0.7500,0.7917,
+                    0.8333,0.8750,0.9167,0.9583,1.0000,1.0833,1.1667,1.2500,1.3333,1.4167,
+                    1.5000,1.5833,1.6667,1.7500,1.8333,1.9167,2.0000,2.5000,3.0000,3.5000,
+                    4.0000,4.5000,5.0000,5.5000,6.0000,6.5000,7.0000,7.5000,8.0000,8.5000,
+                    9.0000,9.5000,10.0000,10.5000,11.0000,11.5000,12.0000,12.5000,13.0000,13.5000,
+                    14.0000,14.5000,15.0000,15.5000,16.0000,16.5000,17.0000,17.5000,18.0000,18.5000,
+                    19.0000,19.5000,20.0000,20.5000,21.0000,0.0000,0.0027,0.0055,0.0082,0.0110,
+                    0.0137,0.0164,0.0192,0.0219,0.0246,0.0274,0.0301,0.0329,0.0356,0.0383,
+                    0.0575,0.0767,0.0958,0.1150,0.1342,0.1533,0.1725,0.1916,0.2108,0.2300,
+                    0.2491,0.2500,0.2917,0.3333,0.3750,0.4167,0.4583,0.5000,0.5417,0.5833,
+                    0.6250,0.6667,0.7083,0.7500,0.7917,0.8333,0.8750,0.9167,0.9583,1.0000,
+                    1.0833,1.1667,1.2500,1.3333,1.4167,1.5000,1.5833,1.6667,1.7500,1.8333,
+                    1.9167,2.0000,2.5000,3.0000,3.5000,4.0000,4.5000,5.0000,5.5000,6.0000,
+                    6.5000,7.0000,7.5000,8.0000,8.5000,9.0000,9.5000,10.0000,10.5000,11.0000,
+                    11.5000,12.0000,12.5000,13.0000,13.5000,14.0000,14.5000,15.0000,15.5000,16.0000,
+                    16.5000,17.0000,17.5000,18.0000,18.5000,19.0000,19.5000,20.0000,20.5000,21.0000)
+  ref.nl5def[,5]<-rep(1,570)
+  ref.nl5def[,6]<-c(NA,NA,NA,NA,NA,NA,NA,NA,52.33,52.44,52.55,
+                    52.65,52.76,52.86,52.96,53.70,54.44,55.17,55.91,56.64,57.36,58.08,
+                    58.79,59.49,60.17,60.84,60.87,62.27,63.59,64.83,65.98,67.04,68.03,
+                    68.96,69.82,70.64,71.41,72.15,72.87,73.55,74.21,74.85,75.47,76.07,
+                    76.65,77.78,78.85,79.89,80.89,81.87,82.83,83.78,84.71,85.64,86.55,
+                    87.47,88.37,93.48,97.83,101.66,105.52,109.46,113.15,116.56,119.89,123.02,
+                    126.18,129.38,132.52,135.57,138.52,141.24,143.71,146.20,148.95,151.98,155.24,
+                    158.54,161.76,165.06,168.53,172.05,175.24,177.59,179.10,180.16,180.99,181.70,
+                    182.43,183.13,183.64,183.83,183.83,183.83,183.83,NA,NA,NA,NA,
+                    51.24,51.35,51.45,51.55,51.66,51.76,51.87,51.97,52.07,52.18,52.28,
+                    53.00,53.71,54.41,55.10,55.78,56.45,57.12,57.77,58.41,59.03,59.64,
+                    59.67,60.95,62.17,63.31,64.40,65.42,66.38,67.30,68.16,68.98,69.76,
+                    70.49,71.21,71.89,72.55,73.20,73.82,74.43,75.03,76.19,77.32,78.41,
+                    79.46,80.49,81.48,82.46,83.41,84.36,85.28,86.19,87.08,92.15,96.97,
+                    101.16,104.94,108.58,112.12,115.52,118.81,122.07,125.26,128.34,131.34,134.33,
+                    137.32,140.36,143.46,146.60,149.74,152.80,155.71,158.39,160.77,162.80,164.47,
+                    165.82,166.88,167.70,168.32,168.80,169.17,169.46,169.69,169.90,170.10,170.28,
+                    170.45,170.61,170.75,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,54.21,54.88,55.54,56.21,
+                    56.87,57.54,58.20,58.85,59.51,60.16,60.79,60.82,62.18,63.47,64.71,
+                    65.89,66.99,68.04,69.02,69.95,70.83,71.67,72.46,73.23,73.97,74.68,
+                    75.36,76.02,76.66,77.29,78.48,79.61,80.66,81.66,82.64,83.58,84.48,
+                    85.35,86.20,87.03,87.84,88.64,93.19,97.38,101.50,105.59,109.35,112.64,
+                    115.57,118.33,120.96,123.65,126.50,129.48,132.44,135.30,138.02,140.61,143.20,
+                    145.92,148.86,151.93,155.09,158.25,161.25,163.97,166.46,168.73,170.84,172.70,
+                    174.21,175.36,176.22,176.77,177.05,177.18,177.29,177.43,177.54,177.53,NA,
+                    NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,53.49,54.13,54.77,55.41,56.05,56.68,57.31,
+                    57.93,58.55,59.15,59.18,60.47,61.69,62.87,63.99,65.07,66.09,67.06,
+                    67.99,68.87,69.72,70.54,71.33,72.09,72.83,73.55,74.26,74.95,75.61,
+                    76.91,78.15,79.33,80.45,81.53,82.55,83.52,84.44,85.33,86.18,87.00,
+                    87.81,92.24,96.29,100.21,103.95,107.58,110.92,113.99,116.96,119.93,122.79,
+                    125.48,128.14,130.93,133.89,137.10,140.54,144.01,147.47,150.61,153.18,155.23,
+                    156.88,158.26,159.41,160.30,160.94,161.43,161.88,162.31,162.56,162.59,162.59,
+                    162.59,162.59,162.59,162.59,162.59,162.59,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    54.43,55.08,55.73,56.38,57.03,57.67,58.30,58.94,59.57,60.18,60.21,
+                    61.52,62.77,63.95,65.06,66.11,67.10,68.04,68.92,69.77,70.58,71.36,
+                    72.11,72.83,73.53,74.20,74.86,75.50,76.11,77.29,78.40,79.45,80.45,
+                    81.43,82.38,83.31,84.22,85.12,86.00,86.86,87.72,92.46,96.77,100.79,
+                    104.45,107.96,111.37,114.56,117.73,120.86,123.91,126.76,129.41,131.92,134.55,
+                    137.35,140.07,142.75,145.38,147.92,150.49,153.26,156.27,159.66,163.30,166.81,
+                    169.89,172.30,174.04,175.25,176.12,176.74,177.17,177.44,177.59,177.70,177.78,
+                    177.83,177.83,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,52.92,53.61,54.29,54.97,
+                    55.64,56.31,56.97,57.62,58.26,58.88,58.91,60.23,61.47,62.65,63.77,
+                    64.83,65.84,66.80,67.72,68.61,69.46,70.27,71.06,71.80,72.50,73.17,
+                    73.81,74.41,74.99,76.06,77.11,78.09,79.03,79.98,80.92,81.84,82.76,
+                    83.69,84.62,85.54,86.47,91.63,96.00,99.83,103.48,106.91,110.15,113.47,
+                    116.76,119.91,122.95,125.76,128.45,131.23,134.13,137.09,140.23,143.59,146.98,
+                    150.11,152.90,155.34,157.35,159.01,160.45,161.63,162.42,162.78,162.78,162.78,
+                    162.78,162.78,162.78,162.78,162.78,162.78,162.78,162.78,162.78)
+  ref.nl5def[,7]<-c(NA,NA,NA,NA,NA,NA,NA,NA,0.0443,0.0442,0.0442,
+                    0.0441,0.0440,0.0440,0.0439,0.0434,0.0430,0.0425,0.0421,0.0417,0.0413,0.0409,
+                    0.0405,0.0402,0.0398,0.0395,0.0395,0.0389,0.0383,0.0379,0.0375,0.0371,0.0368,
+                    0.0366,0.0364,0.0362,0.0361,0.0360,0.0359,0.0358,0.0357,0.0356,0.0356,0.0356,
+                    0.0355,0.0355,0.0355,0.0355,0.0355,0.0355,0.0356,0.0356,0.0357,0.0358,0.0359,
+                    0.0360,0.0361,0.0369,0.0379,0.0389,0.0398,0.0405,0.0411,0.0416,0.0422,0.0427,
+                    0.0431,0.0435,0.0439,0.0443,0.0447,0.0452,0.0457,0.0463,0.0472,0.0481,0.0492,
+                    0.0501,0.0506,0.0504,0.0494,0.0475,0.0452,0.0434,0.0422,0.0414,0.0407,0.0402,
+                    0.0396,0.0391,0.0387,0.0386,0.0386,0.0386,0.0386,NA,NA,NA,NA,
+                    0.0451,0.0451,0.0450,0.0449,0.0449,0.0448,0.0447,0.0447,0.0446,0.0446,0.0445,
+                    0.0441,0.0437,0.0433,0.0429,0.0426,0.0422,0.0419,0.0416,0.0413,0.0410,0.0407,
+                    0.0407,0.0402,0.0397,0.0393,0.0389,0.0386,0.0384,0.0382,0.0380,0.0378,0.0377,
+                    0.0376,0.0375,0.0374,0.0374,0.0373,0.0373,0.0373,0.0373,0.0373,0.0373,0.0373,
+                    0.0373,0.0374,0.0375,0.0375,0.0376,0.0376,0.0377,0.0378,0.0378,0.0382,0.0387,
+                    0.0393,0.0401,0.0409,0.0416,0.0421,0.0426,0.0430,0.0433,0.0437,0.0441,0.0444,
+                    0.0447,0.0449,0.0450,0.0450,0.0449,0.0446,0.0442,0.0436,0.0430,0.0423,0.0417,
+                    0.0411,0.0405,0.0400,0.0396,0.0392,0.0388,0.0385,0.0382,0.0379,0.0377,0.0374,
+                    0.0372,0.0369,0.0367,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,0.0374,0.0373,0.0371,0.0369,
+                    0.0367,0.0366,0.0364,0.0363,0.0361,0.0360,0.0359,0.0359,0.0356,0.0354,0.0352,
+                    0.0351,0.0350,0.0349,0.0348,0.0348,0.0348,0.0348,0.0348,0.0348,0.0348,0.0348,
+                    0.0348,0.0349,0.0349,0.0349,0.0350,0.0351,0.0352,0.0352,0.0353,0.0353,0.0353,
+                    0.0354,0.0354,0.0354,0.0354,0.0354,0.0355,0.0358,0.0362,0.0370,0.0379,0.0388,
+                    0.0397,0.0405,0.0414,0.0424,0.0435,0.0446,0.0457,0.0467,0.0475,0.0482,0.0488,
+                    0.0492,0.0496,0.0497,0.0496,0.0491,0.0482,0.0470,0.0456,0.0442,0.0428,0.0415,
+                    0.0404,0.0396,0.0390,0.0387,0.0385,0.0384,0.0383,0.0382,0.0381,0.0381,NA,
+                    NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,0.0433,0.0429,0.0424,0.0420,0.0415,0.0411,0.0407,
+                    0.0403,0.0400,0.0396,0.0396,0.0388,0.0382,0.0376,0.0371,0.0366,0.0361,0.0358,
+                    0.0354,0.0351,0.0349,0.0346,0.0344,0.0343,0.0341,0.0339,0.0338,0.0337,0.0336,
+                    0.0334,0.0333,0.0332,0.0332,0.0332,0.0332,0.0332,0.0333,0.0334,0.0335,0.0336,
+                    0.0337,0.0346,0.0355,0.0363,0.0371,0.0377,0.0384,0.0390,0.0395,0.0401,0.0407,
+                    0.0413,0.0419,0.0425,0.0432,0.0437,0.0441,0.0441,0.0437,0.0429,0.0419,0.0409,
+                    0.0401,0.0393,0.0386,0.0381,0.0378,0.0375,0.0372,0.0370,0.0369,0.0369,0.0369,
+                    0.0369,0.0369,0.0369,0.0369,0.0369,0.0369,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    0.0368,0.0367,0.0366,0.0365,0.0363,0.0362,0.0361,0.0360,0.0360,0.0359,0.0359,
+                    0.0357,0.0356,0.0356,0.0355,0.0355,0.0356,0.0356,0.0357,0.0358,0.0359,0.0359,
+                    0.0360,0.0362,0.0363,0.0364,0.0365,0.0366,0.0367,0.0369,0.0370,0.0372,0.0373,
+                    0.0375,0.0376,0.0377,0.0378,0.0379,0.0380,0.0381,0.0381,0.0384,0.0384,0.0385,
+                    0.0386,0.0389,0.0393,0.0398,0.0404,0.0412,0.0421,0.0430,0.0438,0.0446,0.0455,
+                    0.0463,0.0473,0.0483,0.0493,0.0503,0.0511,0.0518,0.0521,0.0518,0.0508,0.0495,
+                    0.0479,0.0465,0.0454,0.0446,0.0441,0.0437,0.0434,0.0433,0.0432,0.0431,0.0430,
+                    0.0430,0.0430,NA,NA,NA,NA,NA,NA,NA,NA,NA,
+                    NA,NA,NA,NA,NA,NA,NA,0.0385,0.0382,0.0380,0.0377,
+                    0.0375,0.0372,0.0370,0.0368,0.0366,0.0364,0.0364,0.0361,0.0358,0.0355,0.0353,
+                    0.0352,0.0350,0.0349,0.0349,0.0348,0.0348,0.0347,0.0347,0.0347,0.0348,0.0348,
+                    0.0348,0.0349,0.0349,0.0351,0.0352,0.0354,0.0355,0.0356,0.0358,0.0359,0.0360,
+                    0.0361,0.0362,0.0362,0.0363,0.0368,0.0372,0.0375,0.0377,0.0377,0.0379,0.0382,
+                    0.0386,0.0391,0.0400,0.0410,0.0424,0.0439,0.0457,0.0475,0.0491,0.0501,0.0500,
+                    0.0490,0.0475,0.0458,0.0442,0.0429,0.0417,0.0408,0.0402,0.0399,0.0399,0.0399,
+                    0.0399,0.0399,0.0399,0.0399,0.0399,0.0399,0.0399,0.0399,0.0399)
+
+  dimnames(ref.nl5def)[[2]]<-c("pop","sub","sex","x","L","M","S")
+
+  ref.nl5def$L<-ifelse(is.na(ref.nl5def$M),NA,ref.nl5def$L)
+  ref.nl5def
+}
diff --git a/data-raw/R/create_ref.nl5defSGAgewicht.R b/data-raw/R/create_ref.nl5defSGAgewicht.R
new file mode 100644
index 0000000..d0512d9
--- /dev/null
+++ b/data-raw/R/create_ref.nl5defSGAgewicht.R
@@ -0,0 +1,26 @@
+#' Create table of Dutch birth weight references
+#' @export
+create_ref.nl5defSGAgewicht <- function() {
+  x <- data.frame(matrix(NA, 36L, 8L))
+  x[, 1L] <- rep("pinkeltje", 36L)
+  x[, 2L] <- c(rep("N", 36L))
+  x[, 3L] <- c(rep(c("M", "F"), 18L))
+  x[, 4L] <- rep(0L, 36L)
+  x[, 5L] <- rep(c(1.0885, 0.8089), 18L)
+  x[, 6L] <- c(0.9620, 0.8770, 0.9850, 0.9010, 1.0420, 0.9530,
+               1.1310, 1.0340, 1.2490, 1.1440, 1.3950, 1.2850,
+               1.5650, 1.4520, 1.7600, 1.6410, 1.9800, 1.8510,
+               2.2220, 2.0800, 2.4790, 2.3240, 2.7450, 2.5740,
+               3.0100, 2.8230, 3.2680, 3.0620, 3.5130, 3.2840,
+               3.7400, 3.4840, 3.9470, 3.6600, 3.9470, 3.6600)
+  x[, 7L] <- c(0.2111, 0.1805, 0.2061, 0.1940, 0.2031, 0.2035,
+               0.2014, 0.2086, 0.2003, 0.2089, 0.1984, 0.2051,
+               0.1952, 0.1984, 0.1912, 0.1900, 0.1867, 0.1808,
+               0.1818, 0.1718, 0.1762, 0.1636, 0.1702, 0.1566,
+               0.1646, 0.1511, 0.1602, 0.1470, 0.1577, 0.1441,
+               0.1575, 0.1425, 0.1599, 0.1418, 0.1599, 0.1418)
+  x[, 8L] <- rep(25:42, each = 2L)
+
+  dimnames(x)[[2L]] <- c("pop", "sub", "sex", "x", "L", "M", "S", "ga")
+  x
+}
diff --git a/data-raw/R/create_ref.nl5defSGAlengte.R b/data-raw/R/create_ref.nl5defSGAlengte.R
new file mode 100644
index 0000000..f58e328
--- /dev/null
+++ b/data-raw/R/create_ref.nl5defSGAlengte.R
@@ -0,0 +1,26 @@
+#' Create table of Dutch birth length references
+#' @export
+create_ref.nl5defSGAlengte <- function() {
+  x <- data.frame(matrix(NA, 36L, 8L))
+  x[, 1L] <- rep("pinkeltje", 36L)
+  x[, 2L] <- c(rep("N", 36L))
+  x[, 3L] <- c(rep(c("M", "F"), 18L))
+  x[, 4L] <- rep(0L, 36L)
+  x[, 5L] <- rep(1, 36L)
+  x[, 6L] <- c(28.49, 28.60, 30.68, 30.66, 32.80, 32.66,
+               34.85, 34.58, 36.82, 36.44, 38.71, 38.22,
+               40.52, 39.93, 42.26, 41.57, 43.91, 43.12,
+               45.47, 44.59, 46.94, 45.97, 48.32, 47.27,
+               49.60, 48.47, 50.79, 49.59, 51.88, 50.62,
+               52.87, 51.56, 53.76, 52.42, 53.76, 52.42)
+  x[, 7L] <- c(.1261, .0933, .1134, .0955, .1025, .0951,
+               .0931, .0924, .0846, .0874, .0769, .0811,
+               .0699, .0740, .0638, .0670, .0586, .0605,
+               .0541, .0549, .0502, .0504, .0468, .0470,
+               .0438, .0447, .0413, .0432, .0395, .0424,
+               .0382, .0422, .0377, .0425, .0377, .0425)
+  x[, 8L] <- rep(25:42, each = 2L)
+
+  dimnames(x)[[2L]] <- c("pop", "sub", "sex", "x", "L", "M", "S", "ga")
+  x
+}
diff --git a/data-raw/R/messages.R b/data-raw/R/messages.R
new file mode 100644
index 0000000..239b56a
--- /dev/null
+++ b/data-raw/R/messages.R
@@ -0,0 +1,7 @@
+path <- path.expand("~/Package/growthscreener/growthscreener")
+
+# Added message 20 by hand on 2019-08-12, so we need to work for
+# this file now.
+messages_hgt <- read.table(file = file.path(path, "data-raw", "data", "messages_hgt.txt"),
+                  sep = "\t", header = TRUE, stringsAsFactors = FALSE)
+usethis::use_data(messages_hgt, overwrite = TRUE)
diff --git a/data-raw/R/save_all_data.R b/data-raw/R/save_all_data.R
new file mode 100644
index 0000000..4973c18
--- /dev/null
+++ b/data-raw/R/save_all_data.R
@@ -0,0 +1,23 @@
+library(readxl)
+library(AGD)
+library(usethis)
+
+path <- path.expand("~/Package/growthscreener/growthscreener/data-raw")
+
+source(file.path(path, "R", "create_ref.nl5def.R"))
+source(file.path(path, "R", "create_ref.nl5defSGAgewicht.R"))
+source(file.path(path, "R", "create_ref.nl5defSGAlengte.R"))
+
+ref.nl5defSGAgewicht <- create_ref.nl5defSGAgewicht()
+ref.nl5defSGAlengte <- create_ref.nl5defSGAlengte()
+
+ref.nl5def <- create_ref.nl5def()
+nl45.hgt <- rbind(nl4.hgt, ref.nl5def)
+nl5.hgt <- nl45.hgt[nl45.hgt[, 1] == "nl5", ]
+
+fn <- file.path(path, "data", "Hindo.xlsx")
+Hindo.hgt <- data.frame(read_excel(fn))
+nl5.hgt <- rbind(nl5.hgt, Hindo.hgt)
+
+usethis::use_data(nl5.hgt, ref.nl5defSGAgewicht, ref.nl5defSGAlengte,
+                  overwrite = TRUE)
diff --git a/data-raw/data/Hindo.xlsx b/data-raw/data/Hindo.xlsx
new file mode 100644
index 0000000..decbc2d
Binary files /dev/null and b/data-raw/data/Hindo.xlsx differ
diff --git a/data-raw/data/messages_hgt.txt b/data-raw/data/messages_hgt.txt
new file mode 100644
index 0000000..f6be1c3
--- /dev/null
+++ b/data-raw/data/messages_hgt.txt
@@ -0,0 +1,41 @@
+msgcode	msg
+11	Het advies kan niet worden bepaald. Voer de vorige lengtemeting in.
+12	Het advies kan niet worden bepaald. Voer eerst de lengte bij geboorte in.
+13	Het advies kan niet worden bepaald. Voer eerst het geboortegewicht in.
+14	Het advies kan niet worden bepaald. Voer zowel de lengte van de moeder als de vader in.
+15	Voer de datum van de huidige meting in.
+16	Voer de geboortedatum in.
+17	Voer de gegevens in.
+18	Voer de lengte bij de huidige meting in.
+19	Voer het geslacht in.
+20	Voer etniciteit in.
+21	De JGZ-richtlijn lengtegroei is bedoeld voor kinderen tot 18 jaar.
+22	Het kind is nog te jong om een advies te kunnen geven.
+31	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: In principe geen verwijzing nodig, naar eigen inzicht handelen.
+41	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de huidige lengte en de lengte bij geboorte onder de -2 SD vallen.
+42	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de huidige lengte en het geboortgegewicht onder de -2 SD vallen.
+43	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte < -2 SDS is en de lengte meer dan 1,6 SD afwijkt ten opzichte van de target height.
+44	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte < -2,5 SD is.
+45	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte < -3 SDS is en het geboortegewicht >= 2500 gram is.
+46	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte > 2 SDS is en de lengte meer dan 2 SD afwijkt ten opzichte van de target height.
+47	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte > 2,5 SD is.
+48	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte > 3 SDS is.
+49	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte herhaaldelijk < -2,5 SDS is en het geboorgewicht >= 2500 gram is.
+50	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte herhaaldelijk > 2,5 SDS is.
+51	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte SDS < -3 en het geboortegewicht >= 2500 gram is.
+52	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte SDS > 3 is.
+53	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat de lengte tussen de -2 en -1 SD valt en de lengte meer dan 2 SD afwijkt ten opzichte van de target height.
+54	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat er sprake is van een groeiversnelling van meer dan 2 SD.
+55	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Verwijzen naar huisarts/kinderarts, omdat er sprake is van een lengteafbuiging van meer dan 2 SD.
+71	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Als er naast een grote lengte (>2 SDS) ook sprake is van een laat beginnende puberteitsontwikkeling en/of een traag voortschrijdende puberteitsontwikkeling, dan Verwijzen naar de huisarts/kinderarts. Bij ongerustheid over de verwachte eindlengte, kan een lange adolescent worden verwezen naar de huisarts/kinderarts als de puberteit nog niet is voltooid bij een lengte van minimaal 170 cm. Indien nee, dan is er in principe geen verwijzing nodig. Naar eigen inzicht handelen.
+72	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Als er naast een grote lengte (>2 SDS) ook sprake is van een laat beginnende puberteitsontwikkeling en/of een traag voortschrijdende puberteitsontwikkeling, dan Verwijzen naar de huisarts/kinderarts. Bij ongerustheid over de verwachte eindlengte, kan een lange adolescent worden verwezen naar de huisarts/kinderarts als de puberteit nog niet is voltooid bij een lengte van minimaal 185 cm. Indien nee, dan is er in principe geen verwijzing nodig. Naar eigen inzicht handelen.
+73	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Als er naast een grote lengte (>2 SDS) ook sprake is van een laat beginnende puberteitsontwikkeling en/of een traag voortschrijdende puberteitsontwikkeling, dan Verwijzen naar de huisarts/kinderarts. Indien nee, dan is er in principe geen verwijzing nodig. Naar eigen inzicht handelen.
+74	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Bij ongerustheid over de verwachte eindlengte, kan een lange adolescent worden verwezen naar de huisarts/kinderarts als de puberteit nog niet is voltooid bij een lengte van minimaal 170 cm.
+75	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Bij ongerustheid over de verwachte eindlengte, kan een lange adolescent worden verwezen naar de huisarts/kinderarts als de puberteit nog niet is voltooid bij een lengte van minimaal 185 cm.
+76	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Is de afbuiging minder dan 1,5 SDS? Herhaal dan de meting na 6 maanden en controleer dan of de lengte verder is afgenomen met meer dan 0,5 SDS. Voer ook het gewicht en de lengte bij geboorte in. Is de afbuiging al meer dan 1,5 SDS gedurende drie metingen en is de lengte in het afgelopen half jaar met 0,5 SD verder afgenomen, dan verwijzen naar huisarts/kinderarts.
+77	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke andere symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte? Met name achterstanden in de ontwikkeling, gedragsproblemen en macrocefalie (groot hoofd). Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, in principe geen verwijzing nodig, naar eigen inzicht handelen.
+78	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte, met name achterstand in motorische en/of spraak/taal achterstand, gedragsproblemen en macrocefalie? Het is ook van groot belang om de puberteitskenmerken te onderzoeken, want pubertas praecox is een van de meest voorkomende oorzaken van groeiversnelling van kinderen <8 (meisjes) of <9 (jongens) jaar. Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, dan is er in principe geen verwijzing nodig. Naar eigen inzicht handelen.
+79	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte, met name achterstand in motorische en/of spraak/taal achterstand, gedragsproblemen en macrocefalie? Het is ook van groot belang om de puberteitskenmerken te onderzoeken, want pubertas praecox is een van de meest voorkomende oorzaken van groeiversnelling van kinderen <8 (meisjes) of <9 (jongens) jaar. Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, voer zowel de lengte van de moeder als de vader in om te bepalen of een verwijzing nodig is.
+80	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte, met name achterstand in motorische en/of spraak/taal achterstand, gedragsproblemen en macrocefalie? Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, dan is er in principe geen verwijzing nodig. Naar eigen inzicht handelen.
+81	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte, met name achterstand in motorische en/of spraak/taal achterstand, gedragsproblemen en macrocefalie? Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, voer de eerdere meting in om te bepalen of een verwijzing nodig is.
+82	Het advies volgens de JGZ-richtlijn lengtegroei is als volgt: Zijn er duidelijke symptomen die samenhangen met aandoeningen die gepaard gaan met een grote lengte, met name achterstand in motorische en/of spraak/taal achterstand, gedragsproblemen en macrocefalie? Indien ja, Verwijzen naar huisarts/kinderarts. Indien nee, voer zowel de lengte van de moeder als de vader in om te bepalen of een verwijzing nodig is.
diff --git a/data/messages_hgt.rda b/data/messages_hgt.rda
new file mode 100644
index 0000000..3a9b09b
Binary files /dev/null and b/data/messages_hgt.rda differ
diff --git a/data/nl5.hgt.rda b/data/nl5.hgt.rda
new file mode 100644
index 0000000..6f9d6ce
Binary files /dev/null and b/data/nl5.hgt.rda differ
diff --git a/data/ref.nl5defSGAgewicht.rda b/data/ref.nl5defSGAgewicht.rda
new file mode 100644
index 0000000..75f433a
Binary files /dev/null and b/data/ref.nl5defSGAgewicht.rda differ
diff --git a/data/ref.nl5defSGAlengte.rda b/data/ref.nl5defSGAlengte.rda
new file mode 100644
index 0000000..181d09e
Binary files /dev/null and b/data/ref.nl5defSGAlengte.rda differ
diff --git a/growthscreener.Rproj b/growthscreener.Rproj
new file mode 100644
index 0000000..1440aa0
--- /dev/null
+++ b/growthscreener.Rproj
@@ -0,0 +1,20 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageRoxygenize: rd,collate,namespace
diff --git a/man/advice_hgt.Rd b/man/advice_hgt.Rd
new file mode 100644
index 0000000..de5989e
--- /dev/null
+++ b/man/advice_hgt.Rd
@@ -0,0 +1,95 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_advice_hgt.R,
+%   R/calculate_helpers_hgt.R
+\name{calculate_advice_hgt}
+\alias{calculate_advice_hgt}
+\alias{calculate_helpers_hgt}
+\title{Referral advice for body height}
+\usage{
+calculate_advice_hgt(sex = NA_character_, dob = as.Date(NA), bw = NA,
+  bl = NA, ga = NA, etn = NA_character_, hgtf = NA, hgtm = NA,
+  dom1 = as.Date(NA), y1 = NA, dom0 = as.Date(NA), y0 = NA,
+  d = NULL)
+
+calculate_helpers_hgt(sex = NA_character_, dob = as.Date(NA),
+  bw = NA, bl = NA, ga = NA, etn = NA_character_, hgtf = NA,
+  hgtm = NA, dom1 = as.Date(NA), y1 = NA, dom0 = as.Date(NA),
+  y0 = NA)
+}
+\arguments{
+\item{sex}{Character, either \code{"M"} (male) or \code{"F"} (female)}
+
+\item{dob}{Date of birth (class Date)}
+
+\item{bw}{Birth weight (grammes)}
+
+\item{bl}{Birth length (cm)}
+
+\item{ga}{Gestational age, completed weeks (Integer or character)}
+
+\item{etn}{Etnicity, one of \code{"N"} (dutch), \code{"T"} (turkish),
+\code{"M"} (moroccan) or \code{"H"} (hindustani).}
+
+\item{hgtf}{Height of father (cm)}
+
+\item{hgtm}{Height of mother (cm)}
+
+\item{dom1}{Date of last measurement (Date)}
+
+\item{y1}{Height at last measurement (cm)}
+
+\item{dom0}{Date of previous measurement (Date)}
+
+\item{y0}{Height at previous measurement (cm)}
+
+\item{d}{Optional, list of derived variables, obtained by
+\code{calculate_derived_variables()}}
+}
+\value{
+\code{calculate_advice_hgt} returns an integer, the \code{msgcode}
+
+\code{calculate_helpers_hgt()} returns a \code{list} with
+the following elements:
+\describe{
+\item{\code{bw_z}}{Birth weight SDS}
+\item{\code{bl_z}}{Birth length SDS}
+\item{\code{th}}{Target height (cm)}
+\item{\code{th_z}}{Target height SDS}
+\item{\code{age1}}{Age at last measurement}
+\item{\code{age0}}{Age at previous measurement}
+\item{\code{z1}}{Height SDS at \code{age1}}
+\item{\code{z0}}{Height SDS at \code{age0}}
+}
+}
+\description{
+This function traverses the decision tree of the
+"JGZ-Richtlijn Lengtegroei 2019".
+}
+\details{
+The decision tree assesses both single and paired measurements.
+The last observations (\code{y1}) is generally taken as the
+last measurement, whereas \code{y0} can be one of the previous
+measurements. For more than two measurements, there are many
+pairs possible, and these pairs need not be consecutive.
+The \code{y0} measurement needs to be defined by the user,
+and is informally taken as an earlier measurement that maximumizes
+the referal probability. On the other hand, defining pairs that are
+remote in ages (e.g. period between 1 month and 14 years) is probably
+not that useful. In practice, we may be interested in setting the
+maximum period to, say, five years.
+
+\code{calculate_helpers_hgt()} provides an optional
+pre-calculation for \code{calculate_advice_hgt()}. The user may
+wish to divide up calculations into two steps if intermediate
+results are needed.
+}
+\examples{
+msg(calculate_advice_hgt())
+msgcode <- calculate_advice_hgt(sex = "M", etn = "N",
+                                dob = as.Date("2018-07-31"),
+                                dom1 = as.Date("2018-12-12"), y1 = 64)
+msg(msgcode)
+}
+\author{
+Paula van Dommelen, Stef van Buuren, 2019
+}
diff --git a/man/calculate_birth_z.Rd b/man/calculate_birth_z.Rd
new file mode 100644
index 0000000..1eab7e5
--- /dev/null
+++ b/man/calculate_birth_z.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_birth_z.R
+\name{calculate_birth_z}
+\alias{calculate_birth_z}
+\title{Calculate birth weight SDS relative to Dutch references}
+\usage{
+calculate_birth_z(y, sex, ga, yname = "wgt", dec = 3)
+}
+\arguments{
+\item{y}{Birth weight (grammes) or birth length (cm).
+May be a vector. Converted to numeric.}
+
+\item{sex}{Character, either \code{"M"} (male) or \code{"F"} (female)}
+
+\item{ga}{Gestational age, completed week (Integer or character)}
+
+\item{yname}{Either \code{"wgt"} (for birth weight) or \code{"hgt"}
+(for birth length)}
+
+\item{dec}{Number of decimals for rounding}
+}
+\value{
+\preformatted{ Numeric vector of \code{length(bw)} elements with
+         standard deviation scores relative to Dutch birth
+         weight references
+}
+}
+\description{
+Calculate birth weight SDS relative to Dutch references
+}
+\examples{
+calculate_birth_z(c(2500, 3000), sex = "M", ga = 36)
+}
+\author{
+\preformatted{ Stef van Buuren, 2019
+}
+}
diff --git a/man/calculate_th.Rd b/man/calculate_th.Rd
new file mode 100644
index 0000000..3623294
--- /dev/null
+++ b/man/calculate_th.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_th.R
+\name{calculate_th}
+\alias{calculate_th}
+\title{Calculate target height SDS}
+\usage{
+calculate_th(hgtf, hgtm, sex = NULL, etn = NULL, dec = c(1L, 3L))
+}
+\arguments{
+\item{hgtf}{Length of biological father (cm)}
+
+\item{hgtm}{Length of biological mother (cm)}
+
+\item{sex}{Character, either \code{"M"} (male) or \code{"F"} (female)}
+
+\item{etn}{Etnicity, one of \code{"N"} (dutch), \code{"T"} (turkish),
+\code{"M"} (moroccan) or \code{"H"} (hindustani).}
+
+\item{dec}{Integer vector, length 2, indicating rounding for
+th and th_z, respectively}
+}
+\value{
+\preformatted{ Numeric, length 2: target height (cm) and target height
+         standard deviation score (z), relative to populations
+         living in The Netherlands
+}
+}
+\description{
+Calculate target height SDS
+}
+\examples{
+calculate_th(180, 170, "M", "N")
+}
+\author{
+\preformatted{ Stef van Buuren, 2019
+}
+}
diff --git a/man/datasets.Rd b/man/datasets.Rd
new file mode 100644
index 0000000..f28c313
--- /dev/null
+++ b/man/datasets.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/datasets.R
+\docType{data}
+\name{datasets}
+\alias{datasets}
+\alias{ref.nl5defSGAgewicht}
+\alias{ref.nl5defSGAlengte}
+\alias{nl5.hgt}
+\alias{messages_hgt}
+\title{Table of Dutch birth weight references}
+\format{An object of class \code{data.frame} with 36 rows and 8 columns.}
+\usage{
+ref.nl5defSGAgewicht
+
+ref.nl5defSGAlengte
+
+nl5.hgt
+
+messages_hgt
+}
+\description{
+Table of Dutch birth weight references
+
+Table of Dutch birth length references
+
+Table of Dutch nl5 references
+
+Table of return message and return code
+}
+\keyword{datasets}
diff --git a/man/fold.Rd b/man/fold.Rd
new file mode 100644
index 0000000..9833669
--- /dev/null
+++ b/man/fold.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fold.R
+\name{fold}
+\alias{fold}
+\title{Fold long string across multiple lines}
+\usage{
+fold(s)
+}
+\arguments{
+\item{s}{string}
+}
+\description{
+Fold long string across multiple lines
+}
+\seealso{
+https://stackoverflow.com/questions/2351744/insert-line-breaks-in-long-string-word-wrap
+}
diff --git a/man/msg.Rd b/man/msg.Rd
new file mode 100644
index 0000000..9027057
--- /dev/null
+++ b/man/msg.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/msg.R
+\name{msg}
+\alias{msg}
+\title{Find message string}
+\usage{
+msg(msgcode)
+}
+\arguments{
+\item{msgcode}{Integer with message code}
+}
+\value{
+A string with the message code
+}
+\description{
+Find message string
+}
+\examples{
+msg(31)
+}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..2135931
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(growthscreener)
+
+test_check("growthscreener")
diff --git a/tests/testthat/test-calculate_advice_hgt.R b/tests/testthat/test-calculate_advice_hgt.R
new file mode 100644
index 0000000..0f57de8
--- /dev/null
+++ b/tests/testthat/test-calculate_advice_hgt.R
@@ -0,0 +1,54 @@
+context("calculate_advice_hgt")
+
+kids <- data.frame(
+  dob  = c(as.Date("2018-07-31"), NA, NA, NA, as.Date("2018-07-31"), as.Date("2018-07-31"), as.Date("2018-07-31"), as.Date("2018-07-31"), as.Date("2018-07-31"), as.Date("2018-07-31"),
+           as.Date("2018-07-31"), as.Date("2018-07-31")),
+  dom1 = c(as.Date("2018-12-12"), NA, NA, NA, NA, as.Date("2018-12-12"), as.Date("2018-12-12"), as.Date("2018-12-12"), as.Date("2018-12-12"), as.Date("2018-12-12"),
+           as.Date("2018-12-12"), as.Date("2018-12-12")),
+  dom0 = c(as.Date("2019-03-01"), NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  y1   = c(64, NA, NA, NA, NA, NA, 64, 64, 40, 40,
+           40, 75),
+  y0   = c(60, NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  sex  = c("M", NA_character_, "X", "F", "F", "F", "F", "F", "F", "F",
+           "F", "F"),
+  bw   = c(3000, NA, NA, NA, NA, NA, NA, NA, NA, 3000,
+           2000, 2000),
+  bl   = c(50, NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  ga   = c(40, NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  etn  = c("N", NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, "N", "N", "N",
+           "N", "N"),
+  hgtf = c(180, NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  hgtm = c(170, NA, NA, NA, NA, NA, NA, NA, NA, NA,
+           NA, NA),
+  code = c(31, 19, 19, 16, 15, 18, 20, 31, 13, 45,
+           31, 48),
+  stringsAsFactors = FALSE)
+
+# apply algorithm to kids
+results <- matrix(NA, nrow = nrow(kids), ncol = 3)
+colnames(results) <- c("k", "expected", "found")
+for (k in 1:nrow(kids)) {
+  found <- calculate_advice_hgt(sex  = kids[k, "sex"],
+                                dob  = kids[k, "dob"],
+                                bw   = kids[k, "bw"],
+                                bl   = kids[k, "bl"],
+                                ga   = kids[k, "ga"],
+                                etn  = kids[k, "etn"],
+                                hgtf = kids[k, "hgtf"],
+                                hgtm = kids[k, "hgtm"],
+                                dom1 = kids[k, "dom1"],
+                                y1   = kids[k, "y1"],
+                                dom0 = kids[k, "dom0"],
+                                y0   = kids[k, "y0"])
+  results[k, ] <- c(k, kids$code[k], found)
+}
+
+test_that("expected equals found", {
+  expect_equal(results[, "expected"], results[, "found"])
+})
+
diff --git a/tests/testthat/test-calculate_birth_z.R b/tests/testthat/test-calculate_birth_z.R
new file mode 100644
index 0000000..a89f1c9
--- /dev/null
+++ b/tests/testthat/test-calculate_birth_z.R
@@ -0,0 +1,59 @@
+context("calculate_birth_z")
+
+sex <- "M"
+bw  <- 3000
+ga <- 40
+yname <- "bw"
+
+test_that("functions as intended", {
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga),
+               c(-2.882, -1.245))
+})
+
+test_that("sets extreme bw to NA", {
+  expect_equal(is.na(calculate_birth_z(c(NA, -1000, NaN, NULL), sex, ga)),
+               rep(TRUE, 3))
+})
+
+test_that("takes bw as character", {
+  expect_equal(calculate_birth_z(c("2000", "0", "-1000"), sex, ga),
+               c(-2.882, -5.833, NA))
+})
+
+test_that("sets extreme bw to NA", {
+  expect_equal(is.na(calculate_birth_z(c("NULL", "NA", "", "-", "!"), sex, ga)),
+               rep(TRUE, 5))
+})
+
+test_that("sets invalid sex to NA", {
+  expect_equal(calculate_birth_z(c(2000, 3000), sex = "X", ga),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex = "", ga),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex = TRUE, ga),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex = 0, ga),
+               rep(NA, 2))
+})
+
+test_that("sets invalid ga to NA", {
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = 25),
+               c(5.301, 10.657))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = 25.8),
+               c(5.301, 10.657))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = 24),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = 43),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = "25"),
+               c(5.301, 10.657))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = "25.8"),
+               c(5.301, 10.657))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = "24"),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = "43"),
+               rep(NA, 2))
+  expect_equal(calculate_birth_z(c(2000, 3000), sex, ga = ""),
+               rep(NA, 2))
+})
+