-
Notifications
You must be signed in to change notification settings - Fork 12
Description
I'm trying to use FIESTA to compare post-stratified estimates of forest area and live biomass using custom stratification layers and custom estimation units, and then I'm comparing those to the results from using a simple random sample (no stratification). I'm running into an issue with modGBratio where I think it's using an incorrect denominator (estd) in the latest version of FIESTA (3.7.1.9100). I think the latest official release (3.7.1) has another issue with getting the correct forest area when using a simple random sample, so I don't want to use that version either.
I created an example that calculates live biomass density and forest area with post-stratification and with a simple random sample. It pulls out the denominator in modGBratio and the forest area estimate for comparison.
library(FIESTA)
# Setup stratification
WYspplt <- spMakeSpatialPoints(
xyplt = FIESTA::WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269
)
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package = "FIESTA")
WYbhfn <- system.file("extdata", "sp_data/WYbighorn_adminbnd.shp", package = "FIESTA")
strata_obj <- spGetStrata(
WYspplt,
uniqueid = "CN",
unit_layer = WYbhfn,
strattype = "RASTER",
strat_layer = fornffn
)
# Using post-stratification (strata=TRUE) and SRS (strata=FALSE) for forest area
# and live aboveground biomass density on forest land
for (strata in c(TRUE, FALSE)){
GBpop <- modGBpop(
popTabs = list(cond = FIESTA::WYcond, tree = FIESTA::WYtree),
stratdat = strata_obj,
strata = strata
)
ratio_res <- modGBratio(
GBpopdat = GBpop,
ratiotype = "PERACRE",
landarea = "FOREST",
estvarn = "DRYBIO_AG",
estvarn.filter = "STATUSCD == 1",
sumunits = TRUE
)
area_res <- modGBarea(
GBpopdat = GBpop,
landarea = "FOREST",
sumunits = TRUE
)
if (strata) {
cat("-------------Post-stratified estimates-------------\n")
} else{
cat("-------------Simple random sample------------------\n")
}
agb_est <- ratio_res$est$Estimate
agb_units <- ratio_res$raw$estunitsn
cat("AGB Estimate: ", agb_est, agb_units, "\n")
forest_area <- area_res$est$Estimate[1]
cat("Forest Area (modGBarea): ",
formatC(forest_area, format="d", big.mark=","), "acres\n")
ratio_denom <- ratio_res$raw$totest$estd
cat("Denominator in modGBratio totest: ",
formatC(ratio_denom, format="d", big.mark=","), "acres\n")
# The actual total land area of the estimation units
total_area <- sum(strata_obj$unitarea$ACRES)
cat("Total Unit Area (Unfiltered Frame): ",
formatC(total_area, format="d", big.mark=","), "acres\n")
}
Here are the results for different versions of FIESTA --
FIESTA version 3.7.1.9100 results:
-------------Post-stratified estimates-------------
AGB Estimate: 25 pounds
Forest Area (modGBarea): 665,942 acres
Denominator in modGBratio totest: 1,117,356 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
-------------Simple random sample------------------
AGB Estimate: 25.1 pounds
Forest Area (modGBarea): 668,438 acres
Denominator in modGBratio totest: 1,117,389 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
FIESTA version 3.7.1 results:
-------------Post-stratified estimates-------------
AGB Estimate: 41.9 pounds
Forest Area (modGBarea): 665,942 acres
Denominator in modGBratio totest: 665,942 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
-------------Simple random sample------------------
AGB Estimate: 41.9 pounds
Forest Area (modGBarea): 1,336,877 acres
Denominator in modGBratio totest: 1,336,877 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
FIESTA version 3.6.4 results:
-------------Post-stratified estimates-------------
AGB Estimate: 41.9
Forest Area (modGBarea): 665,942 acres
Denominator in modGBratio totest: 665,942 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
-------------Simple random sample------------------
AGB Estimate: 41.9
Forest Area (modGBarea): 668,438 acres
Denominator in modGBratio totest: 668,438 acres
Total Unit Area (Unfiltered Frame): 1,112,401 acres
I believe that the results from 3.6.4 are the only ones that are correct. I am using 3.6.4 currently but I run into another issue with collapsing strata so I would prefer to use one of the newer versions if this denominator issue can be resolved. Another aspect of these results that confuses me is why the estimate is correct in 3.7.1 with a simple random sample when the denominator is wrong. Also, $raw$estunitsn is saying the units are pounds in 3.7.1.9100 when I think they should be short tons.