Forest Inventory: Application of Horvitz-Thompson estimators in a two-stage sample

Selecting stands

The population of stands could be sampled randomly, but since we are trying to obtain an estimate of total volume for the property, we need to make sure that we sample the stands that a) contribute the most to the total volume, and b) contribute the most variance to the population. Using the Horvitz-Thompson estimators, we can sample the stands however we want so long as we know the probability of each stand being included in the sample! This is a case of unequal probability sampling and the Horvitz-Thompson estimators will be statistically efficient if the probability of inclusion for each primary sampling unit is correlated with the variable of interest. In this case, we are going to weight our sample on age x acres since we believe the oldest stands are going to have the greatest stocking, and that the largest stands have the greatest total volume.

stands$prop <- (stands$age * stands$acres) /
  sum(stands$age * stands$acres)

stands %>%
  ggplot(aes(x = age * acres, y = cuftTot)) +
  geom_point() +
  ylab(bquote("total volume (" ~ ft^{3} ~ ")"))

simulate_sample <- function(sims = 10000, stands, n, weight_var,
                            id_var, cores) {
  weight_var <- enquo(weight_var)
  id_var <- enquo(id_var)
  
  allIds <- stands %>% pull(!! id_var)
  
  simFrame <- bind_rows(
    parallel::mclapply(1:sims,
      function(i, ...) {
        stands %>%
          sample_n(n, weight = !! weight_var)
      },
      mc.cores = cores
    ),
    .id = "sim"
  )
  
  sampledIds <- simFrame %>% pull(!! id_var)
  
  if (!all(allIds %in% sampledIds)) {
    missing <- allIds[!allIds %in% sampledIds]
    warning(
      paste0(
        "these elements that were not included in any sample: ",
        paste(missing, collapse = " ")
      )
    )
  }
  
  list(
    simFrame = simFrame,
    sims = sims,
    n = n,
    ids = allIds
  )
}

estimate_pi_i <- function(simList, id_var) {
  id_var <- enquo(id_var)
  
  simList[["simFrame"]] %>%
    group_by(!! id_var) %>%
    summarize(pi_i = n() / simList[["sims"]]) %>%
    select(!! id_var, pi_i)
}

estimate_joint_probs <- function(simList, id_var, pairs = NULL, cores) {
  id_var <- enquo(id_var)
  
  if(is.null(pairs)) {
    ids <- simList[["ids"]]
    idxs <- 1:length(ids)
    pairs <- combn(ids, m = 2, simplify = FALSE)
  } else {
    ids <- unique(unlist(pairs))
    idxs <- 1:length(ids)
  }

  estimate_joint_prob <- function(pair, simFrame, sims, cores, id_var) {
    simFrame %>%
      filter(!! id_var %in% pair) %>%
      group_by(sim) %>%
      mutate(n = n()) %>%
      filter(n > 1) %>%
      ungroup() %>%
      select(sim) %>%
      distinct() %>%
      tally %>%
      rowwise() %>%
      mutate(
        id1 = pair[1], id2 = pair[2],
        row = idxs[which(ids == id1)],
        col = idxs[which(ids == id2)],
        jointProb = n / sims
      ) %>%
      select(-n)
  }
  
  f <- bind_rows(
    parallel::mclapply(
      pairs,
      estimate_joint_prob,
      simFrame = simList[["simFrame"]],
      id_var = id_var,
      sims = simList[["sims"]],
      mc.cores = cores
    )
  )

  frame <- matrix(NA, nrow = length(ids), ncol = length(ids))
  frame[as.matrix(f[3:4])] <- f$jointProb
  frame[as.matrix(f[4:3])] <- f$jointProb

  output <- data.frame(frame)
  row.names(output) <- ids
  names(output) <- ids

  return(output)
}

simList20 <- simulate_sample(
  sims = 10000,
  stands = stands,
  n = 20,
  id_var = id,
  weight_var = prop,
  cores = 2
)

sampleSim20 <- estimate_pi_i(simList = simList20, id_var = id)

jointProbs <- estimate_joint_probs(
  simList = simList20,
  id_var = id,
  cores = 2
)

