### Vegetation Modeling I (5.3 MB ppt)

Vegetation Modeling
1
Outline
Model types
Predictive models
Predictor data
Predictive model types (parametric/nonparametric)
Model Example (Tree-based/Random Forests)
Modeling dataset
• Response/Predictor data
• Discussion – scale of predictors
• Predictor data extraction
Data exploration
•
•
•
•
Summary statistics/NA values
attaching data frame
Predictor variables
Response variables (binary, continuous)
Model generation
• Tree-based models
• Random Forests (classification/regression trees)
• Variable importance/proximity
Model prediction
• Polygon data
• Predictor data (clipping, stacking)
• Apply model and display map
Modelmap
• Build model
• Model diagnostics
• Make map
2
Vegetation Modeling
## In general, there are 3 reasons we model vegetation:
# 1. Explanatory: to explain why something is happening… find a pattern or a cause.
# 2. Descriptive: to see an association between variables, not aimed at prediction.
# 3. Predictive: to predict an occurrence based on known attributes.
## Examples:
Is the height of a tree related to its diameter?
Did the last 20 years of drought affect tree mortality rates?
Are forest disturbances changing the carbon cycle?
Can we predict the distribution of vegetation 100 years from now based on climate models?
Can we predict the current distribution of vegetation across the landscape from the spectral
signature of remotely-sensed data?
## An article on different reasons for modeling.
Shmueli, G. 2010. To Explain or to Predict? Statistical Science 25(3):289-310.
3
Predictive Vegetation Models
## Predicting the distribution of vegetation across the landscape from available maps of
environmental variables, such as geology, topography, and climate and spectral data
from remotely-sensed imagery or products.
## Statistical models are built by finding relationships between vegetation data and the
available digital data and then applied across given digital landscapes.
Integrate forest inventory data . .
. .with available digital data . .
• Satellite imagery
• DEMs
• Soils
. .and generate maps of forest attributes.
y = f(x1 + x2 + xn)
. .using flexible statistical
models and GIS tools . .
4
Predictor Data
## Digital environmental data and remotely-sensed products are becoming increasingly
more available at many different scales.
Remotely-sensed data and derived products
•
•
Snapshot of what is on the ground
Caution: represents current vegetation distributions based on reflectance patterns, but does not
explain potential occurrences of vegetation.
Topography variables (elevation, aspect, slope)
•
•
No direct physiological influence on vegetation
Caution: these variables are local to the model domain, be careful when extrapolating over space
or time.
Climate variables (temperature, precipitation)
•
•
•
Directly related to physiological responses of vegetation
Surrogates for topographical variables
Caution: these variables are better for extrapolating over space and time, but recognize other
limitations, such as current location and dispersal range, and species competition may also
change through time..
Other surrogates (geology, soils, soil moisture availability, solar radiation)
•
•
Resources directly used by plants for growth
Caution: similar to climate variables
5
Predictive Model Types
Model Objective:
Find a relationship between vegetation data and predictor data for prediction purposes.
The model, in its simplest form, looks like the following:
y ~ f (x1 + x2 + x3 + …) + ε,
where y is a response variable, the x's are predictor variables, and ε is the associated error.
There are many different types of models
Parametric models – make assumptions about the underlying distribution of the data.
• Maximum likelihood classification
• Discriminant analysis
• General linear models (ex. linear regression)
• ...
Nonparametric models – make no assumptions about the underlying data distributions.
• Machine-learning models
• Artificial neural networks
• ...
6
Predictive Model Types
Parametric Models
Parametric models – make assumptions about the underlying distribution of the data
Assumptions:
• The shape of the distribution of the underlying population is bell-shaped (normal).
• Errors are independent.
• Errors are normally distributed with mean of 0 and a constant variance.
• Easy to interpret
• High power with low sample sizes
• If sample data are not from a normally-distributed population, may lead to incorrect
conclusions.
Examples:
• Maximum likelihood classification
• Discriminant Analysis
• Linear Regression
• Multiple linear regression
• Generalized linear models (parametric/nonparametric)
7
Predictive Model Types
Parametric Models – Regression Example
Previous example: Using regression to fill in missing values.
## Import data and subset for sp19 only
path <- "C:/Peru/Rworkshop" # Where workshop material is
setwd(path)
sp19 <- tree[( tree\$SPCD==19 & !is.na(tree\$DIA) ),]
# Start off with a scatter plot
par(mfrow=c(1,2))
plot(sp19\$DIA, sp19\$HT, xlab ="Diameter", ylab ="Height",
main = "Subalpine Fir")
abline(lm(sp19\$HT ~ sp19\$DIA))
# We saw some heteroscedasticity (unequal variance), and transformed data to log scale.
plot(log(sp19\$DIA), log(sp19\$HT), xlab ="Diameter", ylab ="Height",
main = "Subalpine Fir, log scale")
abline(lm(log(sp19\$HT) ~ log(sp19\$DIA)))
# We looked at summary of models and saw lower residual error and higher R2 values.
r.mod <- lm(HT~DIA,data=sp19)
summary(r.mod)
r.mod.log.ht.dia <- lm(log(HT)~log(DIA),data=sp19)
summary(r.mod.log.ht.dia)
8
Predictive Model Types
Parametric Models – Regression Example
Previous example cont..
# Then we looked at the residuals versus the fitted values and normal QQ-plots.
par(mfrow=c(2,2))
plot(r.mod\$fitted,r.mod.log.ht.dia\$residuals, xlab="Fitted",ylab="Residuals",
main ="Fitted versus Residuals")
abline(h=0)
qqnorm(r.mod\$residuals, main="Normal Q-Q Plot")
qqline(r.mod\$residuals)
plot(r.mod.log.ht.dia\$fitted,r.mod.log.ht.dia\$residuals,
xlab="Fitted",ylab="Residuals",
main ="Fitted versus Residuals, log scale")
abline(h=0)
qqnorm(r.mod.log.ht.dia\$residuals, main="Normal Q-Q Plot, log scale")
qqline(r.mod.log.ht.dia\$residuals)
# NOTE:
# Using transformations, we were able to work with data that was nonlinear, with unequal
# variance structure using a parametric model.
9
Predictive Model Types
Nonparametric Models
Nonparametric models – make no assumptions about the underlying distribution of
the data.
Note: Vegetation data typically are not normally distributed across the landscape,
therefore it is most often better to use a nonparametric model.
# If sample data are not from a normally-distributed population, using a parametric may
# Need larger sample sizes to have the same power as parametric statistics
# Harder to interpret
# Can overfit data
## Examples:
# Classification and regression trees (i.e. CART)
# Artificial neural networks (ANN)
# Multivariate adaptive regression splines (MARS)
# Ensemble modeling (i.e. Random Forests)
10
Random Forests
Random Forests (Breiman, 2001)
Generates a series of classification
and regression tree models..
.. sampling, with replacement, from
training data (bootstrap)
.. selecting predictor variables at
random for each node
.. outputting the class that most
frequently results
.. and calculating an out-of-bag error
estimate
.. and measuring variable importance
through permutation
randomForest – Liaw & Wiener
ModelMap – Freeman & Frescino
Modeling Example
12
Model Example
Objective:
Using Random Forests to find relationships between forest inventory data and six
spatial predictor layers. The models will be used to make predictions across a
continuous, pixel-based surface.
80 TM B3
9200’ Elev
95° Aspect
20% Slope
60 TM B3
8000’ Elev
160° Aspect
15% Slope
120 TM B3
10500’ Elev
10° Aspect
12 % Slope
Landsat TM
Elevation
Aspect
Slope
10 %
cover
35 %
cover
Extract data from each
layer at each sample plot
location.
Landsat TM
Elevation
Aspect
Slope
80 %
cover
Prediction
(% Tree crown
cover)
Generate spatially explicit maps of
forest attributes based on cell by cell
predictions.
Model Example
The model form:
y ~ f (x1 + x2 + x3 + x4 + x5 + x6) + ε,
where y is forest inventory data, and the x's are the predictor variables listed below,
including satellite spectral data and topographic variables.
## For this example,
# We will look at 2 responses:
# Presence of aspen
Binary response of 0 and 1 values (1=presence)
# Total carbon
Continuous response
# With 5 predictor variables:
# Landsat Thematic Mapper, band 5
# Landsat Thematic Mapper, NDVI
# Classified forest/nonforest map
# Elevation
# Slope
# Aspect
30-m spectral data, band 5
30-m spectral data, NDVI
250-m classified MODIS, resampled to 30 m
30-m DEM
30-m DEM –derived
30-m DEM - derived
14
Study Area
Utah, USA
15
Study Area
Model data set:
Uinta Mountains, Utah,USA
Highest East-West oriented
mountain range in the
contiguous U.S. - up to
13,528 ft (4,123 m)
Apply model to:
High Uinta Wilderness
## Vegetation
5 different life zones:
1. shrub-montane
2. aspen
3. lodgepole pine
4. spruce-fir
5. alpine
16
Modeling Dataset
17
Data for Modeling
library(rgdal)
library(raster)
library(rpart)
library(car)
library(randomForest)
library(PresenceAbsence)
library(ModelMap)
#
#
#
#
#
#
#
GDAL operations for spatial data
Analyzing gridded spatial data
Recursive partitioning and regression trees
For book (An R Companion to Applied Regression)
Generates Random Forest models
Evaluates results of presence-absence models
Generates and applies Random Forest models
18
Response Data
# Forest Inventory data (Model response)
# We have compiled this data before, let's review.
options(scipen=6)
## The plt table contains plot-level data, where there is 1 record (row) for each plot. This
table has the coordinates (fuzzed) of each plot, which we will need later.
dim(plt)
## Total number of plots
## Display first six plot records.
## The tree table contains tree-level data, where there is 1 record (row) for each tree on a
plot. We need to summarize the tree data to plot-level for modeling.
dim(tree)
## Total number of trees
## Display first six tree records.
# First, let's add the species names to the table.
# Merge species names to table using a reference table of species codes
tree <- merge(tree, ref, by="SPCD")
19
Response Data
# Forest Inventory data (Model response)
## We have 2 responses:
# Presence of aspen
# Total carbon
Binary response of 0 and 1 values (1=presence)
Continuous response
## Let's compile presence of aspen and append to plot table.
# First, create a table of counts by species and plot.
spcnt <- table(tree[,c("PLT_CN", "SPNM")])
# For this exercise, we don't care about how many trees per species, we just want presence
# or absence, therefore, we need to change all values greater than 1 to 1.
spcnt[spcnt > 1] <- 1
# We are only interested in aspen presence, so let's join the aspen column to the plot table.
spcntdf <- data.frame(PLT_CN=row.names(spcnt), ASPEN=spcnt[,"aspen"])
plt2 <- merge(plt, spcntdf, by.x="CN", by.y="PLT_CN")
dim(plt)
dim(plt2)
# Notice there are fewer records (plots with no trees)
plt2 <- merge(plt, spcntdf, by.x="CN", by.y="PLT_CN", all.x=TRUE)
dim(plt2)
20
Response Data
## Forest Inventory data (Model response)
## Now, let's compile total carbon by plot and append to plot table.
# First, create a table of counts by plot.
pcarbon <- aggregate(tree\$CARBON_AG, list(tree\$PLT_CN), sum)
names(pcarbon) <- c("PLT_CN", "CARBON_AG")
# Carbon is stored in FIA database with units of pounds. Let's add a new variable,
CARBON_KG, with conversion to kg.
pcarbon\$CARBON_KG <- round(pcarbon\$CARBON_AG * 0.453592)
# Now we can join this column to the plot table (plt2).
plt2 <- merge(plt2, pcarbon, by.x="CN", by.y="PLT_CN", all.x=TRUE)
dim(plt2)
# Change NA values to 0 values.
plt2[is.na(plt2[,"ASPEN"]), "ASPEN"] <- 0
plt2[is.na(plt2[,"CARBON_KG"]), "CARBON_KG"] <- 0
plt2\$CARBON_AG <- NULL
21
Response Data
# We need to extract data from spatial layers, so let's convert the plot table to a
SpatialPoints object in R.
## We know the projection information, so we can add it to the SpatialPoints object.
prj4str <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"
ptshp <- SpatialPointsDataFrame(plt[,c("LON","LAT")], plt,
proj4string = CRS(prj4str))
## Display the points
plot(ptshp)
22
Predictor Data
## Predictor variables:
# Landsat Thematic Mapper, band 5
# Landsat Thematic Mapper, NDVI
# Classified forest/nonforest map
# Elevation
# Slope
# Aspect
30-m spectral data, band 5, resampled to 90 m
30-m spectral data, NDVI, resampled to 90 m
250-m classified MODIS, resampled to 90 m
30-m DEM, resampled to 90 m
90-m DEM – derived in a following slides
90-m DEM – derived in a following slides
## Set file names
b5fn <- "SpatialData/uintaN_TMB5.img"
ndvifn <- "SpatialData/uintaN_TMndvi.img"
fnffn <- "SpatialData/uintaN_fnfrcl.img"
elevfn <- "SpatialData/uintaN_elevm.img"
#
#
#
#
Landsat TM–Band5
Landsat TM–NDVI
Forest type map (reclassed)
Elevation (meters)
Note: If you don't have uintaN_fnfrcl.img, follow steps on last slide of this presentation (Appendix 1)
## Check rasters
rastfnlst <- c(b5fn, ndvifn, fnffn, elevfn)
rastfnlst
sapply(rastfnlst, raster)
23
Predictor Data
# TM Band 5 (uintaN_TMB5.img)
# TM NDVI (uintaN_TMndvi.img)
24
Predictor Data
# Forest/Nonforest map (uintaN_fnf.img)
# DEM (uintaN_elevm.img)
25
Predictor Data
## Now, let's generate slope from DEM. Save it to your SpatialData folder.
help(terrain)
help(writeRaster)
slpfn <- "SpatialData/uintaN_slp.img"
slope <- terrain(raster(elevfn), opt=c('slope'), unit='degrees',
filename=slpfn, datatype='INT1U', overwrite=TRUE)
plot(slope, col=topo.colors(6))
# Add slope file name to rastfnlst
rastfnlst <- c(rastfnlst, slpfn)
rastfnlst
26
Predictor Data
## We can also generate aspect from DEM. Save it to your SpatialData folder.
help(terrain)
## This is an intermediate step, so we are not going to save it.
aspectr
# Note: Make sure to use radians, not degrees
plot(aspectr, col=terrain.colors(6))
27
Predictor Data
## Aspect is a circular variable. There are a couple ways to deal with this:
## 1. Convert the values to a categorical variable (ex. North, South, West, East)
aspectd <- round(aspectr * 180/pi)
aspectd
## Now, create a look-up table of reclass values.
help(reclassify)
frommat <- matrix(c(0,45, 45,135, 135,225, 225,315, 315,361), 5, 2)
frommat
frommat <- matrix(c(0,45, 45,135, 135,225, 225,315, 315,361), 5, 2, byrow=TRUE)
frommat
tovect <- c(1, 2, 3, 4, 1)
rclmat <- cbind(frommat, tovect)
rclmat
## Reclassify raster to new values.
aspcl <- reclassify(x=aspectd, rclmat, include.lowest=TRUE)
aspcl
unique(aspcl)
bks <- c(0,sort(unique(aspcl)))
# Break points
cols <- c("dark green", "wheat", "yellow", "blue")
# Colors
labs <- c("North", "East", "South", "West")
# Labels
lab.pts <- bks[-1]-diff(bks)/2
# Label position
plot(aspcl, col=cols, axis.args=list(at=lab.pts, labels=labs), breaks=bks)
28
Predictor Data
## 2. Convert to a linear variable (ex. solar radiation index; Roberts and Cooper 1989)
aspval <- (1 + cos(aspectr+30))/2 ## Roberts and Cooper 1989
aspval
plot(aspval)
## Let's multiply by 100 and round so it will be an integer (less memory)
aspval <- round(aspval * 100)
aspval
plot(aspval)
# Save this layer to file
aspvalfn <- "SpatialData/uintaN_aspval.img"
writeRaster(aspval, filename=aspvalfn, datatype='INT1U', overwrite=TRUE)
rastfnlst <- c(rastfnlst, aspvalfn)
## Converts aspect into solar radiation equivalents, with a correction of 30 degrees to reflect
## the relative heat of the atmosphere at the time the peak radiation is received.
## Max value is 1.0, occurring at 30 degrees aspect; min value is 0, at 210 degrees aspect.
## Roberts, D.W., and S. V. Cooper. 1989. Concepts and techniques in vegetation mapping. In Land
classifications based on vegetation: applications for resource management. D. Ferguson, P. Morgan, and F.
D. Johnson, editors. USDA Forest Service General Technical Report INT-257, Ogden, Utah, USA.
29
Discussion - Scale of Predictors
## Tools in R for handling scale issues. See help on functions for further details.
# focal - applies a moving window function across pixels without changing the resolution.
# Note: It works but takes a lot of time on large rasters.
# aggregate - aggregates pixels to lower resolution.
# resample – resamples pixels to match extent or pixel size of another raster.
Predictor Data Extraction
## The next step is to extract the values of each raster at each sample point location.
80 TM B3
9200’ Elev
95° Aspect
20% Slope
60 TM B3
8000’ Elev
160° Aspect
15% Slope
120 TM B3
10500’ Elev
10° Aspect
12 % Slope
Landsat TM
Elevation
Aspect
Slope
10 %
cover
35 %
cover
Extract data from each
layer at each sample plot
location.
80 %
cover
Predictor Data Extraction
## We need to check the projections of the rasters. If the projections are different,
## reproject the points to the projection of the rasters, it is much faster.
## We will use the plt2 table with LON/LAT coordinates and the response data attached.
## We know the LON/LAT coordinates have the following projection:
prj4str <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"
# Check projections of each raster..
sapply(rastfnlst, function(x){ projection(raster(x)) })
## Reproject SpatialPoints object to match raster projections.
help(project)
rast.prj <- projection(raster(rastfnlst[1]))
xy <- cbind(plt\$LON, plt\$LAT)
xyprj <- project(xy, proj=rast.prj)
32
Predictor Data Extraction
## Extract values (raster package)
help(extract)
# Let's extract values from 1 layer.
tmp <- extract(raster(elevfn), xyprj)
# Now, let's create a function to extract, so we can extract from all the rasters at the same time.
extract.fun <- function(rast, xy){ extract(raster(rast), xy) }
# Now, apply this function to the vector list of raster file names.
rastext <- sapply(rastfnlst, extract.fun, xyprj)
# Look at the output and check the class.
class(rastext)
33
Predictor Data Extraction
## Extract values (raster package) cont.. change names
# Let's make the column names shorter.
colnames(rastext)
# Use the rastfnlst vector of file names to get new column names.
# First, get the base name of each raster, without the extension.
cnames <- unlist(strsplit(basename(rastfnlst), ".img"))
cnames
# We could stop here, but let's make the names even shorter and remove
# 'uintaN_' from each name.
cnames2 <- substr(cnames, 8, nchar(cnames))
cnames2
# Now, add names to matrix. Because the output is a matrix, we will use colnames.
colnames(rastext) <- cnames2
34
Predictor Data Extraction
# Now, let's append this matrix to the plot table with the response data (plt2).
# We just want the response variables, so let's extract these columns along with the unique
identifier of the table (CN, ASPEN, CARBON_KG).
modeldat <- cbind(plt2[,c("CN", "ASPEN", "CARBON_KG")], rastext)
# Check if this is a data frame
is.data.frame(modeldat)
dim(modeldat)
# Let's also append the projected xy coordinates for overlaying with raster layers.
modeldat <- cbind(xyprj, modeldat)
colnames(modeldat)[1:2] <- c("X", "Y")
35
Data Exploration
36
Model Data Exploration
## What to look for:
# NA values
# Outliers
# Correlations
# Non-normal distributions
# Changes in variability
# Clustering
# Non-linear data structures
37
Model Data Exploration
Summary statistics
## Summary statistics
str(modeldat)
summary(modeldat)
dim(modeldat)
# We need to convert categorical variables to factors
modeldat\$ASPEN <- factor(modeldat\$ASPEN)
modeldat\$fnfrcl <- factor(modeldat\$fnfrcl)
## Now, display summary statistics again and notice changes for ASPEN and fnfrcl
str(modeldat)
summary(modeldat)
38
Model Data Exploration
NA Values
## Check for NA values.
modeldat[!complete.cases(modeldat),]
modeldat.NA <- modeldat[!complete.cases(modeldat),]
dim(modeldat.NA)
modeldat.NA
# We can overlay plots with NA values on raster.
plot(raster(aspvalfn))
points(modeldat.NA, pch=20)
# Most R model functions will handle or remove NA values, but for this example, let's
remove the plots with NA values from our dataset now.
modeldat <- modeldat[complete.cases(modeldat),]
dim(modeldat)
39
Model Data Exploration
attach data frame
## Attaching a data frame. This is useful if you are exclusively working with 1 data frame.
## Caution: data frame variable names must be unique to data frame.
## Let's save modeldat object and clean up before we attach the data frame.
save(modeldat, file="Outfolder/modeldat.Rda")
# Now, remove all objects except modeldat.
ls()[ls() != "modeldat"]
rm(list = ls()[ls() != "modeldat"])
ls()
# .. and attach modeldat
attach(modeldat)
ASPEN
# Display column vector without using \$
# Notes:
# To load the saved model object
# Make sure to detach data frame when done using it..
# detach(modeldat)
40
Model Data Exploration
Predictors
## Check for outliers and correlations among predictors.
## Let's look at an example using 2 predictors (elevm, TMB5)
preds <- cbind(elevm, TMB5)
## Correlation between predictor variables, to determine strength of relationship
cor(preds)
round(cor(preds),4)
## Scatterplots
plot(preds)
abline(lm(TMB5 ~ elevm))
lines(loess.smooth(elevm, TMB5), col="red")
## Another way
library(car)
scatterplot(TMB5 ~ elevm)
help(scatterplot)
41
Model Data Exploration
Predictors
## Correlation cont..
# We can now see a non-linear trend in the data where, TMB5 spectral values decrease as
elevations rise up to around 3000 meters, then begins to increase in value as elevations
continue to rise.
# This suggests a non-parametric relationship.
help(cor)
# The default uses pearson correlation coefficient, but spearman may be a better choice for
nonlinear data structures.
round(cor(preds, method="spearman"),4)
## Let's look at all predictors at once
names(modeldat)
preds <- modeldat[-c(1:5)]
## Correlation between predictor variables (to minimize number of predictors)
cor(preds)
str(preds)
## Note Error : 'x' must be numeric (remove factor variable from analysis)
round(cor(preds[,-3], method="spearman"),4)
42
Model Data Exploration
Predictors
## Scatterplots
pairs(preds)
help(pairs) #select Scatterplot Matrices from graphics package
## In help doc, scroll down to Examples and copy and paste the 2 functions:
# panel.hist
# panel.cor
pairs(preds, lower.panel = panel.smooth, upper.panel = panel.cor)
## Change cor function to spearman within panel.cor
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{ usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, method="spearman"))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(preds, lower.panel = panel.smooth, upper.panel = panel.cor)
## Another way
scatterplotMatrix(preds[,-3])
43
Model Data Exploration
Predictors
## Check for outliers
plot(TMB5, TMndvi)
identify(TMB5, TMndvi)
# Click on outliers and press esc key to escape
## Display the outliers
modeldat[c(58,60),]
plot(raster("SpatialData/uintaN_TMB5.img"))
points(modeldat[c(58,60),], pch=20)
## Let's remove these outliers from our dataset
modeldat <- modeldat[-c(58,60),]
dim(modeldat)
length(TMB5)
## Detach and reattach data frame
detach(modeldat)
attach(modeldat)
dim(modeldat)
length(TMB5)
save(modeldat, file="Outfolder/modeldat.Rda") ## Resave modeldat object
## Check scatterplot again
plot(TMB5, TMndvi)
## ..and for all predictors
preds <- modeldat[-c(1:5)]
pairs(preds, lower.panel = panel.smooth, upper.panel = panel.cor)
44
Model Data Exploration
Predictors
## We can also look at the distribution of each predictor across its range of values and its
diversion from normality using the cumulative density function we looked at earlier.
## This is only meaningful for continuous predictors
plot(sort(elevm))
lines(c(1, length(elevm)), range(elevm), col=2)
## Now, let's make a little function to look at the rest of the continuous predictors
pdistn <- function(x){
plot(sort(x))
lines(c(1, length(x)), range(x), col=2)
}
# Press enter to go to next display
pdistn(elevm)
pdistn(slp)
pdistn(TMB5)
pdistn(TMndvi)
pdistn(aspval)
## For categorical predictors, let's use table to look at the distribution of samples by value.
table(fnfrcl)
45
Model Data Exploration
Response
## Let's look at the distribution and amount of variability of the sample response data.
## Check for normality (bell-shaped curve)
## Histograms and density functions
par(mfrow=c(2,1))
hist(CARBON_KG)
hist(CARBON_KG, breaks=5)
## Overlay density function (smoothed) on the histogram with continuous response.
hist(CARBON_KG, breaks=10, prob=TRUE)
lines(density(CARBON_KG))
## Data with lots of 0 values tend to have this shape. This is not a normal distribution
## Therefore, using a nonparametric model, such as Random Forests is a good idea.
## Let's look at the distribution if there were no 0 values (plots without trees)
hist(CARBON_KG[CARBON_KG>0], breaks=10, prob=TRUE)
lines(density(CARBON_KG[CARBON_KG>0]))
## The shape of the distribution is similar, with much more higher values than lower values.
par(mfrow=c(1,1))
46
Model Data Exploration
Data distributions – Binomial response
## Let's look at the distribution of the sample response data with a couple predictors.
## We will start with presence/absence of aspen. We know this is a binomial distribution
with 0 and 1 values.
## Let's explore a little more.
## Plot elevation as a function of ASPEN, bar representing median value.
boxplot(elevm ~ ASPEN, main="Aspen Presence/Absence")
boxplot(elevm ~ ASPEN, main="Aspen Presence/Absence", names=c("Absence",
"Presence"))
# Note: the bold line is median value
## Add points with mean values
means <- tapply(elevm, ASPEN, mean)
points(means, col="red", pch=18)
# Note: overall, presence of aspen tends to be at the lower elevations of the sample plots.
# The distributions are slightly skewed with the median elevm values higher than mean
elevm values.
47
Model Data Exploration
Data distributions – Binomial response
## Other ways to explore relationships between the response and predictors.
## We can look at the differences of presence vs absence with 2 predictors.
plot(elevm, ASPEN)
## We can look at the differences of presence vs absence with 2 predictors (elevation and slope).
par(mfrow=c(1,2))
plot(elevm[ASPEN==1], slp[ASPEN==1])
plot(elevm[ASPEN==0], slp[ASPEN==0])
## Now, let's make it more meaningful by adding labels and using the same scale for each.
xlab <- "Elevation"
ylab <- "Slope"
xlim <- c(2000, 4000)
ylim <- c(0,40)
plot(elevm[ASPEN==1], slp[ASPEN==1], xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim,
main="Aspen Present")
plot(elevm[ASPEN==0], slp[ASPEN==0], xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim,
main="No Aspen")
par(mfrow=c(1,1))
48
Model Data Exploration
Data distributions – Binomial response cont.
## Other ways to explore relationships between the response and predictors (1 plot)
## We can also color the points based on a factor (ASPEN)
plot(slp, elevm, col=ASPEN, pch=20, xlab="Slope", ylab="Elevation")
## Add a legend (using default colors)
legend(x=35, y=2200, legend=levels(ASPEN), col=1:length(ASPEN), pch=20)
help(legend)
## Now, do it again using your own color choice
palette(c("blue", "green"))
# Change color palette first
plot(slp, elevm, col=ASPEN, pch=20, xlab="Slope", ylab="Elevation")
legend(x=35, y=2200, legend=levels(ASPEN), col=1:length(ASPEN), pch=20)
palette("default")
# Change color palette back to default colors
## Or..
plot(slp, elevm, col=c("red", "blue"), pch=20, xlab="Slope", ylab="Elevation")
legend(x=35, y=2200, legend=levels(ASPEN), col=c("red", "blue"), pch=20)
## Another way:
scatterplot(slp ~ elevm|ASPEN, data=modeldat)
## For categorical predictors, use table function again
table(ASPEN, fnfrcl)
49
Model Data Exploration
Data distributions – Continuous response
## Other ways to explore relationships between the response and predictors.
## We can look at relationship between elevation and CARBON_KG
plot(elevm, CARBON_KG, xlab="Elevation", ylab="Carbon(kg)")
## Without 0 values
plot(elevm[CARBON_KG>0], CARBON_KG[CARBON_KG>0], xlab="Elevation",
ylab="Carbon (kg)")
## Add regression and smoother lines.
par(mfrow=c(2,1))
## With 0 values
plot(elevm, CARBON_KG, xlab="Elevation", ylab="Carbon(kg)")
line.lm <- lm(CARBON_KG ~ elevm)
abline(line.lm, col="red")
line.sm <- lowess(elevm, CARBON_KG)
lines(line.sm, col="blue")
## Without 0 values
plot(elevm[CARBON_KG>0], CARBON_KG[CARBON_KG>0],xlab="Elevation",ylab="Carbon(kg)")
line.lm.w0 <- lm(CARBON_KG[CARBON_KG>0] ~ elevm[CARBON_KG>0])
abline(line.lm.w0, col="red")
line.sm <- lowess(elevm[CARBON_KG>0], CARBON_KG[CARBON_KG>0])
lines(line.sm, col="blue")
par(mfrow=c(1,1))
50
Model Generation
51
Model Generation
y ~ f (x1 + x2 + x3 + x4 + x5 + x6) + ε
80 TM B3
9200’ Elev
95° Aspect
20% Slope
60 TM B3
8000’ Elev
160° Aspect
15% Slope
120 TM B3
10500’ Elev
10° Aspect
12 % Slope
Landsat TM
Elevation
Aspect
Slope
10 %
cover
35 %
cover
Extract data from each
layer at each sample plot
location.
80 %
cover
Tree-based Models
## Strengths of Tree-based models
# Easy to interpret
# Adaptable for handling missing values
# Handles correlated predictor variables
# Predictor variable interactions are automatically included
# Handles categorical or continuous predictor variables
## Weaknesses of Tree-based models, such as Random Forests
# Optimization is based on each split, not necessary the overall tree model
# Continuous predictors are treated as categorical, thus inefficient
# Nonparametric, thus loses some the power of parametric statistics
# Tendency to overfit data
53
Tree-based Models
## What does this mean, in simple form..
# Using our example dataset,
library(rpart)
## Classification tree
asp.tree <- rpart(ASPEN ~ TMB5 + TMndvi + fnfrcl + elevm + slp + aspval,
data=modeldat, method="class")
plot(asp.tree)
text(asp.tree, cex=0.75)
## Regression tree
carb.tree <- rpart(CARBON_KG ~ TMB5 + TMndvi + fnfrcl + elevm + slp + aspval,
data=modeldat)
plot(carb.tree)
text(carb.tree, cex=0.75)
54
Random Forests
Random Forests (Breiman, 2001)
Generates a series of classification and
regression tree models..
.. sampling, with replacement, from
training data (bootstrap)
.. selecting predictor variables at random
for each node
.. outputting the class that most
frequently results
.. and calculating an out-of-bag error
estimate
.. and measuring variable importance
through permutation
randomForest – Liaw & Wiener
ModelMap – Freeman & Frescino
Breiman, L. (2001). Random forests. Machine Learning J. 45 5- 32.
Liaw, A.; Wiener, M. 2002. Classification and Regression by randomForest. ISSN 2/3:18-22.
Freeman, E.A. et al. 2012. ModelMap: an R package for Model Creation and Map Production. CRAN R vignette.
Random Forests Model
## Strengths Random Forests (Breiman 2001)
# Bootstrap sample – A random selection of plots used to construct one tree.
# Boosting – Successive trees are constructed using information from previous trees, and a
weighted vote is used for prediction.
# Bagging – Each tree is independently constructed, where the majority (or average) vote is
used for prediction. Breiman's Random Forests uses bagging.
# Predictor selection - Each node is split using a randomly selected subset of predictors. In
standard trees, all variables are used. This difference is more robust against overfitting.
# Two main parameters
# 1. Number of trees (bootstrap samples) to generate (ntree)
# 2. Number of variables in the random subset of predictors at each node (mtry)
# Predictions – Aggregate predictions of n trees
# For categorical response (Classification trees) – majority votes
# For continuous response (Regression trees) – average of regression
# Error rate – based on training data ('out-of-bag', or OOB)
# For each tree (bootstrap sample)
# For continuous response (Regression trees) – average of regression
56
Random Forests Model
classification tree
## Now, let's use the randomForests package – Classification tree
library(randomForest)
help(randomForest)
## Let's try with ASPEN binary, categorical response (presence/absence)
set.seed(66)
asp.mod <- randomForest(ASPEN ~ TMB5 + TMndvi + fnfrcl + elevm + slp + aspval,
data=modeldat, importance = TRUE)
## Default parameters:
# ntree = 500
# mtry = sqrt(p)
# nodesize = 1
# replace = TRUE
# Number of trees
# Number of predictors (p) randomly sampled at each node
# Minimum size of terminal nodes
# Bootstrap samples are selected with replacement
## Look at results
asp.mod
summary(asp.mod)
names(asp.mod)
57
Random Forests Model
Classification tree
## Classification tree - Output
names(asp.mod)
err <- asp.mod\$err.rate
tail(err)
mat <- asp.mod\$confusion
mat
# Out-Of-Bag (OOB) error rate (of i-th element)
# Confusion matrix
58
Random Forests Model
Classification tree
## Classification tree - Output
# Plot the number of trees by the error rate
plot(1:500, err[,"OOB"], xlab="Number of trees", ylab="Error rate")
# Note: how many trees needed to stabilize prediction
## Calculate the percent correctly classified from confusion (error) matrix
mat
pcc <- sum(diag(mat[,1:2]))/sum(mat) * 100
pcc
pcc <- round(pcc, 2)
## Round to nearest 2 decimals
pcc
library(PresenceAbsence)
pcc(mat[,1:2], st.dev=TRUE)
Kappa(mat[,1:2], st.dev=TRUE)
## The Kappa statistic summarizes all the available information in the confusion matrix.
## Kappa measures the proportion of correctly classified units after accounting for the probability of
chance agreement.
59
Random Forests Model
## Now, let's use the randomForests package – regression tree
## Now, let's try with the continuous, CARBON_KG response
set.seed(66)
carb.mod <- randomForest(CARBON_KG ~ TMB5 + TMndvi + elevm + slp + aspval,
data=modeldat, importance = TRUE)
## Default parameters:
# ntree = 500
# mtry = p/3
# nodesize = 5
# replace = TRUE
# Number of trees
# Number of predictors (p) randomly sampled at each node
# Minimum size of terminal nodes
# Bootstrap samples are selected with replacement
## Look at results
carb.mod
summary(carb.mod)
names(carb.mod)
60
Random Forests Model
Regression tree
## Regression tree - Output
names(carb.mod)
mse <- carb.mod\$mse
rsq <- carb.mod\$rsq
# Mean square error (of i-th element)
# Pseudo R-squared (1-mse/Var(y))(of i-th element)
tail(mse)
tail(rsq)
61
Random Forests Model
Regression tree
## Regression tree - Output
# Plot the number of trees by the mse (Mean Square Error)
plot(1:500, mse, xlab="Number of trees", ylab="Mean Square Error rate")
# Note: how many trees needed to stabilize prediction
# Similarly, plot the number of trees by the rsq (R-Squared)
plot(1:500, mse, xlab="Number of trees", ylab="R-Squared")
# Again: how many trees needed to stabilize prediction
62
Random Forests Model
Variable Importance
## Other information from RandomForest model (importance=TRUE)
# Variable importance (Breiman 2002)
# Estimated by how much the prediction error increases when the OOB data for that variable is
permuted while all others are left unchanged.
# randomForests computes different measures of variable importance
# 1. Computed from OOB data, averaged over all trees and normalized by the standard
deviation of the differences.
Classification trees – error rate (Mean Decrease Accuracy)
Regression trees – Mean Square Error (%IncMSE)
# 2. The total decrease in node impurities from splitting on the variable, averaged over all trees.
Classification trees – measured by Gini index (Mean Decrease Gini)
Regression trees – measured by residual sum of squares (IncNodePurity)
63
Random Forests Model
Variable Importance - Classification
## Variable importance – Classification tree
## Get importance table
asp.imp <- abs(asp.mod\$importance)
asp.imp
## Get the number of measures (columns) and number of predictors
ncols <- ncol(asp.imp)
numpred <- nrow(asp.imp)
## Number of measures
## Get number of predictors
## Plot the measures of variable importance for ASPEN presence/absence
par(mfrow=c(2,2))
for(i in 1:ncols){
## Loop thru the different importance measures
ivect <- sort(asp.imp[,i], dec=TRUE)
## Get 1st measure, descending order
iname <- colnames(asp.imp)[i]
## Get name of measure
# Generate histogram plot (type='h') with no x axis (xaxt='n')
plot(ivect, type = "h", main = paste("Measure", iname), xaxt="n",
xlab = "Predictors", ylab = "", ylim=c(0,max(ivect)))
# Add x axis with associated labels
axis(1, at=1:numpred, lab=names(ivect))
}
64
Random Forests Model
Variable Importance - Regression
## Let’s make a function and plot importance values for CARBON_KG model.
plotimp <- function(itab){
ncols <- ncol(itab)
numpred <- nrow(itab)
## Number of measures
## Get number of predictors
## Plot the measures of variable importance
par(mfrow = c(ncols/2,2))
for(i in 1:ncols){ ## Loop thru the different importance measures
ivect <- sort(itab[,i], dec=TRUE) ## Get 1st measure, sorted decreasing
iname <- colnames(itab)[i]
## Get name of measure
# Generate histogram plot (type='h') with no x axis (xaxt='n')
plot(ivect, type = "h", main = paste("Measure", iname), xaxt="n",
xlab = "Predictors", ylab = "", ylim=c(0,max(ivect)))
# Add x axis with associated labels
axis(1, at=1:numpred, lab=names(ivect)) }
}
## Check function with ASPEN model
plotimp(asp.imp)
## Now, run funtion with CARBON_KG model
plotimp(carb.mod\$importance)
65
Random Forests Model
Proximity
## Other information from RandomForest model (proximity=TRUE)
# Measure of internal structure (Proximity measure)
# - The fraction of trees in which each plot falls in the same terminal node.
# - Similarity measure - in theory, similar data points will end up in the same terminal node.
## Let's try adding proximity to CARBON_KG model
set.seed(66)
carb.mod <- randomForest(CARBON_KG ~ TMB5 + TMndvi + elevm + slp + aspval,
data=modeldat, importance = TRUE, proximity = TRUE)
names(carb.mod)
carb.prox <- carb.mod\$proximity
66
Model Prediction
67
Study Area
Model data set:
Uinta Mountains, Utah,USA
Highest East-West oriented
mountain range in the
contiguous U.S. - up to
13,528 ft (4,123 m)
Apply model to:
High Uinta Wilderness
## Vegetation
5 different life zones:
1. shrub-montane
2. aspen
3. lodgepole pine
4. spruce-fir
5. alpine
68
Polygon Data
## Let's import and display the 2 polygon layers as well.
## Set dsn and polygon layer names
dsn <- "SpatialData"
aoinm <- "uintaN_aoi"
wildnm <- "uintaN_wild"
# AOI boundary
# Wilderness boundary (mapping AOI)
## Import polygon shapefiles
## Check projections of all 3 layers to see if we can display them together.
sapply(c(bndpoly, wildpoly), projection)
## Now we can display all 3 layers
par(mfrow=c(1,1))
plot(bndpoly, border="black", lwd=3)
lwd=2)
69
Predictor Data
## Now, we need to clip the raster predictor layers to the extent of the wilderness
polygon
## Set file names
b5fn <- "SpatialData/uintaN_TMB5.img"
ndvifn <- "SpatialData/uintaN_TMndvi.img"
fnffn <- "SpatialData/uintaN_fnfrcl.img"
elevfn <- "SpatialData/uintaN_elevm.img"
slpfn <- "SpatialData/uintaN_slp.img"
(degrees)
aspfn <- "SpatialData/uintaN_aspval.img"
#
#
#
#
#
Landsat TM–Band5
Landsat TM–NDVI
Forest type map
Elevation (meters)
Derived slope
# Derived aspect value
## Check rasters
rastfnlst <- c(b5fn, ndvifn, fnffn, elevfn, slpfn, aspfn)
sapply(rastfnlst, raster)
## Compare projections of rasters with projection of wilderness polygon
projection(wildpoly)
70
Predictor Data
Clip layers
## Clip raster layers
## Let's clip the elevm raster using crop function
help(crop)
elevclip <- crop(raster(elevfn), extent(wildpoly))
## Now display the new raster.
plot(elevclip)
## Note: Notice it didn't clip to boundary, it clipped to extent of boundary
## Add polygon layer to display
## We need to crop further by applying a mask with the polygon layer.
plot(elevclip)
71
Predictor Data
Clip layers
## Create a function to clip the all the layers, saving to the working directory. Use
rastfnlst for raster names and build new names from the list. Return the new name
of the raster.
cliprast <- function(rastfn, poly){
## Create new name from rastfn
rastname <- strsplit(basename(rastfn), ".img")[[1]]
newname <- substr(rastname, 8, nchar(rastname))
newname <- paste("Outfolder/",newname, "_clip.img", sep="")
## Crop raster
rastclip <- crop(raster(rastfn), extent(poly))
## Mask raster and save to working directory with newname
rastclip <- mask(rastclip, poly, filename=newname, overwrite=TRUE)
print(paste("finished clipping",rastname))
flush.console()
return(newname)
}
clipfnlst <- {}
for(rastfn in rastfnlst){
clipfn <- cliprast(rastfn, wildpoly)
clipfnlst <- c(clipfnlst, clipfn)
}
72
Predictor Data
Create Raster Stack
## Check clipped rasters.
sapply(clipfnlst, raster)
## For prediction, the predictor layers must be consistent. Check the following:
# dimensions – all layers must have the same number of rows and columns
# resolution – all layers must have the same resolution
# extent – all layers must have the same extent
# projection – all layers must have same projection
## Create a stack of all the predictor layers
clipstack <- stack(clipfnlst)
clipstack
## Add names to stack layers
stacknames <- unlist(strsplit(basename(clipfnlst), "_clip.img"))
names(clipstack) <- stacknames
clipstack
73
Apply Model
y ~ f (x1 + x2 + x3 + x4 + x5 + x6) + ε
80 TM B3
9200’ Elev
95° Aspect
20% Slope
60 TM B3
8000’ Elev
160° Aspect
15% Slope
10 %
cover
120 TM B3
10500’ Elev
10° Aspect
12 % Slope
Landsat TM
Landsat TM
Elevation
Elevation
Aspect
Aspect
Slope
Slope
80 %
cover
35 %
cover
Extract data from each
layer at each sample plot
location.
Prediction
(% Tree crown
cover)
Generate spatially explicit maps of
forest attributes based on cell by cell
predictions.
Apply Model & Display Map
ASPEN
## Predict across stack pixels.
asp.predict <- predict(clipstack, asp.mod)
asp.predict
plot(asp.predict)
# Plot with color breaks
cols <- c("dark grey", "green")
plot(asp.predict, col=cols, breaks=c(0,0.5,1))
colors()
# Or, a little fancier.. create function to color categorical raster classes
colclasses <- function(rast, cols, labs){
nc <- length(cols)
minval <- cellStats(rast, 'min')
maxval <- cellStats(rast, 'max')
bks <- seq(minval,maxval,length.out=nc+1)
#
#
#
#
Number of classes
minimum value of raster
maximum value of raster
break points
lab.pts <- bks[-1]-diff(bks)/2
# label points
plot(rast, col=cols, axis.args=list(at=lab.pts, labels=labs),
breaks=bks)
}
colclasses(asp.predict, cols, labs=c("Absence", "Presence"))
par()
help(par)
par(omi=c(0,0,0,0.5))
# Changes outside margins
colclasses(asp.predict, cols, labs=c("Absence", "Presence"))
75
Apply Model & Display Map
CARBON_KG
## Predict across stack pixels.
carb.predict <- predict(clipstack, carb.mod)
## Plot with heat color ramp
plot(carb.predict, heat.colors(10))
## Plot with grey scale
plot(carb.predict, col=grey(256:0/256))
## Plot with green scale
my.colors <- colorRampPalette(c("white", "dark green"))
plot(carb.predict, col=my.colors(10))
76
ModelMap
Build Model
library(ModelMap)
help(package=ModelMap)
predList <- stacknames
#predList <- c("elevm", "TMndvi")
## Build random forest model
asp.mod2 <- model.build(
model.type = "RF",
qdata.trainfn = modeldat,
folder = "Outfolder",
predList = predList,
predFactor = "fnfrcl",
response.name = "ASPEN",
response.type = "categorical",
unique.rowname = "CN",
seed = 66)
asp.mod2
ModelMap
Model diagnostics
asp.mod2d <- model.diagnostics(
model.obj = asp.mod2,
qdata.trainfn = modeldat,
folder = "Outfolder",
response.name = "ASPEN",
prediction.type = "OOB",
unique.rowname = "CN")
ModelMap
Make Map
# Add full path to beginning of each name in predList
predPath <- sapply(predList,
function(x){paste("SpatialData/uintaN_", x, ".img", sep="")})
predPath
# Generate rastLUTfn. See help(model.mapmake) for details.
numpreds <- length(predList)
rastLUTfn <- data.frame(matrix(c(predPath, predList, rep(1,numpreds)),
numpreds, 3), stringsAsFactors=FALSE)
rastLUTfn
## Check rasters
sapply(rastfnlst, raster)
## Make Map
a <- model.mapmake(
model.obj = asp.mod2,
folder = "Outfolder",
rastLUTfn = rastLUTfn,
make.img = TRUE,
na.action = "na.omit")
Exercise
## 1. Create a map of presence of lodgepole within the Uinta Wilderness area, using a Random
Forests model and the predictors from modeldat. Hint: Load the plt and tree table again and use
code from slide #20 to create a binary response variable for lodgepole presence.
Also, set a seed to 66 so you can recreate the model.
## 2. How many plots in model data set have presence of lodgepole pine?
## 3. Display 1 scatterplot of the relationship of elevation and NDVI values with lodgepole pine
presence and absence as different colors. Hint: see slide #49. Make sure to label plots and add
a legend.
## 4. What is the variable with the highest importance value based on the Gini index?
## 5. What percentage of the total area is lodgepole pine? How does this compare with the
percentage of aspen?
Appendix 1
## Reclass raster layer to 2 categories
fnf <- raster("SpatialData/uintaN_fnf.img")
## Create raster look-up table
fromvect <- c(0,1,2,3)
tovect <- c(2,1,2,2)
rclmat <- matrix(c(fromvect, tovect), 4, 2)
## Generate raster and save to SpatialData folder
fnfrcl <- reclassify(x=fnf, rclmat, datatype='INT2U',
filename="SpatialData/uintaN_fnfrcl.img", overwrite=TRUE)