Simulation of BMI values from Figure 1 of Rizk et al. 10.1016/j.fertnstert.2004.07.960

Rizk-BMI.R

# Simulation of BMI values from Figure 1 of Rizk et al. 10.1016/j.fertnstert.2004.07.960

exact.rnorm <- function(n, m, sd) {
  x <- rnorm(n)
  x <- (x - mean(x)) / sd(x)        # standardize to mean 0, sd 1
  x <- x * sd + m                   # rescale to target mean and sd
  return(x)
}

# Mean and SD of height and weight from Table 1
m.h1.t1 <- 164.1 + 0.05      # round to maximise height
sd.h1.t1 <- 5.31 + 0.005
m.h2.t1 <- 162.5 + 0.05
sd.h2.t1 <- 5.7 + 0.05
m.w1.t1 <- 101.1 - 0.05      # round to minimise weight
sd.w1.t1 <- 12.4 - 0.05
m.w2.t1 <- 99.2 + 0.05
sd.w2.t1 <- 12.7 + 0.05

paper.bmi1 <- (m.w1.t1 / (m.h1.t1 ^ 2)) * 10000
paper.bmi2 <- (m.w2.t1 / (m.h2.t1 ^ 2)) * 10000

n.runs <- 100000
for (i in 1:n.runs) {
  if (i == 1) {
    bmi1.less <- 0
    bmi2.less <- 0
    min.bmi1 <- 999
    min.bmi2 <- 999
    set.seed(0)
  }

  h1 <- exact.rnorm(75, m.h1.t1, sd.h1.t1)
  h2 <- exact.rnorm(75, m.h2.t1, sd.h2.t1)
  w1 <- exact.rnorm(75, m.w1.t1, sd.w1.t1)
  w2 <- exact.rnorm(75, m.w2.t1, sd.w2.t1)

  m.bmi1 <- mean(w1 / (h1 ^ 2)) * 10000
  m.bmi2 <- mean(w2 / (h2 ^ 2)) * 10000
  if (m.bmi1 <= (30.5 + 0.05)) {        # round up comparison value to give authors the best chance
    bmi1.less <- bmi1.less + 1
  }
  if (m.bmi1 < min.bmi1) {
    min.bmi1 <- m.bmi1
  }
  if (m.bmi2 <= (30.1 + 0.05)) {
    bmi2.less <- bmi2.less + 1
  }
  if (m.bmi2 < min.bmi2) {
    min.bmi2 <- m.bmi2
  }
}

cat(
  "BMI1: "
, "count <= paper:", bmi1.less
, "  Minimum:", sprintf("%.3f", min.bmi1)
, "  Paper 'mean BMI':", sprintf("%.3f", paper.bmi1)
, "\n", sep=""
)
cat(
  "BMI2: "
, "count <= paper:", bmi2.less
, "  Minimum:", sprintf("%.3f", min.bmi2)
, "  Paper 'mean BMI':", sprintf("%.3f", paper.bmi2)
, "\n", sep=""
)