0% found this document useful (0 votes)

9 views69 pages

ModelMap: R Package for GIS Modeling

The ModelMap package for R enables user-friendly modeling, validation, mapping of large geographic areas with a single function. It builds predictive models using random forests or stochastic gradient boosting. It validates these models and creates graphs and tables of validation diagnostics. It applies these models to GIS files to create detailed prediction maps, handling large files by reading in sections.

Uploaded by

Carlos Eduardo Melo Martínez

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

0% found this document useful (0 votes)

9 views69 pages

ModelMap: R Package for GIS Modeling

Uploaded by

Carlos Eduardo Melo Martínez

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

ModelMap: an R Package for Model Creation and Map

Production
Elizabeth A. Freeman, Tracey S. Frescino, Gretchen G. Moisen

September 10, 2018

Abstract
The ModelMap package (Freeman, 2009) for R (R Development Core Team, 2008) enables
user-friendly modeling, validation, and mapping over large geographic areas though a single R
function or GUI interface. It constructs predictive models of continuous or discrete responses
using Random Forests or Stochastic Gradient Boosting. It validates these models with an
independent test set, cross-validation, or (in the case of Random Forest Models) with Out
OF Bag (OOB) predictions on the training data. It creates graphs and tables of the model
validation diagnostics. It applies these models to GIS image files of predictors to create
detailed prediction surfaces. It will handle large predictor files for map making, by reading in
the GIS data in sections,thus keeping memory usage reasonable.

1 Introduction
Maps of tree species presence and silvicultural metrics like basal area are needed throughout the
world for a wide variety of forest land management applications. Knowledge of the probable
location of certain key species of interest as well as their spatial patterns and associations to
other species are vital components to any realistic land management activity. Recently developed
modeling techniques such as Random Forest (Breiman, 2001) and Stochastic Gradient Boosting
(Friedman, 2001, 2002) offer great potential for improving models and increasing map accuracy
(Evans and Cushman, 2009; Moisen et al., 2006).
The R software environment offers sophisticated new modeling techniques, but requires advanced
programming skills to take full advantage of these capabilities. In addition, spatial data files can
be too memory intensive to analyze easily with standard R code. The ModelMap package provides
an interface between several existing R packages to automate and simplify the process of model
building and map construction.
While spatial data is typically manipulated within a Geographic Information System (GIS), the
ModelMap package facilitates modeling and mapping extensive spatial data in the R software
environment. ModelMap has simple to use GUI prompts for non-programmers, but still has the
flexibility to be run at the command line or in batch mode, and the power to take full advantage
of sophisticated new modeling techniques. ModelMap uses the raster package to read and predict
over GIS raster data. Large maps are read in by row, to keep memory usage reasonable.
The current implementation of ModelMap builds predictive models using Random Forests, Quan-
tile Regression Forests, and Conditional Inference Forests. Stochastic Gradient Boosting models
are not currently supported. Random Forest models are constructed using the randomForest
package (Liaw and Wiener, 2002). For more information on Quantile Regression Forests and Con-
ditional Inference Forests see the additional vignette, ”‘Pick Your Flavor of Random Forest”. The
ModelMap package models both continuous and binary response variables. For binary response,
the PresenceAbsence package (Freeman, 2007) package is used for model diagnostics.
Random Forest models are built as an ensemble of classification or regression trees (Breiman et al.,
1984). Classification and regression trees are intuitive methods, often described in graphical or

1
biological terms. Typically shown growing upside down, a tree begins at its root. An observation
passes down the tree through a series of splits, or nodes, at which a decision is made as to which
direction to proceed based on the value of one of the explanatory variables. Ultimately, a terminal
node or leaf is reached and predicted response is given.
Trees partition the explanatory variables into a series of boxes (the leaves) that contain the most
homogeneous collection of outcomes possible. Creating splits is analogous to variable selection in
regression. Trees are typically fit via binary recursive partitioning. The term binary refers to the
fact that the parent node will always be split into exactly two child nodes. The term recursive is
used to indicate that each child node will, in turn, become a parent node, unless it is a terminal
node. To start with a single split is made using one explanatory variable. The variable and the
location of the split are chosen to minimize the impurity of the node at that point. There are
many ways to minimizing the impurity of each node. These are known as splitting rules. Each of
the two regions that result from the initial split are then split themselves according to the same
criteria, and the tree continues to grow until it is no longer possible to create additional splits or
the process is stopped by some user-defined criteria. The tree may then be reduced in size using a
process known as pruning. Overviews of classification and regression trees are provided by De’ath
and Fabricius (2000), Vayssieres et al. (2000), and Moisen (2008).
While classification and regression trees are powerful methods in and of themselves, much work
has been done in the data mining and machine learning fields to improve the predictive ability of
these tools by combining separate tree models into what is often called a committee of experts, or
ensemble. Random Forests and Stochastic Gradient Boosting are two of these newer techniques
that use classification and regression trees as building blocks.
Random Forests — In a Random Forests model, a bootstrap sample of the training data is chosen.
At the root node, a small random sample of explanatory variables is selected and the best split
made using that limited set of variables. At each subsequent node, another small random sample
of the explanatory variables is chosen, and the best split made. The tree continues to be grown
in this fashion until it reaches the largest possible size, and is left un-pruned. The whole process,
starting with a new bootstrap sample, is repeated a large number of times. As in committee
models, the final prediction is a (weighted) plurality vote or average from prediction of all the
trees in the collection.
Stochastic Gradient Boosting — Stochastic gradient boosting is another ensemble technique in
which many small classification or regression trees are built sequentially from pseudo-residuals
from the previous tree. At each iteration, a tree is built from a random sub-sample of the dataset
(selected without replacement) producing an incremental improvement in the model. Ultimately,
all the small trees are stacked together as a weighted sum of terms. The overall model accuracy
gets progressively better with each additional term.

2 Package Overview
The ModelMap package for R enables user-friendly modeling, diagnostics, and mapping over large
geographic areas though simple R function calls: [Link](), [Link](), and
[Link](). The function [Link]() constructs predictive models of continuous or
discrete responses using Random Forests. The function [Link]() validates these
models with an independent test set, cross-validation, or (in the case of Random Forest Models)
with Out OF Bag (OOB) predictions on the training data. This function also creates graphs and
tables of the basic model validation diagnostics. The functions [Link]() and
[Link] provide additional graphical tools to examine the relationships between
the predictor variable. The function [Link]() applies the models to GIS image files of
predictors to create detailed prediction surfaces. This function will handle large predictor files for
map making, by reading in the GIS data in sections, thus keeping memory usage reasonable. The
raster package is used to read and write to the GIS image files.

2
2.1 Interactive Model Creation
The ModelMap package can be run in a traditional R command line mode, where all arguments
are specified in the function call. However, in a Windows environment, ModelMap can also be used
in an interactive, pushbutton mode. If the functions [Link](), [Link](), and
[Link]() are called without argument lists, pop up windows ask questions about the type
of model, the file locations of the data, response variable, predictors, etc . . .
To provide a record of the options chosen for a particular model and map, a text file is generated
each time these functions are called, containing a list of the selected arguments.
This paper concentrates on the traditional command line function calls, but does contain some
tips on using the GUI prompts.

2.2 File Names

File names in the argument lists for the functions can be provided either as the full path, or as
the base name, with the path specified by the folder argument. However, file names in the Raster
Look Up Table (the rastLUTfn, described in section 2.8) must include the full path.

2.3 Training Data

Training and test data can be supplied in two forms. The argument [Link] can be either
an R data frame containing the training data, or the file name (full path or base name) of the
comma separated values (CSV) training data file. If a filename is given, the file must be a comma-
delimited text file with column headings. The data frame or CSV file should include columns for
both response and predictor variables.
In a Windows environment, if [Link] = NULL (the default), a GUI interface prompts the
user to browse to the training data file.
Note: If [Link] = "binary", any response with a value greater than 0 is treated as a
presence. If there is a cutoff value where anything below that value is called trace, and treated as
an absence, the response variable must be transformed before calling the functions.

2.4 Independent Test Set for Model Validation

The argument [Link] is the file name (full path or base name) of the independent data set
for testing (validating) the model’s predictions, or alternatively, the R data frame containing the
test data. The column headings must be the same as those in the training data (qdatatrainfn).
If no test set is desired (for example, cross-validation will be performed, or RF models with out-
of-bag estimation), set [Link] = FALSE.
In a Windows environment, if [Link] = NULL (default), a prompt will ask the a test set is
available, and if so asks the user to browse to the test data file. If no test set is available, the a
prompt asks if a proportion of the data should be set aside as an independent test set. If this is
desired, the user will be prompted to specify the proportion to set aside as test data, and two new
data files will be generated in the output folder. The new file names will be the original data file
name with "_train" and "_test" pasted on the end.

2.5 Missing predictor values

There are three circumstances that can lead to having missing predictor values. First, there are
true missing values for predictors within the test set or study area. Second, there are categorical
predictors with categories that are present in the test or mapping data but not in the training

3
data. And finally, portions of the mapping rectangle lie outside of the study area. Each of the
three cases is handled slightly differently by ModelMap.
In the first instance, true mising data values in the test set or within the study area for production
mapping could be caused by data collection errors. These are data points or pixels for which you
may still need be interested in a prediction based on the other remaining predictors. These missing
values should be coded as NA in the training or test data. In Imagine image files, pixels of the
specified NODATA value will be read into R as NA. The argument [Link] will determine how these
NA pixels will be treated. For model diagnostics, there are 2 options: (1) [Link] = "[Link]"
(the default) where any data point or pixel with any NA predictors is omitted from the model
building process and the diagnostic predictions, or returned as -9999 in the map predictions; (2)
[Link] = "[Link]" where before making predictions, a missing categorical predictor is
replaced with the most common category for that predictor, and a missing continuous predictor
is replaced with the median for that predictor. Currently, for map making only one option is
available: [Link] = "[Link]".
The second type of missing value occurs when using categorical predictors. There may be cases
where a category is found in the validation test set or in the map region that was not present
in the training data. This is a particularly common occurrence when using cross-validation on a
small dataset. Again, the argument [Link] will determine how these data points or pixels are
treated. If [Link] = "[Link]", no prediction will be made for these locations. For model
diagnostics, with [Link] = "[Link]" the most common category will be substituted
for the unknown category. Again, for map making [Link] = "[Link]" is the only available
option. In either instance, a warning will be generated with a list of the categories that were
missing from the training data. After examining these categories, you may decide that rather
than omitting these locations or substituting the most common category, a better option would
be to collapse similar categories into larger groupings. In this case you would need to pre-process
your data and run the models and predictions again.
The final type of missing predictor occurs when creating maps of non-rectangular study regions.
There may be large portions of the rectangle where you have no predictors, and are uninterested
in making predictions. The suggested value for the pixels outside the study area is -9999. These
pixels will be ignored, thus saving computing time, and will be exported as NA.
Note: in Imagine image files, if the specified NODATA is set as -9999, any -9999 pixels will be read
into R as NA.

2.6 The Model Object

The models built by the ModelMap package are stochastic models. If a seed is not specified
(with argument seed) each function call will result in a slightly different model. The function
[Link]() returns the model object. To keep this particular model for use in later R sessions,
assign the function output to an R object, then use the functions save() and load().
Random Forest is implemented through the randomForest package within R. Random Forest has
relatively few user set parameters and is not very sensitive to tuning of these parameters. The
number of predictors used to select the splits (the mtry argument) is the primary user specified
parameter that can affect model performance, and the default for ModelMap is to automatically
optimize this parameter with the tuneRF() function from the randomForest package. In most
circumstance, Random Forest is less likely to over fit data. For an in depth discussion of the
possible penalties of increasing the number of trees (the ntree argument) see Lin and Jeon (2002).
The randomForest package provides two measures to evaluate variable importance. The first is
the percent increase in Mean Standard Error (MSE) as each variable is randomly permuted. The
second is the increase in node purity from all the splits in the forest based on a particular variable,
as measured by the gini criterion (Breiman, 2001). These importance measures should be used
with caution when predictors vary in scale or number of categories (Strobl et al., 2007).
Stochastic gradient boosting is currently disabled in ModelMap.

4
2.7 Spatial Raster Layers
The ModelMap uses the raster package to read spatial rasters into R. The data for predictive
mapping in ModelMap should be in the form of pixel-based raster layers representing the predictors
in the model. The layers must be a file type recognizeable by the raster package, for example
ERDAS Imagine image (single or multi-band) raster data formats, having continuous or categorical
data values. For effective model development and accuracy, if there is more than one raster layer,
the layers must have the same extent, projection, and pixel size.
To speed up processing of predictor layers, [Link]() builds a single raster brick object
containing one layer for each predictor in the model. By default, this brick is stored as a temp
file and deleted once the map is complete. If the argument [Link] = TRUE then
the brick will be saved in a native raster package format file, with the file name constructed by
appending ’_brick.grd’ to the OUTPUTfn.
The function [Link]() by default outputs an ERDAS Imagine image file of map informa-
tion suitable to be imported into a GIS. (Note: The file extension of the OUTPPUTfn argument can
be used to specify other file types, see the help file for the writeFormats() function in the raster
package for a list of possible file types and extensions.) Maps can then be imported back into R
and view graphically using the raster package.
The supplementary materials in Elith et al. (2008) also contain R code to predict to grids imported
from a GIS program, including large grids that need to be imported in pieces. However this code
requires pre-processing of the raster data in the GIS software to produce ASCII grids for each layer
of data before they can be imported into R. ModelMap simplifies and automates this process, by
reading Imagine image files directly, (including multi band images). ModelMap also will verify
that the extent of all rasters is identical and will produce informative error messages if this is not
true. ModelMap also simplifies working with masked values and missing predictors.

2.8 Raster Look Up Table

The Raster Look Up Table (rastLUTfn) provides the link between the spatial rasters for map
production and the column names of the Training and Test datasets. The Raster Look Up Table
can be given as an R data frame specified by the argument rastLUTfn or read in from a CSV file
specified by rastLUTfn.
The rastLUTfn must include 3 columns: (1) the full path and base names of the raster file or files;
(2) the column headers from the Training and Test datasets for each predictor; (3) the layer (band)
number for each predictor. The names (column 2) must match not only the column headers in
Training and Test data sets ([Link] and [Link]), but also the predictor names in
the arguments predList and predFactor, and the predictor names in [Link].
In a windows environment, the function [Link]() may be used to help build the look-up-
table with the aid of GUI prompts.

3 Examples
These examples demonstrate some of the capabilities of the ModelMap package by building three
Random Forest models: continuous response, binary response and categorical response. The
continuous response variables are percent cover for two species of interest: Pinyon and Sage. The
binary response variables are Presence/Absence of these same [Link] categorical response is
vegetation category.
Next, model validation diagnostics are performed with three techniques: an independent test set,
Out Of Bag estimation, and cross-validation. Note: in an actual model comparison study, rather
than a package demonstration, the models would be compared with the same validation technique,
rather than mixing techniques.

5
Name Type Description
ELEV250 Continuous 90m NED elevation (ft)
resampled to 250m, average of 49 points
NLCD01 250 Categorical National Land Cover Dataset 2001
resampled to 250m - min. value of 49 points
EVI2005097 Continuous MODIS Enhanced vegetation index
NDV2005097 Continuous MODIS Normalized difference vegetation index
NIR2005097 Continuous MODIS Band 2 (Near Infrared)
RED2005097 Continuous MODIS Band 1 (Red)

Table 1: Predictor variables

Finally, spatial maps are produced by applying these models to remote sensing raster layers.

3.1 Example dataset

The dataset is from a pilot study in Nevada launched in 2004 involving acquisition and photo-
interpretation of large-scale aerial photography, the Nevada Photo-Based Inventory Pilot (NPIP)
(Frescino et al., 2009). The data files for these examples are included in the ModelMap package
instalation in the R library dirrectory. The datasets are under the ’external’ then under ’vignette-
examples’.
The predictor data set consists of 6 predictor variables: 5 continuous variables, and 1 categor-
ical variable (Table 1). The predictor layers are 250-meter resolution, pixel-based raster layers
including Moderate Resolution Imaging Spectro-radiometer (MODIS) satellite imagery (Justice
et al., 2002), a Landsat Thematic Mapper-based, thematic layer of predicted land cover, National
Land Cover Dataset (NLCD) (Homer et al., 2004), and a topographic layer of elevation from the
National Elevation Dataset (Gesch et al., 2002).
The continuous response variables are percent cover of Pinyon and Sage. The binary response
variables are presence of Pinyon and Sage. The categorical response variable is the vegetation
category: TREE, SHRUB, OTHERVEG, and NONVEG.
The MODIS data included 250-meter, 16-day, cloud-free, composites of MODIS imagery for April
6, 2005: visible-red (RED) and near-infrared (NIR) bands and 2 vegetation indices, normalized
difference vegetation index (NDVI) and enhanced vegetation index (EVI) (Huete et al., 2002). The
land cover and topographic layers were 30-meter products re-sampled to 250 meter using majority
and mean summaries, respectively.
The rectangular subset of Nevada chosen for these maps contains a small mountain range sur-
rounded by plains, and was deliberately selected to lie along the edge of the study region to
illustrate how ModelMap handles unsampled regions of a rectangle (Figure 13).

3.2 Example 1 - Random Forest - Continuous Response

Example 1 builds Random Forest models for two continuous response variables: Percent Cover for
Pinyon and Percent Cover for Sage. An independent test set is used for model validation.

3.2.1 Set up

After installing the ModelMap package, find the sample datasets from the R istallation and copy
them to your working directory. The data consists of five files and is located in the vignette
directory of ModelMap, for example, in C:\R\R-2.15.0\library\ModelMap\vignettes.
There are 5 files:

6
[Link]
VModelMapData [Link]
VModelMapData dem ELEVM [Link]
VModelMapData modis [Link]
VModelMapData nlcd NLCD01 [Link]
Load the ModelMap package.

R> library("ModelMap")

Next define some of the arguments for the models.

Specify model type. The choices are "RF" for Random Forest models, "QRF" for quantile regression
forest models, and "CF" for conditional inference forest models. See the ”Pick your flavor of
Random Forest” vignette for further information on QRF and CF models.

R> [Link] <- "RF"

Define training and test data file names. Note that the arguments [Link] and [Link]
will accept either character strings giving the file names of CSV files of data, or the data itself in
the form of a data frame.

R> qdatafn <- "[Link]"

R> [Link] <- "VModelMapData_TRAIN.csv"
R> [Link] <- "VModelMapData_TEST.csv"

Define the output folder.

R> folder <- getwd()

Split the data into training and test sets. In example 1, an independent test set is used for
model validation diagnostics. The function [Link]() randomly divides the original data into
training and test sets. This function writes the training and test sets to the folder specified by
folder, under the file names specified by [Link] and [Link]. If the arguments
[Link] and [Link] are not included, filenames will be generated by appending
"_train" and "_test" to qdatafn.

R> [Link]( [Link]=0.2,

qdatafn=qdatafn,
seed=42,
folder=folder,
[Link]=[Link],
[Link]=[Link])

Define file names to store model output. This filename will be used for saving the model itself. In
addition, since we are not defining other output filenames, the names for other output files will be
generated based on MODELfn.

R> MODELfn.a <- "VModelMapEx1a"

R> MODELfn.b <- "VModelMapEx1b"

Define the predictors and define which predictors are categorical. Example 1 uses five continuous
predictors: the four predictor layers from the MODIS imagery plus the topographic elevation layer.
As none of the chosen predictors are categorical set predFactor to FALSE.

7
R> predList <- c( "ELEV250",
"EVI2005097",
"NDV2005097",
"NIR2005097",
"RED2005097")
R> predFactor <- FALSE

Define the response variable, and whether it is continuous, binary or categorical.

R> [Link].a <- "PINYON"

R> [Link].b <- "SAGE"
R> [Link] <- "continuous"

Define the seeds for each model.

R> seed.a <- 38

R> seed.b <- 39

Define the column that contains unique identifiers for each data point. These identifiers will be
used to label the output file of observed and predicted values when running model validation.

R> [Link] <- "ID"

3.2.2 Model creation

Now create the models. The [Link]() function returns the model object itself. The func-
tion also saves a text file listing the values of the arguments from the function call. This file is
particularly useful when using the GUI prompts, as otherwise there would be no record of the
options used for each model.

8
3.2.3 Model Diagnostics

Next make model predictions on an independent test set and run the diagnostics on these predic-
tions. Model predictions on an independent test set are not stochastic, it is not necessary to set
the seed.
The [Link]() function returns a data frame of observed and predicted values. This
data frame is also saved as a CSV file. This function also runs model diagnostics, and creates
graphs and tables of the results. The graphics are saved as files of the file type specified by
[Link].
For a continuous response model, the model validation diagnostics graphs are the variable impor-
tance plot (Figure 1 and Figure 2), and a scatter plot of observed verses predicted values, labeled
with the Pearson’s and Spearman’s correlation coefficients and the slope and intercept of the linear
regression line (Figure 3 and Figure 4).
For Random forest models, the model diagnostic graphs also include the out of bag model error
as a function of the number of trees(Figure 5 and Figure 6)
In example 1, the diagnostic plots are saved as PDF files.
These diagnostics show that while the most important predictor variables are similar for both
models, the correlation coefficients are considerably higher for the Pinyon percent cover model as
compared to the Sage model.

3.2.4 Comparing Variable Importance

The [Link]() function uses a back-to-back barplot to compare variable impor-

tance between two models built from the same predictor variables (Figure 7). Variable Importance
is calculated in various ways depending on the model type, the response type and the importance
type. Importance measures are summarized in (Table 2).

R> [Link]( [Link].1=[Link].ex1a,

[Link].2=[Link].ex1b,
[Link].1="Pinyon",
[Link].2="Sage",
[Link]="predList",
predList=predList,

9
Relative Influence
VModelMapEx1a_pred

ELEV250 ● ELEV250 ●

NDV2005097 ● NDV2005097 ●

NIR2005097 ● EVI2005097 ●

EVI2005097 ● RED2005097 ●

RED2005097 ● NIR2005097 ●

20 60 100 0 40000
%IncMSE IncNodePurity

Figure 1: Example 1 - Variable importance graph for Pinyon percent cover (RF model).

Model Response Importance Type Measured By

RF continuous 1 - permutation Percent Increase in MSE
RF binary 1 - permutation Mean Decrease in Accuracy
RF categorical 1 - permutation Mean Decrease in Accuracy
RF continuous 2 - node impurity Residual Sum of Squares
RF binary 2 - node impurity Mean Decrease in Gini
RF categorical 2 - node impurity Mean Decrease in Gini

Table 2: Importance Measures

10
Relative Influence
VModelMapEx1b_pred

ELEV250 ● ELEV250 ●

NDV2005097 ● EVI2005097 ●

EVI2005097 ● NDV2005097 ●

RED2005097 ● NIR2005097 ●

NIR2005097 ● RED2005097 ●

20 30 40 50 0 20000 40000
%IncMSE IncNodePurity

Figure 2: Example 1 - Variable importance graph for Sage percent cover (RF model).

11
VModelMapEx1a_pred

●
50

●
●
●
40

● ●
● ●
● ● ● ●
●
● ● ● ●
●
30
observed

●●●●●● ● ● ● ● ●
●
● ●● ● ●
●
● ●● ●
● ●
●●● ●
20

●● ● ●●
● ● ●●
●● ●● ●
● ● ●
● ● ● ●
●
● ●●● ●
●
●● ● ●● ● ●
10

● ● ●● ●
● RMSD:9.28
●● ●
●● ● ● ● pearson's cor: 0.69
● ● ● ● ● ●
spearman's cor: 0.76
●● ●● ●● ●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●●
●●
●●
●
●● ●●
●●●●
●●●●●● ●●●●● ● ● ● obs
● = 0.89(pred) + 0.06
0

0 10 20 30 40 50

predicted

Figure 3: Example 1 - Observed verses predicted values for Pinyon percent cover (RF model).

12
VModelMapEx1b_pred

●
70

●
●
60

●●

●
●
50

●
●
●
● ●
●● ● ●
observed

●
●
●
●● ●
●
30

●
● ● ●
● ●● ●
● ● ●
●
20

●● ●● ●
●
● ●
●
●
● ●
● ●● ● ●
● ● RMSD:13.19
●●
● ● ● ●
10

●●●●●● ●● ● ●
●●● ● ●● pearson's cor: 0.43
●
●● ●●● ●● ●
●●●●●
● ●●
●●● ●● ● spearman's cor: 0.41
●
●●
●●●●
●●●
●●● ● ●● ● ●●
●●
●
●
●●
●
●●
●
●●
●●
●
●
●●
●
●
●●
●
●●●
●
●●
●●
●●
●
●●
●
●●●●
●
●●●
●
● ● ● obs = 0.94(pred) + 1.15
0

0 10 20 30 40 50 60 70

predicted

Figure 4: Example 1 - Observed verses predicted values for Sage percent cover (RF model).

13
VModelMapEx1a_pred
OOB − 953 plots
160
140
MSE

120
100

0 100 200 300 400 500

ntree

Figure 5: Example 1 - Out of Bag error as a function of number of trees for Pinyon (RF model).

14
VModelMapEx1b_pred
OOB − 953 plots
240
220
MSE

200
180
160

0 100 200 300 400 500

ntree

Figure 6: Example 1 - Out of Bag error as a function of number of trees for Sage (RF model).

15
Variable Importance

ELEV250

EVI2005097

NDV2005097

NIR2005097

RED2005097

%IncMSE %IncMSE
Pinyon Sage

Figure 7: Example 1 - Variable Importances for Pinyon verses Sage percent cover models.

[Link]="sum",
main="Variable Importance",
[Link]="pdf",
PLOTfn="VModelMapEx1CompareImportance",
folder=folder)
R>

The [Link]() function can also be used to compare the two types of variable
importance (Figure 8).

R> opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))

