Blake et al. (2016) re-analysis

get conception risk and and arrange by FCD

To do this we subtract 9 (as the 14 days are 8 days preceding ovulation, the day of ovulation, and 5 days after) and add 14 (to obtain days relative to menstrual onset). We then add days with 0 fertility for the remaining 40-14 days of the cycle.

conc_risks_rel_to_ovu = conc %>% 
  select(Index1, BCPred) %>% unique() %>% 
  mutate(Index1 = Index1 - 9)
qplot(Index1, BCPred, data = conc_risks_rel_to_ovu) + geom_vline(xintercept = 0)

conc_risks_rel_to_ovu_interpolated = data.frame(Index1 = seq(-8, 5, 0.1)) %>% 
  left_join(conc_risks_rel_to_ovu)

library(mgcv)
m = gam(BCPred ~ s(Index1, k = 14), data = conc_risks_rel_to_ovu)
# plot(m)
conc_risks_rel_to_ovu_interpolated$BCPred2 = as.numeric(predict(m, newdata = conc_risks_rel_to_ovu_interpolated))
cor.test(conc_risks_rel_to_ovu_interpolated$BCPred,conc_risks_rel_to_ovu_interpolated$BCPred2)

normal_curve = function(x, sd, mean) {
  exp(- ( (x - mean)^2/(2 * sd^2) )) / sqrt( 2*sd^2 * pi )
}
conc_risks_rel_to_ovu_interpolated$BCPred_normal = 1.09*normal_curve(conc_risks_rel_to_ovu_interpolated$Index1, sd = 2, mean = -1.3)


conc_risks_rel_to_ovu_interpolated %>% gather(key,value,-Index1) %>% 
  na.omit() %>% 
ggplot(aes(Index1, value, color = key)) + geom_line()

fertile window probabilities from gangestad ea. 2006

days = data.frame(
    RCD_for_merge = c(29:1, 30:41),
    FCD = c(1:41),
    prc_stirn_b = c(.01, .01, .02, .03, .05, .09, .16, .27, .38, .48, .56, .58, .55, .48, .38, .28, .20, .14, .10, .07, .06, .04, .03, .02, .01, .01, .01, .01, .01, rep(.01, times = 12)),
#                   rep(.01, times = 70)), # gangestad uses .01 here, but I think such cases are better thrown than kept, since we might simply have missed a mens
    prc_wcx_b = c(.000, .000, .001, .002, .004, .009, .018, .032, .050, .069, .085, .094, .093, .085, .073, .059, .047, .036, .028, .021, .016, .013, .010, .008, .007, .006, .005, .005, .005, rep(.005, times = 12))
)
days$RCD = -1 * days$RCD_for_merge + 15

m = gam(prc_stirn_b ~ s(RCD, k = 14), data = days)
# plot(m)
days$prc_stirn_b2 = as.numeric(predict(m, newdata = days))
cor.test(days$prc_stirn_b,days$prc_stirn_b2)

days$prc_stirn_b = days$prc_stirn_b2

expand the main dataset into long format

Essentially we’re trying to reproduce the conception risk dataset.

Best guesses as to what variables mean: - LengthTestingCycle: The length of the cycle in which women actually used the LH strips. - BCActual: The day the LH surge occurred, counted from the end of the cycle - FCDay: The day the LH surge occurred, counted from the beginning of the cycle - LengthRep: The average cycle length reported in the screening survey.

conc_risks_rel_to_ovu = conc_risks_rel_to_ovu_interpolated %>% select(Index1, BCPred2) %>% rename(BCPred = BCPred2)
long = main %>% select(Subject, Age, RelStat, LengthTestingCycle, FCDay, BCActual, LengthRep) %>%   filter(!is.na(LengthTestingCycle)) %>%  # can't use those without complete cycle
  group_by(Subject, Age, RelStat, LengthTestingCycle, FCDay, BCActual, LengthRep) %>%
  # group to retain after expanding
  expand(FCD = full_seq(1:LengthTestingCycle, period = 1)) %>% 
  mutate(
         day_of_ovu_surge = FCDay + 1, # ovulation occurs one day after LH surge
         day_of_ovu_bc_est_m14 = LengthTestingCycle - 1 - 14, # counting 14 back from the next menstrual onset, minus 1 cause we're translating a length into days from 0
         day_of_ovu_fc_est_m14 = LengthRep - 1 - 14, # counting 14 back from the next menstrual onset, minus 1 cause we're translating a length into days from 0
         day_of_ovu_fc = 14,
         day_of_ovu_bc_est_half = (LengthTestingCycle - 1)/2 + 3, # three days after the middle of the cycle
         day_of_ovu_rep_half = (LengthRep - 1)/2 + 3, # three days after the middle of the cycle
         
         FCD0 = FCD - 1, # if we want to subtract the day of ovu, need to start counting from 0
         days_until_ovu_bc_half = FCD0 - day_of_ovu_bc_est_half ,
         days_until_ovu_rep_half = FCD0 - day_of_ovu_rep_half,
         days_until_ovu_bc_m14 = FCD0 - day_of_ovu_bc_est_m14,
         days_until_ovu_fc_m14 = FCD0 - day_of_ovu_fc_est_m14,
         days_until_ovu_fc = FCD0 - day_of_ovu_fc,
         days_until_ovu_surge = FCD0 - day_of_ovu_surge,
         days_until_ovu_bc_m14_squished = if_else(
          LengthTestingCycle - FCD < 14,
          29 - (LengthTestingCycle - FCD),
          round((FCD/ (LengthTestingCycle - 14) ) * 15)) - 15
         ) %>% 
  left_join(conc_risks_rel_to_ovu %>% rename(conc_bc_m14 = BCPred), by = c("days_until_ovu_bc_m14" = "Index1")) %>% 
  left_join(conc_risks_rel_to_ovu %>% rename(conc_bc_m14_sq = BCPred), by = c("days_until_ovu_bc_m14_squished" = "Index1")) %>% 
  left_join(days %>% select(RCD, prc_stirn_b) %>% rename(fert_bc_m14_sq = prc_stirn_b), by = c("days_until_ovu_bc_m14_squished" = "RCD")) %>% 
  left_join(days %>% select(RCD, prc_stirn_b) %>% rename(fert_bc_m14 = prc_stirn_b), by = c("days_until_ovu_bc_m14" = "RCD")) %>% 
  left_join(conc_risks_rel_to_ovu %>% rename(conc_fc_m14 = BCPred), by = c("days_until_ovu_fc_m14" = "Index1")) %>% 
  left_join(days %>% select(RCD, prc_stirn_b) %>% rename(fert_fc_m14 = prc_stirn_b), by = c("days_until_ovu_fc_m14" = "RCD")) %>% 
  left_join(conc_risks_rel_to_ovu %>% rename(conc_ovu_fc = BCPred), by = c("days_until_ovu_fc" = "Index1")) %>% 
  left_join(conc_risks_rel_to_ovu %>% rename(conc_ovu_surge = BCPred), by = c("days_until_ovu_surge" = "Index1")) %>% 
  left_join(days %>% select(RCD, prc_stirn_b) %>% rename(fert_ovu_surge = prc_stirn_b), by = c("days_until_ovu_surge" = "RCD")) %>% 
  mutate(
         conc_bc_m14 = if_else(is.na(conc_bc_m14), 0, conc_bc_m14),
         conc_fc_m14 = if_else(is.na(conc_fc_m14), 0, conc_fc_m14),
         conc_ovu_fc = if_else(is.na(conc_ovu_fc), 0, conc_ovu_fc),
         conc_ovu_surge = if_else(is.na(conc_ovu_surge), 0, conc_ovu_surge)
         ) %>% 
  ungroup()

# View(long %>% select(Subject, FCD, RCD, ovu_day_fcd, ovu_day_rcd) %>% filter(ovu_day_rcd|ovu_day_fcd))
# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_rep_half, conc_rep_half))+geom_point() + geom_vline(xintercept=0)
# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_bc_m14, conc_bc_m14))+geom_point() + geom_vline(xintercept=0)
# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_fc, conc_ovu_fc))+geom_point() + geom_vline(xintercept=0)
# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_surge, conc_ovu_surge))+geom_point() + geom_vline(xintercept=0)

# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_bc_half, days_until_ovu_surge))+geom_point() + geom_abline(slope = 1)
# ggplot(long %>% filter(Subject==2), aes(days_until_ovu_bc_m14, days_until_ovu_surge))+geom_point() + geom_abline(slope = 1)

conception risk correlations

long %>% select(starts_with("conc_")) %>% cor(use='pairwise.complete.obs') %>% round(2)# %>% .[,5]

mean fertility on day of ovulation

long %>% filter(days_until_ovu_surge %in% -2:0) %>% 
  ungroup() %>% 
  select(starts_with("conc")) %>% 
  summarise_each(funs(mean)) %>% pander()

days until ovulation

long %>% select(starts_with("days_until_ovu")) %>% cor(use='na.or.complete') %>% round(2)

long %>% filter(!duplicated(Subject)) %>% select(starts_with("day_of_ovu"))  %>% ungroup() %>% 
  gather(key, value, -day_of_ovu_surge) %>% 
  group_by(key) %>% 
  do(cbind(mean_v =mean(.$value), describe(.$value - .$day_of_ovu_surge))) %>% 
  select(key, mean_v, mean, sd, median, min, max) %>% 
  rename(mean_diff = mean, mean = mean_v) %>% 
  pander()

xsect = long %>% filter(!duplicated(Subject)) %>% ungroup()
qplot(day_of_ovu_surge - day_of_ovu_bc_est_m14, data=xsect, binwidth = 0.5)

qplot(day_of_ovu_surge - day_of_ovu_rep_half, data=xsect, binwidth = 0.5)

qplot(day_of_ovu_surge - day_of_ovu_bc_est_half, data=xsect, binwidth = 0.5)