---
title: "Green-book Estimators"
output: rmarkdown::html_vignette
description: >
  Learn how calculate population estimates and their sampling errors using the sample-based estimators from Bechtold and Patterson's (2005), 'Green-Book'.
vignette: >
  %\VignetteIndexEntry{Green-book Estimators}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include = F}
library(knitr)
knitr::opts_chunk$set(message = F, warning = F)
```

```{r, include=FALSE}
# Sets up output folding
hooks = knitr::knit_hooks$get()
hook_foldable = function(type) {
  force(type)
  function(x, options) {
    res = hooks[[type]](x, options)
    
    if (isFALSE(options[[paste0("fold.", type)]])) return(res)
    
    paste0(
      "<details><summary>", type, "</summary>\n\n",
      res,
      "\n\n</details>"
    )
  }
}
knitr::knit_hooks$set(
  output = hook_foldable("output"),
  plot = hook_foldable("plot")
)
```

```{r, echo=-1}
data.table::setDTthreads(2)
```


## Green-Book (GB) module overview

`FIESTA`'s Green-Book (GB) module calculates population estimates and their sampling errors based on Bechtold and Patterson's (2005), 'Green-Book' for FIA's nationally-consistent, systematic annual sample design, chapter 4 (Scott et al. 2005). FIA's sample design is based on 2 phases: the first phase uses remotely-sensed data to stratify the land area to increase precision of estimates; while the 2nd phase obtains photo and ground observations and measurements for a suite of information across a hexagonal grid, each approximately 6000 acres in size. The associated estimators and variance estimators are used for area and tree attribute totals with the assumption of a simple random, stratified design and double sampling for stratification. Adjustment factors are calculated by estimation unit and strata to account for nonsampled (nonresponse) conditions.   

Functions include non-ratio estimators for area and tree estimates by domain and ratio-of-means estimators for per-acre and per-tree estimates within domain. In addition, `FIESTA` adjusts for nonsampled conditions, supports post-stratification for reducing variance, and reports by estimation unit or a summed combination of estimation units. Output from the Green-Book module was tested and compared to output from FIA's publicly-available online tool ([EVALIDator](https://apps.fs.usda.gov/fiadb-api/evalidator)) for state-level population estimates and associated sampling errors generated from the FIA Database (FIADB).

## Objective of tutorial

The Green-Book estimators can be used with FIA's standard state-level population data (i.e, Evaluation) from the FIA database (FIADB) and also population data from a custom boundary. The population data includes a set of FIA plot data and summarized auxiliary information for post-stratification, including a table of area by estimation unit within the population, and a table of strata proportions by estimation unit. This tutorial steps through several examples using FIESTA's Green Book module, for four different populations: (POP1) an FIA standard Evaluation, Wyoming 561301; (POP2) a custom boundary with one population, Bighorn National Forest; (POP3) a custom boundary with sub-populations, Bighorn National Forest Districts; and (POP4) an FIA standard Evaluation, Rhode Island 441901, stored in a SQLite database. All examples can be used with any population, standard or custom.



## GB Examples
<a name="GBfun"></a>

### GB Example Data
<details>
  <summary>View GB Example Data</summary>

*Example FIA plot data from FIADB*

The examples use FIA plot data from FIA Evaluation 561301, including three inventory years of field measurements in the state of Wyoming, from FIADB_1.7.2.00, last updated June 20, 2018, downloaded on June 25, 2018, and stored as internal data objects in `FIESTA`. 



**Wyoming (WY), Inventory Years 2011-2013 (Evaluation 561301)**

Data Frame  | Description
:-----------| :-----------------------------------------------------------------------------------
WYplt       | WY plot-level data
WYcond      | WY condition-level data
WYtree      | WY tree-level data


*Example Auxiliary data*

Auxiliary data for state-level estimates, including plot-level estimation unit and stratum assignments; area by estimation unit; and pixel counts by strata class and estimation unit, were downloaded from FIADB at the same time, from the same FIA Evaluation (i.e., 561301), and stored as internal data objects in `FIESTA`. Estimates using auxiliary data from FIADB can be compared with EVALIDator estimates, using the 2013 evaluation (https://apps.fs.usda.gov/fiadb-api/evalidator).

Auxiliary data for the custom boundaries are summarized from spatial layers stored as external objects in FIESTA, originating from the USDA Forest Service, Automated Lands Program (ALP; 2018) and from a 250m resolution, Moderate Resolution Imaging Spectroradiometer (MODIS), classified map, reclassified from 3 to 2 classes: 1:forest; 2:nonforest (Ruefenacht et al. 2008)



**Wyoming (WY), Auxiliary data from FIADB (Evaluation 561301)**

Data Frame  | Description
:-----------| :-----------------------------------------------------------------------------------
WYpltassgn  | WY plot-level data with strata and estimation unit assignments
WYunitarea  | WY estimation unit look-up table with total acres by estimation unit (ESTUNIT)
WYstratalut | WY strata look-up table with pixel counts (P1POINTCNT) by strata and estimation unit



**Wyoming (WY), Auxiliary data from other sources**

External data             | Description
:-------------------------| :------------------------------------------------------------------------
WYbighorn_adminbnd.shp    | Polygon shapefile of WY Bighorn National Forest Administrative boundary^1^
WYbighorn_districtbnd.shp | Polygon shapefile of WY Bighorn National Forest District boundaries^2^
WYbighorn_forest_nonforest_250m.tif| GeoTIFF raster of predicted forest/nonforest (1/0)^3^ 

^1^USDA Forest Service, Automated Lands Program (ALP). 2018. S_USA.AdministrativeForest (http://data.fs.usda.gov/geodata/edw). Description: An area encompassing all the National Forest System lands administered by an administrative unit. The area encompasses private lands, other governmental agency lands, and may contain National Forest System lands within the proclaimed boundaries of another administrative unit. All National Forest System lands fall within one and only one Administrative Forest Area.

^2^USDA Forest Service, Automated Lands Program (ALP). 2018. S_USA.RangerDistrict (http://data.fs.usda.gov/geodata/edw). Description: A depiction of the boundary that encompareasses a Ranger District.

^3^Based on 250m resolution, Moderate Resolution Imaging Spectroradiometer (MODIS), classified map, reclassified from 3 to 2 classes: 1:forest; 2:nonforest. Projected in Albers Conical Equal Area, Datum NAD27 (Ruefenacht et al. 2008).

</details>


### Set up 

First, you'll need to load the `FIESTA` library:

```{r, warning = F, message = F}
library(FIESTA)
```

Next, you'll need to set up an "outfolder". This is just a file path to a folder where you'd like `FIESTA` to send your data output. For this vignette, we have set our outfolder file path as a temporary directory. 

```{r}
outfolder <- tempdir()
```


### Get auxiliary data for custom examples

Now, we need to get the auxiliary data for the custom boundaries. The FIESTA *spGetStrata* function is a spatial wrapper function to facilitate extraction and summary of user-defined spatial data used for post-stratification. The function uses the FIESTA *spExtractPoly* and *spExtractRast* functions to subset (i.e., clip) plots to the boundary and extract values from estimation unit (i.e., polygon) and strata values (i.e., raster) to plot center locations, respectively. Other internal spatial functions calculate stratum pixel counts and area by estimation unit. If a polygon strata layer is given, the FIESTA *spPoly2Rast* function converts the polygon layer to raster before calculating strata weights. 

Our custom examples demonstrate how to get data for one area of interest, or population (e.g, Bighorn National Forest) and for one area of interest, with multiple estimation units, or subpopulations (e.g., Bighorn National Forest Districts). 



**Bighorn National Forest**


<details>
  <summary>View Getting Strata Data</summary>
  
```{r}

# File names for spatial layers, stored as external data objects in FIESTA. 
WYbhfn <- system.file("extdata", "sp_data/WYbighorn_adminbnd.shp", package = "FIESTA")
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package = "FIESTA")


# Get estimation unit and strata information for Bighorn National Forest.
stratdat.bh <- spGetStrata(xyplt = WYplt,
                           uniqueid = "CN", 
                           unit_layer = WYbhfn, 
                           strat_layer = fornffn,
                           spMakeSpatial_opts = list(xvar = "LON_PUBLIC", 
                                                     yvar = "LAT_PUBLIC", 
                                                     xy.crs = 4269))

# Get names of output components
names(stratdat.bh)

# Plot assignment of strata and estimation unit (ONEUNIT, STRATUMCD)
head(stratdat.bh$pltassgn)

