Censuses miss young children everywhere, and no smoothing repairs that, because the children were never counted in the first place. This part rebuilds them from a different source — the births that must have happened, survived forward with a life table — then joins that reconstruction to the rest of the distribution.

You do not need to have done Parts 1 and 2. Sheets 01 and 02 are what those parts produce. Fill in the decision log (sheet 06) as you go.

Running it. Put UGA_C2_lab3_data.xlsx in the same folder as this file. Then either run the chunks one at a time in RStudio, or press Knit to produce this report — from the console that is rmarkdown::render("C2_lab3_children.Rmd"). Nothing else is needed: the protocol functions travel with the file, in an appendix at the end.

library(readxl)     # reads the workbook
library(dplyr)      # data handling the workflow functions use
library(tibble)     # tibble(), used by one of the protocol functions
library(DemoTools)  # the demographic methods underneath everything
library(knitr)      # kable(), for readable tables

DATA <- "UGA_C2_lab3_data.xlsx"

## Stop early, with directions, if the workbook cannot be opened: a missing
## file, the wrong working directory, or a OneDrive file that has not been
## downloaded to this computer all land here.
if (!isTRUE(tryCatch({ excel_sheets(DATA); TRUE }, error = function(e) FALSE)))
  stop("Cannot open ", DATA, " from ", getwd(),
       ". Put the workbook in the same folder as this file, or give DATA ",
       "its full path. If it sits in OneDrive, right-click the file, choose ",
       "'Always keep on this device', close it in Excel, and try again.")

1 Read the workbook

smoothed   <- read_excel(DATA, sheet = "01_start_smoothed",   skip = 7, n_max = 106)
unsmoothed <- read_excel(DATA, sheet = "02_start_unsmoothed", skip = 3, n_max = 106)
lt         <- read_excel(DATA, sheet = "03_lifetables",       skip = 3, n_max = 37)
asfr       <- read_excel(DATA, sheet = "04_asfr",             skip = 3, n_max = 14)
srb        <- read_excel(DATA, sheet = "05_srb",              skip = 2, n_max = 1)

details <- read_excel(DATA, sheet = "01_start_smoothed", range = "A4:B6",
                      col_names = c("field", "value"))
val <- function(f) details$value[details$field == f]
REFDATE <- as.numeric(val("census_reference_date"))
LOCID   <- as.numeric(val("LocID"))
SMOOTH  <- val("smooth_method")

Reverse survival needs its inputs as small tables — ages down the rows, one column per reference date. Sheet 03 holds both life tables stacked, so the two dates are separated here and put side by side. The earlier date comes first; that ordering is what the method expects.

dates <- sort(unique(lt$RefDate))
side_by_side <- function(column, up_to) {
  ages <- sort(unique(lt$Age[lt$Age <= up_to]))
  m <- sapply(dates, function(d) lt[[column]][lt$RefDate == d][match(ages, lt$Age[lt$RefDate == d])])
  rownames(m) <- ages
  m
}
nLxF <- side_by_side("nLx_F", 55)   # mothers: female survivorship to age 55
nLxM <- side_by_side("nLx_M", 5)    # the children themselves

Fertility the same way. The workbook publishes rates per 1,000 women, which is how national reports print them; the method works in births per woman, so divide by 1,000 here rather than storing a converted number in the data.

adates  <- sort(unique(asfr$RefDate))
AsfrMat <- sapply(adates, function(d) asfr$ASFR_per1000[asfr$RefDate == d] / 1000)
rownames(AsfrMat) <- seq(15, 45, 5)

2 Step F — rebuild the young children, then splice

The reconstruction asks a different question from the census: given the women of childbearing age counted in this census, the fertility rates of the years before it, and the chance of surviving from birth to the census, how many children under 10 should there be? Where that number exceeds the count, the census missed children, and the protocol replaces the count.

It does not replace everything. The splice keeps the reconstruction only up to the age where it stops exceeding the count, then returns to the census’s own numbers. The function finds that age itself and reports it.