R> [Link]( [Link].1=[Link].ex1a,
[Link].2=[Link].ex1a,
[Link].1="",
[Link].2="",
[Link].1=1,
[Link].2=2,
[Link]="predList",
predList=predList,
[Link]="sum",
main="Pinyon",

16
[Link]="none",
cex=0.9)
R> [Link]( [Link].1=[Link].ex1b,
[Link].2=[Link].ex1b,
[Link].1="",
[Link].2="",
[Link].1=1,
[Link].2=2,
[Link]="predList",
predList=predList,
[Link]="sum",
main="Sage",
[Link]="none",
cex=0.9)
R> mtext("Comparison of Importance Types",side=3,line=0,cex=1.8,outer=TRUE)
R> par(opar)

3.2.5 Interaction Plots

The [Link]() function provides a diagnostic plot useful in visualizing two-way

interactions between predictor variables. Two of the predictor variables from the model are used to
produce a grid of possible combinations over the range of both variables. The remaining predictor
variables are fixed at either their means (for continuous predictors) or their most common value
(for categorical predictors). Model predictions are generated over this grid and plotted as the z
axis. The [Link]() function was developed from the [Link] function
from the tutorial provided as appendix S3 in Elith et al. (2008).
The [Link]() function provides two graphical options: an image plot, and a
3-D perspective plot. These options are selected by setting [Link] = "image" or [Link]
= "persp"
The x and y arguments are used to specify the predictor variables for the X and Y axis. The
predictors can be specified by name, or by number, with the numbers referring to the order the
variables appear in predList.

17
Comparison of Importance Types
Pinyon

ELEV250

EVI2005097

NDV2005097

NIR2005097

RED2005097

%IncMSE IncNodePurity

Sage

ELEV250

EVI2005097

NDV2005097

NIR2005097

RED2005097

%IncMSE IncNodePurity

Figure 8: Example 1 - Variable Importances Types for continuous response models - the relative
importance of predictor variables as measured by the mean decrease in accuracy from randomly
permuting each predictor as compared to the decrease in node impurities from splitting on the
variable.

18
PINYON

8000

8
6000

7
RED2005097

4000

5
2000

2000 3000 4000 5000 6000 7000

NIR2005097

Figure 9: Example 1 - Interaction plot for Pinyon percent cover (RF model), showing interactions
between two of the satellite based predictors (NIR2005097 and RED2005097). Image plot, with
darker green indicating higher percent cover. Here we can see that predicted Pinyon cover is
highest at low values of either NIR or RED. However low values of both predictors does not
further raise the predicted cover.

main=[Link].a,
[Link]="image",
[Link]="pdf",
MODELfn=MODELfn.a,
folder=folder)
R> [Link]( [Link].ex1b,
x=1,
y=3,
main=[Link].b,
[Link]="image",
[Link]="pdf",
MODELfn=MODELfn.b,
folder=folder)

19
SAGE
8000

20
6000

18
RED2005097

16
4000

12
2000

2000 3000 4000 5000 6000 7000

NIR2005097

Figure 10: Example 1 - Interaction plot for Sage percent cover (RF model), showing interactions
between two of the satellite based predictors (NIR2005097 and RED2005097). Image plot, with
darker green indicating higher percent cover. Here we can see that predicted Sage cover is lowest
when high values of NIR are combined with RED lying between 2000 and 5000. When NIR is
lower than 2900, RED still has an effect on the predicted cover, but the effect is not as strong.

20
PINYON
5000

20
4000
3000

15
NDV2005097

2000

10
1000

5
0

500 1000 1500 2000 2500 3000

ELEV250

Figure 11: Example 1 - Interaction plot for Pinyon percent cover (RF model), showing interactions
between elevation and a satellite based predictor (ELEV250 and NDV2005097). Image plot, with
darker green indicating higher percent cover. Here we can see that predicted Pinyon cover is
highest at elevation greater than 2000m. In addition, high values of NDV slightly increase the
predicted cover, but there seems to be little interaction between the two predictors.

21
SAGE
5000

25
4000

20
3000
NDV2005097

15
2000

10
1000

5
0

500 1000 1500 2000 2500 3000

ELEV250

Figure 12: Example 1 - Interaction plot for Sage percent cover (RF model), showing interactions
between elevation and a satellite based predictor (ELEV250 and NDV2005097). Image plot, with
darker green indicating higher percent cover. Here we do see an interaction between the two
predictors. At low elevations, predicted Sage cover is low throughout the range of NDV, and
particularly low at mid-values. At mid elevations, predicted Sage cover is high throughout the
range of NDV. At high elevations NDV has a strong influence on predicted Sage cover with high
cover tied to low to mid values of NDV.

22
3.2.6 Map production

Before building maps of the responses, examine the predictor variable for elevation (Figure 13):

R> elevfn <- paste(folder,"/VModelMapData_dem_ELEVM_250.img",sep="")

R> mapgrid <- raster(elevfn)

R> opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))

R> [Link]<-[Link](101)
R> zlim <- c(1500,maxValue(mapgrid))
R> [Link]<-rev(pretty(zlim,n=5))
R> [Link]<-[Link][trunc(([Link]/max([Link]))*100)+1]
R> [Link]<-paste([Link],"m",sep="")
R> [Link]<-paste((7:3)*500,"m")
R> [Link]<-[Link][c(100,75,50,25,1)]
R> image( mapgrid,
col = [Link],
xlab="", ylab="",
zlim=zlim,
asp=1, bty="n", main="")
R> legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=[Link],
fill=[Link],
bty="n",
cex=1.2)
R> mtext("Elevation of Study Region",side=3,line=1,cex=1.5)
R> par(opar)

Run the function [Link]() to map the response variable over the study area.
The [Link]() function can extract information about the model from the [Link], so
it is not necessary to re-supply the arguments that define the model, such as the type of model, the
predictors, etc . . . (Note: If model was created outside of ModelMap, it may be necessary to supply
the [Link] argument) Also, unlike model creation, map production is not stochastic, so
it is not necessary to set the seed.
The [Link]() uses a look up table to associate the predictor variables with the rasters.
The function argument rastLUTfn will accept either a file name of the CSV file containing the
table, or the data frame itself.
Although in typical user applications the raster look up table must include the full path for
predictor rasters, the table provided for the examples will be incomplete when initially downloaded,
as the working directory of the user is unknown and will be different on every computer. This
needs to be corrected by pasting the full paths to the user’s working directory to the first column,
using the value from folder defined above.

R> rastLUTfn <- "VModelMapData_LUT.csv"