# Area by estimation unit
stratdat.bh$unitarea

#  Pixel counts and strata weights (strwt) by strata and estimation unit
stratdat.bh$stratalut

# Variable names
stratdat.bh$unitvar        # Estimation unit variable
stratdat.bh$strvar         # Strata variable
stratdat.bh$areavar        # Area variable


```

</details>

**Bighorn National Forest Districts**

<details>
  <summary>View Getting Strata Data (Districts)</summary>
  
```{r}
# File names for external spatial data 
WYbhdistfn <- system.file("extdata", "sp_data/WYbighorn_districtbnd.shp", package = "FIESTA")
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif", package = "FIESTA")

# Get estimation unit and strata information for Bighorn National Forest Districts
stratdat.bhdist <- spGetStrata(xyplt = WYplt,
                               uniqueid = "CN", 
                               unit_layer = WYbhdistfn, 
                               unitvar = "DISTRICTNA",
                               strat_layer = fornffn,
                               spMakeSpatial_opts = list(xvar = "LON_PUBLIC", 
                                                         yvar = "LAT_PUBLIC", 
                                                         xy.crs = 4269))

# Get names of output list components
names(stratdat.bhdist)

# Plot assignment of strata and estimation unit (DISTRICTNA, STRATUMCD)
head(stratdat.bhdist$pltassgn)

# Area by estimation units (Districts)
stratdat.bhdist$unitarea

# Pixel counts and strata weights (strwt) by strata and estimation unit
stratdat.bhdist$stratalut


# Variable names
stratdat.bhdist$unitvar        # Estimation unit variable
stratdat.bhdist$strvar         # Strata variable
stratdat.bhdist$areavar        # Area variable

```

</details>

### GB Module
* [modGBpop()](#modGBpop)
* [modGBarea()](#modGBarea)
* [modGBtree()](#modGBtree)
* [modGBratio()](#modGBratio)


<a name="modGBpop"></a>

### modGBpop()

`FIESTA`'s population functions (`mod*pop`) check input data and perform population-level calculations, such as: summing number of sampled plots; adjusting for partial nonresponse; and standardizing auxiliary data. These functions are specific to each `FIESTA` module and are run prior to or within a module for any population of interest. 

For `FIESTA`'s GB Module, the `modGBpop` function calculates and outputs: number of plots, adjustment factors, and an expansion factor by strata. The outputs are similar to data found in FIADB's pop_stratum table. The output from `modGBpop` can be used for one or more estimates from `modGBarea`, `modGBtree`, or `modGBratio` functions.


#### POP1: FIADB POPULATION - Get population data for post-stratified area and tree estimates for Wyoming

<details>
  <summary>View Example</summary>
  
In this example, we use the sample Wyoming data (2013 Evaluation) stored in `FIESTA` to generate population data for the GB module. We check this output with the FIADB pop_stratum table from FIA DataMart for 561301 Evalid, using the `FIESTA::DBqryCSV` function.


```{r, results = FALSE, message = FALSE}
GBpopdat <- modGBpop(popTabs = list(cond = FIESTA::WYcond,          # FIA plot/condition data
                                    tree = FIESTA::WYtree,          # FIA tree data
                                    seed = FIESTA::WYseed),         # FIA seedling data
                     popTabIDs = list(cond = "PLT_CN"),             # unique ID of plot in cond
                     pltassgn = FIESTA::WYpltassgn,                 # plot assignments
                     pltassgnid = "CN",                             # unique ID of plot in pltassgn
                     pjoinid = "PLT_CN",                            # plot id to join to pltassgn
                     unitarea = FIESTA::WYunitarea,                 # area by estimation units
                     unitvar = "ESTN_UNIT",                         # name of estimation unit variable
                     strata = TRUE,                                 # if using post-stratification
                     stratalut = FIESTA::WYstratalut,               # strata classes and pixels counts
                     strata_opts = strata_options(getwt = TRUE))    # strata options

```

To get the names of the list components associated with the output of our call of `modGBpop`, we run the following code:

```{r, results = T}
names(GBpopdat)
```

From this list outputted by `GBpopdat` we can access many things. Some examples include the number of plots by plot status (`plotsampcnt`), the number of conditions by condition status (`condsampcnt`), the number of plots and adjustment factors by strata (`stratalut`), and the adjustment factors added to the plot-level, tree-level, and seedling data (`pltidsadj`, `treex`, and `seedx`). These objects can be seen below:

```{r}
# Number of plots by plot status
GBpopdat$plotsampcnt	

# Number of conditions by condition status
GBpopdat$condsampcnt

# Number of plots and adjustment factors by strata
GBpopdat$stratalut  

# Adjustment factors added to plot-level data
head(GBpopdat$pltidsadj)

# Adjustment factors added to tree data
head(GBpopdat$treex)

# Adjustment factors added to seedling data
head(GBpopdat$seedx)
```

One may also want to compare `FIESTA` output with FIADB pop_stratum table for WY in the 2013 evaluation to check for consistency. The can be done as follows: 

```{r, results = TRUE, messages = FALSE}
qry <- "SELECT ESTN_UNIT, STRATUMCD, P1POINTCNT, P2POINTCNT,
               ADJ_FACTOR_MACR, ADJ_FACTOR_SUBP, ADJ_FACTOR_MICR
        FROM POP_STRATUM
        WHERE EVALID = 561301
        ORDER BY ESTN_UNIT, STRATUMCD"

pop_stratum <- DBqryCSV(qry, 
                        states="Wyoming",
                        sqltables="POP_STRATUM")


head(pop_stratum)

head(GBpopdat$stratalut)
```

</details>


#### POP2: CUSTOM POPULATION - Get population data for post-stratified area and tree estimates for the Bighorn National Forest

<details>
  <summary>View Example</summary>
  
In this example, we use the sample WY plot data (2013 Evaluation) in FIESTA and output from `spGetStrata` to generate population data for the Bighorn National Forest. Here, we have only one estimation unit within the population of interest (Bighorn National Forest), therefore strata and pixel counts are summarized to the population.

If the `FIESTA::spGetStrata` function is used to obtain stratification data, the output list object can be input directly into `modGBpop` through the `stratdat` parameter.

```{r, message = FALSE}
GBpopdat.bh <- modGBpop(popTabs = popTables(plt = WYplt,
                                            cond = WYcond,
                                            tree = WYtree,
                                            seed = WYseed),
                        stratdat = stratdat.bh)

# Get names of output list components
names(GBpopdat.bh)
```

The utilization of the `stratdat` parameter shown above is simply a shortcut that is available when `spGetStrata()` was used. Alternatively, the data can be input through individual parameters as shown below.

```{r, message=FALSE}
## Using output as individual parameter inputs
GBpopdat.bh <- modGBpop(popTabs = popTables(plt = WYplt,
                                            cond = WYcond,
                                            tree = WYtree,
                                            seed = WYseed),
                        popTabIDs = popTableIDs(plt = "CN"),
                        pltassgn = stratdat.bh$pltassgn, 
                        pltassgnid = "CN", 
                        unitvar = stratdat.bh$unitvar, 
                        unitarea = stratdat.bh$unitarea, 
                        areavar = stratdat.bh$areavar, 
                        strata = TRUE,
                        stratalut  =stratdat.bh$stratalut, 
                        strvar = stratdat.bh$strvar)

## Get names of output list components
names(GBpopdat.bh)
```



</details>


</details>


#### POP3: CUSTOM SUB-POPULATIONS - Get sub-population data for area and tree estimates for the Bighorn National Forest Districts, using post-stratification

<details>
  <summary>View Example</summary>

In this example, we use the sample Wyoming plot data (2013 Evaluation) stored in FIESTA and output from `spGetStrata` to generate sub-population data for Bighorn National Forest Districts. Here, we have more than one estimation unit (i.e., sub-population) within the population of interest (i.e., Bighorn National Forest Districts), therefore strata and pixel counts are summarized by each District within the population.


```{r, message = FALSE}
# Bighorn National Forest District

# Using output list from spGetStrata()
GBpopdat.bhdist <- modGBpop(popTabs = popTables(plt = WYplt,
                                                cond = WYcond,
                                                tree = WYtree,
                                                seed = WYseed), 
                            stratdat = stratdat.bhdist)

