library(knitr)
opts_chunk$set(fig.width=10,
               fig.height=6,
               fig.path='Figs/',
               #echo=FALSE,         # Toggle this to show the code
               warning=FALSE,
               message=FALSE,
               cache=TRUE,
               results='asis',
               root.dir = '~/janus/')

Reproduce BEIR 10-2 (fail)

Last update: May 2014

Reproduce the BEIR estimates of dose reponse for atomic bomb survivors.

# Common
setwd('~/janus')
source('scripts/util/util.R') # http://goo.gl/VYzkAs

# Data
data <- read.csv('~/janus/data/lss14/lss14_spruced_up.csv')

# Constants
graphable_threshold = 2.0
modelable_threshold = 1.5

# Filter the data
data <- data %>%
  select(city, sex, gd3, ahs, ctime, solid,
         subjects, agex, age, dose, pyr, death) %>%
  filter(dose <= graphable_threshold)

# Model
model <- function(new_terms) {
  formula <- solid ~ city +
                     sex*I(log(age + 1)) +
                     sex*I(log(agex + 1)) +
                     offset(log(pyr)) -
                     1

  formula <- update(formula, new_terms)

  glm(formula,
      family='poisson',
      data = data)
}

factor_fit <- model(~ . + factor(dose))

# Data points to predict
to_predict <- expand.grid(
  city = unique(data$city),
  sex = unique(data$sex),
  agex = 30,
  age = 60,
  pyr = 1,
  dose = unique(data$dose))

# Make predictions
err <- function(model, data) {
  baseline <- data
  baseline$dose <- min(data$dose)

  risk <- exp(predict(model, newdata=data))
  control_risk <- exp(predict(model, newdata=baseline))
  err <- risk / control_risk - 1

  err
}
to_predict$err <- err(factor_fit, to_predict)

# Average by Dose and city
aggregate <- to_predict %>%
  group_by(dose) %>%
  summarize(err=mean(err))

Show

10-2

This is the original 10-2 figure from the beir VII report.

:

Reproduce 10-2

This is my attempt to reproduce 10-2 to prove that I am faithfully applying their methodology.

# Show data
lm_smooth <- function(...) geom_smooth(method='lm',
                                       se=FALSE,
                                       color='black',
                                       fullrange=TRUE,
                                       ...)

data <- aggregate[aggregate$dose <= modelable_threshold,]
ggplot(data, aes(dose, err)) +
  lm_smooth(formula='y ~ x - 1') +
  lm_smooth(formula='y ~ I(x + 0.3*x^2) - 1') +
  lm_smooth(formula='y ~ I(x + 0.7*x^2) - 1') +
  geom_point(data=aggregate, aes(dose, err))