R> rastLUTfn <- [Link]( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
R> rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")

To produce a map from a raster larger than the memory limits of R, predictions are made one row
at a time.
Since this is a Random Forest model of a continuous response, the prediction at each pixel is the
mean of all the trees. Therefore these individual tree predictions can also be used to map measures

23
Elevation of Study Region

3500 m
1960000

3000 m
2500 m
2000 m
1500 m
1955000
1950000
1945000

−1544000 −1538000 −1532000

Figure 13: Elevation of study region. Projection: Universal Transverse Mercator (UTM) Zone 11,
Datum: NAD83

24
of uncertainty such as standard deviation and coefficient of variation for each pixel. To do so, set
[Link] = "TRUE". To calculate these pixel uncertainty measures, [Link]() must keep all the
individual trees in memory, so [Link] = "TRUE" is much more memory intensive.

The function [Link]() creates an Imagine image file of map information suitable to be
imported into a GIS. As this sample dataset is relatively small, we can also import it into R for
display.
We need to define a color ramp. For this response variable, zero values will display as white,
shading to dark green for high values.

R> l <- seq(100,0,[Link]=101)

R> c <- seq(0,100,[Link]=101)
R> [Link] <- hcl(h = 120, c = c, l = l)

Next, we import the data and create the map (Figure 14). From the map, we can see that Pinyon
percent cover is higher in the mountains, while Sage percent cover is higher in the foothills at the
edges of the mountains.
Note that the sample map data was taken from the South Eastern edge of our study region, to
illustrate how ModelMap deals with portions of the rectangle that fall outside of the study region.
The empty wedge at lower right in the maps is the portion outside the study area. ModelMap
uses -9999 for unsampled data. When viewing maps in a GIS, a mask file can be used to hide
unsampled regions, or other commands can be used to set the color for -9999 values.
Since we know that percent cover can not be negative, we will set zlim to range from zero to the
maximum value found in our map.

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> mapgrid.a <- raster(paste(MODELfn.a,"_map.img",sep=""))
R> mapgrid.b <- raster(paste(MODELfn.b,"_map.img",sep=""))
R> zlim <- c(0,max(maxValue(mapgrid.a),maxValue(mapgrid.b)))
R> [Link]<-rev(pretty(zlim,n=5))
R> [Link]<-[Link][trunc(([Link]/max([Link]))*100)+1]
R> [Link]<-paste([Link],"%",sep="")
R> image( mapgrid.a,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,

25
Percent Cover
PINYON SAGE

40%
30%
20%
10%
0%

Figure 14: Example 1 - Maps of percent cover for Pinyon and Sage (RF models).

asp=1,bty="n",main="")
R> mtext([Link].a,side=3,line=1,cex=1.2)
R> image( mapgrid.b,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
R> mtext([Link].b,side=3,line=1,cex=1.2)
R> legend( x=xmax(mapgrid.b),y=ymax(mapgrid.b),
legend=[Link],
fill=[Link],
bty="n",
cex=1.2)
R> mtext("Percent Cover",side=3,line=1,cex=1.5,outer=T)
R> par(opar)

Next, we will define color ramps for the standard deviation and the coefficient of variation, and
map these uncertainty measures. Often, as the mean increases, so does the standard deviation
(Zar, 1996), therefore, a map of the standard deviation of the pixels (Figure 15) will look to the
naked eye much like the map of the mean. However, mapping the coefficient of variation (dividing
the standard deviation of each pixel by the mean of the pixel), can provide a better visualization
of spatial regions of higher uncertainty (Figure 16). In this case, for Pinyon the coefficient of
variation is interesting as it is higher in the plains on the upper left portion of the map, where
percent cover of Pinyon is lower.

R> [Link] <- hcl(h = 15, c = c, l = l)

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> mapgrid.a <- raster(paste(MODELfn.a,"_map_stdev.img",sep=""))

26
Standard Deviation of Percent Cover
PINYON SAGE

25%
20%
15%
10%
5%
0%

Figure 15: Example 1 - Map of standard deviation of Random Forest trees at each pixel for Pinyon
and Sage (RF models).

R> mapgrid.b <- raster(paste(MODELfn.b,"_map_stdev.img",sep=""))

R> zlim <- c(0,max(maxValue(mapgrid.a),maxValue(mapgrid.b)))
R> [Link]<-rev(pretty(zlim,n=5))
R> [Link]<-[Link][trunc(([Link]/max([Link]))*100)+1]
R> [Link]<-paste([Link],"%",sep="")
R> image( mapgrid.a,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
R> mtext([Link].a,side=3,line=1,cex=1.2)
R> image( mapgrid.b,
col=[Link],xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
R> mtext([Link].b,side=3,line=1,cex=1.2)
R> legend( x=xmax(mapgrid.b),y=ymax(mapgrid.b),
legend=[Link],
fill=[Link],
bty="n",
cex=1.2)
R> mtext("Standard Deviation of Percent Cover",side=3,line=1,cex=1.5,outer=T)
R> par(opar)

R> [Link] <- hcl(h = 70, c = c, l = l)

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

27
Coefficient of Variation of Percent Cover
PINYON SAGE

25
20
15
10
5
0

Figure 16: Example 1 - Map of coefficient of variation of Random Forest trees at each pixel for
Pinyon and Sage (RF models).

R> mapgrid.a <-raster(paste(MODELfn.a,"_map_coefv.img",sep=""),[Link]=TRUE)

R> mapgrid.b <- raster(paste(MODELfn.b,"_map_coefv.img",sep=""),[Link]=TRUE)
R> zlim <- c(0,max(maxValue(mapgrid.a),maxValue(mapgrid.b)))
R> [Link]<-rev(pretty(zlim,n=5))
R> [Link]<-[Link][trunc(([Link]/max([Link]))*100)+1]
R> image( mapgrid.a,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",zlim=zlim,
asp=1,bty="n",main="")
R> mtext([Link].a,side=3,line=1,cex=1.2)
R> image( mapgrid.b,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
R> mtext([Link].b,side=3,line=1,cex=1.2)
R> legend( x=xmax(mapgrid.b),y=ymax(mapgrid.b),
legend=[Link],
fill=[Link],
bty="n",
cex=1.2)
R> mtext("Coefficient of Variation of Percent Cover",side=3,line=1,cex=1.5,outer=T)
R> par(opar)

28
3.3 Example 2 - Random Forest - Binary Response
Example 2 builds a binary response model for presence of Pinyon and Sage. A catagorical predictor
is added to the model. Out-of-bag estimates are used for model validation.

3.3.1 Set up

Define model type.

R> [Link] <- "RF"

Define data.

R> qdatafn <- "[Link]"

Define folder.

R> folder <- getwd()

Define model filenames.

R> MODELfn.a <- "VModelMapEx2a"

R> MODELfn.b <- "VModelMapEx2b"

Define the predictors. These are the five continuous predictors from the first example, plus one
categorical predictor layer, the thematic layer of predicted land cover classes from the National
Land Cover Dataset. The argument predFactor is used to specify the categorical predictor.

R> predList <- c( "ELEV250",

"NLCD01_250",
"EVI2005097",
"NDV2005097",
"NIR2005097",
"RED2005097")
R> predFactor <- c("NLCD01_250")

Define the data column to use as the response, and if it is continuous, binary or categorical.
Since [Link] = "binary" this variable will be automatically translated so that zeros are
treated as Absent and any value greater than zero is treated as Present.

R> [Link].a <- "PINYON"

R> [Link].b <- "SAGE"
R> [Link] <- "binary"

Define the seeds for each model.

R> seed.a <- 40

R> seed.b <- 41

Define the column that contains unique identifiers.

R> [Link] <- "ID"

29
Define raster look up table.

R> rastLUTfn <- "VModelMapData_LUT.csv"

R> rastLUTfn <- [Link]( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
R> rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")

3.3.2 Model creation

Create the model. Because Out-Of-Bag predictions will be used for model diagnostics, the full
dataset can be used as training data. To do this, set [Link] <- qdatafn, [Link]
<- FALSE and [Link] = FALSE.

3.3.3 Model Diagnostics

Make Out-Of-Bag model predictions on the training data and run the diagnostics on these pre-
dictions. This time, save JPEG, PDF, and PS versions of the diagnostic plots.
Out of Bag model predictions for a Random Forest model are not stochastic, so it is not necessary
to set the seed.
Since this is a binary response model model diagnostics include ROC plots and other thresh-
old selection plots generated by PresenceAbsence (Freeman, 2007; Freeman and Moisen, 2008a)
(Figure 17 and Figure 18) in addition to the variable importance graph (Figure 19 and Figure 20).
For binary response models, there are also CSV files of presence-absence thresholds optimized
by 12 possible criteria, along with their associated error statistics. For more details on these
12 optimization criteria see Freeman and Moisen (2008a). Some of these criteria are dependent
on user selected parameters. In this example, two of these parameters are specified: required
sensitivity ([Link]) and required specificity ([Link]). Other user defined parameters, such
as False Positive Cost (FPC) and False Negative Cost (FNC) are left at the default values. When
default values are used for these parameters, [Link]() will give a warning. In this
case:

30
1: In [Link](PRED, [Link] = [Link](), ... :
costs assumed to be equal

The variable importance graphs show NLCD was a very important predictor for Pinyon presence,
but not an important variable when predicting Sage presence.

Take a closer look at the text file of thresholds optimized by multiple criteria. These thresholds
are used later to display the mapped predictions, so read this file into R now.

R> [Link].a <- [Link]( paste(MODELfn.a,"_pred_optthresholds.csv",sep=""),

header=TRUE,
sep=",",
stringsAsFactors=FALSE)
R> [Link].a[,-1]<-signif([Link].a[,-1],2)
R> [Link].b <- [Link]( paste(MODELfn.b,"_pred_optthresholds.csv",sep=""),
header=TRUE,
sep=",",
stringsAsFactors=FALSE)
R> [Link].b[,-1]<-signif([Link].b[,-1],2)
R> [Link].a <- [Link]( paste(MODELfn.a,"_pred_prevalence.csv",sep=""),
header=TRUE,
sep=",",
stringsAsFactors=FALSE)
R> [Link].a[,-1]<-signif([Link].a[,-1],2)
R> [Link].b <- [Link]( paste(MODELfn.b,"_pred_prevalence.csv",sep=""),
header=TRUE,
sep=",",
stringsAsFactors=FALSE)
R> [Link].b[,-1]<-signif([Link].b[,-1],2)

Optimized thresholds for Pinyon:

R> [Link].a

[Link] threshold PCC sensitivity specificity

1 Default 0.50 0.92 0.92 0.92

31
2 Sens=Spec 0.52 0.92 0.92 0.92
3 MaxSens+Spec 0.62 0.92 0.91 0.93
4 MaxKappa 0.62 0.92 0.91 0.93
5 MaxPCC 0.64 0.92 0.90 0.94
6 PredPrev=Obs 0.57 0.92 0.91 0.92
7 ObsPrev 0.46 0.92 0.92 0.91
8 MeanProb 0.47 0.92 0.92 0.91
9 MinROCdist 0.62 0.92 0.91 0.93
10 ReqSens 0.76 0.90 0.85 0.95
11 ReqSpec 0.22 0.90 0.96 0.85
12 Cost 0.64 0.92 0.90 0.94
Kappa
1 0.83
2 0.83
3 0.84
4 0.84
5 0.84
6 0.84
7 0.83
8 0.83
9 0.84
10 0.81
11 0.81
12 0.84
And for Sage:
R> [Link].b
[Link] threshold PCC sensitivity specificity
1 Default 0.50 0.66 0.81 0.48
2 Sens=Spec 0.61 0.67 0.67 0.66
3 MaxSens+Spec 0.61 0.67 0.67 0.66
4 MaxKappa 0.61 0.67 0.67 0.66
5 MaxPCC 0.60 0.67 0.69 0.65
6 PredPrev=Obs 0.59 0.66 0.70 0.61
7 ObsPrev 0.56 0.66 0.74 0.57
8 MeanProb 0.60 0.67 0.69 0.64
9 MinROCdist 0.61 0.67 0.67 0.66
10 ReqSens 0.44 0.66 0.86 0.40
11 ReqSpec 0.79 0.56 0.33 0.86
12 Cost 0.60 0.67 0.69 0.65
Kappa
1 0.29
2 0.33
3 0.33
4 0.33
5 0.33
6 0.32
7 0.31
8 0.33
9 0.33
10 0.27
11 0.18
12 0.33

32
Observed and predicted prevalence for Pinyon:

R> [Link].a

[Link] threshold [Link] pred

1 Default 0.50 0.46 0.47
2 Sens=Spec 0.52 0.46 0.47
3 MaxSens+Spec 0.62 0.46 0.46
4 MaxKappa 0.62 0.46 0.46
5 MaxPCC 0.64 0.46 0.45
6 PredPrev=Obs 0.57 0.46 0.46
7 ObsPrev 0.46 0.46 0.47
8 MeanProb 0.47 0.46 0.47
9 MinROCdist 0.62 0.46 0.46
10 ReqSens 0.76 0.46 0.42
11 ReqSpec 0.22 0.46 0.52
12 Cost 0.64 0.46 0.45

And for Sage:

R> [Link].b

[Link] threshold [Link] pred

1 Default 0.50 0.56 0.69
2 Sens=Spec 0.61 0.56 0.53
3 MaxSens+Spec 0.61 0.56 0.53
4 MaxKappa 0.61 0.56 0.53
5 MaxPCC 0.60 0.56 0.54
6 PredPrev=Obs 0.59 0.56 0.57
7 ObsPrev 0.56 0.56 0.60
8 MeanProb 0.60 0.56 0.54
9 MinROCdist 0.61 0.56 0.53
10 ReqSens 0.44 0.56 0.74
11 ReqSpec 0.79 0.56 0.25
12 Cost 0.60 0.56 0.54

The model quality graphs show that the model of Pinyon presence is much higher quality than
the Sage model. This is illustrated with four plots: a histogram plot, a calibration plot, a ROC
plot with it’s associated Area Under the Curve (AUC), and an error rate verses threshold plot
Pinyon has a double humped histogram plot, with most of the observed presences and absences
neatly divided into the two humps. Therefor the optimized threshold values fall between the two
humps and neatly divide the data into absences and presences. For Sage, on the other hand, the
observed presences and absences are scattered throughout the range of predicted probabilities,
and so there is no single threshold that will neatly divide the data into present and absent groups.
In this case, the different optimization criteria tend to be widely separated, each representing a
different compromise between the error statistics (Freeman and Moisen, 2008b).
Calibration plots provide a goodness-of-fit plot for presence-absence models, as described by Pearce
and Ferrier (2000), Vaughan and Ormerod (2005), and Reineking and Schröder (2006). In a
Calibration plot the predicted values are divided into bins, and the observed proportion of each
bin is plotted against the predicted value of the bin. For Pinyon, the standard errors for the bins
overlap the diagonal, and the bins do not show a bias. For Sage, however, the error bars for the
highest and lowest bins do not overlap the diagonal, and there is a bias where low probabilities
tend to be over predicted, and high probabilities tend to be under predicted.

33
VModelMapEx2a_pred

Observed vs. Predicted

600

484

1.0
present ●

observed as proportion of bin

500

absent 64

0.8
400
number of plots

●
28

0.6
57
300

0.4
●
200

0.2
469
100

0.0
0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

predicted probability Predicted Probability of Occurrence

ROC Plot Error Rate verses Threshold

1.0

1.0
Sensitivity (true positives)

0.8

0.8
Accuracy Measures
0.6

0.6
0.4

0.4

sensitivity
0.2

0.2

AUC: specificity
0.97 pred Kappa
0.0

0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

1−Specificity (false positives) Threshold

Figure 17: Example 2 - Model quality and threshold selection graphs for Pinyon presence (RF
model).

The ROC plot from a good model will rise steeply to the upper left corner then level off quickly,
resulting in an AUC near 1.0. A poor model (i.e. a model that is no better than random assign-
ment) will have a ROC plot lying along the diagonal, with an AUC near 0.5. The Area Under
the Curve (AUC) is equivalent to the chance that a randomly chosen plot with an observed value
of present will have a predicted probability higher than that of a randomly chosen plot with an
observed value of absent. The PresenceAbsence package used to create the model quality graphs
for binary response models uses the method from DeLong et al. (1988) to calculate Area Under
the Curve (AUC). For these two models, the Area Under the Curve (AUC) for Pinyon is 0.97 and
the ROC plot rises steeply, while the AUC for Sage is only 0.70, and the ROC plot is much closer
to the diagonal.
In the Error Rate verses Threshold plot sensitivity, specificity and Kappa are plotted against all
possible values of the threshold (Fielding and Bell, 1997). In the graph of Pinyon error rates,
sensitivity and specificity cross at a higher value, and also, the error statistics show good values
across a broader range of thresholds. The Kappa curve is high and flat topped, indicating that
for this model, Kappa will be high across a wide range of thresholds. For Sage, sensitivity and
specificity cross at a lower value, and the Kappa curve is so low that it is nearly hidden behind
the graph legend. For this model even the most optimal threshold selection will still result in a
relatively low Kappa value.

34
VModelMapEx2b_pred

Observed vs. Predicted

1.0
200

present
absent observed as proportion of bin 280

0.8
364
●
number of plots

150

0.6 ●
289
177
●
100

0.4

●
81
0.2
50

●
0.0
0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

predicted probability Predicted Probability of Occurrence

ROC Plot Error Rate verses Threshold

1.0

1.0
Sensitivity (true positives)

0.8

0.8
Accuracy Measures
0.6

0.6
0.4

0.4

sensitivity
0.2

0.2

AUC: specificity
0.71 pred Kappa
0.0

0.0

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

1−Specificity (false positives) Threshold

Figure 18: Example 2 - Model quality and threshold selection graphs for Sage presence (RF model).

35
Relative Influence
VModelMapEx2a_pred

NLCD01_250 ● NLCD01_250 ●

ELEV250 ● ELEV250 ●

EVI2005097 ● NIR2005097 ●

NIR2005097 ● EVI2005097 ●

NDV2005097 ● NDV2005097 ●

RED2005097 ● RED2005097 ●

25 35 45 0 50 150
MeanDecreaseAccuracy MeanDecreaseGini

Figure 19: Example 2 - Variable importance graph for Pinyon presence (RF model).

36
Relative Influence
VModelMapEx2b_pred

ELEV250 ● ELEV250 ●

RED2005097 ● NIR2005097 ●

EVI2005097 ● EVI2005097 ●

NDV2005097 ● NDV2005097 ●

NIR2005097 ● RED2005097 ●

NLCD01_250 ● NLCD01_250 ●

15 25 35 0 40 80
MeanDecreaseAccuracy MeanDecreaseGini

Figure 20: Example 2 - Variable importance graph for Sage presence (RF model).

37
Variable Importance

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

MeanDecAccuracy MeanDecAccuracy
Pinyon Sage

Figure 21: Example 2 - Variable Importances for Pinyon verses Sage presence models.

3.3.4 Comparing Variable Importance

The [Link]() function uses a back-to-back barplot to compare variable impor-

tance between two models built from the same predictor variables (Figure 21).

Because a binary response model is a two-class example of a categorical response model, we can
use categorical tools to investigate the class specific variable importances. (Figure 22) compares
the relative importance of the predictor variables inr predicting Presences to their importances in
predicting Absences.

38
R> opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))
R> [Link]( [Link].1=[Link].ex2a,
[Link].2=[Link].ex2a,
[Link].1="Absence",
[Link].2="Presence",
class.1="0",
class.2="1",
[Link]="predList",
predList=predList,
[Link]="sum",
main="Pinyon Variable Importance",
[Link]="none",
cex=0.9)
R> [Link]( [Link].1=[Link].ex2b,
[Link].2=[Link].ex2b,
[Link].1="Absence",
[Link].2="Presence",
class.1="0",
class.2="1",
[Link]="predList",
predList=predList,
[Link]="sum",
main="Sage Variable Importance",
[Link]="none",
cex=0.9)
R> mtext("Presence-Absence Variable Importance Comparison",side=3,line=0,cex=1.8,outer=TRUE)
R> par(opar)

3.3.5 Interaction Plots

Here we will look at how the [Link]() function works with a factored predictor
variable.
In image plots the levels of the factored predictor are shown as vertical or horizontal bars across
the plot region.
In 3-D perspective plots the levels are represented by ribbons across the prediction surface.

39
Presence−Absence Variable Importance Comparison
Pinyon Variable Importance

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

MeanDecAccuracy MeanDecAccuracy
Absence Presence

Sage Variable Importance

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

MeanDecAccuracy MeanDecAccuracy
Absence Presence

Figure 22: Example 2 - Relative variable importances (permutation based) for predicting Presences
verses predicting Absences for Pinyon and Sage.

40
R> [Link]( [Link].ex2a,
x="ELEV250",
y="NLCD01_250",
main=[Link].a,
[Link]="persp",
[Link]="pdf",
MODELfn=MODELfn.a,
folder=folder,
theta=300,
phi=55)
R> [Link]( [Link].ex2b,
x="ELEV250",
y="NLCD01_250",
main=[Link].b,
[Link]="persp",
[Link]="pdf",
MODELfn=MODELfn.b,
folder=folder,
theta=300,
phi=55)

3.3.6 Map production

The function [Link]() creates ascii text files of map predictions.

When working with categorical predictors, sometimes there are categories in the prediction data
(either the test set, or the map data) not found in the training data. In this case, there were three
classes for the predictor NLCD01_250 that were not present in the training data. With the default
[Link] = "[Link]" the [Link]() function generated the following warnings, and
these pixels will show up as blank pixels in the maps.

2: In [Link]([Link] = [Link], rastLUTfn = rastLUTfn, :

categorical factored predictor NLCD01_250 contains levels 41, 43, 20 not
found in training data
3: In [Link]([Link] = [Link], rastLUTfn = rastLUTfn, :
Returning -9999 for data points with levels not found in the training
data

Begin by mapping the probability surface, in other words, the probability that the species is
present at each grid point (Figure 27).

41
PINYON

1.0
80

0.8
70

0.6
NLCD01_250
(categorical)

0.4

0.2

30
0.0

500 1000 1500 2000 2500 3000

ELEV250

Figure 23: Example 2 - Interaction plot for Pinyon presence-absence (RF model), showing interac-
tions between elevation and National Land Cover Dataset classes (ELEV250 and NLCD01 250).
Image plot, with darker green indicating higher probability of presence. Here we see that in all
NLCD classes predicted Pinyon presence is strongly tied to elevation, with low presence below
2000m, and moderate presence at higher elevations.

42
SAGE

1.0
80

0.8
70

0.6
NLCD01_250
(categorical)

0.4

0.2

30
0.0

500 1000 1500 2000 2500 3000

ELEV250

Figure 24: Example 2 - Interaction plot for Sage presence-absence (RF model), showing interac-
tions between elevation and National Land Cover Dataset classes (ELEV250 and NLCD01 250).
Image plot, with darker green indicating higher probability of presence. Here we see that predicted
Sage presence is influenced by both elevation and NLCD class. In most NLCD classes predicted
presence is highest between 1400m and 2000m, with very low presence at lower elevations and
low presence at higher elevations. In contrast, in NLCD class 50 while predicted presence drops
slightly as elevation increases, it remains quite high all the way to 3000m.

43
PINYON

0.8
value

0.6
0.4 3000
fitted

0.2 2500

2000
0
25

80 1500
EV

70
EL

50 1000
NL
CD 42
(ca 01 500
teg _2
ori 50 30
ca
l)

Figure 25: Example 2 - Interaction plot for Pinyon presence-absence (RF model), showing interac-
tions between elevation and National Land Cover Dataset classes (ELEV250 and NLCD01 250).
Perspective plot, with probability or presence shown on the Z axis. In the perspective plot (as
compared to the image plot) it is easier to see that while predicted Pinyon presence is influenced
by both elevation and NLCD class, the shape of the relationship between elevation and presence is
similar in all NLCD classes, and while the overall probability is higher in some classes, the curves
relating probability to elevation are generally parallel from class to class. Therefore there appears
to be little 2-way interaction between these predictors.

44
SAGE

0.8
value

0.6
0.4 3000
fitted

0.2 2500

2000
0
25

80 1500
EV

70
EL

50 1000
NL
CD 42
(ca 01 500
teg _2
ori 50 30
ca
l)

Figure 26: Example 2 - Interaction plot for Sage presence-absence (RF model), showing interac-
tions between elevation and National Land Cover Dataset classes (ELEV250 and NLCD01 250).
Perspective plot, with probability or presence shown on the Z axis. Here we can see that while all
NLCD classes have low predicted Sage presence at low elevation, in NLCD class 50 at mid eleva-
tions predicted presence shoots higher than the other classes, and then does not drop as far as the
other classes at high elevations. Resulting in a different shape of the elevation verses probability
curve for class 50.

45
First Define a color ramp. For this map, pixels with a high probability of presence will display as
green, low probability will display as brown, and model uncertainty (probabilities near 50%) will
display as yellow. Notice that the map for Pinyon, is mostly dark green and dark brown, with
a thin dividing line of yellow. With a high quality model, most of the pixels are assigned high
or low probabilities. The map for Sage, however, is mostly yellow, with only occasional areas of
green and brown. With poor quality models, many of the pixels are inderminate, and assigned
probabilities near 50%.

R> h=c( seq(10,30,[Link]=10),

seq(31,40,[Link]=10),
seq(41,90,[Link]=60),
seq(91,100,[Link]=10),
seq(101,110,[Link]=10))
R> l =c( seq(25,40,[Link]=10),
seq(40,90,[Link]=35),
seq(90,90,[Link]=10),
seq(90,40,[Link]=35),
seq(40,10,[Link]=10))
R> [Link] <- hcl(h = h, c = 80, l = l)

Import the data and create the map. Since we know that probability of presence can range from
zero to one, we will use those values for zlim.

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> mapgrid.a <- raster(paste(MODELfn.a,"_map.img",sep=""))
R> mapgrid.b <- raster(paste(MODELfn.b,"_map.img",sep=""))
R> [Link]<-c(100,80,60,40,20,1)
R> [Link]<-[Link][[Link]]
R> [Link]<-c("100%"," 80%"," 60%"," 40%"," 20%"," 0%")
R> image( mapgrid.a,
col=[Link],
xlab="",ylab="",yaxt="n",main="",zlim=c(0,1),
asp=1,bty="n",xaxt="n")
R> mtext([Link].a,side=3,line=1,cex=1.2)
R> image( mapgrid.b,
col=[Link],
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=c(0,1),
asp=1,bty="n",main="")
R> mtext([Link].b,side=3,line=1,cex=1.2)
R> legend( x=xmax(mapgrid.b),y=ymax(mapgrid.b),
legend=[Link],
fill=[Link],
bty="n",
cex=1.2)
R> mtext("Probability of Presence",side=3,line=1,cex=1.5,outer=T)
R> par(opar)

To translate the probability surface into a Presence-Absence map it is necessary to select a cutoff
threshold. Probabilities below the selected threshold are mapped as absent while probabilities
above the threshold are mapped as present. Many criteria that can be used for threshold selection,
ranging from the traditional default of 50 percent, to thresholds optimized to maximize Kappa, to
thresholds picked to meet certain management criteria. The choice of threshold criteria can have

46
Probability of Presence
PINYON SAGE

100%
80%
60%
40%
20%
0%

Figure 27: Example 2 - Probability surface map for presence of Pinyon and Sage (RF models).

a dramatic effect on the final map. For further discussion on this topic see Freeman and Moisen
(2008b).
Here are examples of Presence-Absence maps for Pinyon and Sage produced by four different
threshold optimization criteria (Figures 28 and 29). For a high quality model, such as Pinyon, the
various threshold optimization criteria tend to result in similar thresholds, and the models tend to
be less sensitive to threshold choice, therefore the Presence Absence maps from the four criteria
are very similar. Poor quality models, such as this model for Sage, tend to have no single good
threshold, as each criteria is represents a different compromise between errors of omission and
errors of commission. It is therefore particularly important to carefully match threshold criteria
to the intended use of the map.

R> opar <- par(mfrow=c(2,2),mar=c(2.5,3,4,1),oma=c(0,0,4,6),xpd=NA)

R> mapgrid <- raster(paste(MODELfn.a,"_map.img",sep=""))
R> criteria <- c("Default","MaxKappa","ReqSens","ReqSpec")
R> [Link]<-c("Default","MaxKappa","ReqSens = 0.9","ReqSpec = 0.9")
R> for(i in 1:4){
thresh <- [Link].a$threshold[[Link].a$[Link]==criteria[i]]
presencegrid <- mapgrid
v <- getValues(presencegrid)
v <- ifelse(v > thresh,1,0)
presencegrid <- setValues(presencegrid, v)

image( presencegrid,
col=c("white","forestgreen"),
zlim=c(0,1),
asp=1,
bty="n",
xaxt="n", yaxt="n",

47
main="",xlab="",ylab="")
if(i==2){
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=c("Present","Absent"),
fill=c("forestgreen","white"),
bty="n",
cex=1.2)}
mtext([Link][i],side=3,line=2,cex=1.2)
mtext(paste("threshold =",thresh),side=3,line=.5,cex=1)
}
R> mtext(MODELfn.a,side=3,line=0,cex=1.2,outer=TRUE)
R> mtext([Link].a,side=3,line=2,cex=1.5,outer=TRUE)
R> par(opar)

R> opar <- par(mfrow=c(2,2),mar=c(2.5,3,4,1),oma=c(0,0,4,6),xpd=NA)

R> mapgrid <- raster(paste(MODELfn.b,"_map.img",sep=""))
R> criteria <- c("Default","MaxKappa","ReqSens","ReqSpec")
R> [Link]<-c("Default","MaxKappa","ReqSens = 0.9","ReqSpec = 0.9")
R> for(i in 1:4){
thresh <- [Link].b$threshold[[Link].b$[Link]==criteria[i]]
presencegrid <- mapgrid
v <- getValues(presencegrid)
v <- ifelse(v > thresh,1,0)
presencegrid <- setValues(presencegrid, v)

image( presencegrid,
col=c("white","forestgreen"),
xlab="",ylab="",xaxt="n", yaxt="n",
zlim=c(0,1),
asp=1,bty="n",main="")
if(i==2){
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=c("Present","Absent"),
fill=c("forestgreen","white"),
bty="n",
cex=1.2)}
mtext([Link][i],side=3,line=2,cex=1.2)
mtext(paste("threshold =",thresh),side=3,line=.5,cex=1)
}
R> mtext(MODELfn.b,side=3,line=0,cex=1.2,outer=TRUE)
R> mtext([Link].b,side=3,line=2,cex=1.5,outer=TRUE)
R> par(opar)

3.4 Example 3 - Random Forest - Categorical Response

Example 3 builds a categorical response model for vegetation category. The response variable
consists of four categories: TREE, SHRUB, OTHERVEG, and NONVEG. This model will use the
same predictors as Model 2. Out-of-bag estimates are used for model validation.

3.4.1 Set up

Define model type.

48
PINYON
VModelMapEx2a
Default MaxKappa
threshold = 0.5 threshold = 0.62

Present
Absent

ReqSens = 0.9 ReqSpec = 0.9

threshold = 0.76 threshold = 0.22

Figure 28: Example 2 - Presence-Absence maps by four different threshold selection criteria for
Pinyon (RF model).

49
SAGE
VModelMapEx2b
Default MaxKappa
threshold = 0.5 threshold = 0.61

Present
Absent

ReqSens = 0.9 ReqSpec = 0.9

threshold = 0.44 threshold = 0.79

Figure 29: Example 2 - Presence-Absence maps by four different threshold selection criteria for
Sage (RF model).

50
R> [Link] <- "RF"

Define data.

R> qdatafn <- "[Link]"

Define folder.

R> folder <- getwd()

Define model filenames.

R> MODELfn <- "VModelMapEx3"

R> predList <- c( "ELEV250",

"NLCD01_250",
"EVI2005097",
"NDV2005097",
"NIR2005097",
"RED2005097")
R> predFactor <- c("NLCD01_250")

Define the data column to use as the response, and if it is continuous, binary or categorical.

R> [Link] <- "VEGCAT"

R> [Link] <- "categorical"

Define the seeds for each model.

R> seed <- 44

Define the column that contains unique identifiers.

R> [Link] <- "ID"

Define raster look up table.

R> rastLUTfn <- "VModelMapData_LUT.csv"

R> rastLUTfn <- [Link]( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
R> rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")

51
3.4.2 Model creation

3.4.3 Model Diagnostics

Make Out-Of-Bag model predictions on the training data and run the diagnostics on these pre-
dictions. Save PDF versions of the diagnostic plots.
Out of Bag model predictions for a Random Forest model are not stochastic, so it is not necessary
to set the seed.
Since this is a categorical response model model diagnostics include a CSV file the observed and
predicted values, as well as a CSV file of the confusion matrix and its associated Kappa value and
MAUC.

Take a closer look at the text file output for the confusion matrix. Read this file into R now.

R> [Link] <- [Link]( paste( MODELfn,"_pred_cmx.csv",sep=""),

header=FALSE,
sep=",",
stringsAsFactors=FALSE)
R> [Link]

V1 V2 V3
1 Kappa [Link] observed
2 0.279154 0.0238329 NONVEG
3 predicted NONVEG 492
4 predicted OTHERVEG 7
5 predicted SHRUB 80

52
6 predicted TREE 90
7 total total 669
8 Omission Omission 0.26457399103139
9 MAUC 0.794949436880685 cmx
V4 V5 V6
1 observed observed observed
2 OTHERVEG SHRUB TREE
3 37 101 143
4 17 8 1
5 19 114 1
6 3 2 76
7 76 225 221
8 0.776315789473684 0.493333333333333 0.656108597285068
9 cmx cmx cmx
V7 V8
1 total Commission
2 total Commission
3 773 0.363518758085382
4 33 0.484848484848485
5 214 0.467289719626168
6 171 0.555555555555556
7 1191 PCC
8 PCC 0.586901763224181
9 cmx cmx

The PresenceAbsence package function Kappa() is used to calculate Kappa for the confusion
matrix. Note that while most of the functions in the PresenceAbsence package are only applicable
to binary confusion matrices, the Kappa() function will work on any size confusion matrix.
The HandTill2001 package is used to calculate the Multiple class Area under the Curve (MAUC)
and described by Hand and Till (2001).
The text file output of the confusion matrix is designed to be easily interpreted in Excel, but is
not very workable for carrying out analysis in R. However, it is relativly easy to calculate the
confusion matrix from the CSV file of the observed and predicted values.

R> PRED <-[Link]( paste( MODELfn,"_pred.csv",sep=""),

header=TRUE,
sep=",",
stringsAsFactors=TRUE)
R> head(PRED)

ID obs pred NONVEG OTHERVEG SHRUB

1 1 NONVEG NONVEG 0.4730539 0.029940120 0.32934132
2 2 NONVEG TREE 0.4406780 0.000000000 0.01129944
3 3 NONVEG NONVEG 0.7807487 0.000000000 0.02673797
4 4 TREE NONVEG 0.5606936 0.000000000 0.06936416
5 5 NONVEG NONVEG 0.5164835 0.005494505 0.12637363
6 6 SHRUB NONVEG 0.4480874 0.000000000 0.37704918
TREE
1 0.1676647
2 0.5480226
3 0.1925134
4 0.3699422
5 0.3516484
6 0.1748634

53
For categorical models, this file contains the observed category for each location, the category
predicted by majority vote, as well as one column for each category observed in the data, giving
the proportion of trees that voted for that category.
To calculate the confusion matrix from the file we will use the observed and predicted columns.
The [Link]() function will convert columns containing character strings to factors. If the
categories had been numerical, the [Link]() function can be used to convert the columns to
factors. Because there may be catergories present in the observed data that are missing from the
predictions (and vice versa), to get a symetric confusion matrix it is important to make sure all
levels are present in both factors.
The following code will work for both numerical and character categories:

R> #
R> #these lines are needed for numeric categories, redundant for character categories
R> #
R> PRED$pred<-[Link](PRED$pred)
R> PRED$obs<-[Link](PRED$obs)
R> #
R> #adjust levels so all values are included in both observed and predicted
R> #
R> LEVELS<-unique(c(levels(PRED$pred),levels(PRED$obs)))
R> PRED$pred<-factor(PRED$pred,levels=LEVELS)
R> PRED$obs<- factor(PRED$obs, levels=LEVELS)
R> #
R> #calculate confusion matrix
R> #
R> CMX<-table( predicted=PRED$pred, observed= PRED$obs)
R> CMX

observed
predicted NONVEG OTHERVEG SHRUB TREE
NONVEG 492 37 101 143
OTHERVEG 7 17 8 1
SHRUB 80 19 114 1
TREE 90 3 2 76

To calculate the errors of Omission and Comission:

R> [Link] <- diag(CMX)

R> [Link] <- 1-([Link]/apply(CMX,2,sum))
R> [Link] <- 1-([Link]/apply(CMX,1,sum))
R> [Link]

NONVEG OTHERVEG SHRUB TREE

0.2645740 0.7763158 0.4933333 0.6561086

R> [Link]

NONVEG OTHERVEG SHRUB TREE

0.3635188 0.4848485 0.4672897 0.5555556

To calculate PCC:

R> [Link] <- sum([Link])/sum(CMX)

R> [Link]

54
[1] 0.5869018

To calculate Kappa:

R> [Link] <- PresenceAbsence::Kappa(CMX)

R> [Link]

Kappa [Link]
NONVEG 0.2791541 0.02383285

The MAUC is calculated from the category specific predictions (the percent of trees that voted
for each category):

R> VOTE <- HandTill2001::multcap( response = PRED$obs,

predicted= [Link](PRED[,-c(1,2,3)]) )
R> MAUC <- HandTill2001::auc(VOTE)
R> MAUC

[1] 0.7949483

Note, both the PresenceAbsence package and the HandTill2001 package have functions named
auc(). The :: operator is used to specify that we are calling the auc() function from the
HandTill2001 package.
As in continuous and binary response models, the [Link]() function creates a vari-
able importance graph (Figure 30).
With categorical response models, the [Link]() function alse creates category spe-
cific variable importance graphs, for example (Figure 31).

3.4.4 Comparing Variable Importance

In example 1 the [Link]() function was used to compare the importance be-
tween two continuous Random Forest models, for percent cover of Pinyon and of Sage. Here we
will compare the variable importances of the binary models from Example 2 with the categorical
model we have just created in example 3 (Figure 32). We are examining the question “Are the
same predictors important for determining vegetaion category as were important for determining
species presence?” Keep in mind that to use [Link]() the two models must be
built from the same predictor variables. For example, we could not use it to compare the models
from Example 1 and Example 3, because in Example 1 "NLCD01_250" was not included in the
predictors.

R> opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))

R> [Link]( [Link].1=[Link].ex2a,
[Link].2=[Link].ex3,
[Link].1="Pinyon",
[Link].2="VEGCAT",
[Link]=FALSE,
[Link]="predList",
predList=predList,
[Link]="sum",
main="Pinyon Presence vs VEGCAT",
[Link]="none",
cex=0.9)
R> [Link]( [Link].1=[Link].ex2b,

55
Relative Influence
VModelMapEx3_pred

NLCD01_250 ● ELEV250 ●

ELEV250 ● EVI2005097 ●

EVI2005097 ● NDV2005097 ●

NDV2005097 ● RED2005097 ●

RED2005097 ● NIR2005097 ●

NIR2005097 ● NLCD01_250 ●

20 40 60 0 50 100
MeanDecreaseAccuracy MeanDecreaseGini

Figure 30: Example 3 - Overall variable importance graph for predicting vegetation category.

56
VModelMapEx3_pred
Relative Influence − TREE − 221 plots

NLCD01_250 ●

ELEV250 ●

NDV2005097 ●

EVI2005097 ●

NIR2005097 ●

RED2005097 ●

10 20 30 40 50 60

TREE

Figure 31: Example 3 - Category specific variable importance graph for vegetation category ”Tree”.

57
Variable Importance Comparison
Pinyon Presence vs VEGCAT

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

Pinyon VEGCAT

Sage Presence vs VEGCAT

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

Sage VEGCAT

Figure 32: Example 3 - Comparison of variable importances between categorical model of vegeta-
tion category (VEGCAT) and binary models for Pinyon presence and Sage presence.

[Link].2=[Link].ex3,
[Link].1="Sage",
[Link].2="VEGCAT",
[Link]=FALSE,
[Link]="predList",
predList=predList,
[Link]="sum",
main="Sage Presence vs VEGCAT",
[Link]="none",
cex=0.9)
R> mtext("Variable Importance Comparison",side=3,line=0,cex=1.8,outer=TRUE)
R> par(opar)

With categorical models the [Link]() function also can be usde to compare
the variable importance between categories of the same model. The importance measure used
for category specific importance is the relative influence of each variable, calculated by randomly
permuting each predictor variable, and looking at the decrease in model accuracy associated with
each predictor.

R> opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))

58
R> [Link]( [Link].1=[Link].ex3,
[Link].2=[Link].ex3,
[Link].1="SHRUB",
[Link].2="TREE",
class.1="SHRUB",
class.2="TREE",
[Link]="predList",
predList=predList,
[Link]="sum",
main="VEGCAT - SHRUB vs. TREE",
[Link]="none",
cex=0.9)
R> [Link]( [Link].1=[Link].ex3,
[Link].2=[Link].ex3,
[Link].1="OTHERVEG",
[Link].2="NONVEG",
class.1="OTHERVEG",
class.2="NONVEG",
[Link]="predList",
predList=predList,
[Link]="sum",
main="VEGCAT - OTHERVEG vs. NONVEG",
[Link]="none",
cex=0.9)
R> mtext("Category Specific Variable Importance",side=3,line=0,cex=1.8,outer=TRUE)
R> par(opar)

3.4.5 Interaction Plots

We will look at how the [Link]() function behaves with a categorical response
variable. With categorical models, interactions can affect one prediction category, without influ-
encing other categories. For example, if modelling disturbance type, it is possible that landslides
might be influenced by an interaction between soil type and slope, while fires might be influenced
by both variables individually, but without any interaction.
Therefore when calling [Link]() it is neccessary to specify a particular cate-
gory. The function then will graph how the probability of that category varies as a function of the
two specified predictor variables.
Here we look at the interaction between elevation and land cover class for two of our response
categories (Figure 34, Figure 35 ).

59
Category Specific Variable Importance
VEGCAT − SHRUB vs. TREE

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

MeanDecAccuracy MeanDecAccuracy
SHRUB TREE

VEGCAT − OTHERVEG vs. NONVEG

ELEV250

NLCD01_250

EVI2005097

NDV2005097

NIR2005097

RED2005097

MeanDecAccuracy MeanDecAccuracy
OTHERVEG NONVEG

Figure 33: Example 3 - Category specific comparison of variable importances for the model of
vegetation category (VEGCAT). National Land Cover DATA (NLCD) and elevation (ELEV)
are the most important predictor variables for both TREE and SHRUB categories, though the
relative importance of the remote sensing bands differs between these two categories. ELEV is
relativly less important for the NONVEG and OTHERVEG categories, while ELEV is important
for classifying OTHERVEG but relativly unimportant for the classifying NONVEG. In other
words, if the model lost the information contained in NLCD and ELEV, the predictions for TREE
and SHRUB categories would suffer, but there would be less of an effect on the prediction accuracy
for NONVEG. The predictions for OTHERVEG would suffer if NLCD were removed from the
model, but would lose relativly little accuracy if ELEV were removed.

60
main=[Link],
[Link]="image",
[Link]="pdf",
MODELfn=MODELfn,
folder=folder,
[Link]="NONVEG")
R>

3.4.6 Map production

The function [Link]() creates an ascii text files and an imagine image file of predictions
for each map pixel.

R> [Link]( [Link]=[Link].ex3,

folder=folder,
MODELfn=MODELfn,
rastLUTfn=rastLUTfn,
[Link]="[Link]")

With categorical models, the [Link]() function outputs a map file, using integer codes
for each category, along with a table relating these codes to the original categories. In this example
[Link] was set to "[Link]", therefore pixels with factored predictors with values not found
in the training data will be omited from the map.
Take a look at the codes:

R> [Link]<-[Link]( paste(MODELfn,"_map_key.csv",sep=""),

header=TRUE,
sep=",",
stringsAsFactors=FALSE)
R> [Link]

row category integercode

1 1 NONVEG 1
2 2 OTHERVEG 2
3 3 SHRUB 3
4 4 TREE 4

Column one gives the row number for each code. Column two gives the category names. Column
three gives the integer codes used to represent each category in the map output. In this example the
categories in the training data are character strings, and [Link]() assigned the integers
1 through the Number of categories. If the training data categories were already numeric codes,
for example, c(30,42,50,80), then [Link]() would keep the original values in the map
output, and columns two and three would contain the same values.
Next we define a color for each category. The colors() function will generate a list of the possible
color names. Some are quite poetical.

R> [Link]$colors<-c("bisque3","springgreen2","paleturquoise1","green4")
R> [Link]

Import the map output and transform the values of the map output to intergers from 1 to n (the
number of map categories).

61
VEGCAT

1.0
80

0.8
70

probability of SHRUB
0.6
NLCD01_250
(categorical)

0.4

0.2

30
0.0

500 1000 1500 2000 2500 3000

ELEV250

Figure 34: Example 3 - Interaction plot for SHRUB vegetation category (RF model), showing inter-
actions between elevation and National Land Cover Dataset classes (ELEV250 and NLCD01 250).
Image plot, with darker green indicating higher probability of being assigned to the specified cate-
gory. Here we see the direct effect of NLCD class: SHRUB has a low probability of being assigned
to landcover class 42. We also see a direct effect of elevation, with SHRUB having a slightly
higher probability of being assigned to middle elevations. There is little evidence of interaction
between these two predictors, since relationship of probability to elevation is similar for all land
cover classes.

62
VEGCAT

1.0
80

0.8
70

probability of NONVEG
0.6
NLCD01_250
(categorical)

0.4

0.2

30
0.0

500 1000 1500 2000 2500 3000

ELEV250

Figure 35: Example 3 - Interaction plot for NONVEG vegetation category (RF model), show-
ing interactions between elevation and National Land Cover Dataset classes (ELEV250 and
NLCD01 250). Image plot, with darker green indicating higher probability of being assigned
to the specified category. NONVEG has a much higher chance of being assigned than SHRUB
in all combinations of the two predictor variables, reflecting its higher prevalence in the training
data (56% as opposed to 19%). NONVEG does show some interaction between landcover class
and elevation. In all land cover classes NONVEG has a higher chance of being predicted at low
elevations, but while in most land cover classes the probability goes down at elevations above
1400m, in land cover class 42 NONVEG matains a high chance of being predicted to over 2000m.

63
Note that here in example 3, where the categorical responses were character strings, the [Link]()
function generated integer codes from 1 to n for the raster output. So for this example this step is
not actually neccessary. However, if the categories in the training data had been unevenly spaced
numeric codes, for example, c(30,42,50,80), then [Link]() would keep these original
numeric codes in the raster output. In such cases, to produce a map in R with the image()
function (as opposed to viewing the image in a GIS environment) creating a new raster where
the numeric codes are replaced with the numbers 1 to n makes assigning specific colors to each
category simpler.

R> mapgrid <- raster(paste(MODELfn,"_map.img",sep=""))

R> integergrid <- mapgrid
R> v <- getValues(mapgrid)
R> v <- [Link]$row[match(v,[Link]$integercode)]
R> integergrid <- setValues(integergrid, v)

Produce the map (Figure 36).

R> opar <- par(mfrow=c(1,1),mar=c(3,3,2,1),oma=c(0,0,3,8),xpd=NA)

R> image( integergrid,
col = [Link]$colors,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=c(1,nrow([Link])),
main="",asp=1,bty="n")
R> mtext([Link],side=3,line=1,cex=1.2)
R> legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=[Link]$category,
fill=[Link]$colors,
bty="n",
cex=1.2)
R> par(opar)