## Get names of output list components
names(GBpopdat.bhdist)
```

</details>


#### POP4: FIADB POPULATION - Get population data for area and tree estimates for Rhode Island, using post-stratification, with data stored in a SQLite database

<details>
  <summary>View Example</summary>
  
In this example, we use the sample Rhode Island data (441901 Evaluation) stored in a SQLite database as external data in `FIESTA`. Data were extracted from the FIA database on June 6, 2022. All output can be compared with output from other FIA tools.

First, let's look at the SQLite database. Use the DBI package explore the contents.

```{r, results = FALSE, message = F}

SQLitefn <- system.file("extdata", "FIA_data/RIdat_eval2019.db", package="FIESTA")

conn <- DBI::dbConnect(RSQLite::SQLite(), SQLitefn)
DBI::dbListTables(conn)
DBI::dbDisconnect(conn)

```


```{r, results = FALSE, message = F}
GBpopdat.RI <- modGBpop(popTabs = popTables(plt = "plot",
                                            cond = "cond",
                                            tree = "tree",
                                            seed = "seed"),
                        dsn = SQLitefn,
                        pltassgn = "pop_plot_stratum_assgn",
                        stratalut = "pop_stratum",
                        unitarea = "pop_estn_unit",
                        unitvar = "ESTN_UNIT",
                        areavar = "AREA_USED",
                        strata_opts = strata_options(getwt = TRUE,
                                                     getwtvar = "P1POINTCNT"))

names(GBpopdat.RI)

# Strata-level population data, including number of plots and adjustment factors
GBpopdat.RI$stratalut  
```

</details>


<a name="modGBarea()"></a>

### modGBarea()

`FIESTA`'s `modGBarea` function generates acre estimates. Calculations are based on Scott et al. 2015 ('Green-Book') for mapped forest inventory plots. The non-ratio estimator for estimating area by stratum and domain is used. Plots that are totally nonsampled are excluded from the estimation dataset. Next, an adjustment factor is calculated by strata to adjust for nonsampled (nonresponse) conditions that have proportion less than 1. The attribute is the proportion of the plot which is divided by the adjustment factor, and averaged by stratum. Strata means are combined using the strata weights and then expanded to acres using the total land area in the population.

If there are more than one estimation unit (i.e., subpopulation) within the population, estimates are generated by estimation unit. If `sumunits=TRUE`, the estimates and percent standard errors returned are a sum combination of all estimation units. If `rawdata=TRUE`, the raw data returned will include estimates by estimation unit.  

Parameters defined in the following examples are organized by category: population data (pop); estimation information (est); and output details (out).

#### POP1: 1.1 Area of forest land, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

Using the `modGBarea` function we generate estimates by estimation unit (i.e., ESTN_UNIT) and sum to population (i.e., WY). `FIESTA` then returns raw data for area of forest land, Wyoming, 2011-2013 (sum estimation units). 

The following estimates match output from [EVALIDator](https://apps.fs.usda.gov/fiadb-api/evalidator) using the WY 2013 Evaluation. 

```{r}
area1.1 <- modGBarea(GBpopdat = GBpopdat,      
                     landarea = "FOREST",
                     sumunits = TRUE) 

```

To get the names of the list components associated with the output of our call of `modGBarea`, we run the following code:

```{r, results = T}
names(area1.1)
```

To easily access our estimate and percent sampling error of estimate we can just grab the `est` object from out outputted list: 

```{r, results = T}
area1.1$est
```

We can also look at raw data and estimates, as shown below:

```{r, results = TRUE}
raw1.1 <- area1.1$raw        # extract raw data list object from output
names(raw1.1)

head(raw1.1$unit_totest)    # estimates by estimation unit (i.e., ESTN_UNIT)
raw1.1$totest               # estimates for population (i.e., WY)
```

</details>

#### POP1: 1.2 Area by forest type on forest land, Wyoming, 2011-2013 
<details>
  <summary>View Example</summary>

In this example, we look at adding rows to the output and include `returntitle=TRUE` to return title information. We use the variable `rowvar` to specify a domain whose levels we want to calculate estimates by.

```{r}
## Area of forest land by forest type, Wyoming, 2011-2013
area1.2 <- modGBarea(GBpopdat = GBpopdat,        
                     landarea = "FOREST",         # est - forest land filter
                     rowvar = "FORTYPCD",         # est - row domain
                     sumunits = TRUE,             # est - sum estimation units to population
                     returntitle = TRUE)          # out - return title information
```

Again, we can look at the contents of the output list. The output now includes titlelst, a list of associated titles.

```{r}
names(area1.2)
```

And the estimates: 

```{r}
## Estimate and percent sampling error of estimate
area1.2$est
```

Along with raw data and titles:

```{r}
raw1.2 <- area1.2$raw      # extract raw data list object from output
names(raw1.2)