bp <- census_workflow_adjust_basepop(
  popM1 = smoothed$Male, popF1 = smoothed$Female,
  popM_unsmoothed = unsmoothed$Male, popF_unsmoothed = unsmoothed$Female,
  Age_unsmoothed  = unsmoothed$AgeStart,
  smooth_method = SMOOTH, LocID = LOCID, census_reference_date = REFDATE,
  nLxMatFemale = nLxF, nLxMatMale = nLxM, nLxMatDatesIn = dates,
  AsfrMat = AsfrMat, AsfrDatesIn = adates,
  ## basepop_five uses a supplied sex ratio at birth as-is at three internal
  ## dates, so a single value would leave two of them empty; repeat it.
  SRB = rep(srb$SRB[1], 3), SRBDatesIn = NULL,
  radix = 100000)

take <- function(label, sex) bp$DataValue[bp$BPLabel == label & bp$SexID == sex]
base <- data.frame(Age = 0:105, Male = take("BP4", 1), Female = take("BP4", 2))

Below, the count the census gave and the population the reconstruction produced, at the ages where they differ. The last column is the share of children the census missed at each age.

kable(data.frame(
  Age        = 0:9,
  `Census, smoothed` = format(round((smoothed$Male + smoothed$Female)[1:10]), big.mark = ","),
  `Base population`  = format(round((base$Male + base$Female)[1:10]), big.mark = ","),
  `Missed (%)` = round(100 * ((base$Male + base$Female)[1:10] /
                              (smoothed$Male + smoothed$Female)[1:10] - 1), 1),
  check.names = FALSE))
Age Census, smoothed Base population Missed (%)
0 1,399,725 1,660,243 18.6
1 1,378,362 1,633,255 18.5
2 1,395,906 1,608,845 15.3
3 1,434,508 1,588,140 10.7
4 1,476,026 1,570,737 6.4
5 1,468,939 1,553,271 5.7
6 1,411,668 1,530,135 8.4
7 1,363,184 1,497,963 9.9
8 1,320,891 1,456,195 10.2
9 1,296,724 1,405,177 8.4 
pstart()
plot(0:14, (smoothed$Male + smoothed$Female)[1:15] / 1000, type = "l", lwd = 2, col = GREY,
     bty = "l", xlab = "Age", ylab = "Thousands",
     ylim = range(c((smoothed$Male + smoothed$Female)[1:15],
                    (base$Male + base$Female)[1:15])) / 1000 * c(0.98, 1.04))
lines(0:14, (base$Male + base$Female)[1:15] / 1000, lwd = 2, col = BLUE)
legend("topright", c("Census, smoothed", "Base population, children rebuilt"),
       lwd = 2, col = c(GREY, BLUE), bty = "n")

The gap closes with age, which is the shape this correction should have: the youngest children are the most likely to be left off a census form, and by school age the count is nearly right.

Your turn. Record in the decision log the age at which the splice stopped using the reconstruction, and how many children the step added in total. Then say in one sentence why a census would miss an infant.

3 The base population

Reported at 100+, as the protocol asks: the internal 105+ is a working convenience, not a published age group.

base100 <- rbind(base[base$Age < 100, ],
                 data.frame(Age = 100,
                            Male   = sum(base$Male[base$Age >= 100]),
                            Female = sum(base$Female[base$Age >= 100])))

cat("population before this part ",
    format(round(sum(unsmoothed$Male + unsmoothed$Female)), big.mark = ","), "\n")
## population before this part  47,968,041
cat("base population             ",
    format(round(sum(base100$Male + base100$Female)), big.mark = ","), "\n")
## base population              49,564,228
cat("children added by step F    ",
    format(round(sum(base100$Male + base100$Female) -
                 sum(unsmoothed$Male + unsmoothed$Female)), big.mark = ","), "\n")
## children added by step F     1,596,187
pstart()
plot(base$Age, (base$Male + base$Female) / 1000, type = "l", lwd = 2, col = BLUE,
     bty = "l", xlab = "Age", ylab = "Thousands per single year of age")
lines(unsmoothed$AgeStart, (unsmoothed$Male + unsmoothed$Female) / 1000,
      lwd = 1.2, lty = 3, col = GREY)
legend("topright", c("Base population", "Before this part"),
       lwd = c(2, 1.2), lty = c(1, 3), col = c(BLUE, GREY), bty = "n")

dir.create("out3", showWarnings = FALSE)
write.csv(base100, "out3/uganda_base_population.csv", row.names = FALSE)

4 Where this part ends

That table is a base population: evaluated, adjusted, single years of age and sex, with every choice recorded in the log beside it.

5 Bring your own census

Copy the provided input file, clear the green cells, and enter your own values with their sources and reference dates.