4 Conclusion
In summary, the ModelMap software package for R creates sophisticated models from training data
and validates the models with an independent test set, cross-validation, or in the case of Random
Forest Models, with out-of-bag (OOB) predictions on the training data. It creates graphs and
tables of the model diagnostics. It applies these models to GIS image files of predictors to create
detailed prediction surfaces. It will handle large predictor files for map making, by reading in the
GIS data in sections, and output the prediction for each of these sections, before reading the next
section.

64
VEGCAT

NONVEG
OTHERVEG
SHRUB
TREE

Figure 36: Example 3 - Map of predicted vegetation category (RF model). The white pixels found
just below the center of this map idicate pixels with factored predictor variables that have values
not found in the training data.

65
Appendices
Arguments for [Link]()
[Link] Model type: "RF" or "QRF" or "CF".
[Link] Filename of the training data file for building model.
folder Folder for all output.
MODELfn Filename to save model object.
predList Predictor short names used to build the model.
predFactor Predictors from predList that are factors (i.e categorical).
[Link] Response variable used to build the model.
[Link] Response type: "binary" or "continuous".
[Link] Unique identifier for each row in the training data.
seed Seed to initialize randomization to build stochastic models.
[Link] Specifies the action to take if there are NA values in the prediction data
[Link] Should a copy of the predictor data be included in the model object. Useful
if [Link] will be used later.
Random Forest Models:
ntree Number of random forest trees.
mtry Number of variables to try at each node of Random Forest trees.
replace Should sampling be done with or without replacement.
strata A (factor) variable that is used for stratified sampling.
sampsize For classification, if strata provided, sampling is stratified by strata. For
binary response models, if argument strata is not provided then sampling
is stratified by presence/absence.

