how to caclulate HFD's Bongaarts-Feeney adjusted TFR using dplyr

HFD_TFRadj.R

library(HMDHFDplus)
library(tidyverse)
library(data.table)
# read and reshape births by year, age, order
Bxi <- readHFDweb("SWE","birthsRRbo", 
                  username = Sys.getenv("us"), 
                  password = Sys.getenv("pw"))|> 
  select(-OpenInterval, -Total) |> 
  pivot_longer(B1:B5p, names_to = "order", values_to = "Bxi") |> 
  mutate(order = parse_number(order))

# read exposures
Ex <- readHFDweb("SWE","exposRR", username =Sys.getenv("us"), password = Sys.getenv("pw"))|> 
  select(-OpenInterval)


TFRadj <-
  # merge births and exposures
  Bxi |> 
  left_join(Ex, by = join_by(Year, Age)) |> 
  # calculate order-specific rates 
  # (using UNCONDITIONAL exposures)
  mutate(Fxi = Bxi / Exposure,
         # age mid is x bar in HFDMP
         age_mid = Age + .5) |> 
  group_by(Year, order) |> 
  # calc ingredients for eq 6.9b
  summarize(TFRi = sum(Fxi),
            MABi = sum(Fxi * age_mid) / TFRi,
            .groups = "drop") |> 
  arrange(order, Year) |> 
  group_by(order) |> 
  mutate(ri = (shift(MABi,-1) - shift(MABi,1)) / 2, # eq 6.10
         TFRadji = TFRi / (1 - ri)) |> # eq 6.9b
  ungroup() |> 
  group_by(Year) |> 
  summarize(TFRadj = sum(TFRadji)) # 6.9a

# calculate TFR to compare
# also repeat TFRadj ignoring birth order:
TFRadj2 <-
  Bxi |> 
  left_join(Ex, by = join_by(Year, Age)) |> 
  group_by(Year, Age) |> 
  summarize(Fx = sum(Bxi) / Exposure[1],
            .groups = "drop") |> 
  group_by(Year) |> 
  summarize(TFR = sum(Fx),
            MAB = sum(Fx * (Age + .5)) / TFR,
            .groups = "drop") |> 
  mutate(r = (shift(MAB,-1) - shift(MAB,1))/2,
         TFRadj2 = TFR / (1 - r)) 

# plot and compare, adjusted is indeed more tame,
# but still maybe smooth it?
# See how taking birth order into account gives a less
# erratic signal?
left_join(TFRadj, TFRadj2, by = join_by(Year)) |> 
  ggplot(aes(x = Year, y = TFRadj)) +
  geom_line(color = "red") +
  geom_line(mapping = aes(y = TFR), color = "blue") +
  geom_line(mapping = aes(y = TFRadj2), color = "red", linetype = 2) +
  theme_minimal()