Expect that to take much longer than the lab does, and to happen outside the session. Many offices may find that some of the inputs are not published for their country, and finding out which ones is worth as much as the run itself. Record each one in the decision log.

6 Appendix — the protocol functions this part uses

These are the Population Division’s own censusAdjust files.Underneath they call DemoTools; what they add is the protocol’s decision rules.

## ==== census_workflow_adjust_basepop.R ======================
## census_workflow_adjust_basepop
## 
## This function carries out the basepop adjustment for children missing from census counts
## as described in protocol flowchart xx
## 
## @param popM1 numeric. vector of male population counts by single year of age
## @param popF1 numeric. vector of female population counts by single year of age
## @return A list of objects describing the steps implemented in adjusting for missing children
## @export

census_workflow_adjust_basepop <- function(popM1, # male pop counts by single year of age
                                           popF1, # female pop counts by single year of age
                                           popM_unsmoothed, # male pop counts from before smoothing step of workflow (can be single, abridged or five-year)
                                           popF_unsmoothed, # female pop counts from before smoothing step of workflow
                                           Age_unsmoothed, # starting age of age groups for unsmoothed series
                                           smooth_method = NA, # smoothing method used
                                           LocID,
                                           census_reference_date, # decimal year
                                           nLxMatFemale = NULL, # matrix of nLx life table values for females. If NULL then values from DemoToolsData will be used.
                                           nLxMatMale = NULL, # matrix of nLx life table values for males
                                           nLxMatDatesIn = NULL, # dates associated with nLx matrices
                                           AsfrMat = NULL, # matrix of age-specific fertility rates. If NULL then DemoToolsData values are used.
                                           AsfrDatesIn = NULL, # dates associated with ASFR matrix
                                           SRB = NULL, # vector of sex ratio at birth
                                           SRBDatesIn = NULL, # dates associated with SRB vector
                                           radix = NULL)  { # radix associated with nLx values

  Age1 <- 1:length(popM1)-1
  # group to abridged age groups
  popM_abr <- DemoTools::single2abridged(popM1)
  popF_abr <- DemoTools::single2abridged(popF1)
  Age_abr  <- as.numeric(row.names(popM_abr))

  # run basepop_five()
  BP1 <- DemoTools::basepop_five(location = LocID,
                                 refDate = census_reference_date,
                                 Age = Age_abr,
                                 Females_five = popF_abr,
                                 Males_five = popM_abr,
                                 nLxFemale = nLxMatFemale,
                                 nLxMale   = nLxMatMale,
                                 nLxDatesIn = nLxMatDatesIn,
                                 AsfrMat = AsfrMat,
                                 AsfrDatesIn = AsfrDatesIn,
                                 SRB = SRB,
                                 SRBDatesIn = SRBDatesIn,
                                 radix = radix,
                                 verbose = FALSE)

  # graduate result to single year of age
  popM_BP1 <- DemoTools::graduate_mono(Value = BP1[[2]], Age = Age_abr, AgeInt = DemoTools::age2int(Age_abr), OAG = TRUE)
  popF_BP1 <- DemoTools::graduate_mono(Value = BP1[[1]], Age = Age_abr, AgeInt = DemoTools::age2int(Age_abr), OAG = TRUE)

  # what is the minimum age at which BP1 is not higher than input population for both males and females
  BP1_higher <- popM_BP1 > popM1 & popF_BP1 > popF1
  minLastBPage1 <- min(Age1[!BP1_higher],10) - 1

  if (minLastBPage1 >= 0) {

  # graduate the unsmoothed series to single age if necessary
  AgeInt_unsmoothed <- DemoTools::age2int(Age_unsmoothed)
  if (!(max(AgeInt_unsmoothed, na.rm=TRUE)==1)) {
    popM_unsmoothed <- DemoTools::graduate_mono(Value = popM_unsmoothed,
                                                Age = Age_unsmoothed,
                                                AgeInt = AgeInt_unsmoothed,
                                                OAG = TRUE)
    popF_unsmoothed <- DemoTools::graduate_mono(Value = popF_unsmoothed,
                                                Age = Age_unsmoothed,
                                                AgeInt = AgeInt_unsmoothed,
                                                OAG = TRUE)
  }
  Age1_unsmoothed <- 1:length(popM_unsmoothed) - 1

  # splice the BP1 series for ages at or below minLastBPage1 with unsmoothed single age series to age 15 and smoothed series thereafter
  popM_BP2 <- c(popM_BP1[Age1 <= minLastBPage1], popM_unsmoothed[Age1_unsmoothed > minLastBPage1 & Age1_unsmoothed < 15], popM1[Age1 >= 15])
  popF_BP2 <- c(popF_BP1[Age1 <= minLastBPage1], popF_unsmoothed[Age1_unsmoothed > minLastBPage1 & Age1_unsmoothed < 15], popF1[Age1 >= 15])

    # if we are smoothing, then smooth BP2 using the best method from child smoothing before
    if (!is.na(smooth_method)) {

      if (substr(smooth_method, 1, 8) == "bestMavN") { # if the best smoothing was mav on one year data

        mavN <- as.numeric(substr(smooth_method, nchar(smooth_method)-1, nchar(smooth_method)))
        popM_BP3_mav <- mavPop1(popM_BP2, Age1)
        popM_BP3     <- unlist(select(popM_BP3_mav$MavPopDF, !!paste0("Pop", mavN)))
        popF_BP3_mav <- mavPop1(popF_BP2, Age1)
        popF_BP3     <- unlist(select(popF_BP3_mav$MavPopDF, !!paste0("Pop", mavN)))

      } else { # if the best smoothing was on five-year data

        popM5_BP2 <- DemoTools::groupAges(popM_BP2, N=5)
        popF5_BP2 <- DemoTools::groupAges(popF_BP2, N=5)
        Age5      <- seq(0,max(Age_abr),5)

        bestGrad5 <- as.numeric(substr(smooth_method, nchar(smooth_method), nchar(smooth_method)))

        if (bestGrad5 == 1) {
          popM_BP3 <- DemoTools::graduate_mono(popM5_BP2, AgeInt = DemoTools::age2int(Age5), Age = Age5, OAG = TRUE)
          popF_BP3 <- DemoTools::graduate_mono(popF5_BP2, AgeInt = DemoTools::age2int(Age5), Age = Age5, OAG = TRUE)
        }
        if (bestGrad5 == 2) {
          popM5_BP2_mav2 <- DemoTools::smooth_age_5(popM5_BP2, Age5, method = "MAV", n = 2)
          popF5_BP2_mav2 <- DemoTools::smooth_age_5(popF5_BP2, Age5, method = "MAV", n = 2)
          # splice infants from BP1 to 1-4 year olds from smoothed BP2 and remaining smoothed BP2 thereafter
          popMabr_BP2_mav2 <- c(popM_BP1[1], popM5_BP2_mav2[1]-popM_BP1[1], popM5_BP2_mav2[2:length(popM5_BP2_mav2)])
          popFabr_BP2_mav2 <- c(popF_BP1[1], popF5_BP2_mav2[1]-popF_BP1[1], popF5_BP2_mav2[2:length(popF5_BP2_mav2)])

          popM_BP3 <- DemoTools::graduate_mono(popMabr_BP2_mav2, AgeInt = DemoTools::age2int(Age_abr), Age = Age_abr, OAG = TRUE)
          popF_BP3 <- DemoTools::graduate_mono(popFabr_BP2_mav2, AgeInt = DemoTools::age2int(Age_abr), Age = Age_abr, OAG = TRUE)
        }

      }
      popM_BP3 <- c(popM_BP3[Age1 < 15], popM1[Age1 >=15])
      popF_BP3 <- c(popF_BP3[Age1 < 15], popF1[Age1 >=15])

    } else { # if no smoothing then BP3 = BP2
      popM_BP3 <- popM_BP2
      popF_BP3 <- popF_BP2

    }

    # what is the minimum age at which BP1 is higher than BP3 for both males and females
    BP1_higher <- popM_BP1 >= popM_BP3 & popF_BP1 >= popF_BP3
    minLastBPage3 <- min(Age1[!BP1_higher],10) - 1

    # splice the BP1 up to age minLastBPage3 with the BP3
    popM_BP4 <- c(popM_BP1[Age1 <= minLastBPage3], popM_BP3[Age1 > minLastBPage3 & Age1 < 15], popM1[Age1 >= 15])
    popF_BP4 <- c(popF_BP1[Age1 <= minLastBPage3], popF_BP3[Age1 > minLastBPage3 & Age1 < 15], popF1[Age1 >= 15])

  } else { # IF BP1 CAME IN LOWER THAN ORIGINAL AT ALL CHILD AGES, THEN WE SKIP THE STEPS AND JUST RETURN ORIGINAL

    popM_BP2 <- popM1
    popM_BP3 <- popM1
    popM_BP4 <- popM1
    popF_BP2 <- popF1
    popF_BP3 <- popF1
    popF_BP4 <- popF1

  }

nAge <- length(Age1)
pop_basepop <- data.frame(SexID = rep(c(rep(1,nAge),rep(2,nAge)),4),
                          AgeStart = rep(Age1,8),
                          BPLabel = c(rep("BP1",nAge*2),
                                      rep("BP2",nAge*2),
                                      rep("BP3",nAge*2),
                                      rep("BP4",nAge*2)),
                          DataValue = c(popM_BP1, popF_BP1,
                                        popM_BP2, popF_BP2,
                                        popM_BP3, popF_BP3,
                                        popM_BP4, popF_BP4))

return(pop_basepop)

}