head(raw1.2$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
raw1.2$totest            # estimates for population (i.e., WY)
head(raw1.2$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest)      # estimates by row for population (i.e., WY)


## Titles (list object) for estimate
titlelst1.2 <- area1.2$titlelst
names(titlelst1.2)
titlelst1.2
```

</details>

#### POP1: 1.3 Area by forest type and stand-size class on forest land, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

In this example, we look at generating estimates by multiple domains. In order to enhance the readability of the output we use FIA names as labels. We also output estimates and percent standard error in the same cell with the `allin1` argument in `table_options` and save data to an outfolder with the `outfolder` argument in `savedata_options`. 

```{r, message = FALSE}
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013

area1.3 <- modGBarea(GBpopdat = GBpopdat,                                     # pop - population calculations for WY, post-stratification
                     landarea = "FOREST",                                     # est - forest land filter
                     rowvar = "FORTYPCD",                                     # est - row domain
                     colvar = "STDSZCD",                                      # est - column domain
                     sumunits = TRUE,                                         # est - sum estimation units to population
                     savedata = TRUE,                                         # out - save data to outfolder
                     returntitle = TRUE,                                      # out - return title information
                     table_opts = table_options(row.FIAname = TRUE,           # table - row domain names
                                                col.FIAname = TRUE,           # table - column domain names
                                                allin1 = TRUE),               # table - return output with est(pse)
                     savedata_opts = savedata_options(outfolder = outfolder,  # save - outfolder for saving data
                                                      outfn.pre = "WY"))      # save - prefix for output files
```

We can again look at the output list, estimates, raw data, and titles:

```{r}
# Look at output list
names(area1.3)

# Estimate and percent sampling error of estimate
head(area1.3$est)


# Raw data (list object) for estimate
raw1.3 <- area1.3$raw      # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest)      # estimates for population (i.e., WY)
head(raw1.3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest)      # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest)      # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest)      # estimates by row and column for population (i.e., WY)


## Titles (list object) for estimate
titlelst1.3 <- area1.3$titlelst
names(titlelst1.3)
titlelst1.3


# List output files in outfolder
list.files(outfolder, pattern = "WY_area")
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_area")
```


#### POP2: 2.1 Area by forest type and stand-size class, Bighorn National Forest
<details>
  <summary>View Example</summary>

We now shift to producing estimates for the Bighorn National Forest. In this example, we utilize the `title_opts` argument to customize the title outputs.

```{r, messages = FALSE}
area2.1 <- modGBarea(GBpopdat = GBpopdat.bh,       # pop - population calculations for Bighorn NF, post-stratification
                     landarea = "FOREST",          # est - forest land filter
                     rowvar = "FORTYPCD",          # est - row domain
                     colvar = "STDSZCD",           # est - column domain
                     returntitle = TRUE,           # out - return title information
                     title_opts = title_options(title.ref = "Bighorn National Forest, 2011-2013"),  # title - customize title reference
                     table_opts = table_options(row.FIAname = TRUE,         # table - return FIA row names
                                                col.FIAname = TRUE))        # table - return FIA column names
```

To get the names of the list components associated with the output of our call of `modGBarea`, we run the following code:

```{r, results = T}
names(area2.1)
```

To easily access our estimate and percent sampling error of estimate we can just grab the `est` object from our ouputted list: 

```{r, results = T}
area2.1$est
```

We can also look at raw data and titles for estimate, as shown below. Note the change in titles.

```{r, results = TRUE}
raw2.1 <- area2.1$raw      # extract raw data list object from output
names(raw2.1)
head(raw2.1$unit_grpest)  # estimates by row and group domains


# Titles (list object) for estimate
titlelst2.1 <- area2.1$titlelst
names(titlelst2.1)
titlelst2.1
```

</details>


#### POP2: 2.2 Area by forest type group and primary disturbance class, Bighorn National Forest
<details>
  <summary>View Example</summary>

Note: Since we only have one estimation unit within the population of interest, we set sumunits=FALSE. Let's also add a few more table options to control the forest type groups displayed in the table (rowlut) and fill the NULL values with 0s (estnull).

```{r}
area2.2 <- modGBarea(GBpopdat = GBpopdat.bh,        # pop - population calculations for Bighorn NF, post-stratification
                     landarea = "FOREST",           # est - forest land filter
                     sumunits = TRUE,               # est - sum estimation units to population
                     rowvar = "FORTYPGRPCD",        # est - row domain
                     colvar = "DSTRBCD1",           # est - column domain
                     returntitle = TRUE,            # out - return title information
                     title_opts = title_options(
                        title.ref = "Bighorn National Forest, 2011-2013"  # title - customize title reference
                     ),
                    table_opts = table_options(   
                      row.FIAname = TRUE,           # table - return FIA row names
                      col.FIAname = TRUE,           # table - return FIA column names
                      estnull = 0,
                      rowlut = c(180, 200, 220, 260, 280, 900, 999),
                      raw.keep0 = TRUE
                    ))

```

To get the names of the list components associated with the output of our call of `modGBarea`, we run the following code:

```{r, results = T}
names(area2.2)
```

To easily access our estimate and percent sampling error of estimate we can just grab the `est` object from our ouputted list: 

```{r, results = T}
area2.2$est
```

We can also look at raw data and titles for estimate, as shown below. Note the change in titles.

```{r, results = TRUE}
## Raw data (list object) for estimate
raw2.2 <- area2.2$raw      # extract raw data list object from output
names(raw2.2)
head(raw2.2$unit_grpest)  # estimates by row and group domains


## Titles (list object) for estimate
titlelst2.2 <- area2.2$titlelst
names(titlelst2.2)
titlelst2.2

```

</details>



#### POP3: 3.1 Area by forest type group and primary disturbance class, Bighorn National Forest Districts
<details>
  <summary>View Example</summary>

In this example, we add a filter to remove instances where there is no visible disturbance from the table output. This filter does not change the population data set the estimates are derived from, it only changes the output.

```{r}
area3.1 <- modGBarea(GBpopdat = GBpopdat.bhdist,    # pop - population calculations for Bighorn NF, post-stratification
                     landarea = "FOREST",           # est - forest land filter
                     sumunits = TRUE,               # est - sum estimation units to population
                     pcfilter = "DSTRBCD1 > 0",     # est - condition filter for table output
                     rowvar = "FORTYPGRPCD",        # est - row domain
                     colvar = "DSTRBCD1",           # est - column domain
                     returntitle = TRUE,            # out - return title information
                     title_opts = title_options(
                       title.ref = "Bighorn National Forest, 2011-2013"  # title - customize title reference
                     ),
                     table_opts = table_options(   
                       row.FIAname = TRUE,          # table - return FIA row names
                       col.FIAname = TRUE           # table - return FIA column names
                     ))

```

To get the names of the list components associated with the output of our call of `modGBarea`, we run the following code:

```{r, results = T}
names(area3.1)
```

To easily access our estimate and percent sampling error of estimate we can just grab the `est` object from our ouputted list: 

```{r, results = T}
area3.1$est
```

We can also look at raw data and titles for estimate, as shown below:

```{r, results = TRUE}
raw3.1 <- area3.1$raw       
names(raw3.1)
head(raw3.1$unit_rowest)   # estimates by estimation unit for row domains
raw3.1$rowest              # estimates for population for row domains

head(raw3.1$unit_colest)   # estimates by estimation unit for column domains
raw3.1$colest              # estimates for population for column domains


## Titles (list object) for estimate
titlelst3.1 <- area3.1$titlelst
names(titlelst3.1)
titlelst3.1
```

</details>


#### POP4: 4.1 Area by forest type group and stand-size class, Rhode Island, 2019
<details>
  <summary>View Example</summary>

Note: estimates should match other FIA tools.
```{r}
area4.1 <- modGBarea(GBpopdat = GBpopdat.RI,        # pop - population calculations for Bighorn NF, post-stratification
                     landarea = "FOREST",           # est - forest land filter
                     sumunits = TRUE,               # est - sum estimation units to population
                     rowvar = "FORTYPCD",           # est - row domain
                     colvar = "STDSZCD",            # est - column domain
                     returntitle = TRUE,            # out - return title information
                     table_opts = table_options(   
                       row.FIAname = TRUE,          # table - return FIA row names
                       col.FIAname = TRUE           # table - return FIA column names
                     ))

```

To get the names of the list components associated with the output of our call of `modGBarea`, we run the following code:

```{r, results = T}
names(area4.1)
```

To easily access our estimate and percent sampling error of estimate we can just grab the `est` object from our ouputted list: 

```{r, results = T}
area4.1$est
```

We can also look at raw data and titles for estimate, as shown below. Note the change in titles.

```{r, results = TRUE}
raw4.1 <- area4.1$raw 
names(raw4.1)
head(raw4.1$unit_grpest)  # estimates by row and group domains

## Titles (list object) for estimate
titlelst4.1 <- area4.1$titlelst
names(titlelst4.1)
titlelst4.1
```

</details>


</details>

<a name="modGBtree"></a>

### modGBtree

`FIESTA`'s `modGBtree` function generates tree estimates based Scott et al. 2005 ('the green-book') for mapped forest inventory plots. The non-ratio estimator for estimating tree attributes by stratum and domain is used. Plots that are totally nonsampled are excluded from estimation dataset. Next, an adjustment factor is calculated by strata to adjust for nonsampled (nonresponse) conditions that have proportion less than 1. Attributes adjusted to a per-acre value are summed by plot, divided by the adjustment factor, and averaged by stratum. Strata means are combined using the strata weights and then expanded to using the total land area in the population. 

If there are more than one estimation unit (i.e., subpopulation) within the population, estimates are generated by estimation unit. If `sumunits = TRUE`, the estimates and percent standard errors returned are a sum combination of all estimation units. If `rawdata = TRUE`, the raw data returned will include estimates by estimation unit.  

Parameters defined in the following examples are organized by category: population data (pop); estimation information (est); and output details (out).

The following reference table can be used for defining `estvar` and `estvar.filter`:

```{r}
FIESTAutils::ref_estvar[, c("ESTTITLE", "ESTVAR", "ESTFILTER", "ESTUNITS")] 
```


#### POP1: 1.1 Net cubic-foot volume of live trees, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

Now, we can generate estimates by estimation unit (i.e., ESTN_UNIT) and sum to population (i.e., WY) with `modGBtree`: 

```{r}
tree1.1 <- modGBtree(GBpopdat = GBpopdat,               # pop - population calculations
                     landarea = "FOREST",               # est - forest land filter
                     sumunits = TRUE,                   # est - sum estimation units to population
                     estvar = "VOLCFNET",               # est - net cubic-foot volume
                     estvar.filter = "STATUSCD == 1",   # est - live trees only
                     returntitle = TRUE)                # out - return title information
```

We can now take a look at the output list, estimates and percent sampling errors, raw data, and titles:

```{r}
# Look at output list
names(tree1.1)

# Estimate and percent sampling error of estimate
tree1.1$est

# Raw data (list object) for estimate
raw1.1 <- tree1.1$raw      # extract raw data list object from output
names(raw1.1)
head(raw1.1$unit_totest)   # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.1$totest)        # estimates for population (i.e., WY)


# Titles (list object) for estimate
titlelst1.1 <- tree1.1$titlelst
names(titlelst1.1)
titlelst1.1
```

</details>

#### POP1: 1.2 Net cubic-foot volume of live trees by forest type, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

This example adds rows to the output for net cubic-foot volume of live trees (at least 5 inches diameter) by forest type, Wyoming, 2011-2013. We also choose to return titles with `returntitle = TRUE`. 

```{r}
tree1.2 <- modGBtree(GBpopdat = GBpopdat,               # pop - population calculations
                     landarea = "FOREST",               # est - forest land filter
                     sumunits = TRUE,                   # est - sum estimation units to population
                     estvar = "VOLCFNET",               # est - net cubic-foot volume
                     estvar.filter = "STATUSCD == 1",   # est - live trees only
                     rowvar = "FORTYPCD",               # est - row domain 
                     returntitle = TRUE)                # out - return title information
```

Again, we investigate the output of the returned list:

```{r}
# Look at output list
names(tree1.2)

# Estimate and percent sampling error of estimate
tree1.2$est

# Raw data (list object) for estimate
raw1.2 <- tree1.2$raw      # extract raw data list object from output
names(raw1.2)
head(raw1.2$unit_totest)   # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$totest)        # estimates for population (i.e., WY)
head(raw1.2$unit_rowest)   # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest)        # estimates by row for population (i.e., WY)

# Titles (list object) for estimate
titlelst1.2 <- tree1.2$titlelst
names(titlelst1.2)
titlelst1.2
```


We can also create a simple barplot from the output:

```{r}
datBarplot(raw1.2$unit_rowest, 
           xvar = titlelst1.2$title.rowvar, 
           yvar = "est")
```


</details>

#### POP1: 1.3 Net cubic-foot volume of live trees by forest type and stand-size class, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

This examples adds rows and columns to the output, including FIA names, for net cubic-foot volume of live trees (at least 5 inches diameter) by forest type and stand-size class, Wyoming, 2011-2013. We also use the `*_options` functions to return output with estimates (est) and percent standard error (pse) in same cell - est(pse) with `allin1 = TRUE` and save data to an outfolder with `savedata = TRUE` and `outfolder = outfolder`. 

```{r}
tree1.3 <- modGBtree(GBpopdat = GBpopdat,                # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "VOLCFNET",                # est - net cubic-foot volume
                     estvar.filter = "STATUSCD  == 1",   # est - live trees only
                     rowvar = "FORTYPCD",                # est - row domain
                     colvar = "STDSZCD",                 # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     savedata = TRUE,                    # out - save data to outfolder
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       col.FIAname = TRUE,               # est - column domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ),
                     savedata_opts = savedata_options(
                       outfolder = outfolder,            # out - outfolder for saving data
                       outfn.pre = "WY"                  # out - prefix for output files
                     ))
```

Again, we investigate the output of the returned list:

```{r}
# Look at output list from modGBarea()
names(tree1.3)

# Estimate and percent sampling error of estimate
tree1.3$est


# Raw data (list object) for estimate
raw1.3 <- tree1.3$raw      # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest)   # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest)        # estimates for population (i.e., WY)
head(raw1.3$unit_rowest)   # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest)        # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest)   # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest)        # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest)   # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest)        # estimates by row and column for population (i.e., WY)


# Titles (list object) for estimate
titlelst1.3 <- tree1.3$titlelst
names(titlelst1.3)
titlelst1.3


# List output files in outfolder
list.files(outfolder, pattern = "WY_tree")
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_tree")

```

</details>

#### POP1: 1.4 Number of live trees by species, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>
  

```{r}
tree1.4 <- modGBtree(GBpopdat = GBpopdat,                # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "SPCD",                    # est - row domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(    
                       row.FIAname = TRUE,               # est - row domain names
                       allin1 = FALSE                    # out - return output with est and pse
                     ))