sampleSim20 %>%
  left_join(stands) %>%
  mutate(pi_i_approx = 20 * prop) %>%
  ggplot(aes(x = pi_i_approx, y = pi_i)) +
  geom_point() +
  geom_abline() +
  xlim(0, 1.2 * max(sampleSim20$pi_i)) +
  ylim(0, 1.2 * max(sampleSim20$pi_i)) +
  xlab(bquote("approximate" ~ pi[i])) +
  ylab(bquote("simulated" ~ pi[i]))

sampleSim20 %>%
  left_join(stands) %>%
  ggplot(
    aes(
      x = acres,
      y = cuftPerAc,
      color = pi_i)
    ) +
  geom_point() +
  ylab(bquote("volume (" ~ ft^{3} ~ "/ac)")) +
  scale_color_continuous(
    name = bquote("simulated" ~ pi[i])
  )

sampleSim20 %>%
  left_join(stands) %>%
  ggplot(
    aes(
      x = pi_i,
      y = cuftPerAc * acres
    )
  ) +
  geom_point() +
  ylab(bquote("total volume (" ~ ft^{3} ~ ")")) +
  xlab(bquote("simulated" ~ pi[i]))

Measure plots

In each sample stand we are installing a systematic random sample of 1/15 acre circular plots. On each plot we are estimating the cubic foot volume of the living trees. This is the second stage of our two-stage sample. We are simulating plots using a normal distribution based on the true mean and coefficient of variation of each stand.

plotTab <- sampStands %>%
  rowwise() %>%
  mutate(
    cuftObs = list(
      rnorm(
        n = nPlots,
        mean = cuftPerAc,
        sd = (cv / 100) * cuftPerAc
      )
    )
  ) %>%
  select(id, cuftObs) %>%
  unnest() %>%
  ungroup()

Equations

The Horvitz-Thompson estimator for the population total is: \[ \hat{t}_{HT} = \sum_{i \in S} \displaystyle \frac{\hat{t}_{i}}{\pi_{i}} \] where \(S\) is the sample of primary sample units from population sized \(N\)

The Sen-Yates-Grundy estimator for variance of the population total is:

\[ \hat{V}(\hat{t}_{HT}) = \sum_{i \in S} \sum_{\substack{k \in S \\ k > i}} \displaystyle \frac{\pi_{ik} - \pi_{i}\pi_{k}}{\pi_{ik}} \left(\frac{\hat{t}_{i}}{\pi_{i}} - \frac{\hat{t}_{k}}{\pi_{k}}\right)^{2}+ \sum_{i \in S} \displaystyle \frac{\hat{V}(\hat{t_{i}})}{\pi_{i}} \]

where \(\pi_{ik}\) is the joint inclusion probability for primary sampling units \(i\) and \(k\) and \(\hat{V}(\hat{t_{i}})\) is the variance of the estimate for the \(i\)th sample unit:

\[ \hat{V}(\hat{t_{i}}) = \displaystyle \frac{\sum_{j=1}^{n_{j}}(\hat{t}_{ij} - \hat{t}_{i})^2}{n_{i}(n_{i} - 1)} \]

where \(n_{i}\) is the number of sample observations in primary sample unit \(i\) and \(\hat{t}_{i}\) is the estimate of the total for the \(i\)th sample unit.

\(\hat{t}_{ij}\) is the \(j\)th estimate of the total on the \(i\)th primary sample unit. In this case, we observed the cubic feet per acre on each plot then multiplied by stand acres to get total cubic feet.

The Sen-Yates-Grundy formulation of the variance calculation requires each observation to be compared to the rest of the observations in the sample. There may be a more elegant way to code that part, but for now this is what I have come up with:

# Sen-Yates-Grundy formulation
var_syg <- function(row, data, y_var, pi_i_var, id_var, jointProbs) {
  
  data$idx <- 1:nrow(data)
  i <- data[row, ]
  b <- data[data$idx > row, ]
  n <- nrow(data)
  
  if(row == max(data$idx)) {
    return(0)
  }
  
  r <- as.character(i %>% pull(!! id_var))
  pi_i <- i %>% pull(!! pi_i_var)
  y_i <- i %>% pull(!! y_var)
  
  x <- c()
  for (l in 1:nrow(b)) {
    k <- b[l, ]
    pi_k <- k %>% pull(!! pi_i_var)
    y_k <- k %>% pull(!! y_var)
    
    c <- as.character(k %>% pull(!! id_var))
    jp <- jointProbs[row.names(jointProbs) == r, c]

    x[l] <- ((pi_i * pi_k - jp) / jp) *
      (y_i / pi_i - y_k / pi_k)^2
  }
  
  sum(x)
}

var_syg_wrap <- function(data, y_var, pi_i_var, id_var, jointProbs) {
  unlist(
    purrr::map(
      1:nrow(data),
      var_syg,
      data = data,
      y_var = enquo(y_var),
      pi_i_var = enquo(pi_i_var),
      id_var = enquo(id_var),
      jointProbs = jointProbs
    )
  )
}

Stand-level estimates

We can start by summarizing the stand-level estimates. For each stand we will estimate the total volume using the mean volume per acre from the sample and the total acres of each stand. We will also compute the sample variance and estimate the total variance for each stand.

standDat <- plotTab %>%
  left_join(
    sampStands %>% select(id, acres, cuftPerAc, nPlots, pi_i),
    by = "id"
  ) %>%
  group_by(id) %>%
  summarize(
    tHat = mean(cuftObs * acres),
    varTHat = var(cuftObs * acres) / n(),
    nPlots = n()
  )

standDat %>%
  left_join(stands) %>%
  mutate(se = sqrt(varTHat)) %>%
  ggplot(aes(x = cuftTot, y = tHat)) +
  geom_point() +
  geom_errorbar(aes(ymin = tHat - 1.7 * se, ymax = tHat + 1.7 * se)) +
  geom_abline() +
  xlab(bquote("true total volume (" ~ ft^{3} ~ ")"))+
  ylab(bquote("estimated total volume (" ~ ft^{3} ~ ")"))

Population-level estimate

We can now combine the estimates from the primary sampling units to obtain an estimate of the population. To do this we will use the inclusion and selection probabilities to weight each sampled stand’s contribution to the population total.

Using the equations listed above we will compute estimates of the population total and sampling variance of that estimate, then convert those into estimates of cubic foot volume on a per-acre basis.

sampleSummary <- standDat %>%
  left_join(
    sampStands %>%
      select(id, pi_i, prop),
    by = "id"
  ) %>%
  ungroup() %>%
  nest() %>%
  mutate(
    var_I_syg = map(
      .x = .$data,
      .f = var_syg_wrap,
      y_var = tHat,
      pi_i_var = pi_i,
      id_var = id,
      jointProbs = jointProbs
    )
  ) %>%
  unnest() %>%
  summarize(
    totalCuft = sum(tHat / pi_i),
    fpcI = 1,
    fpcII = 1,
    varI = sum(var_I_syg),
    varII = sum(varTHat / pi_i),
    varTot = fpcI * varI + fpcII * varII,
    nObs = n()
  ) %>%
  mutate(
    seTot = sqrt(varTot),
    ci90Tot = qt(
      1 - 0.1 / 2,
      df = nObs - 1
      ) * seTot,
    meanCuft = totalCuft / totalAcres,
    seMeanCuft = seTot / totalAcres,
    ci90MeanCuft = ci90Tot / totalAcres,
    source = "estimate"
  ) %>%
  select(
    source, totalCuft, seTot, ci90Tot,
    meanCuft, seMeanCuft, ci90MeanCuft,
    nObs
  )

results <- bind_rows(sampleSummary, truth)

results %>%
  ggplot(aes(x = source, y = meanCuft)) +
  geom_point() +
  geom_errorbar(
    aes(
      ymax = meanCuft + ci90MeanCuft,
      ymin = meanCuft - ci90MeanCuft
    )
  ) +
  ylim(0, max((results$meanCuft + results$ci90MeanCuft)) * 1.5) +
  ylab(bquote("volume (" ~ ft^{3} ~ "/ac)"))

## Warning: Removed 1 rows containing missing values (geom_errorbar).