66
Arguments for [Link]
[Link] The model object to use for prediction, if the model has been previously
created.
[Link] Filename of the training data file for building model.
[Link] Filename of independent data set for testing (validating) model.
folder Folder for all output.
MODELfn Filename to save model object.
[Link] Response variable used to build the model.
[Link] Name of column in training and test that uniquely identifies each row .
[Link] Name of column in training that indicates subset of data to use for diag-
nostics.
seed Seed to initialize randomization to build stochastic models.
[Link] Type of prediction to use for model validation: "TEST", "CV", "OOB" or
"TRAIN"
MODELpredfn Filename for output of validation prediction *.csv file.
[Link] Specifies the action to take if there are NA values in the prediction data or if
there is a level or class of a categorical predictor variable in the validation
test set or the mapping data set, but not in the training data set.
[Link] The number of cross-validation folds.
[Link] Vector of one or more device types for graphical output: "default",
"jpeg", "pdf", "postscript", "[Link]". "default" refers to the
default graphics device for your computer
DIAGNOSTICfn Filename for output files from model validation diagnostics.
[Link] Pixels per inch for jpeg output.
[Link] Device width for diagnostic plots in inches.
[Link] Device height for diagnostic plots in inches.
cex Cex for diagnostic plots.
[Link] Required sensitivity for threshold optimization for binary response model.
[Link] Required specificity for threshold optimization for binary response model.
FPC False Positive Cost for threshold optimization for binary response model.
FNC False Negative Cost for threshold optimization for binary response model.