```

We can also look at the output list and estimates again:

```{r}
# Look at output list
names(tree1.4)

# Estimate and percent sampling error of estimate
tree1.4$est
```

</details>

#### POP1: 1.5 Number of live trees (including seedlings) by species, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

We can also add seedlings. 

Note: seedling data are only available for estimating number of trees (estvar = TPA_UNADJ).

Note: must include seedling data in population data calculations.

```{r}
tree1.5 <- modGBtree(GBpopdat = GBpopdat,                # pop - population calculations
                     estseed = "add",                    # est - add seedling data
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "SPCD",                    # est - row domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       allin1 = FALSE,                   # out - return output with est and pse
                     ))
```

And again we can look at our outputs and compare estimates:

```{r}
# Look at output list
names(tree1.5)

# Estimate and percent sampling error of estimate
tree1.5$est

# Compare estimates with and without seedlings
head(tree1.4$est)
head(tree1.5$est)
```

</details>

#### POP1: 1.6 Number of seedlings by species, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

Of course, we can also look at *only* seedlings.

```{r}
tree1.6 <- modGBtree(GBpopdat = GBpopdat,           # pop - population calculations
                     estseed = "only",              # est - add seedling data
                     landarea = "FOREST",           # est - forest land filter
                     sumunits = TRUE,               # est - sum estimation units to population
                     estvar = "TPA_UNADJ",          # est - number of trees per acre 
                     rowvar = "SPCD",               # est - row domain
                     returntitle = TRUE,            # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,          # est - row domain names
                       allin1 = FALSE               # out - return output with est and pse
                     ))
```

And again we can look at our outputs and compare estimates:

```{r}
# Look at output list
names(tree1.6)

# Estimate and percent sampling error of estimate
tree1.6$est

# Compare estimates with, without, and only seedlings
head(tree1.4$est)
head(tree1.5$est)
head(tree1.6$est)
```

</details>

#### POP2: 2.1 Number of live trees by forest type and species on forest land, Bighorn National Forest
<details>
  <summary>View Example</summary>


```{r}
tree2.1 <- modGBtree(GBpopdat = GBpopdat.bh,             # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "FORTYPCD",                # est - row domain
                     colvar = "SPCD",                    # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       col.FIAname = TRUE,               # est - column domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ))
```

And we can see our output:

```{r}
# Look at output list
names(tree2.1)

# Estimate and percent sampling error of estimate
head(tree2.1$est)
```

</details>

#### POP2: 2.2 Net cubic-foot volume of standing dead trees by species and cause of death on forest land, Bighorn National Forest
<details>
  <summary>View Example</summary>

We can also examine dead trees with the filter `estvar.filter = "STATUSCD == 2 & STANDING_DEAD_CD == 1"`. 

```{r}
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"
tree2.2 <- modGBtree(GBpopdat = GBpopdat.bh,                  # pop - population calculations
                     landarea = "FOREST",                     # est - forest land filter
                     sumunits = TRUE,                         # est - sum estimation units to population
                     estvar = "VOLCFNET",                     # est - number of trees per acre 
                     estvar.filter = deadtree.filter,         # est - standing dead trees only
                     rowvar = "SPCD",                         # est - row domain
                     colvar = "AGENTCD",                      # est - column domain
                     returntitle = TRUE,                      # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,                    # est - row domain names
                       col.FIAname = TRUE,                    # est - column domain names
                       allin1 = TRUE                          # out - return output with est(pse)
                     ))
```

And we can see our output of dead trees:

```{r}
# Look at output list
names(tree2.2)

# Estimate and percent sampling error of estimate
head(tree2.2$est)
```

</details>

#### Adding diameter classes to population tree data (for `modGBtree` and `modGBratio` examples 7-9)
<details>
  <summary>View Example</summary>

```{r}
# Look at tree data in GBpopdat.bh
head(GBpopdat.bh$treex)

# Use reference data frame stored as an R object in FIESTA
head(FIESTAutils::ref_diacl2in)

# Appends a new column to GBpopdat$treex classifying the DIA variable based on MIN and MAX columns in ref_diacl2in
dat <- datLUTclass(x = GBpopdat.bh$treex, 
                   xvar = "DIA", 
                   LUT = FIESTAutils::ref_diacl2in, 
                   LUTclassnm = "DIACL2IN")

GBpopdat.bh$treex <- dat$xLUT

# Look at tree data, with new column (DIACL2IN)  
head(GBpopdat.bh$treex)

# Look at table of new diameter classes (DIACL2IN)
table(GBpopdat.bh$treex$DIACL2IN)



# Another way to append diameter classes
# First, create a new variable using cut function to define 4 diameter classes
dat <- datLUTclass(x = GBpopdat.bh$treex, 
                   xvar = "DIA", 
                   cutbreaks = c(0, 5, 10, 20, 100))

GBpopdat.bh$treex <- dat$xLUT

# Look at tree data, with new column (DIACL2IN)  
head(GBpopdat.bh$treex)

