Unexpected Behavior in Neonate Simulations

[schnebeni]

I have the following code block, which simulates 2000 neonates and then ages them to 2 years of age. I plotted the growth curve of a sample of the simulated individuals using ggplot2. Attached is the resulting growth curves.

library(SimKid)
simpop<-sim_kid(2000,agemin=0,agemax=0,prob_female=0.5,age0isbirth=T,masterseed=123)
Y     <-grow_kid(data=simpop,grow_time=24,tstep=1,age0isbirth=T)

w<-unique(Y$ID)[1:20]
ggplot(Y[Y$ID%in%w,],aes(x=AGE,y=WTKG,col=factor(ID)))+
  geom_line()+
  theme(legend.position="none")+
  facet_wrap(~SEXF)+xlab("AGE")+ylab("WEIGHT")

This doesn't look correct to me for simulations from the CDC growth curves. If I use the following code, which uses the wtage.xls and wtageinf.xls charts available at https://www.cdc.gov/growthcharts/cdc-data-files.htm, I get much smoother curves.

df <- dplyr::bind_rows(
  readxl::read_xls(here::here("data","source","wtage.xls")), 
  readxl::read_xls(here::here("data","source","wtageinf.xls"))
)
df <- df[order(df$Agemos,df$Sex),]
df <- df[!duplicated(df[,c("Agemos","Sex")]) & df$Agemos<=24,]

seed   = 123

N      = 1000

X   = rep(1:nrow(df), each=N)
SEX = df[X,][["Sex"]]
AGE = df[X,][["Agemos"]]
L   = df[X,][["L"]]
M   = df[X,][["M"]]
S   = df[X,][["S"]]

set.seed(seed)
Z   = rep(rnorm(n=N), nrow(df))

WEIGHT = M*((1+ L*S*Z)**(1/L))
Y      = data.frame(Z, SEX, AGE, WEIGHT)

library(ggplot2)

w<-unique(Y$Z)[1:9]
ggplot(Y[Y$Z%in%w,],aes(x=AGE,y=WEIGHT,col=factor(Z)))+
  geom_line()+
  theme(legend.position="none")+
  facet_wrap(~SEX)

Do you know what is causing the discrepancy between the smoothness of these two curves?

[Andy00000000000]

Hi schnebeni,

What you are observing is a function of I believe two things:

1. For 0-2 yr old, SimKid uses the height-for-age and weight-for-height growth charts. Therefore, the assumption that grow_kid() implements is that a given subject will grow at the same respective percentiles of height-for-age and weight-for-height from birth to 24 months.
2. For 0-2 yr old, SimKid will round HTCM in order to allow for an internal merge which causes some jaggedness.

The reason SimKid used weight-for-height growth charts for 0-2 yr old instead of weight-for-age was to preserve the correlation between weight and height by sex and age. This becomes impactful when looking at BMI or BSA rather than just WTKG.

I've attached two examples below highlighting the differences. The first is looking at a single subject at the median and comparing the two methods. The second is looking at a population similar to what you attached above, but expanding to include HTCM and BMI.

I'll try to document these methods more clearly in a future release and perhaps add a note or warning to bring it to the user's attention. Maybe another function or option to use weight-for-age instead of weight-for-height for 0-2 year would be helpful as well when height-weight correlation, BMI, and BSA are not of interest. Let me know if you have any suggestions.

Best,
Andy

Examples:

`

rm(list = ls())

library(SimKid)
library(dplyr)
library(ggplot2)

Example 1. sim_kid and grow_kid for 0-2 yr

simpop <- sim_kid(1, agemin=0, agemax=0, prob_female=0, age0isbirth=T, masterseed=123, htwt_percentile_min = 0.49999, htwt_percentile_max = 0.50001) # Z ~= 0
Y <- grow_kid(data=simpop, grow_time=24, tstep=1, age0isbirth=T)

Example 1. CDC height-for-age and weight-for-age

df <- SimKid::cdc0 %>% # dataframe of CDC growth charts
  mutate(AGEMO = gsub(",.*", "", AGEGRP))%>%
  mutate(AGEMO = readr::parse_number(AGEMO))%>%
  distinct(CHART, VAR, SEXF, AGEMO, .keep_all = T)%>% # keep age 0 as birth, not 0-1 month
  filter(AGEMO <= 24, SEXF == 0) # match to simpop example

df_htage <- df %>% filter(VAR == "HTCM")
df_wtage <- df %>% filter(VAR == "WTKG")

Y2 <- df_htage %>%
  mutate(ID = 2)%>%
  select(ID, AGEMO)

Y2$HTCM <- SimKid::lms_calc(z = 0, l = df_htage$L, m = df_htage$M, s = df_htage$S)
Y2$WTKG <- SimKid::lms_calc(z = 0, l = df_wtage$L, m = df_wtage$M, s = df_wtage$S)

Y2 <- SimKid::calc_bmi_bsa(Y2)

Combine Example 1.

YALL <- Y[,which(colnames(Y) %in% colnames(Y2))] %>%
  bind_rows(Y2)%>%
  mutate(ID = ifelse(ID == 1, "SimKid", "LMS vs. Age"))%>%
  mutate(ID = as.factor(ID))

Plot Example 1.

p_htage <- ggplot(YALL, aes(x = AGEMO, y = HTCM, color = ID))+ 
  geom_line()+
  theme_bw()

p_wtage <- ggplot(YALL, aes(x = AGEMO, y = WTKG, color = ID))+
  geom_line()+
  theme_bw()

p_wtht <- ggplot(YALL, aes(x = HTCM, y = WTKG, color = ID))+
  geom_line()+
  theme_bw()

p_bmiage <- ggplot(YALL, aes(x = AGEMO, y = BMI, color = ID))+
  geom_line()+
  theme_bw()

Example 2. sim_kid and grow_kid for 0-2 yr

rm(list=ls())

simpop <-sim_kid(2000, agemin=0, agemax=0, prob_female=0.5, age0isbirth=T, masterseed=123)
Y <- grow_kid(data=simpop, grow_time=24, tstep=1, age0isbirth=T)

Example 2. CDC height-for-age and weight-for-age

df <- SimKid::cdc0 %>% # dataframe of CDC growth charts
  mutate(AGEMO = gsub(",.*", "", AGEGRP))%>%
  mutate(AGEMO = readr::parse_number(AGEMO))%>%
  distinct(CHART, VAR, SEXF, AGEMO, .keep_all = T)%>% # keep age 0 as birth, not 0-1 month
  filter(AGEMO <= 24) # match to simpop example

df_htage <- df %>% filter(VAR == "HTCM")
df_wtage <- df %>% filter(VAR == "WTKG")

N <- 2000

Y2 <- data.frame(
  ID = sort(rep(1:N, nrow(df_htage)/2)),
  AGEMO = rep(df_htage$AGEMO, N/2),
  SEXF = rep(df_htage$SEXF, N/2)
)%>%
  left_join(df_htage %>% select(AGEMO, SEXF, L, M, S) %>% rename_at(c("L", "M", "S"), ~paste0(., "_HT")))%>%
  left_join(df_wtage %>% select(AGEMO, SEXF, L, M, S) %>% rename_at(c("L", "M", "S"), ~paste0(., "_WT")))

set.seed(123)
z_ht <- rnorm(n = N)

set.seed(124)
z_wt <- rnorm(n = N)

Y2 <- Y2 %>%
  left_join(data.frame(ID = 1:N, ZHTCM = z_ht, ZWTKG = z_wt))

Y2$PHTCM <- 100*pnorm(Y2$ZHTCM)
Y2$PWTKG <- 100*pnorm(Y2$ZWTKG)

Y2$HTCM <- SimKid::lms_calc(z = Y2$ZHTCM, l = Y2$L_HT, m = Y2$M_HT, s = Y2$S_HT)
Y2$WTKG <- SimKid::lms_calc(z = Y2$ZWTKG, l = Y2$L_WT, m = Y2$M_WT, s = Y2$S_WT)

Y2 <- SimKid::calc_bmi_bsa(Y2)

Y2 <- Y2 %>%
  mutate(GAWK = 40, CHART = "CDC", MONTH = AGEMO, AGE = AGEMO/12)

Y2 <- Y2[,colnames(Y)] # match format to SimKid dataframe

Example 2. Compare Validation

p <- validate_kid(Y, age0isbirth = T, overlay_percentile = 0.50)
p2 <- validate_kid(Y2 %>% mutate(HTCM = round(HTCM)), age0isbirth = T, overlay_percentile = 0.50)

p[[1]] + coord_cartesian(xlim = c(0, 24), ylim = c(40, 120)) + labs(title = "SimKid")
p2[[1]] + coord_cartesian(xlim = c(0, 24), ylim = c(40, 120))+ labs(title = "LMS vs. Age")

p[[2]] + coord_cartesian(xlim = c(0, 24), ylim = c(0, 25)) + labs(title = "SimKid")
p2[[2]] + coord_cartesian(xlim = c(0, 24), ylim = c(0, 25))+ labs(title = "LMS vs. Age")

p[[3]] + coord_cartesian(xlim = c(35, 100), ylim = c(0, 20)) + labs(title = "SimKid")
p2[[3]] + coord_cartesian(xlim = c(35, 100), ylim = c(0, 20))+ labs(title = "LMS vs. Age")

p[[4]] + coord_cartesian(xlim = c(0, 25), ylim = c(0, 25)) + labs(title = "SimKid")
p2[[4]] + coord_cartesian(xlim = c(0, 25), ylim = c(0, 25))+ labs(title = "LMS vs. Age")

End

`