Skip to content

Incorrect denominator in modGBratio #57

@FevetS

Description

@FevetS

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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions