FIESTA - Model-Assisted Estimators
FIESTA_tutorial_MA.RmdModel-Assisted (MA) module overview
FIESTA’s Model-Assisted (MA) module calculates
population estimates and their sampling errors by taking advantage of
available model-assisted survey estimators from the mase R
package (McConville, et al. 2018). These estimators can use a variety of
auxiliary data to build models and predict over a response variable of
interest, while using a bias-correction term so that the bias of the
model does not depend on model mis-specification.
Functions in FIESTA used for fitting model-assisted
estimators include the modMAarea function for area
estimates and modMAtree for tree estimates. The
modMApop function is used to get population data needed for
model-assisted estimation. Below is a description and table of contents
for the sections related to these functions:
| FUNCTION | DESCRIPTION |
|---|---|
| modMApop | Creates population data for model-assisted estimation. |
| modMAarea | Produces area level estimates through model-assisted estimation. |
| modMAtree | Produces tree level estimates through model-assisted estimation. |
Objective of tutorial
The main objective of this tutorial is to demonstrate how to use
FIESTA for generating estimates using estimators from
mase. The model-assisted 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 following examples are for generating estimates and estimated
variances using standard FIA Evaluation data from FIA’s National
database, with custom Estimation unit and Stratification information.
The examples use data from 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.
Example data - Wyoming (WY), Inventory Years 2011-2012
View MA Example Data
| Data Frame | Description |
|---|---|
| WYplt | WY plot-level data |
| WYcond | WY condition-level data |
| WYtree | WY tree-level data |
| External data | Description |
|---|---|
| WYbighorn_adminbnd.shp | Polygon shapefile of WY Bighorn National Forest Administrative boundary* |
| WYbighorn_districtbnd.shp | Polygon shapefile of WY Bighorn National Forest District boundaries** |
| WYbighorn_forest_nonforest_250m.tif | GeoTIFF raster of predicted forest/nonforest (1/0) for stratification*** |
| WYbighorn_dem_250m.img | Erdas Imagine raster of elevation change, in meters**** |
*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.
**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 encompasses a Ranger District.
***Based on MODIS-based classified map resampled from 250m to 500m resolution and reclassified from 3 to 2 classes: 1:forest; 2:nonforest. Projected in Albers Conical Equal Area, Datum NAD27 (Ruefenacht et al. 2008). Clipped to extent of WYbighorn_adminbnd.shp.
****USGS National Elevation Dataset (NED), resampled from 30m resolution to 250m. Projected in Albers Conical Equal Area, Datum NAD27 (U.S. Geological Survey 2017). Clipped to boundary of WYbighorn_adminbnd.shp.
Set up
First, you’ll need to load the FIESTA library:
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 our purposes in this vignette, we have saved our outfolder
file path as the outfolder object in a temporary directory.
We also set a few default options preferred for this vignette.
outfolder <- tempdir()Get data for examples
View Getting Data
Now that we’ve loaded FIESTA and setup our outfolder, we
can retrieve the data needed to run the examples. First, we point to
some external data and predictor layers stored in FIESTA
and derive new predictor layers using the terra
package.
# File names for external spatial data
WYbhfn <- system.file("extdata", "sp_data/WYbighorn_adminbnd.shp",
package = "FIESTA")
WYbhdistfn <- system.file("extdata", "sp_data/WYbighorn_districtbnd.shp",
package = "FIESTA")
## predictor variables
fornffn <- system.file("extdata", "sp_data/WYbighorn_forest_nonforest_250m.tif",
package = "FIESTA")
demfn <- system.file("extdata", "sp_data/WYbighorn_dem_250m.img",
package = "FIESTA")
# Derive new predictor layers from dem
library(terra)
dem <- rast(demfn)
slpfn <- paste0(outfolder, "/WYbh_slp.img")
slp <- terra::terrain(dem,
v = "slope",
unit = "degrees",
filename = slpfn,
overwrite = TRUE)
aspfn <- paste0(outfolder, "/WYbh_asp.img")
asp <- terra::terrain(dem,
v = "aspect",
unit = "degrees",
filename = aspfn,
overwrite = TRUE)Next, we can get our FIA plot data and set up our auxiliary data. We
can get our FIA plot data with the spMakeSpatialPoints
function from FIESTA. For more information on how to use
this function, please see the sp vignette included with
FIESTA (link).
WYspplt <- spMakeSpatialPoints(
xyplt = WYplt,
xy.uniqueid = "CN",
xvar = "LON_PUBLIC",
yvar = "LAT_PUBLIC",
xy.crs = 4269
)
rastlst.cont <- c(demfn, slpfn, aspfn)
rastlst.cont.name <- c("dem", "slp", "asp")
rastlst.cat <- fornffn
rastlst.cat.name <- "fornf"Next, we must generate dataset for model-assisted estimation. We can
do this with the spGetAuxiliary function from
FIESTA. Again, see the sp vignette for further
information on this function.
modeldat <- spGetAuxiliary(
xyplt = WYspplt,
uniqueid = "CN",
unit_layer = WYbhfn,
unitvar = NULL,
rastlst.cont = rastlst.cont,
rastlst.cont.name = rastlst.cont.name,
rastlst.cat = rastlst.cat,
rastlst.cat.name = rastlst.cat.name,
rastlst.cont.stat = "mean",
asptransform = TRUE,
rast.asp = aspfn,
keepNA = FALSE,
showext = FALSE,
savedata = FALSE)
str(modeldat, max.level = 1)output
## List of 12
## $ unitvar : chr "ONEUNIT"
## $ pltassgn :'data.frame': 56 obs. of 25 variables:
## $ pltassgnid : chr "CN"
## $ unitarea :'data.frame': 1 obs. of 2 variables:
## $ areavar : chr "ACRES_GIS"
## $ unitzonal :'data.frame': 1 obs. of 8 variables:
## $ inputdf :Classes 'data.table' and 'data.frame': 5 obs. of 7 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## $ prednames : chr [1:5] "dem" "slp" "asp_cos" "asp_sin" ...
## $ zonalnames : chr [1:7] "dem" "slp" "asp_cos" "asp_sin" ...
## $ predfac : chr "fornf"
## $ npixelvar : chr "npixels"
## $ predfac.levels:List of 1Examples
modMApop
Example 1: Creating our population dataset with
modMApop
View Example
We can create our population data for model-assisted estimation. To
do so, we use the modMApop function in FIESTA.
We must assign our population tables with the popTabs
argument (and unique identifiers for these tables with the
popTabIDs argument if they are not the default), the plot
assignment with the pltassgn argument, and in auxiliary
dataset we just created with the auxdat argument. The
spGetAuxiliary function has done much of the hard work for
us so far, so we can just run a simple call to
modMApop:
MApopdat <- modMApop(popTabs = list(tree = WYtree, cond = WYcond),
pltassgn = WYpltassgn,
auxdat = modeldat)Note that the modMApop function returns a list with lots
of information and data for us to use. For a quick look at what this
list includes we can use the str function:
str(MApopdat, max.level = 1)output
## List of 26
## $ module : chr "MA"
## $ popType : chr "VOL"
## $ condx :Classes 'data.table' and 'data.frame': 66 obs. of 17 variables:
## ..- attr(*, "sorted")= chr [1:2] "PLT_CN" "CONDID"
## ..- attr(*, ".internal.selfref")=<externalptr>
## $ pltcondx :Classes 'data.table' and 'data.frame': 66 obs. of 29 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr [1:2] "PLT_CN" "CONDID"
## $ cuniqueid : chr "PLT_CN"
## $ condid : chr "CONDID"
## $ ACI.filter : chr "COND_STATUS_CD == 1"
## $ unitarea :Classes 'data.table' and 'data.frame': 1 obs. of 2 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "ONEUNIT"
## $ areavar : chr "ACRES_GIS"
## $ areaunits : chr "acres"
## $ unitvar : chr "ONEUNIT"
## $ unitvars : chr "ONEUNIT"
## $ unitlut :Classes 'data.table' and 'data.frame': 1 obs. of 7 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "ONEUNIT"
## $ plotsampcnt:'data.frame': 2 obs. of 3 variables:
## $ condsampcnt:'data.frame': 4 obs. of 3 variables:
## $ npixels :Classes 'data.table' and 'data.frame': 1 obs. of 2 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## $ npixelvar : chr "npixels"
## $ states : chr "Wyoming"
## $ invyrs :List of 1
## $ estvar.area: chr "CONDPROP_ADJ"
## $ adj : chr "plot"
## $ treex :Classes 'data.table' and 'data.frame': 1691 obs. of 21 variables:
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "PLT_CN"
## $ tuniqueid : chr "PLT_CN"
## $ adjtree : logi FALSE
## $ prednames : chr [1:5] "dem" "slp" "asp_cos" "asp_sin" ...
## $ predfac : chr "fornf"Now that we’ve created our population dataset, we can move on to estimation.
modMAarea
Example 2: Area of forest land, Wyoming, 2011-2013
View Example
In this example, we look at estimating the area of forest land in
Wyoming from 2011 to 2013 summed to the population unit
(sumunit = TRUE) with the generalized regression estimator
(MAmethod = "greg"). FIESTA returns raw data
for area of forest land, Wyoming, 2011-2013 (sum estimation units).
area1 <- modMAarea(
MApopdat = MApopdat, # pop - population calculations for WY, post-stratification
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
)We can look at the structure of this output with str and
the estimates below. Note that again FIESTA outputs a
list.
str(area1, max.level = 2)output
## List of 4
## $ est :Classes 'data.table' and 'data.frame': 1 obs. of 3 variables:
## ..$ ONEUNIT : num 1
## ..$ Estimate : num 660395
## ..$ Percent Sampling Error: num 9.42
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "ONEUNIT"
## $ raw :List of 9
## ..$ unit_totest :'data.frame': 1 obs. of 18 variables:
## ..$ domdat :'data.frame': 66 obs. of 19 variables:
## ..$ module : chr "MA"
## ..$ esttype : chr "AREA"
## ..$ MAmethod : chr "greg"
## ..$ predselectlst:List of 1
## ..$ rowvar : chr "TOTAL"
## ..$ colvar : chr "NONE"
## ..$ areaunits : chr "acres"
## $ statecd: int 56
## $ invyr : int [1:3] 2011 2012 2013
area1$estoutput
## ONEUNIT Estimate Percent Sampling Error
## 1: 1 660395.1 9.42Example 3: Area of forest land, Wyoming, 2011-2013, using the Elastic Net for variable selection
View Example
Here, we fit the same model as the above example, but rather than
using "greg" are our model-assisted method, we can use
"gregEN" where the EN stands for “elastic net”. The elastic
net performs variable selection for us, grabbing predictors it finds to
be most useful in the model.
area2 <- modMAarea(
MApopdat = MApopdat, # pop - population calculations for WY, post-stratification
MAmethod = "gregEN", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
)We can again see that the structure of the list is very similar to that in the above example:
str(area2, max.level = 2)output
## List of 4
## $ est :Classes 'data.table' and 'data.frame': 1 obs. of 3 variables:
## ..$ ONEUNIT : num 1
## ..$ Estimate : num 666568
## ..$ Percent Sampling Error: num 10.2
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "ONEUNIT"
## $ raw :List of 9
## ..$ unit_totest :'data.frame': 1 obs. of 18 variables:
## ..$ domdat :'data.frame': 66 obs. of 19 variables:
## ..$ module : chr "MA"
## ..$ esttype : chr "AREA"
## ..$ MAmethod : chr "gregEN"
## ..$ predselectlst:List of 1
## ..$ rowvar : chr "TOTAL"
## ..$ colvar : chr "NONE"
## ..$ areaunits : chr "acres"
## $ statecd: int 56
## $ invyr : int [1:3] 2011 2012 2013However now the raw list has an item call
predselectlst. We can look at this item now:
area2$raw$predselectlst$totestoutput
## ONEUNIT TOTAL dem slp asp_cos asp_sin fornf2
## 1: 1 1 -0.0002711995 0 0 0 0Notably, we can see that dem, slp,
asp_cos, and asp_sin were removed from the
model.
Example 4: Area by forest type on forest land, Wyoming, 2011-2013
View Example
In this example, we look at adding rows to the output and include returntitle=TRUE to return title information.
area3 <- modMAarea(
MApopdat = MApopdat, # pop - population calculations for WY, post-stratification
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
rowvar = "FORTYPCD", # est - row domain
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.
str(area3, max.level = 1)output
## List of 5
## $ est :'data.frame': 8 obs. of 3 variables:
## $ titlelst:List of 9
## $ raw :List of 10
## $ statecd : int 56
## $ invyr : int [1:3] 2011 2012 2013And the estimates:
## Estimate and percent sampling error of estimate
area3$estoutput
## Forest type Estimate Percent Sampling Error
## 1 201 42653.6 60.56
## 2 265 47917.2 67.78
## 3 266 71039.9 43.05
## 4 268 39520.3 68.8
## 5 281 407596.6 16.05
## 6 901 24905.9 78.03
## 7 999 26761.6 63.54
## 8 Total 660395.1 9.42Along with raw data and titles:
## Raw data (list object) for estimate
raw3 <- area3$raw # extract raw data list object from output
names(raw3)output
## [1] "unit_totest" "unit_rowest" "domdat" "module"
## [5] "esttype" "MAmethod" "predselectlst" "rowvar"
## [9] "colvar" "areaunits"
head(raw3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS AREAUSED est
## 1 1 0.5936603 0.003129767 56 37 1112412 1112412 660395.1
## est.var est.se est.cv pse CI99left CI99right CI95left CI95right
## 1 3872964796 62233.15 0.09423624 9.423624 500093.1 820697 538420.4 782369.8
## CI68left CI68right
## 1 598506.8 722283.3
raw3$totest # estimates for population (i.e., WY)output
## NULL
head(raw3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Forest type nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS
## 1 1 201 0.03834336 0.0005392280 56 2 1112412
## 2 1 265 0.04307505 0.0008523855 56 4 1112412
## 3 1 266 0.06386110 0.0007557639 56 5 1112412
## 4 1 268 0.03552668 0.0005973938 56 3 1112412
## 5 1 281 0.36640782 0.0034566717 56 23 1112412
## 6 1 901 0.02238906 0.0003052203 56 1 1112412
## AREAUSED est est.var est.se est.cv pse CI99left CI99right
## 1 1112412 42653.63 667273775 25831.64 0.6056142 60.56142 0.0 109191.53
## 2 1112412 47917.21 1054793996 32477.59 0.6777855 67.77855 0.0 131573.95
## 3 1112412 71039.88 935228490 30581.51 0.4304836 43.04836 0.0 149812.62
## 4 1112412 39520.32 739251639 27189.18 0.6879799 68.79799 0.0 109555.01
## 5 1112412 407596.59 4277497216 65402.58 0.1604591 16.04591 239130.7 576062.46
## 6 1112412 24905.87 377698351 19434.46 0.7803167 78.03167 0.0 74965.72
## CI95left CI95right CI68left CI68right
## 1 0.00 93282.72 16965.145 68342.11
## 2 0.00 111572.13 15619.617 80214.81
## 3 11101.23 130978.53 40627.860 101451.90
## 4 0.00 92810.14 12481.820 66558.81
## 5 279409.89 535783.28 342556.478 472636.70
## 6 0.00 62996.71 5579.111 44232.62
head(raw3$rowest) # estimates by row for population (i.e., WY)output
## NULL
## Titles (list object) for estimate
titlelst3 <- area3$titlelst
names(titlelst3)output
## [1] "title.estpse" "title.unitvar" "title.ref" "outfn.estpse"
## [5] "outfn.rawdat" "outfn.param" "title.rowvar" "title.row"
## [9] "title.unitsn"
titlelst3output
## $title.estpse
## [1] "Area, in acres, and percent sampling error on forest land by forest type"
##
## $title.unitvar
## [1] "ONEUNIT"
##
## $title.ref
## [1] ", 2011-2013"
##
## $outfn.estpse
## [1] "area_FORTYPCD_forestland"
##
## $outfn.rawdat
## [1] "area_FORTYPCD_forestland_rawdata"
##
## $outfn.param
## [1] "area_FORTYPCD_forestland_parameters"
##
## $title.rowvar
## [1] "Forest type"
##
## $title.row
## [1] "Area, in acres, on forest land by forest type; , 2011-2013"
##
## $title.unitsn
## [1] "acres"Example 5: Area by forest type and stand-size class on forest land, Wyoming, 2011-2013
View Example
In this example, we look at adding rows and columns to output,
including FIA names. 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.
## Area of forest land by forest type and stand-size class, Wyoming, 2011-2013
area4 <- modMAarea(
MApopdat = MApopdat, # pop - population calculations for WY, post-stratification
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
rowvar = "FORTYPCD", # est - row domain
colvar = "STDSZCD", # est - column domain
savedata = TRUE, # out - save data to outfolder
returntitle = TRUE, # out - return title information
table_opts = list(
row.FIAname = TRUE, # table - row domain names
col.FIAname = TRUE, # table - column domain names
allin1 = TRUE # table - return output with est(pse)
),
savedata_opts = list(
outfolder = outfolder, # save - outfolder for saving data
outfn.pre = "WY" # save - prefix for output files
)
)
area4$estoutput
## Forest type Large diameter Medium diameter
## 1: Douglas-fir 23,330.7 ( 79.00) -- ( --)
## 2: Engelmann spruce 36,892.1 ( 74.16) 11,025.2 (175.15)
## 3: Engelmann spruce / subalpine fir 40,287.9 ( 68.99) 14,443.2 ( 66.30)
## 4: Subalpine fir 26,255.3 ( 75.29) -- ( --)
## 5: Lodgepole pine 140,798.2 ( 33.20) 230,288.3 ( 20.42)
## 6: Aspen -- ( --) -- ( --)
## 7: Nonstocked -- ( --) -- ( --)
## 8: Total 267,564.1 ( 21.17) 255,756.7 ( 19.38)
## Small diameter Nonstocked Total
## 1: 19,322.9 ( 99.12) -- ( --) 42,653.6 ( 60.56)
## 2: -- ( --) -- ( --) 47,917.2 ( 67.78)
## 3: 16,308.7 ( 71.46) -- ( --) 71,039.9 ( 43.05)
## 4: 13,265.0 (141.72) -- ( --) 39,520.3 ( 68.80)
## 5: 36,510.1 ( 74.46) -- ( --) 407,596.6 ( 16.05)
## 6: 24,905.9 ( 78.03) -- ( --) 24,905.9 ( 78.03)
## 7: -- ( --) 26,761.6 ( 63.54) 26,761.6 ( 63.54)
## 8: 110,312.6 ( 37.41) 26,761.6 ( 63.54) 660,395.1 ( 9.42)We can again look at the output list, estimates, raw data, and titles:
## Look at output list
names(area4)output
## [1] "est" "pse" "titlelst" "raw" "statecd" "invyr"
## Estimate and percent sampling error of estimate
head(area4$est)output
## Forest type Large diameter Medium diameter
## 1: Douglas-fir 23,330.7 ( 79.00) -- ( --)
## 2: Engelmann spruce 36,892.1 ( 74.16) 11,025.2 (175.15)
## 3: Engelmann spruce / subalpine fir 40,287.9 ( 68.99) 14,443.2 ( 66.30)
## 4: Subalpine fir 26,255.3 ( 75.29) -- ( --)
## 5: Lodgepole pine 140,798.2 ( 33.20) 230,288.3 ( 20.42)
## 6: Aspen -- ( --) -- ( --)
## Small diameter Nonstocked Total
## 1: 19,322.9 ( 99.12) -- ( --) 42,653.6 ( 60.56)
## 2: -- ( --) -- ( --) 47,917.2 ( 67.78)
## 3: 16,308.7 ( 71.46) -- ( --) 71,039.9 ( 43.05)
## 4: 13,265.0 (141.72) -- ( --) 39,520.3 ( 68.80)
## 5: 36,510.1 ( 74.46) -- ( --) 407,596.6 ( 16.05)
## 6: 24,905.9 ( 78.03) -- ( --) 24,905.9 ( 78.03)
## Raw data (list object) for estimate
raw4 <- area4$raw # extract raw data list object from output
names(raw4)output
## [1] "unit_totest" "unit_rowest" "unit_colest" "unit_grpest"
## [5] "domdat" "module" "esttype" "MAmethod"
## [9] "predselectlst" "rowvar" "colvar" "areaunits"
head(raw4$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS AREAUSED est
## 1 1 0.5936603 0.003129767 56 37 1112412 1112412 660395.1
## est.var est.se est.cv pse CI99left CI99right CI95left CI95right
## 1 3872964796 62233.15 0.09423624 9.423624 500093.1 820697 538420.4 782369.8
## CI68left CI68right
## 1 598506.8 722283.3
head(raw4$totest) # estimates for population (i.e., WY)output
## NULL
head(raw4$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Forest type nhat nhat.var NBRPLT
## 1 1 Douglas-fir 0.03834336 0.0005392280 56
## 2 1 Engelmann spruce 0.04307505 0.0008523855 56
## 3 1 Engelmann spruce / subalpine fir 0.06386110 0.0007557639 56
## 4 1 Subalpine fir 0.03552668 0.0005973938 56
## 5 1 Lodgepole pine 0.36640782 0.0034566717 56
## 6 1 Aspen 0.02238906 0.0003052203 56
## NBRPLT.gt0 FORTYPCD ACRES_GIS AREAUSED est est.var est.se
## 1 2 201 1112412 1112412 42653.63 667273775 25831.64
## 2 4 265 1112412 1112412 47917.21 1054793996 32477.59
## 3 5 266 1112412 1112412 71039.88 935228490 30581.51
## 4 3 268 1112412 1112412 39520.32 739251639 27189.18
## 5 23 281 1112412 1112412 407596.59 4277497216 65402.58
## 6 1 901 1112412 1112412 24905.87 377698351 19434.46
## est.cv pse CI99left CI99right CI95left CI95right CI68left
## 1 0.6056142 60.56142 0.0 109191.53 0.00 93282.72 16965.145
## 2 0.6777855 67.77855 0.0 131573.95 0.00 111572.13 15619.617
## 3 0.4304836 43.04836 0.0 149812.62 11101.23 130978.53 40627.860
## 4 0.6879799 68.79799 0.0 109555.01 0.00 92810.14 12481.820
## 5 0.1604591 16.04591 239130.7 576062.46 279409.89 535783.28 342556.478
## 6 0.7803167 78.03167 0.0 74965.72 0.00 62996.71 5579.111
## CI68right
## 1 68342.11
## 2 80214.81
## 3 101451.90
## 4 66558.81
## 5 472636.70
## 6 44232.62
head(raw4$rowest) # estimates by row for population (i.e., WY)output
## NULL
head(raw4$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Stand-size class nhat nhat.var NBRPLT NBRPLT.gt0 STDSZCD
## 1 1 Large diameter 0.24052605 0.0025935705 56 18 1
## 2 1 Medium diameter 0.22991178 0.0019845421 56 14 2
## 3 1 Small diameter 0.09916524 0.0013761149 56 6 3
## 4 1 Nonstocked 0.02405726 0.0002336838 56 1 5
## ACRES_GIS AREAUSED est est.var est.se est.cv pse CI99left
## 1 1112412 1112412 267564.15 3209443017 56651.95 0.2117322 21.17322 121638.40
## 2 1112412 1112412 255756.70 2455793975 49555.97 0.1937621 19.37621 128108.99
## 3 1112412 1112412 110312.64 1702888801 41266.07 0.3740829 37.40829 4018.28
## 4 1112412 1112412 26761.59 289174599 17005.13 0.6354306 63.54306 0.00
## CI99right CI95left CI95right CI68left CI68right
## 1 413489.89 156528.37 378599.92 211226.172 323902.12
## 2 383404.42 158628.79 352884.62 206475.381 305038.03
## 3 216607.00 29432.62 191192.66 69275.269 151350.01
## 4 70563.91 0.00 60091.04 9850.701 43672.48
head(raw4$colest) # estimates by column for population (i.e., WY)output
## NULL
head(raw4$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT grpvar nhat
## 1 1 Douglas-fir#Large diameter 0.020973074
## 2 1 Douglas-fir#Small diameter 0.017370285
## 3 1 Engelmann spruce#Large diameter 0.033164010
## 4 1 Engelmann spruce#Medium diameter 0.009911038
## 5 1 Engelmann spruce / subalpine fir#Large diameter 0.036216716
## 6 1 Engelmann spruce / subalpine fir#Medium diameter 0.012983698
## nhat.var NBRPLT NBRPLT.gt0 Forest type
## 1 2.745298e-04 56 1 Douglas-fir
## 2 2.964436e-04 56 1 Douglas-fir
## 3 6.048455e-04 56 3 Engelmann spruce
## 4 3.013453e-04 56 1 Engelmann spruce
## 5 6.243626e-04 56 3 Engelmann spruce / subalpine fir
## 6 7.410764e-05 56 1 Engelmann spruce / subalpine fir
## Stand-size class STDSZCD FORTYPCD ACRES_GIS AREAUSED est est.var
## 1 Large diameter 1 201 1112412 1112412 23330.71 339719955
## 2 Small diameter 3 201 1112412 1112412 19322.92 366837446
## 3 Large diameter 1 265 1112412 1112412 36892.05 748472910
## 4 Medium diameter 2 265 1112412 1112412 11025.16 372903083
## 5 Large diameter 1 266 1112412 1112412 40287.92 772624533
## 6 Medium diameter 2 266 1112412 1112412 14443.23 91705337
## est.se est.cv pse CI99left CI99right CI95left CI95right CI68left
## 1 18431.49 0.7900101 79.00101 0 70807.09 0 59455.77 5001.3621
## 2 19153.00 0.9912064 99.12064 0 68657.78 0 56862.11 276.0668
## 3 27358.23 0.7415752 74.15752 0 107362.19 0 90513.21 9685.4439
## 4 19310.70 1.7515118 175.15118 0 60766.22 0 48873.43 0.0000
## 5 27796.12 0.6899369 68.99369 0 111885.99 0 94767.33 12645.8475
## 6 9576.29 0.6630299 66.30299 0 39110.12 0 33212.41 4920.0086
## CI68right
## 1 41660.05
## 2 38369.77
## 3 64098.66
## 4 30228.84
## 5 67930.00
## 6 23966.44
head(raw4$grpest) # estimates by row and column for population (i.e., WY)output
## NULL
## Titles (list object) for estimate
titlelst4 <- area4$titlelst
names(titlelst4)output
## [1] "title.estpse" "title.unitvar" "title.ref" "outfn.estpse"
## [5] "outfn.rawdat" "outfn.param" "title.rowvar" "title.row"
## [9] "title.colvar" "title.col" "title.unitsn"
titlelst4output
## $title.estpse
## [1] "Area, in acres (percent sampling error), by forest type and stand-size class on forest land"
##
## $title.unitvar
## [1] "ONEUNIT"
##
## $title.ref
## [1] ", 2011-2013"
##
## $outfn.estpse
## [1] "WY_area_FORTYPNM_STDSZNM_forestland"
##
## $outfn.rawdat
## [1] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata"
##
## $outfn.param
## [1] "WY_area_FORTYPNM_STDSZNM_forestland_parameters"
##
## $title.rowvar
## [1] "Forest type"
##
## $title.row
## [1] "Area, in acres (percent sampling error), by forest type on forest land; , 2011-2013"
##
## $title.colvar
## [1] "Stand-size class"
##
## $title.col
## [1] "Area, in acres (percent sampling error), by stand-size class on forest land; , 2011-2013"
##
## $title.unitsn
## [1] "acres"
## List output files in outfolder
list.files(outfolder, pattern = "WY_area")output
## [1] "WY_area_FORTYPNM_STDSZNM_forestland_modMA_mase_greg.csv"
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_area")output
## [1] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_domdat.csv"
## [2] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_colest.csv"
## [3] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_grpest.csv"
## [4] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_rowest.csv"
## [5] "WY_area_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_totest.csv"
modMAtree
We will set our estimate variable and filter now. We set
estvar to "VOLCFNET" for net cubic foot
volume, and filter with estvar.filter set to
"STATUSCD == 1" so we only consider live trees in our
estimation.
estvar <- "VOLCFNET"
live_trees <- "STATUSCD == 1"Example 6: Net cubic-foot volume of live trees, Wyoming, 2011-2013
View Example
We now will generate estimates by estimation unit (i.e., ESTN_UNIT)
and sum to population (i.e., WY) with modMAtree.
## Return raw data and titles
## Total net cubic-foot volume of live trees (at least 5 inches diameter), Wyoming, 2011-2013
tree1 <- modMAtree(
MApopdat = MApopdat, # pop - population calculations
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
estvar = estvar, # est - net cubic-foot volume
estvar.filter = live_trees, # est - live trees only
returntitle = TRUE # out - return title information
)
names(tree1)output
## [1] "est" "titlelst" "raw" "statecd" "invyr"
tree1$raw$unit_totestoutput
## ONEUNIT nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS AREAUSED est
## 1 1 1165.825 30923.7 56 34 1112412 1112412 1296877648
## est.var est.se est.cv pse CI99left CI99right CI95left
## 1 3.826688e+16 195619213 0.1508386 15.08386 792995947 1800759350 913471036
## CI95right CI68left CI68right
## 1 1680284261 1102342580 1491412717Example 7: Net cubic-foot volume of live trees, Wyoming, 2011-2013, using the Elastic Net for variable selection
View Example
Here, we fit the same model as the above example, but rather than
using "greg" are our model-assisted method, we can use
"gregEN" where the EN stands for “elastic net”. The elastic
net performs variable selection for us, grabbing predictors it finds to
be most useful in the model.
## Return raw data and titles
## Total net cubic-foot volume of live trees (at least 5 inches diameter), Wyoming, 2011-2013
tree2 <- modMAtree(
MApopdat = MApopdat, # pop - population calculations
MAmethod = "gregEN", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
estvar = estvar, # est - net cubic-foot volume
estvar.filter = live_trees, # est - live trees only
returntitle = TRUE # out - return title information
)We can again see that the structure of the list is very similar to that in the above example:
str(tree2, max.level = 2)output
## List of 5
## $ est :Classes 'data.table' and 'data.frame': 1 obs. of 3 variables:
## ..$ ONEUNIT : num 1
## ..$ Estimate : num 1.31e+09
## ..$ Percent Sampling Error: num 14.9
## ..- attr(*, ".internal.selfref")=<externalptr>
## ..- attr(*, "sorted")= chr "ONEUNIT"
## $ titlelst:List of 10
## ..$ title.estpse : chr "Net volume of live trees (at least 5 in dia), in cubic feet, and percent sampling error on forest land live"
## ..$ title.yvar : chr "Net volume, in cubic feet"
## ..$ title.estvar : chr "Net volume of live trees (at least 5 in dia)"
## ..$ title.unitvar: chr "ONEUNIT"
## ..$ title.ref : chr "Wyoming, 2011-2013"
## ..$ outfn.estpse : chr "tree_VOLCFNET_live_forestland"
## ..$ outfn.rawdat : chr "tree_VOLCFNET_live_forestland_rawdata"
## ..$ outfn.param : chr "tree_VOLCFNET_live_forestland_parameters"
## ..$ title.tot : chr "Net volume of live trees (at least 5 in dia), in cubic feet, on forest land live; Wyoming, 2011-2013"
## ..$ title.unitsn : chr "cubic feet"
## $ raw :List of 13
## ..$ unit_totest :'data.frame': 1 obs. of 18 variables:
## ..$ domdat :'data.frame': 66 obs. of 19 variables:
## ..$ plotweights :List of 1
## ..$ estvar : chr "VOLCFNET"
## ..$ estvar.filter: chr "STATUSCD == 1"
## ..$ module : chr "MA"
## ..$ esttype : chr "TREE"
## ..$ MAmethod : chr "gregEN"
## ..$ predselectlst:List of 1
## ..$ rowvar : chr "TOTAL"
## ..$ colvar : chr "NONE"
## ..$ areaunits : chr "acres"
## ..$ estunits : num 1
## $ statecd : int 56
## $ invyr : int [1:3] 2011 2012 2013However now the raw list has an item call
predselectlst. We can look at this item now:
tree2$raw$predselectlstoutput
## $totest
## ONEUNIT TOTAL dem slp asp_cos asp_sin fornf2
## 1: 1 1 0.6497733 -63.48846 32.55884 -536.5366 -1369.907Notably, we can see that [INSERT CORRECT PREDS] dem,
slp, asp_cos, and asp_sin were
removed from the model.
Example 8: Net cubic-foot volume of live trees by forest type, Wyoming, 2011-2013
View Example
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.
tree3 <- modMAtree(
MApopdat = MApopdat, # pop - population calculations
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
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:
## Look at output list
names(tree3)output
## [1] "est" "titlelst" "raw" "statecd" "invyr"
## Estimate and percent sampling error of estimate
tree3$estoutput
## Forest type Estimate Percent Sampling Error
## 1 201 15329084.7 63.95
## 2 265 191806950.6 70.43
## 3 266 135152582.5 66.37
## 4 268 49388166.9 73.19
## 5 281 905200863.5 19.75
## 6 901 -- --
## 7 999 -- --
## 8 Total 1296877648.2 15.08
## Raw data (list object) for estimate
raw3 <- tree3$raw # extract raw data list object from output
names(raw3)output
## [1] "unit_totest" "unit_rowest" "domdat" "plotweights"
## [5] "estvar" "estvar.filter" "module" "esttype"
## [9] "MAmethod" "predselectlst" "rowvar" "colvar"
## [13] "areaunits" "estunits"
head(raw3$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS AREAUSED est
## 1 1 1165.825 30923.7 56 34 1112412 1112412 1296877648
## est.var est.se est.cv pse CI99left CI99right CI95left
## 1 3.826688e+16 195619213 0.1508386 15.08386 792995947 1800759350 913471036
## CI95right CI68left CI68right
## 1 1680284261 1102342580 1491412717
head(raw3$totest) # estimates for population (i.e., WY)output
## NULL
head(raw3$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Forest type nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS
## 1 1 201 13.78004 77.64629 56 2 1112412
## 2 1 265 172.42433 14746.63678 56 4 1112412
## 3 1 266 121.49504 6501.78676 56 5 1112412
## 4 1 268 44.39736 1055.90466 56 2 1112412
## 5 1 281 813.72781 25826.51205 56 23 1112412
## 6 1 901 0.00000 0.00000 56 0 1112412
## AREAUSED est est.var est.se est.cv pse CI99left
## 1 1112412 15329085 9.608428e+13 9802259 0.6394549 63.94549 0
## 2 1112412 191806951 1.824839e+16 135086607 0.7042842 70.42842 0
## 3 1112412 135152582 8.045709e+15 89697877 0.6636786 66.36786 0
## 4 1112412 49388167 1.306641e+15 36147491 0.7319059 73.19059 0
## 5 1112412 905200864 3.195931e+16 178771664 0.1974939 19.74939 444715574
## 6 1112412 0 0.000000e+00 0 NaN NaN 0
## CI99right CI95left CI95right CI68left CI68right
## 1 40578030 0 34541159 5581151 25077018
## 2 539766992 0 456571836 57469009 326144892
## 3 366199002 0 310957191 45951822 224353343
## 4 142497933 0 120235947 13441010 85335324
## 5 1365686153 554814841 1255586886 727419973 1082981754
## 6 0 0 0 0 0
head(raw3$rowest) # estimates by row for population (i.e., WY)output
## NULL
## Titles (list object) for estimate
titlelst3 <- tree3$titlelst
names(titlelst3)output
## [1] "title.estpse" "title.yvar" "title.estvar" "title.unitvar"
## [5] "title.ref" "outfn.estpse" "outfn.rawdat" "outfn.param"
## [9] "title.rowvar" "title.row" "title.unitsn"
titlelst3output
## $title.estpse
## [1] "Net volume of live trees (at least 5 in dia), in cubic feet, and percent sampling error on forest land by forest type live"
##
## $title.yvar
## [1] "Net volume, in cubic feet"
##
## $title.estvar
## [1] "Net volume of live trees (at least 5 in dia)"
##
## $title.unitvar
## [1] "ONEUNIT"
##
## $title.ref
## [1] "Wyoming, 2011-2013"
##
## $outfn.estpse
## [1] "tree_VOLCFNET_live_FORTYPCD_forestland"
##
## $outfn.rawdat
## [1] "tree_VOLCFNET_live_FORTYPCD_forestland_rawdata"
##
## $outfn.param
## [1] "tree_VOLCFNET_live_FORTYPCD_forestland_parameters"
##
## $title.rowvar
## [1] "Forest type"
##
## $title.row
## [1] "Net volume of live trees (at least 5 in dia), in cubic feet, on forest land by forest type live; Wyoming, 2011-2013"
##
## $title.unitsn
## [1] "cubic feet"We can also create a simple barplot from the output:
## Create barplot
datBarplot(
raw3$unit_rowest,
xvar = titlelst3$title.rowvar,
yvar = "est"
)plot
And a fancier barplot:
## Create fancier barplot
datBarplot(
raw3$unit_rowest,
xvar = titlelst3$title.rowvar,
yvar = "est",
errbars = TRUE,
sevar = "est.se",
main = FIESTAutils::wraptitle(titlelst3$title.row, 75),
ylabel = titlelst3$title.yvar,
divideby = "million"
)plot
Example 9: Net cubic-foot volume of live trees by forest type and stand-size class, Wyoming, 2011-2013
View Example
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.
tree4 <- modMAtree(
MApopdat = MApopdat, # pop - population calculations
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
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:
## Look at output list from modGBarea()
names(tree4)output
## [1] "est" "pse" "titlelst" "raw" "statecd" "invyr"
## Estimate and percent sampling error of estimate
tree4$estoutput
## Forest type Large diameter
## 1: Douglas-fir 11,866,612.7 ( 79.00)
## 2: Engelmann spruce 149,865,021.5 ( 78.94)
## 3: Engelmann spruce / subalpine fir 113,586,516.8 ( 78.74)
## 4: Subalpine fir 49,388,166.9 ( 73.19)
## 5: Lodgepole pine 456,907,508.9 ( 34.84)
## 6: Aspen -- ( --)
## 7: Nonstocked -- ( --)
## 8: Total 781,613,826.8 ( 23.93)
## Medium diameter Small diameter Nonstocked
## 1: -- ( --) 3,462,472.0 ( 99.12) -- ( --)
## 2: 41,941,929.1 (175.15) -- ( --) -- ( --)
## 3: 16,395,905.1 ( 66.30) 5,170,160.5 ( 71.46) -- ( --)
## 4: -- ( --) -- ( --) -- ( --)
## 5: 442,719,166.9 ( 21.99) 5,574,187.7 ( 77.89) -- ( --)
## 6: -- ( --) -- ( --) -- ( --)
## 7: -- ( --) -- ( --) -- ( --)
## 8: 501,057,001.1 ( 23.52) 14,206,820.3 ( 43.06) -- ( --)
## Total
## 1: 15,329,084.7 ( 63.95)
## 2: 191,806,950.6 ( 70.43)
## 3: 135,152,582.5 ( 66.37)
## 4: 49,388,166.9 ( 73.19)
## 5: 905,200,863.5 ( 19.75)
## 6: -- ( --)
## 7: -- ( --)
## 8: 1,296,877,648.2 ( 15.08)
## Raw data (list object) for estimate
raw4 <- tree4$raw # extract raw data list object from output
names(raw4)output
## [1] "unit_totest" "unit_rowest" "unit_colest" "unit_grpest"
## [5] "domdat" "plotweights" "estvar" "estvar.filter"
## [9] "module" "esttype" "MAmethod" "predselectlst"
## [13] "rowvar" "colvar" "areaunits" "estunits"
head(raw4$unit_totest) # estimates by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT nhat nhat.var NBRPLT NBRPLT.gt0 ACRES_GIS AREAUSED est
## 1 1 1165.825 30923.7 56 34 1112412 1112412 1296877648
## est.var est.se est.cv pse CI99left CI99right CI95left
## 1 3.826688e+16 195619213 0.1508386 15.08386 792995947 1800759350 913471036
## CI95right CI68left CI68right
## 1 1680284261 1102342580 1491412717
head(raw4$totest) # estimates for population (i.e., WY)output
## NULL
head(raw4$unit_rowest) # estimates by row, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Forest type nhat nhat.var NBRPLT
## 1 1 Douglas-fir 13.78004 77.64629 56
## 2 1 Engelmann spruce 172.42433 14746.63678 56
## 3 1 Engelmann spruce / subalpine fir 121.49504 6501.78676 56
## 4 1 Subalpine fir 44.39736 1055.90466 56
## 5 1 Lodgepole pine 813.72781 25826.51205 56
## 6 1 Aspen 0.00000 0.00000 56
## NBRPLT.gt0 FORTYPCD ACRES_GIS AREAUSED est est.var est.se
## 1 2 201 1112412 1112412 15329085 9.608428e+13 9802259
## 2 4 265 1112412 1112412 191806951 1.824839e+16 135086607
## 3 5 266 1112412 1112412 135152582 8.045709e+15 89697877
## 4 2 268 1112412 1112412 49388167 1.306641e+15 36147491
## 5 23 281 1112412 1112412 905200864 3.195931e+16 178771664
## 6 0 901 1112412 1112412 0 0.000000e+00 0
## est.cv pse CI99left CI99right CI95left CI95right CI68left
## 1 0.6394549 63.94549 0 40578030 0 34541159 5581151
## 2 0.7042842 70.42842 0 539766992 0 456571836 57469009
## 3 0.6636786 66.36786 0 366199002 0 310957191 45951822
## 4 0.7319059 73.19059 0 142497933 0 120235947 13441010
## 5 0.1974939 19.74939 444715574 1365686153 554814841 1255586886 727419973
## 6 NaN NaN 0 0 0 0 0
## CI68right
## 1 25077018
## 2 326144892
## 3 224353343
## 4 85335324
## 5 1082981754
## 6 0
head(raw4$rowest) # estimates by row for population (i.e., WY)output
## NULL
head(raw4$unit_colest) # estimates by column, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT Stand-size class nhat nhat.var NBRPLT NBRPLT.gt0 STDSZCD
## 1 1 Large diameter 702.62958 28262.84167 56 18 1
## 2 1 Medium diameter 450.42380 11221.16423 56 14 2
## 3 1 Small diameter 12.77118 30.23649 56 4 3
## 4 1 Nonstocked 0.00000 0.00000 56 0 5
## ACRES_GIS AREAUSED est est.var est.se est.cv pse
## 1 1112412 1112412 781613827 3.497417e+16 187013826 0.2392663 23.92663
## 2 1112412 1112412 501057001 1.388576e+16 117837837 0.2351785 23.51785
## 3 1112412 1112412 14206820 3.741649e+13 6116902 0.4305609 43.05609
## 4 1112412 1112412 0 0.000000e+00 0 NaN NaN
## CI99left CI99right CI95left CI95right CI68left CI68right
## 1 299898133 1263329520 415073463 1148154191 595636453 967591201
## 2 197526847 804587155 270099085 732014918 383872235 618241767
## 3 0 29962915 2217913 26195727 8123819 20289821
## 4 0 0 0 0 0 0
head(raw4$colest) # estimates by column for population (i.e., WY)output
## NULL
head(raw4$unit_grpest) # estimates by row and column, by estimation unit (i.e., ESTN_UNIT)output
## ONEUNIT grpvar nhat
## 1 1 Douglas-fir#Large diameter 10.66746
## 2 1 Douglas-fir#Small diameter 3.11258
## 3 1 Engelmann spruce#Large diameter 134.72075
## 4 1 Engelmann spruce#Medium diameter 37.70358
## 5 1 Engelmann spruce / subalpine fir#Large diameter 102.10828
## 6 1 Engelmann spruce / subalpine fir#Medium diameter 14.73905
## nhat.var NBRPLT NBRPLT.gt0 Forest type
## 1 71.021071 56 1 Douglas-fir
## 2 9.518513 56 1 Douglas-fir
## 3 11310.014570 56 3 Engelmann spruce
## 4 4361.052579 56 1 Engelmann spruce
## 5 6464.613685 56 3 Engelmann spruce / subalpine fir
## 6 95.500443 56 1 Engelmann spruce / subalpine fir
## Stand-size class STDSZCD FORTYPCD ACRES_GIS AREAUSED est est.var
## 1 Large diameter 1 201 1112412 1112412 11866613 8.788582e+13
## 2 Small diameter 3 201 1112412 1112412 3462472 1.177879e+13
## 3 Large diameter 1 265 1112412 1112412 149865022 1.399570e+16
## 4 Medium diameter 2 265 1112412 1112412 41941929 5.396634e+15
## 5 Large diameter 1 266 1112412 1112412 113586517 7.999709e+15
## 6 Medium diameter 2 266 1112412 1112412 16395905 1.181781e+14
## est.se est.cv pse CI99left CI99right CI95left CI95right
## 1 9374744 0.7900101 79.00101 0 36014353 0 30240773
## 2 3432024 0.9912064 99.12064 0 12302781 0 10189116
## 3 118303443 0.7894000 78.94000 0 454594497 0 381735510
## 4 73461783 1.7515118 175.15118 0 231166942 0 185924378
## 5 89441092 0.7874270 78.74270 0 343971502 0 288887835
## 6 10870975 0.6630299 66.30299 0 44397681 0 37702625
## CI68left CI68right
## 1 2543824.72 21189401
## 2 49468.39 6875476
## 3 32217229.75 267512813
## 4 0.00 114996578
## 5 24641118.08 202531916
## 6 5585178.40 27206632
head(raw4$grpest) # estimates by row and column for population (i.e., WY)output
## NULL
## Titles (list object) for estimate
titlelst4 <- tree4$titlelst
names(titlelst4)output
## [1] "title.estpse" "title.yvar" "title.estvar" "title.unitvar"
## [5] "title.ref" "outfn.estpse" "outfn.rawdat" "outfn.param"
## [9] "title.rowvar" "title.row" "title.colvar" "title.col"
## [13] "title.unitsn"
titlelst4output
## $title.estpse
## [1] "Net volume of live trees (at least 5 in dia), in cubic feet (percent sampling error), by forest type and stand-size class on forest land live"
##
## $title.yvar
## [1] "Net volume, in cubic feet"
##
## $title.estvar
## [1] "Net volume of live trees (at least 5 in dia)"
##
## $title.unitvar
## [1] "ONEUNIT"
##
## $title.ref
## [1] "Wyoming, 2011-2013"
##
## $outfn.estpse
## [1] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland"
##
## $outfn.rawdat
## [1] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata"
##
## $outfn.param
## [1] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_parameters"
##
## $title.rowvar
## [1] "Forest type"
##
## $title.row
## [1] "Net volume of live trees (at least 5 in dia), in cubic feet (percent sampling error), by forest type on forest land live; Wyoming, 2011-2013"
##
## $title.colvar
## [1] "Stand-size class"
##
## $title.col
## [1] "Net volume of live trees (at least 5 in dia), in cubic feet (percent sampling error), by stand-size class on forest land live; Wyoming, 2011-2013"
##
## $title.unitsn
## [1] "cubic feet"
## List output files in outfolder
list.files(outfolder, pattern = "WY_tree")output
## [1] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_modMA_mase_greg.csv"
list.files(paste0(outfolder, "/rawdata"), pattern = "WY_tree")output
## [1] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_domdat.csv"
## [2] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_colest.csv"
## [3] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_grpest.csv"
## [4] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_rowest.csv"
## [5] "WY_tree_VOLCFNET_live_FORTYPNM_STDSZNM_forestland_rawdata_modMA_mase_greg_unit_totest.csv"Example 10: Number of live trees by species, Wyoming, 2011-2013
View Example
We can use tree domain in estimation output rows:
## Number of live trees (at least 1 inch diameter) by species
tree5 <- modMAtree(
MApopdat = MApopdat, # pop - population calculations
MAmethod = "greg", # est - model-assisted method
landarea = "FOREST", # est - forest land filter
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:
## Look at output list
names(tree5)output
## [1] "est" "titlelst" "raw" "statecd" "invyr"
## Estimate and percent sampling error of estimate
tree5$estoutput
## Species Estimate Percent Sampling Error
## 1 quaking aspen 11153476.9 73.73
## 2 limber pine 4461866.7 90.07
## 3 lodgepole pine 184898322.4 18.59
## 4 subalpine fir 68658860.9 29.14
## 5 Engelmann spruce 50806944.2 29.77
## 6 Douglas-fir 24655901 73.84
## 7 Total 344635372.3 13.07Example 11: Number of live trees (plus seedlings) by species, Wyoming, 2011-2013
View Example
We can also add seedlings.
Note: seedling data are only available for number of trees (estvar = TPA_UNADJ).
Note: must include seedling data in population data calculations.
MApopdat_seed <- modMApop(popTabs = list(tree = WYtree,
cond = WYcond,
seed = WYseed),
pltassgn = WYpltassgn,
auxdat = modeldat)
## Number of live trees by species, including seedlings
tree6 <- modMAtree(
MApopdat = MApopdat_seed, # pop - population calculations
MAmethod = "greg", # est - model-assisted method
estseed = "add", # est - add seedling data
landarea = "FOREST", # est - forest land filter
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:
## Look at output list
names(tree6)output
## [1] "est" "titlelst" "raw" "statecd" "invyr"
## Estimate and percent sampling error of estimate
tree6$estoutput
## Species Estimate Percent Sampling Error
## 1 quaking aspen 105867206.7 60.84
## 2 limber pine 38477186.7 92.27
## 3 lodgepole pine 251790363.6 18.03
## 4 subalpine fir 396440112.5 23.94
## 5 Engelmann spruce 107965157.9 33.92
## 6 Douglas-fir 49218941.8 80.24
## 7 Total 949758969.2 16.07
## Compare estimates with and without seedlings
head(tree5$est)output
## Species Estimate Percent Sampling Error
## 1 quaking aspen 11153476.9 73.73
## 2 limber pine 4461866.7 90.07
## 3 lodgepole pine 184898322.4 18.59
## 4 subalpine fir 68658860.9 29.14
## 5 Engelmann spruce 50806944.2 29.77
## 6 Douglas-fir 24655901 73.84
head(tree6$est)output
## Species Estimate Percent Sampling Error
## 1 quaking aspen 105867206.7 60.84
## 2 limber pine 38477186.7 92.27
## 3 lodgepole pine 251790363.6 18.03
## 4 subalpine fir 396440112.5 23.94
## 5 Engelmann spruce 107965157.9 33.92
## 6 Douglas-fir 49218941.8 80.24