## ==== getmavPop1.R ==========================================
# Test of functions to create moving average estimates of population by single years of age
# and birth cohort estimates based on single year data and reference date
# note Bachi defaults ageMin=23 and ageMax=77
# 11/19/20  rc code to fill in youngest oldest ages using seq of mavs
# 11/25/20 fixed mutate
# 12/17/20 more fixes 
# 2021-01-20 remove BirthCohort1 (switch to DT function )
#            No longer do birth cohorts here
# 2021-01-29 remove sex, CenDate, and CenYr entries from mavPop1

#

SCcount <- function(Pop) {
  # count number of sign changes age to age
  # get Pop(a) - Pop(a-1) only look at ages 23-67
  PopCh <- Pop[25:68] - Pop[24:67]
  PopCHsign <- sign(PopCh)
  PopCHsignChange <- PopCHsign[1:43] != PopCHsign[2:44]
  sum(PopCHsignChange)
}

#

# new version of function with help from rc 11/19/20
mavPop1 <- function(Pop, Age) {
  
  # store pop data for sending back
  MavPopDF <- tibble(Age = Age, Pop1 = Pop)
  Nage <- length(Age) - 2

  # get results for original data
  BachiN <- check_heaping_bachi(Pop, Age, method = "pasex")
  MySCcount <- SCcount(Pop)
  MavSummaryDF <- tibble(MavN = 1,BachiN, SCcountN = MySCcount)

  for  (MavN in 2:10)  {
    MyMavNew <- mav(Pop, n = MavN, Age = Age)
    MavPopDF <- cbind(MavPopDF, MyMavNew)
    colnames(MavPopDF)[ncol(MavPopDF)]  <- paste("Pop", MavN, sep = "")
  }
  
  OpenAge <- Nage + 1

  MavPopDF <- mutate( MavPopDF,  across(3:11, ~(
    function(x) {
      for(i in 1:length(x))
      {
        if(is.na(x[i]) & (Age[i] == 0 | Age[i] >= (OpenAge - 1)))
          x[i] <- Pop1[i]
        
        else if(is.na(x[i]) & (Age[i] == 1 | Age[i] == (OpenAge - 2)))
          x[i] <- Pop2[i]
        
        else if(is.na(x[i]) & (Age[i] == 2 | Age[i] == (OpenAge - 3)))
          x[i] <- Pop4[i]
        
        else if(is.na(x[i]) & (Age[i] == 3 | Age[i] == (OpenAge - 4)))
          x[i] <- Pop6[i]
        
        else if(is.na(x[i]) & (Age[i] == 4 | Age[i] == (OpenAge - 5)))
          x[i] <- Pop8[i]
        
      }
      return(x)
    } ) # end of function spec
    (.) ) # end of across spec  I still don't know the function of (.)
  ) # end of mutate
  
  
  for (i in 3:ncol(MavPopDF)){ 
    BachiN <- check_heaping_bachi(MavPopDF[[i]], Age, method = "pasex")
    SCcountN <- SCcount(MavPopDF[[i]])
    
    MavSummaryDF <- rbind(MavSummaryDF,c(MavN=i, BachiN, SCcountN))
  } 

  rownames(MavPopDF) <- MavPopDF$Age
  
  list(MavPopDF = MavPopDF,
       MavSummaryDF = MavSummaryDF)
  
} # end of mavPop1 function