How can we use the interpolated data produced by forestTIME to get population-level (i.e. state-level) per-area estimates?
Warning
The interpolated values produced by
forestTIMEare inferences and not samples, so it may not be appropriate to use (modified) design-based estimators, as we do in this vignette, which treat the data, including interpolated values, as a probability sample.
Example data
We’ll use RI as an example state because it is small. We’ll compare our estimates with rFIA, adapting examples in the FIA demystified vignette to our annualized data. We’ll use the built-in example data that comes with forestTIME, but if you were to download your own and you wanted to use it with rFIA also, make sure to set extract = "rFIA" to extract all the tables needed for both forestTIME and rFIA.
fia_download(states = "RI", download_dir = "fia", extract = "rFIA")Important
rFIAproduces estimates of carbon from 33 – 40.7 tons/acre using design-based estimators. “Correct” estimates should be in this ballpark.rfia_RI <- readFIA( dir = system.file("exdata", package = "forestTIME"), states = "RI" ) agc_rfia_annual <- biomass( rfia_RI, totals = TRUE, method = "annual", treeType = "live", landType = 'forest', component = "AG", areaDomain = COND_STATUS_CD == 1) |> mutate(method = "rFIA annual") |> select(method, YEAR, carbon_ton_acre = CARB_ACRE, carbon_total = CARB_TOTAL) agc_rfia_ti <- biomass( rfia_RI, totals = TRUE, method = "TI", treeType = "live", landType = 'forest', component = "AG", areaDomain = COND_STATUS_CD == 1) |> mutate(method = "rFIA TI") |> select(method, YEAR, carbon_ton_acre = CARB_ACRE, carbon_total = CARB_TOTAL) mean(agc_rfia_annual$carbon_ton_acre) mean(agc_rfia_ti$carbon_ton_acre)
We’ll use the standard basic workflow to get estimated aboveground carbon for each tree in each year.
# Data prep
db <- fia_load(
"RI",
dir = system.file("exdata", package = "forestTIME")
)
data <- fia_tidy(db) #single tibble
# Expand to include all years between surveys, interpolate/extrapolate, and adjust for mortality
# Assign FIA design-based evaluation information to annualized data
# Estimate carbon
data_midpt <- data |>
fia_annualize(use_mortyr = FALSE) |>
fia_design(db) |>
fia_allometry()Domain Indicator
We’ll add domain indicator columns as is done in the rFIA demystified vignette following equation 4.1, pg. 47 of Bechtold and Patterson(2005). These will be used to select which trees/plots/conditions fall within our population of interest. In this example we set the domain indicator for live trees per acre of forested land.
Reason:
We build separate domain indicators for estimating tree totals and area totals, because we can specify different domains of interest for both. For example, if we used our tree domain (live trees on forest land) to estimate area, then we would not actually be estimating the full forested area in RI. Instead we would estimate the forested area ONLY where live trees are currently present.
So we can’t just filter(STATUSCD == 1 & COND_STATUSCD == 1) to estimate carbon tons/acre.
Without Sampling Error - Using Expansion Factors
Expansion factors are weights that scale tree-level measurements to population-level totals. They were widely used in the older periodic FIA sampling design because they made estimation simple. However, the expansion factors provided in FIADB (POP_STRATUM.EXPNS) do not support valid variance estimation, and in the annual panel system they become unstable because each panel combination produces different weights. They are still widely used so we will demonstrate how to use them with our annualized data set. See Bechtold and Patterson (2005) section 4.4 for more information.
The FIA-supplied EXPNS values are no longer appropriate when using forestTIME because annualization changes the number of plots with “measurements” in each individual year. We therefore compute a new expansion factor (EXPNS_ft) in fia_design using FIA’s post stratified evaluation information with our updated plot counts.
We think these EXPNS values can be used much in the same way as the ones in the “raw” FIA database.
Using these expansion factors, we can follow the methods in the FIA demystified vignette.
tree_totals <- data_midpt |>
group_by(plot_ID, YEAR) |>
summarize(
# purposefully omits adjustment factor `aAdj` because it is assumed to be 1
carbPlot = sum(CARBON_AG_ft * TPA_UNADJ * EXPNS_ft * tDI / 2000, na.rm = TRUE), #tons/plot
)
area_totals <- data_midpt |>
group_by(plot_ID, YEAR) |>
# Keep only one row for each condition in each plot and year
distinct(CONDID, COND_STATUS_CD, CONDPROP_UNADJ, EXPNS_ft, aDI) |>
summarize(
# purposefully omits adjustment factor `aAdj` because it is assumed to be 1
forArea = sum(CONDPROP_UNADJ * EXPNS_ft * aDI, na.rm = TRUE) #acres/plot
)
agc_pop <- inner_join(tree_totals, area_totals) |>
group_by(YEAR) |>
summarize(
CARB_AG_TOTAL = sum(carbPlot, na.rm = TRUE), # tons/plot
AREA_TOTAL = sum(forArea, na.rm = TRUE) # acres/plot
) |>
# the units work out to be tons(live carbon)/acre(forested land)
mutate(method = "forestTIME", carbon_ton_acre = CARB_AG_TOTAL / AREA_TOTAL) |>
select(
method,
YEAR,
carbon_ton_acre,
carbon_total = CARB_AG_TOTAL,
AREA_TOTAL
)
agc_popWe have some ballpark similar estimates to the ones rFIA produces.
With Sampling Error
Note: Small strata merging has not yet been implemented in the forestTIME workflow so some strata may contain too few plots to compute reliable variance estimates. Small strata merging is currently being developed for forestTIME. Disregard if only interested in totals and ratio estimates.
If we want to compute sampling error using FIA’s design-based variance estimator, we cannot use EXPNS. Instead, we need to use post-stratification where the landscape is divided into estimation units (for example, private land, public land, water) and then into smaller strata (for example, tree canopy cover) using remote-sensing data. Each pixel is assigned to a stratum, and the proportion of area in each stratum becomes its stratum weight. Phase 2 ground plots are then linked to their stratum (P2POINTCNT), and plot-level measurements are averaged within each stratum. The stratum means are then combined using the stratum weights first at the strata level and then at the estimation unit level to produced design-based population estimates. This stratified information changes from year to year as the land is resurveyed and this is reflected in the EVALID.
The fia_design() function attaches all of the necessary design information to the annualized data so that population estimates can be produced using FIA’s post-stratified estimators. forestTIME attempts to match the data year to the same evaluation year however, some plots belong to non-matching evaluation years if a). no evaluation was done that year or b). the plot never belonged to that evaluation.
# Estimate Tree attributes -- plot level
tree_pop <- data_midpt |>
filter(EXPVOL == TRUE) |>
distinct(ESTN_UNIT_CN, STRATUM_CN, CONDID, plot_ID, tree_ID, YEAR, .keep_all = TRUE) |>
group_by(YEAR, ESTN_UNIT_CN, STRATUM_CN, plot_ID, EVALID) |>
summarize(
carbPlot = sum(CARBON_AG_ft * TPA_UNADJ * tDI / 2000, na.rm = TRUE),
.groups = "drop"
)
# Estimate Area attributes -- plot level
area_pop <- data_midpt |>
filter(EXPCURR == TRUE) |> # use EXPCURR flag
distinct(ESTN_UNIT_CN, STRATUM_CN, CONDID, plot_ID, YEAR, .keep_all = TRUE) |>
group_by(YEAR, P1POINTCNT, P1PNTCNT_EU, P2POINTCNT_ft, AREA_USED,
ESTN_UNIT_CN, STRATUM_CN, plot_ID, EVALID) |>
summarize(
forArea = sum(CONDPROP_UNADJ * aDI, na.rm = TRUE),
.groups = "drop"
)
# Joining the two tables
carb_pop_plot <- left_join(tree_pop, area_pop)
# Strata level estimates
carb_pop_strat <- carb_pop_plot |>
group_by(YEAR, ESTN_UNIT_CN, STRATUM_CN, EVALID) |>
summarize(aStrat = mean(forArea, na.rm = TRUE), # Area mean
carbStrat = mean(carbPlot, na.rm = TRUE), # Carbon mean
# We don't want a vector of these values, since they are repeated
P2POINTCNT = first(P2POINTCNT_ft),
AREA_USED = first(AREA_USED),
# Strata weight, relative to estimation unit
w = first(P1POINTCNT) / first(P1PNTCNT_EU),
# Strata level variances
aVar = (sum(forArea^2) - sum(P2POINTCNT * aStrat^2)) /
(P2POINTCNT * (P2POINTCNT-1)),
carbVar = (sum(carbPlot^2) - sum(P2POINTCNT * carbStrat^2)) /
(P2POINTCNT * (P2POINTCNT-1)),
# Strata level co-varainces
carbCV = (sum(forArea*carbPlot) - sum(P2POINTCNT * aStrat *carbStrat)) /
(P2POINTCNT * (P2POINTCNT-1)),
.groups = 'drop')
# Moving on to the estimation unit estimates
carb_pop_eu <- carb_pop_strat |>
group_by(YEAR, ESTN_UNIT_CN, EVALID) |>
summarize(aEst = sum(aStrat * w, na.rm = TRUE) * first(AREA_USED), # Mean Area
carbEst = sum(carbStrat * w, na.rm = TRUE) * first(AREA_USED), # Mean carbon
# Estimation unit variances
aVar = (first(AREA_USED)^2 / sum(P2POINTCNT)) *
(sum(P2POINTCNT*w*aVar) + sum((1-w)*(P2POINTCNT/sum(P2POINTCNT))*aVar)),
carbVar = (first(AREA_USED)^2 / sum(P2POINTCNT)) *
(sum(P2POINTCNT*w*carbVar) + sum((1-w)*(P2POINTCNT/sum(P2POINTCNT))*carbVar)),
# Estimation unit covariances
carbCV = (first(AREA_USED)^2 / sum(P2POINTCNT)) *
(sum(P2POINTCNT*w*carbCV) + sum((1-w)*(P2POINTCNT/sum(P2POINTCNT))*carbCV)),
.groups = 'drop')
# sum across Estimation Units for totals
carb_pop <- carb_pop_eu |>
group_by(YEAR) |>
summarize(AREA_TOTAL = sum(aEst, na.rm = TRUE),
CARB_TOTAL = sum(carbEst, na.rm = TRUE),
# Ratios
CARB_ACRE = CARB_TOTAL / AREA_TOTAL,
# Total samping errors
AREA_TOTAL_SE = sqrt(sum(aVar, na.rm = TRUE)) / AREA_TOTAL * 100,
CARB_TOTAL_SE = sqrt(sum(carbVar, na.rm = TRUE)) / CARB_TOTAL * 100,
# Ratio variances
carbAcreVar = (1 / AREA_TOTAL^2) * (sum(carbVar) +
(CARB_ACRE^2) * sum(aVar) -
2 * CARB_ACRE * sum(carbCV)),
CARB_ACRE_SE = sqrt(sum(carbAcreVar, na.rm = TRUE)) /
CARB_ACRE * 100) |>
mutate(method = "forestTIME") |>
# Re ordering
select(YEAR,
carbon_ton_acre = CARB_ACRE,
carbon_total = CARB_TOTAL,
AREA_TOTAL,
method)
mean(carb_pop$carbon_ton_acre)