Arguments for [Link]()

[Link] The model object to use for prediction, if the model has been pre-
viously created.
folder Folder for all output.
MODELfn Filename to save model object.
rastLUTfn Filename of .csv file for a raster look up table.
[Link] Specifies the action to take if there are NA values in the prediction
data or if there is a level or class of a categorical predictor variable
in the validation test set or the mapping data set, but not in the
training data set.
[Link] If TRUE then the raster brick containing the predictors from the
model object is saved as a native raster package format file.
[Link] Should maps of mean, standard deviation, and coefficient of varia-
tion of the predictions be produced: TRUE or FALSE. Only used if
[Link] = "continuous".
OUTPUTfn Filename for output raster file for map production. If NULL,
"modelfn _map.txt".

67
References
L. Breiman. Random forests. Machine Learning, 45(1):5–32, 2001.
L. Breiman, R. A. Friedman, R. A. Olshen, and C. G. Stone. Classification and Regression Trees.
Wadsworth, 1984.
G. De’ath and K. E. Fabricius. Classification and regression trees: a powerful yet simple technique
for ecological data analysis. Ecology, 81:3178–3192, 2000.
E. R. DeLong, D. M. Delong, and D. L. Clarke-Pearson. Comparing areas under two or more
correlated receiver operating characteristic curves: A nonparametric approach. Biometrics,
44(3):387–394, 1988.

J. Elith, J. R. Leathwick, and T. Hastie. A working guide to boosted regression trees. Journal of
Animal Ecology, 77:802–813, 2008.
J. S. Evans and S. A. Cushman. Gradient modeling of conifer species using random forests.
Landscape Ecology, 24(5):673–683, 2009.
A. H. Fielding and J. F. Bell. A review of methods for the assessment of prediction errors in
conservation presence/absence models. Environmental Conservation, 24(1):38–49, 1997.
E. Freeman. PresenceAbsence: An R Package for Presence-Absence Model Evaluation. USDA
Forest Service, Rocky Mountain Research Station, 507 25th street, Ogden, UT, USA, 2007.
URL [Link] eafreeman@[Link].

E. Freeman. ModelMap: An R Package for Modeling and Map production using Random Forest
and Stochastic Gradient Boosting. USDA Forest Service, Rocky Mountain Research Station, 507
25th street, Ogden, UT, USA, 2009. URL [Link] eafreeman@[Link].
E. A. Freeman and G. Moisen. PresenceAbsence: An R package for presence absence analysis.
Journal of Statistical Software, 23(11):1–31, 2008a. URL [Link]

E. A. Freeman and G. G. Moisen. A comparison of the performance of threshold criteria for

binary classification in terms of predicted prevalence and kappa. Ecological Modelling, 217:
48–58, 2008b.
T. S. Frescino, G. G. Moisen, K. A. Megown, V. J. Nelson, Elizabeth, Freeman, P. L. Patterson,
M. Finco, K. Brewer, and J. Menlove. Nevada photo-based inventory pilot(npip) photo sampling
procedures. Gen. Tech. Rep. RMRSGTR-222, U.S. Departmentof Agriculture, Forest Service,
Rocky Mountain Research Station., Fort Collins, CO, 2009.
J. H. Friedman. Greedy function approximation: a gradient boosting machine. Annals of Statistics,
29(5):1189–1232, 2001.

J. H. Friedman. Stochastic gradient boosting. Computational Statistics & Data Analysis, 38(4):
367–378, 2002.
D. Gesch, M. Oimoen, S. Greenlee, C. Nelson, M. Steuck, and D. Tyler. The national elevation
dataset. photogrammetric engineering and remote sensing. Photogrammetric Engineering and
Remote Sensing, 68:5–11, 2002.

D. J. Hand and R. J. Till. A simple generalisation of the area under the roc curve for multiple
class classification problems. Machine Learning, 45(2):171–186, 2001.
C. Homer, C. Huang, L. Yang, B. Wylie, and M. Coan. Development of a 2001 national land-cover
database for the united states. Photogrammetric Engineering and Remote Sensing, 70:829–840,
2004.

68
A. Huete, K. Didan, T. Miura, E. P. Rodriguez, X. Gao, and L. G. Ferreira. Overview of the
radiometric and biophysical performance of the modis vegetation indices. Remote Sensing of
Environment, 83:195–213, 2002.
C. O. Justice, J. R. G. Townshend, E. F. Vermote, E. Masuoka, R. E. Wolfe, N. Saleous, D. P. Roy,
and J. T. Morisette. An overview of modis land data processing and product status. Remote
Sensing of Environment, 83:3–15, 2002.
A. Liaw and M. Wiener. Classification and regression by randomForest. R News, 2(3):18–22, 2002.
URL [Link]
Y. Lin and Y. Jeon. Random forest and adaptive nearest neighbors. Technical Report 1055,
Department of Statistics, University of Wisconsin, 1210 West Dayton St., Madison, WI 53706,
2002.
G. G. Moisen. Classification and regression trees. In S. E. Jørgensen and B. D. Fath, editors,
Encyclopedia of Ecology, volume 1, pages 582–588. Elsevier, 2008.
G. G. Moisen, E. A. Freeman, J. A. Blackard, T. S. Frescino, N. E. Zimmermann, and T. C.
Edwards, Jr. Predicting tree species presence in utah: a comparison of stochastic gradient
boosting, generalized additive models, and tree-based methods. Ecological Modelling, 199:176–
187, 2006.
J. Pearce and S. Ferrier. Evaluating the predicting performance of habitat models developed using
logistic regression. Ecological Modelling, 133:225–245, 2000.

R Development Core Team. R: A Language and Environment for Statistical Computing. R Foun-
dation for Statistical Computing, Vienna, Austria, 2008. URL [Link]
ISBN 3-900051-07-0.
B. Reineking and B. Schröder. Constrain to perform: Regularization of habitat models. Ecological
Modelling, 193:675–690, 2006.
C. Strobl, A.-L. Boulesteix, A. Zeileis, and T. Hothorn. Bias in random forest variable importance
measures: Illustrations, sources and a solution. Bioinformatics, 8:25, 2007.
I. P. Vaughan and S. J. Ormerod. The continuing challenges of testing species distribution models.
Journal of Applied Ecology, 42:720–730, 2005.

M. P. Vayssieres, R. P. Plant, and B. H. Allen-Diaz. Classification trees: An alternative non-

parametric approach for predicting species distributions. Journal of vegetation science, 11:
679–694, 2000.
J. H. Zar. Biostatistical Analysis. Prentice Hall, 1996.