# Look at table of new diameter classes (DIACL)
table(GBpopdat.bh$treex$DIACL)
```

</details>

#### POP2: 2.3 Number of Live Trees by Species Groups and Diameter Class on forest land, Bighorn National Forest
<details>
  <summary>View Example</summary>

We can also look at the number of live trees by species group and diameter class (DIACL2IN):

```{r}
tree2.3 <- modGBtree(GBpopdat = GBpopdat.bh,             # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "SPGRPCD",                 # est - row domain
                     colvar = "DIACL2IN",                # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ))
```

Outputs:

```{r}
# Look at output list
names(tree2.3)

# Estimate and percent sampling error of estimate
head(tree2.3$est)
```

</details>

#### POP2: 2.4 Number of Live Trees by Species Groups and a Different Diameter Class on forest land, Bighorn National Forest
<details>
  <summary>View Example</summary>

Next, we can look at number of live trees by species group and diameter class (DIACL):

```{r}
tree2.4 <- modGBtree(GBpopdat = GBpopdat.bh,             # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "SPGRPCD",                 # est - row domain
                     colvar = "DIACL",                   # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ))
```

Outputs:

```{r}
# Look at output list
names(tree2.4)

# Estimate and percent sampling error of estimate
head(tree2.4$est)
```

</details>

#### POP2: 2.5 Number of Live Trees (+ seedlings) by Species Groups and a Different Diameter Class on forest land, Bighorn National Forest
<details>
  <summary>View Example</summary>

Finally, we add seedlings to Example 8:

```{r}
tree2.5 <- modGBtree(GBpopdat = GBpopdat.bh,             # pop - population calculations
                     estseed = "add",                    # est - add seedling data
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "TPA_UNADJ",               # est - number of trees per acre 
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "SPGRPCD",                 # est - row domain
                     colvar = "DIACL",                   # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ))
```

And look at the outputs:

```{r}
# Look at output list
names(tree2.5)

# Estimate and percent sampling error of estimate
head(tree2.5$est)
```


</details>

#### POP3: 3.1 Volume of Dead Trees by Forest Type Group and Primary Disturbance on forest land, Bighorn National Forest Districts
<details>
  <summary>View Example</summary>

```{r}
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"

tree3.1 <- modGBtree(GBpopdat = GBpopdat.bhdist,             # pop - population calculations
                     landarea = "FOREST",                    # est - forest land filter
                     sumunits = FALSE,                       # est - sum estimation units to population
                     estvar = "VOLCFNET",                    # est - number of trees per acre 
                     estvar.filter = deadtree.filter,        # est - only dead trees 
                     rowvar = "FORTYPGRPCD",                 # est - row domain
                     returntitle = TRUE,                     # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,                   # est - row domain names
                       allin1 = TRUE                         # out - return output with est(pse)
                     ))
```

Outputs:

```{r}
# Look at output list
names(tree3.1)

# Estimate and percent sampling error of estimate
tree3.1$est

# Estimate and percent sampling error by district
tree3.1$raw$unit_rowest
```

</details>


#### POP4: 4.1 Net cubic-foot volume of live trees by forest type group and stand-size class, Rhode Island, 2019
<details>
  <summary>View Example</summary>
  

```{r}
tree4.1 <- modGBtree(GBpopdat = GBpopdat.RI,             # pop - population calculations
                     landarea = "FOREST",                # est - forest land filter
                     sumunits = TRUE,                    # est - sum estimation units to population
                     estvar = "VOLCFNET",                # est - net cubic-foot volume estimate
                     estvar.filter = "STATUSCD == 1",    # est - live trees only
                     rowvar = "FORTYPCD",                # est - row domain
                     colvar = "STDSZCD",                 # est - column domain
                     returntitle = TRUE,                 # out - return title information
                     table_opts = table_options(
                       row.FIAname = TRUE,               # est - row domain names
                       col.FIAname = TRUE,               # est - column domain names
                       allin1 = TRUE                     # out - return output with est(pse)
                     ))
```

Outputs:

```{r}
# Look at output list
names(tree4.1)

# Estimate and percent sampling error of estimate
head(tree4.1$est)
```

</details>

<a name="modGBratio"></a>

### modGBratio

`FIESTA`'s `modGBratio` function generates per-acre and per-tree estimates based on Scott et al. 2005 ('the green-book') for mapped forest inventory plots. The ratio estimator for estimating per-acre or per-tree by stratum and domain is used, referred to as Ratio of Means (ROM). 

If there are more than one estimation unit (i.e., subpopulation) within the population, estimates are generated by estimation unit. If sumunits = TRUE, the estimates and percent standard errors returned are a sum combination of all estimation units. If rawdata = TRUE, the raw data returned will include estimates by estimation unit.  

Parameters defined in the following examples are organized by category: population data (pop); estimation information (est); and output details (out).

#### POP1: 1.1 Net cubic-foot volume per acre of live trees on timberland, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>

We generate estimates by estimation unit (i.e., ESTN_UNIT) and sum to population (i.e., WY):

```{r}
ratio1.1 <- modGBratio(GBpopdat = GBpopdat,                # pop - population calculations
                       landarea = "TIMBERLAND",            # est - forest land filter
                       sumunits = TRUE,                    # est - sum estimation units to population
                       estvarn = "VOLCFNET",               # est - net cubic-foot volume, numerator
                       estvarn.filter = "STATUSCD == 1",   # est - live trees only, numerator
                       returntitle = TRUE)                 # out - return title information
```

And we can look at our output:

```{r}
# Look at output list
names(ratio1.1)

# Estimate and percent sampling error of estimate
head(ratio1.1$est)

# Raw data (list object) for estimate
raw1.1 <- ratio1.1$raw      # extract raw data list object from output
names(raw1.1)
head(raw1.1$unit_totest)    # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.1$totest)         # estimates for population (i.e., WY)

# Titles (list object) for estimate
titlelst <- ratio1.1$titlelst
names(titlelst)
titlelst
```

</details>

#### POP1: 1.2 Net cubic-foot volume per acre of live trees by forest type on timberland, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>


```{r}
ratio1.2 <- modGBratio(GBpopdat = GBpopdat,                # pop - population calculations
                       landarea = "TIMBERLAND",            # est - forest land filter
                       sumunits = TRUE,                    # est - sum estimation units to population
                       estvarn = "VOLCFNET",               # est - net cubic-foot volume
                       estvarn.filter = "STATUSCD == 1",   # est - live trees only
                       rowvar = "FORTYPCD",                # est - row domain 
                       returntitle = TRUE)                 # out - return title information
```

And of course view our outputs:

```{r}
# Look at output list
names(ratio1.2)

# Estimate and percent sampling error of estimate
head(ratio1.2$est)

# Raw data (list object) for estimate
raw1.2 <- ratio1.2$raw      # extract raw data list object from output
names(raw1.2)
head(raw1.2$unit_totest)    # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$totest)         # estimates for population (i.e., WY)
head(raw1.2$unit_rowest)    # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.2$rowest)         # estimates by row for population (i.e., WY)

# Titles (list object) for estimate
titlelst <- ratio1.2$titlelst
names(titlelst)
titlelst
```

</details>

#### POP1: 1.3 Net cubic-foot volume per acre of live trees by forest type and stand-size class on timberland, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>


```{r}
ratio1.3 <- modGBratio(GBpopdat = GBpopdat,                # pop - population calculations
                       landarea = "TIMBERLAND",            # est - forest land filter
                       sumunits = TRUE,                    # est - sum estimation units to population
                       estvarn = "VOLCFNET",               # est - net cubic-foot volume, numerator
                       estvarn.filter = "STATUSCD == 1",   # est - live trees only, numerator
                       rowvar = "FORTYPCD",                # est - row domain
                       colvar = "STDSZCD",                 # est - column domain
                       returntitle = TRUE,                 # out - return title information
                       savedata = TRUE,                    # out - save data to outfolder
                       table_opts = table_options(
                         row.FIAname = TRUE,               # est - row domain names
                         col.FIAname = TRUE,               # est - column domain names
                         allin1 = TRUE                     # out - return output with est(pse)
                       ),
                       savedata_opts = savedata_options(
                         outfolder = outfolder,            # out - outfolder for saving data
                         outfn.pre = "WY"                  # out - prefix for output files
                       ))
```

And look at our output again:

```{r}
# Look at output list from modGBarea()
names(ratio1.3)

# Estimate and percent sampling error of estimate
head(ratio1.3$est)

# Raw data (list object) for estimate
raw1.3 <- ratio1.3$raw      # extract raw data list object from output
names(raw1.3)
head(raw1.3$unit_totest)    # estimates by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$totest)         # estimates for population (i.e., WY)
head(raw1.3$unit_rowest)    # estimates by row, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$rowest)         # estimates by row for population (i.e., WY)
head(raw1.3$unit_colest)    # estimates by column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$colest)         # estimates by column for population (i.e., WY)
head(raw1.3$unit_grpest)    # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)
head(raw1.3$grpest)         # estimates by row and column for population (i.e., WY)

# Titles (list object) for estimate
titlelst <- ratio1.3$titlelst
names(titlelst)
titlelst
```

</details>

#### POP2: 2.1 Number of live trees per acre by species, Bighorn National Forest
<details>
  <summary>View Example</summary>


```{r}
ratio2.1 <- modGBratio(GBpopdat = GBpopdat.bh,              # pop - population calculations
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                       rowvar = "SPCD",                     # est - row domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         allin1 = FALSE,                    # out - return output with est and pse
                       ),
                       title_opts = title_options(
                         title.ref = "Bighorn National Forest"
                       ))
```

And we can look at our output:

```{r}
# Look at output list
names(ratio2.1)

# Estimate and percent sampling error of estimate
head(ratio2.1$est)
```

</details>

#### POP2: 2.2 Number of live trees (including seedlings) per acre by species, Bighorn National Forest
<details>
  <summary>View Example</summary>

```{r}
ratio2.2 <- modGBratio(GBpopdat = GBpopdat.bh,              # pop - population calculations
                       estseed = "add",                     # est - add seedling data
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                       rowvar = "SPCD",                     # est - row domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         allin1 = FALSE                     # out - return output with est and pse
                       ),
                       title_opts = title_options(
                         title.ref = "Bighorn National Forest"
                       ))
```

Output and comparison:

```{r}
# Look at output list
names(ratio2.2)

# Estimate and percent sampling error of estimate
head(ratio2.2$est)

# Compare estimates with and without seedlings
head(ratio2.1$est)
head(ratio2.2$est)
```

</details>

#### POP2: 2.3 Number of live seedlings per acre by species, Bighorn National Forest
<details>
  <summary>View Example</summary>


```{r}
ratio2.3 <- modGBratio(GBpopdat = GBpopdat.bh,         # pop - population calculations
                       estseed = "only",               # est - add seedling data
                       landarea = "FOREST",            # est - forest land filter
                       sumunits = TRUE,                # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",          # est - number of trees per acre, numerator 
                       rowvar = "SPCD",                # est - row domain
                       returntitle = TRUE,             # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,           # est - row domain names
                         allin1 = FALSE                # out - return output with est and pse
                       ))
```

Output and comparisons:

```{r}
# Look at output list
names(ratio2.3)

# Estimate and percent sampling error of estimate
head(ratio2.3$est)

# Compare estimates with, without, and only seedlings
head(ratio2.1$est)
head(ratio2.2$est)
head(ratio2.3$est)
```

</details>

#### POP2: 2.4 Number of live trees by forest type and species, Bighorn National Forest
<details>
  <summary>View Example</summary>


```{r}
ratio2.4 <- modGBratio(GBpopdat = GBpopdat.bh,              # pop - population calculations
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                       rowvar = "FORTYPCD",                 # est - row domain
                       colvar = "SPCD",                     # est - column domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         col.FIAname = TRUE,                # est - column domain names
                         allin1 = TRUE                      # out - return output with est(pse)
                       ))

```

And view our output:

```{r}
# Look at output list
names(ratio2.4)

# Estimate and percent sampling error of estimate
head(ratio2.4$est)
```

</details>

#### POP2: 2.5 Number of standing dead trees by species and cause of death, Bighorn National Forest
<details>
  <summary>View Example</summary>

Next, we look at dead trees: 

```{r}
deadtree.filter <- "STATUSCD == 2 & STANDING_DEAD_CD == 1"

ratio2.5 <- modGBratio(GBpopdat = GBpopdat.bh,              # pop - population calculations
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "VOLCFNET",                # est - number of trees per acre, numerator 
                       estvarn.filter = deadtree.filter,    # est - standing dead trees only, numerator
                       rowvar = "SPCD",                     # est - row domain
                       colvar = "AGENTCD",                  # est - column domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         col.FIAname = TRUE,                # est - column domain names
                         allin1 = TRUE                      # out - return output with est(pse)
                       ))
```

And we can see our output:

```{r}
# Look at output list
names(ratio2.5)

# Estimate and percent sampling error of estimate
head(ratio2.5$est)
```


</details>

#### POP3.1: 3.1 Number of live trees by species group and diameter class (DIACL2IN), Bighorn National Forest 
<details>
  <summary>View Example</summary>


```{r}
ratio2.6 <- modGBratio(GBpopdat = GBpopdat.bh,              # pop - population calculations
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                       rowvar = "SPGRPCD",                  # est - row domain
                       colvar = "DIACL2IN",                 # est - column domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = list(
                         row.FIAname = TRUE,                # est - row domain names
                         allin1 = TRUE                      # out - return output with est(pse)
                       ),
                       title_opts = list(
                         title.ref = "Bighorn National Forest"
                       ))
```

And examine our output:

```{r}
# Look at output list
names(ratio2.6)

# Estimate and percent sampling error of estimate
head(ratio2.6$est)
```

</details>

#### POP3: 3.2 Number of Live Trees per acre by Species Group and Diameter Class, Bighorn National Forest Districts
<details>
  <summary>View Example</summary>


```{r}
ratio3.2 <- modGBratio(GBpopdat = GBpopdat.bhdist,          # pop - population calculations
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                       rowvar = "SPGRPCD",                  # est - row domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         allin1 = TRUE,                     # out - return output with est(pse)
                       ),
                       title_opts = title_options(
                         title.ref="Bighorn National Forest Districts"
                       ))
```

And of course examine our output:

```{r}
# Look at output list
names(ratio3.2)

# Estimate and percent sampling error of estimate
head(ratio3.2$est)
```


</details>

#### POP1: 1.4 Net cubic-foot volume of live trees (at least 5 inches diameter) divided by net cubic-foot volume of all trees by forest type, Wyoming, 2011-2013
<details>
  <summary>View Example</summary>


```{r}
ratio1.4 <- modGBratio(GBpopdat = GBpopdat,                 # pop - population calculations for WY, post-stratification
                       landarea = "FOREST",                 # est - forest land filter
                       sumunits = TRUE,                     # est - sum estimation units to population
                       estvarn = "VOLCFNET",                # est - net cubic-foot volume, numerator
                       estvarn.filter = "STATUSCD == 1",    # est - live trees only
                       estvard = "VOLCFNET",                # est - net cubic-foot volume, numerator
                       rowvar = "FORTYPCD",                 # est - row domain
                       returntitle = TRUE,                  # out - return title information
                       table_opts = table_options(
                         row.FIAname = TRUE,                # est - row domain names
                         allin1 = TRUE,                     # out - return output with est(pse)
                       ))
```

And examine the output:

```{r}
# Look at output list
names(ratio1.4)

# Estimate and percent sampling error of estimate
head(ratio1.4$est)
```

</details>

### No strata

If you want to exclude post-stratification (i.e., simple random sample, Horvitz-Thompson), you must generate a new population dataset using the `modGBpop` function by setting `strata = FALSE`.

<details>
  <summary>View Code</summary>

```{r}
## Get population data for Wyoming estimates, with post-stratification (for comparison)
GBpopdat.strat <- modGBpop(popTabs = popTables(
                                       cond = WYcond,               # FIA plot/condition data
                                       tree = WYtree,               # FIA tree data
                                       seed = WYseed                # FIA seedling data
                                     ),
                           pltassgn = WYpltassgn,                   # plot assignments
                           pltassgnid = "CN",                       # uniqueid of plots
                           unitarea = WYunitarea,                   # area by estimation units
                           unitvar = "ESTN_UNIT",                   # name of estimation unit
                           strata = TRUE,                           # if using post-stratification
                           stratalut = WYstratalut,                 # strata classes and pixels counts
                           strata_opts = strata_options(
                             getwt=TRUE,                            # calculate strata weights
                             getwtvar="P1POINTCNT"                  # use P1POINTCNT in stratalut to calculate weights
                           ))

## Get population data for Wyoming estimates, with no post-stratification
GBpopdat.nostrat <- modGBpop(popTabs = popTables(
                                          cond = WYcond,            # FIA plot/condition data
                                          tree = WYtree,            # FIA tree data
                                          seed = WYseed             # FIA seedling data
                                       ),
                             pltassgn = WYpltassgn,                 # plot assignments
                             pltassgnid = "CN",                     # uniqueid of plots
                             unitarea = WYunitarea,                 # area by estimation units
                             unitvar = "ESTN_UNIT",                 # name of estimation unit
                             strata = FALSE)                        # if using post-stratification
```


```{r}
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013, with post-stratification
area.strat <- modGBarea(GBpopdat = GBpopdat.strat,        # pop - population calculations for WY, post-stratification
                        landarea = "FOREST",              # est - forest land filter
                        sumunits = TRUE,                  # est - sum estimation units to population
                        rowvar = "FORTYPCD",              # est - row domain
                        table_opts = table_options(
                          row.FIAname = TRUE,             # est - row domain names
                          allin1 = FALSE                  # out - return output with est(pse)
                        ))   

## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013, no post-stratification
area.nostrat <- modGBarea(GBpopdat = GBpopdat.nostrat,    # pop - population calculations for WY, no post-stratification
                          landarea = "FOREST",            # est - forest land filter
                          sumunits = TRUE,                # est - sum estimation units to population
                          rowvar = "FORTYPCD",            # est - row domain
                          table_opts = table_options(
                            row.FIAname = TRUE,           # est - row domain names
                            allin1 = FALSE                # out - return output with est(pse)
                          ))   
```

```{r}
# Compare estimates and percent standard errors with and without post-stratification
head(area.strat$est)
head(area.nostrat$est)
```

```{r}
## Number of live trees by species, Wyoming, 2011-2013, with post-stratification
tree.strat <- modGBtree(GBpopdat = GBpopdat.strat,          # pop - population calculations for WY, post-stratification
                        landarea = "FOREST",                # est - forest land filter
                        sumunits = TRUE,                    # est - sum estimation units to population
                        estvar = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                        estvar.filter = "STATUSCD == 1",    # est - live trees only, numerator
                        rowvar = "FORTYPCD",                # est - row domain
                        table_opts = table_options(
                          row.FIAname = TRUE,               # est - row domain names
                          allin1 = FALSE                    # out - return output with est(pse)
                        ))   

## Number of live trees by species, Wyoming, 2011-2013, no post-stratification
tree.nostrat <- modGBtree(GBpopdat = GBpopdat.nostrat,      # pop - population calculations for WY, no post-stratification
                          landarea = "FOREST",              # est - forest land filter
                          sumunits = TRUE,                  # est - sum estimation units to population
                          estvar = "TPA_UNADJ",             # est - number of trees per acre, numerator 
                          estvar.filter = "STATUSCD == 1",  # est - live trees only, numerator
                          rowvar = "FORTYPCD",              # est - row domain
                          table_opts = table_options(
                            row.FIAname = TRUE,             # est - row domain names
                            allin1 = FALSE                  # out - return output with est(pse)
                          ))   
```


```{r}
# Compare estimates and percent standard errors with and without post-stratification
head(tree.strat$est)
head(tree.nostrat$est)
```

```{r}
## Number of live trees per acre by species, Wyoming, 2011-2013, with post-stratification
ratio.strat <- modGBratio(GBpopdat = GBpopdat.strat,           # pop - population calculations for WY, post-stratification
                          landarea = "FOREST",                 # est - forest land filter
                          sumunits = TRUE,                     # est - sum estimation units to population
                          estvarn = "TPA_UNADJ",               # est - number of trees per acre, numerator 
                          estvarn.filter = "STATUSCD == 1",    # est - live trees only, numerator
                          rowvar = "FORTYPCD",                 # est - row domain
                          table_opts = table_options(
                            row.FIAname = TRUE,                # est - row domain names
                            allin1 = FALSE                     # out - return output with est(pse)
                          ))   

## Number of live trees per acre by species, Wyoming, 2011-2013, no post-stratification
ratio.nostrat <- modGBratio(GBpopdat = GBpopdat.nostrat,       # pop - population calculations for WY, no post-stratification
                            landarea = "FOREST",               # est - forest land filter
                            sumunits = TRUE,                   # est - sum estimation units to population
                            estvarn = "TPA_UNADJ",             # est - number of trees per acre, numerator 
                            estvarn.filter = "STATUSCD == 1",  # est - live trees only, numerator
                            rowvar = "FORTYPCD",               # est - row domain
                            table_opts = table_options(
                              row.FIAname = TRUE,              # est - row domain names
                              allin1 = FALSE                   # out - return output with est(pse)
                            ))  
```


```{r}
# Compare estimates and percent standard errors with and without post-stratification
head(ratio.strat$est)
head(ratio.nostrat$est)
```

</details>

### By estimation unit (`sumunits = FALSE`)

If you just wanted estimates by estimation unit and did not want to sum them, set `sumunits = FALSE`. If `sumunits = FALSE`, the estimates and percent standard errors returned are by estimation unit, with an attribute, named 'unit' appended to data frame, with the unit value. The raw data returned will look the same as if `sumunits = TRUE`. 

<details>
  <summary>View Code</summary>

```{r}
area.unit <- modGBarea(GBpopdat = GBpopdat,         # pop - population calculations for WY, post-stratification
                       landarea = "FOREST",         # est - forest land filter
                       sumunits = FALSE,            # est - sum estimation units to population
                       rowvar = "FORTYPCD",         # est - row domain
                       colvar = "STDSZCD",          # est - column domain
                       returntitle = TRUE,          # out - return title information
                       table_opts = table_options(
                         allin1 = TRUE              # out - return output with est(pse)
                       ))

# Estimate and percent sampling error of estimate (first 6 rows)
head(area.unit$est)

# Unique estimation units
unique(area.unit$est$unit)
```

</details>

<!-- ### Population data included in estimation module -->

<!-- If you were only interested in a couple estimates for a population, you could calculate population data within a module by including all information directly into the estimation module. The `modGBpop` function is called if the GBpopdat list object is not included as a parameter. -->

<!-- <details> -->
<!--   <summary>View Code</summary> -->

<!-- ```{r} -->
<!-- ## Get population data for Wyoming, with post-stratification -->
<!-- ## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013 -->
<!-- area.strat <- modGBarea(  -->
<!--     cond = WYcond,               # pop - FIA plot/condition data -->
<!--     tree = WYtree,               # pop - FIA tree data -->
<!--     pltassgn = WYpltassgn,       # pop - plot assignments -->
<!--     pltassgnid = "CN",           # pop - uniqueid of plots -->
<!--     unitarea = WYunitarea,       # pop - area by estimation units -->
<!--     unitvar = "ESTN_UNIT",       # pop - name of estimation unit -->
<!--     strata = TRUE,               # pop - if using post-stratification -->
<!--     stratalut = WYstratalut,     # strata classes and pixels counts -->
<!--     landarea = "FOREST",         # est - forest land filter -->
<!--     sumunits = TRUE,             # est - sum estimation units to population -->
<!--     rowvar = "FORTYPCD",         # est - row domain -->
<!--     colvar = "STDSZCD",          # est - column domain -->
<!--     table_opts = table_options(allin1 = TRUE)               # out - return output with est(pse) -->
<!--     )    -->


<!-- ## Get population data for Wyoming, with no post-stratification -->
<!-- ## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013 -->
<!-- area.nostrat <- modGBarea(  -->
<!--     cond = WYcond,               # pop - FIA plot/condition data -->
<!--     tree = WYtree,               # pop - FIA tree data -->
<!--     pltassgn = WYpltassgn,       # pop - plot assignments -->
<!--     pltassgnid = "CN",           # pop - uniqueid of plots -->
<!--     unitarea = WYunitarea,       # pop - area by estimation units -->
<!--     unitvar = "ESTN_UNIT",       # pop - name of estimation unit -->
<!--     strata = FALSE,              # pop - if using post-stratification -->
<!--     landarea = "FOREST",         # est - forest land filter -->
<!--     sumunits = TRUE,             # est - sum estimation units to population -->
<!--     rowvar = "FORTYPCD",         # est - row domain -->
<!--     colvar = "STDSZCD",          # est - column domain -->
<!--     rawdata = TRUE,              # out - return raw data  -->
<!--     allin1 = TRUE                # out - return output with est(pse) -->
<!--     )    -->


<!-- ## Compare estimates and percent standard errors with and without post-stratification -->
<!-- head(area.strat$est) -->
<!-- head(area.nostrat$est) -->
<!-- ``` -->

</details>

## References

Bechtold, William A.; Patterson, Paul L., Editors. 2005. The enhanced Forest Inventory and Analysis program national sampling design and estimation procedures. Gen. Tech. Rep. SRS-80. Asheville, NC: U.S. Department of Agriculture, Forest Service, Southern Research Station. 85 p.

Patterson, Paul L. 2012. Photo-based estimators for the Nevada photo-based inventory. Res. Pap. RMRS-RP-92. Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. 14 p.

Westfall, James A.; Patterson, Paul L.; Coulston, John W. 2011. Post-stratified estimation: with-in strata and total sample size recommendations. Canadian Journal of Forest Research. 41: 1130-1139.