Modified Assignment Algorithm for Layout
No ratings yet
Modified Assignment Algorithm for Layout
14 pages
Inclusion Measure for Interval Neutrosophic Sets
No ratings yet
Inclusion Measure for Interval Neutrosophic Sets
3 pages
BUC Algorithm For Iceberg Cubes: Implementation and Sensitivity Analysis
No ratings yet
BUC Algorithm For Iceberg Cubes: Implementation and Sensitivity Analysis
5 pages
STRIPS, Subconjunto de PDDL
No ratings yet
STRIPS, Subconjunto de PDDL
7 pages
Data Structures and Algorithms Overview
No ratings yet
Data Structures and Algorithms Overview
7 pages
Stack vs Queue: Data Structure Analysis
No ratings yet
Stack vs Queue: Data Structure Analysis
5 pages
Amphioxus Development: Fertilization to Gastrulation
No ratings yet
Amphioxus Development: Fertilization to Gastrulation
6 pages
C++ Template Functions and Set Operations
No ratings yet
C++ Template Functions and Set Operations
3 pages
Haskell Solution for Snake Cube Puzzle
100% (1)
Haskell Solution for Snake Cube Puzzle
15 pages
Feature Selection in High Dimensional Data
No ratings yet
Feature Selection in High Dimensional Data
6 pages
Boruta Feature Selection in R
No ratings yet
Boruta Feature Selection in R
1 page
Clustering - Fuzzy C-Means
No ratings yet
Clustering - Fuzzy C-Means
4 pages
Depth vs. Breadth in Menu Design
No ratings yet
Depth vs. Breadth in Menu Design
18 pages
3D Object Recognition from Drawings
No ratings yet
3D Object Recognition from Drawings
6 pages
Diagrams in Non-Routine Problem Solving
No ratings yet
Diagrams in Non-Routine Problem Solving
8 pages
Automatic Road Extraction from Satellite Images
No ratings yet
Automatic Road Extraction from Satellite Images
8 pages
Hierarchical Grid Structures for GIS Data
No ratings yet
Hierarchical Grid Structures for GIS Data
42 pages
Comparative Study of Feature Selection Algorithms
No ratings yet
Comparative Study of Feature Selection Algorithms
10 pages
Understanding Software Modifiability
No ratings yet
Understanding Software Modifiability
16 pages
Software Design Principles and Practices
No ratings yet
Software Design Principles and Practices
11 pages
Tree Models of Software Processes
No ratings yet
Tree Models of Software Processes
5 pages
Unit 3: Software Design Concepts
No ratings yet
Unit 3: Software Design Concepts
29 pages
Chapter Four
No ratings yet
Chapter Four
7 pages
Algorithmic Design of Islamic Ornaments
No ratings yet
Algorithmic Design of Islamic Ornaments
11 pages
Property Module
No ratings yet
Property Module
15 pages
Morphological Process in Architecture
No ratings yet
Morphological Process in Architecture
11 pages
DS Module-1
No ratings yet
DS Module-1
44 pages
Group L - Interactive System Development Final With Audio
No ratings yet
Group L - Interactive System Development Final With Audio
17 pages
Efficient Hidden Markov Model System
No ratings yet
Efficient Hidden Markov Model System
3 pages
BRDK Coding Guidelines External
No ratings yet
BRDK Coding Guidelines External
17 pages
FuzzyToolkit: Clustering Software Toolkit
No ratings yet
FuzzyToolkit: Clustering Software Toolkit
6 pages
Algorithm and Flowchart Design Guide
No ratings yet
Algorithm and Flowchart Design Guide
14 pages
Scheduling 20r
No ratings yet
Scheduling 20r
6 pages
Gapless Edges of 2d Topological Orders and Enriched Monoidal Categories
No ratings yet
Gapless Edges of 2d Topological Orders and Enriched Monoidal Categories
26 pages
Soft Computing Assignment on Fuzzy Systems
No ratings yet
Soft Computing Assignment on Fuzzy Systems
4 pages
Geometrical Description of Crystals
No ratings yet
Geometrical Description of Crystals
13 pages
Understanding Web Mining Techniques
No ratings yet
Understanding Web Mining Techniques
16 pages
Revisedddddd Pinaka Latest Shit
No ratings yet
Revisedddddd Pinaka Latest Shit
28 pages
RIA Patterns: Auto-Completion & Refresh
No ratings yet
RIA Patterns: Auto-Completion & Refresh
3 pages
SOftware Requirements2
No ratings yet
SOftware Requirements2
10 pages
COBWEB Algorithm in Conceptual Clustering
No ratings yet
COBWEB Algorithm in Conceptual Clustering
5 pages
Handwritten Document
No ratings yet
Handwritten Document
30 pages
Search Algorithms for Conceptual Graphs
No ratings yet
Search Algorithms for Conceptual Graphs
10 pages
Hexahedral Mesh Tutorial for HM-3200
No ratings yet
Hexahedral Mesh Tutorial for HM-3200
27 pages
Containerizing Server Applications with Docker
No ratings yet
Containerizing Server Applications with Docker
11 pages
Yard Space Allocation Algorithm For Unloading Cont
No ratings yet
Yard Space Allocation Algorithm For Unloading Cont
22 pages
Informatics Tournament Tasks Solutions
No ratings yet
Informatics Tournament Tasks Solutions
3 pages
Understanding Design as Relation
No ratings yet
Understanding Design as Relation
7 pages
Understanding Data Structures in C
No ratings yet
Understanding Data Structures in C
10 pages
Lbo Homology in Semi-Groups and Knots
No ratings yet
Lbo Homology in Semi-Groups and Knots
27 pages
Polymaths Chronoscope
No ratings yet
Polymaths Chronoscope
4 pages
Understanding Support Vector Machines and Algorithms
No ratings yet
Understanding Support Vector Machines and Algorithms
6 pages
Enhancing Security in MapReduce Systems
No ratings yet
Enhancing Security in MapReduce Systems
8 pages
K-means vs Hierarchical Clustering
No ratings yet
K-means vs Hierarchical Clustering
4 pages
White Paper Best Practices For Artifactory Naming Conventions
No ratings yet
White Paper Best Practices For Artifactory Naming Conventions
20 pages
Arch Linux: Embracing Simplicity and Freedom
No ratings yet
Arch Linux: Embracing Simplicity and Freedom
2 pages
Lattice Generator for 3D Printing STL
No ratings yet
Lattice Generator for 3D Printing STL
7 pages
Programmatic Lattice Generation For Additive Manufacture: Sciencedirect
No ratings yet
Programmatic Lattice Generation For Additive Manufacture: Sciencedirect
7 pages
This Guide Provides A Brief Introduction To Random Forests®
No ratings yet
This Guide Provides A Brief Introduction To Random Forests®
49 pages
Estimating Random Utility Models in R
No ratings yet
Estimating Random Utility Models in R
40 pages
Estimating Multinomial Logit in R
No ratings yet
Estimating Multinomial Logit in R
73 pages
Random Forest Variants in R Models
No ratings yet
Random Forest Variants in R Models
19 pages
Mlogit Exercises for Heating System Choice
No ratings yet
Mlogit Exercises for Heating System Choice
33 pages
Exploratory Data Analysis with ModelMap
No ratings yet
Exploratory Data Analysis with ModelMap
20 pages
Combining Multiple Imputations in R
No ratings yet
Combining Multiple Imputations in R
5 pages
adephylo: Phylogenetic Data Analysis
No ratings yet
adephylo: Phylogenetic Data Analysis
20 pages
Journal of Statistical Software: Spacetime: Spatio-Temporal Data in R
No ratings yet
Journal of Statistical Software: Spacetime: Spatio-Temporal Data in R
30 pages
Applied Econometrics With R: Package Vignette and Errata: Christian Kleiber Achim Zeileis
No ratings yet
Applied Econometrics With R: Package Vignette and Errata: Christian Kleiber Achim Zeileis
6 pages
Spatial Statistics Exercises - Fall 2003
No ratings yet
Spatial Statistics Exercises - Fall 2003
2 pages
Introduction to Spatial Statistics Course
No ratings yet
Introduction to Spatial Statistics Course
4 pages
Object-Oriented Programming Course Overview
No ratings yet
Object-Oriented Programming Course Overview
2 pages
Persistent Homology in Data Analysis
No ratings yet
Persistent Homology in Data Analysis
56 pages
Best States for Sunglasses Sales
No ratings yet
Best States for Sunglasses Sales
10 pages
Big Data Hadoop Security Challenges
No ratings yet
Big Data Hadoop Security Challenges
6 pages
Divi - The Ultimate WordPress Theme & Visual Page Builder
No ratings yet
Divi - The Ultimate WordPress Theme & Visual Page Builder
64 pages
Cisco UCS C220 M5 Server Installation and Service Guide: Americas Headquarters
No ratings yet
Cisco UCS C220 M5 Server Installation and Service Guide: Americas Headquarters
142 pages
Overview of Enterprise Computing
No ratings yet
Overview of Enterprise Computing
4 pages
Design vs Coding in Software Engineering
No ratings yet
Design vs Coding in Software Engineering
20 pages
CPGET 2022 M.A. Public Administration Exam
No ratings yet
CPGET 2022 M.A. Public Administration Exam
73 pages
Grade 6 ICT Past Exam Paper
No ratings yet
Grade 6 ICT Past Exam Paper
4 pages
LM238WF2 23.8" Full HD LCD Specs
No ratings yet
LM238WF2 23.8" Full HD LCD Specs
33 pages
Mobile X-Ray Unit Specifications
No ratings yet
Mobile X-Ray Unit Specifications
14 pages
Decision Tree Exercise Overview
No ratings yet
Decision Tree Exercise Overview
42 pages
IBM Power Systems Hardware Planning Guide
No ratings yet
IBM Power Systems Hardware Planning Guide
224 pages
Infosys Management Study Report
No ratings yet
Infosys Management Study Report
56 pages
Quantum Computing's Role in Net-Zero Energy
No ratings yet
Quantum Computing's Role in Net-Zero Energy
22 pages
Shimadzu RF-5301PC Spectrofluorophotometer
No ratings yet
Shimadzu RF-5301PC Spectrofluorophotometer
12 pages
CCVS Training for VSX Specialists
No ratings yet
CCVS Training for VSX Specialists
1 page
Sivakasi Fireworks Price List 2024
No ratings yet
Sivakasi Fireworks Price List 2024
6 pages
Full Stack Developer Profile: Prasun Chakraborty
No ratings yet
Full Stack Developer Profile: Prasun Chakraborty
1 page
Data Science Lab Experiments Overview
No ratings yet
Data Science Lab Experiments Overview
22 pages
Network Basics Revision Worksheet
No ratings yet
Network Basics Revision Worksheet
2 pages
Ilovepdf Merged
No ratings yet
Ilovepdf Merged
44 pages
SONOCHEK - The Ultrasonic Testing Device
No ratings yet
SONOCHEK - The Ultrasonic Testing Device
6 pages
Master VB.NET in 2 Weeks Guide
No ratings yet
Master VB.NET in 2 Weeks Guide
236 pages
Career Path: Information Security Analyst
No ratings yet
Career Path: Information Security Analyst
5 pages
Understanding SHA-1 and SHA-512 Hashing
No ratings yet
Understanding SHA-1 and SHA-512 Hashing
34 pages
SIWES Report: NTA Lagos Internship
No ratings yet
SIWES Report: NTA Lagos Internship
22 pages
FS100 Emergency Stop Instructions
No ratings yet
FS100 Emergency Stop Instructions
25 pages
Optimizing Network Traffic with MARL
No ratings yet
Optimizing Network Traffic with MARL
2 pages

ModelMap: R Package for GIS Modeling

Uploaded by

ModelMap: R Package for GIS Modeling

Uploaded by

ModelMap: an R Package for Model Creation and Map

September 10, 2018

2.2 File Names

2.3 Training Data

2.4 Independent Test Set for Model Validation

2.5 Missing predictor values

2.6 The Model Object

2.8 Raster Look Up Table

Table 1: Predictor variables

3.1 Example dataset

3.2 Example 1 - Random Forest - Continuous Response

Next define some of the arguments for the models.

R> [Link] <- "RF"

R> qdatafn <- "[Link]"

Define the output folder.

R> folder <- getwd()

R> [Link]( [Link]=0.2,

R> MODELfn.a <- "VModelMapEx1a"

Define the response variable, and whether it is continuous, binary or categorical.

R> [Link].a <- "PINYON"

Define the seeds for each model.

R> seed.a <- 38

R> [Link] <- "ID"

3.2.2 Model creation

R> [Link].ex1a <- [Link]( [Link]=[Link],

R> [Link].ex1a <- [Link]( [Link]=[Link].ex1a,

3.2.4 Comparing Variable Importance

The [Link]() function uses a back-to-back barplot to compare variable impor-

R> [Link]( [Link].1=[Link].ex1a,

Model Response Importance Type Measured By

Table 2: Importance Measures

0 100 200 300 400 500

0 100 200 300 400 500

R> opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))

3.2.5 Interaction Plots

The [Link]() function provides a diagnostic plot useful in visualizing two-way

R> [Link]( [Link].ex1a,

2000 3000 4000 5000 6000 7000

2000 3000 4000 5000 6000 7000

500 1000 1500 2000 2500 3000

500 1000 1500 2000 2500 3000

R> elevfn <- paste(folder,"/VModelMapData_dem_ELEVM_250.img",sep="")

R> opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))

R> rastLUTfn <- "VModelMapData_LUT.csv"

−1544000 −1538000 −1532000

R> [Link]( [Link]=[Link].ex1a,

R> l <- seq(100,0,[Link]=101)

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> [Link] <- hcl(h = 15, c = c, l = l)

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> mapgrid.b <- raster(paste(MODELfn.b,"_map_stdev.img",sep=""))

R> [Link] <- hcl(h = 70, c = c, l = l)

R> opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)

R> mapgrid.a <-raster(paste(MODELfn.a,"_map_coefv.img",sep=""),[Link]=TRUE)

Define model type.

R> [Link] <- "RF"

R> qdatafn <- "[Link]"

R> folder <- getwd()

Define model filenames.

R> MODELfn.a <- "VModelMapEx2a"

R> predList <- c( "ELEV250",

R> [Link].a <- "PINYON"

Define the seeds for each model.

R> seed.a <- 40

Define the column that contains unique identifiers.

R> [Link] <- "ID"

R> rastLUTfn <- "VModelMapData_LUT.csv"

3.3.2 Model creation

R> [Link].ex2a <- [Link]( [Link]=[Link],

3.3.3 Model Diagnostics

R> [Link].ex2a <- [Link]( [Link]=[Link].ex2a,

R> [Link].a <- [Link]( paste(MODELfn.a,"_pred_optthresholds.csv",sep=""),

Optimized thresholds for Pinyon:

[Link] threshold PCC sensitivity specificity

[Link] threshold [Link] pred

And for Sage: