Movatterモバイル変換

Modeling and Map Production using Random Forest and Related Stochastic Models

Description

Creates sophisticated models of training data and validates the models with an independent test set, cross validation, or with Out Of Bag (OOB) predictions on the training data. Create graphs and tables of the model validation results. Applies these models to GIS .img files of predictors to create detailed prediction surfaces. Handles large predictor files for map making, by reading in the .img files in chunks, and output to the .txt file the prediction for each data chunk, before reading the next chunk of data.

Details

This package provides a push button approach to complex model building and production mapping. It contains three main functions:model.build,model.diagnostics, andmodel.mapmake.

In addition it contains a simple functionget.test that can be used to randomly divide a training dataset into training and test/validation sets;build.rastLUT that uses GUI prompts to walk a user through the process of setting up a Raster look up table to link predictors from the training data with the rasters used for map contruction;model.explore, for preliminary data exploration; and,model.importance.plot andmodel.interaction.plot for interpreting the effects of individual model predictors.

ModelMap can be run in a traditional R command mode, where all arguments are specified in the function call. However it can also be used in a full push button mode, where you type in the simple command such asmodel.build, and GUI pop-up windows ask questions about the type of model, the file locations of the data, etc...

Random Forest is implemented through therandomForest package withinR. Random Forest is more user friendly than Stochastic Gradient Boosting, as it has fewer parameters to be set by the user, and is less sensitive to tuning of these parameters. A Random Forest model consists of multiple trees that vote on predictions. For each tree a random subset of the training data is used to construct the tree, with the remaining data points used to construct out-of-bag (OOB) error estimates. At each node of the tree a random selection of predictors is chosen to determine the split. The number of predictors used to select the splits is the primary user specified parameter that can affect model performance, and this parameter can be automatically optimized using therandomForest functiontuneRF(). Random Forest will not over fit data, therefore the only penalty of increasing the number of trees is computation time. Random Forest can compute variable importance, an advantage over some "black box" modeling techniques if it is important to understand the ecological relationships underlying a model (Brieman, 2001).

Quantile Regression Forests is implemented through thequantregForest package.

Conditional Forests is implemented with thecforest() function in theparty package. As stated in theparty package, ensembles of conditional inference trees have not yet been extensively tested, so this routine is meant for the expert user only and its current state is rather experimental.

For Presence-Absence data, the packagePresenceAbsence is used for model validation.

For model diagnostics the packagecorrplot is used to plot the correlation between predictor variables.

For map making, theraster is used to read and write.img files.

For interaction plots, thefields package is used to produce image plots.

Author(s)

Author: Elizabeth Freeman and Tracey Frescino

Maintainer: Elizabeth Freeman <elizabeth.a.freeman@usda.gov>

References

Breiman, L. (2001) Random Forests. Machine Learning, 45:5-32.

Elith, J., Leathwick, J. R. and Hastie, T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology. 77:802-813.

Friedman, J.H. (2001). Greedy function approximation: a gradient boosting machine. Ann. Stat., 29(5):1189-1232.

Friedman, J.H. (2002). Stochastic gradient boosting. Comput. Stat. Data An., 38(4):367-378.

Liaw, A. and Wiener, M. (2002). Classification and Regression by randomForest. R News 2(3), 18–22.

N. Meinshausen (2006) "Quantile Regression Forests", Journal of Machine Learning Research 7, 983-999 http://jmlr.csail.mit.edu/papers/v7/

Ridgeway, G., (1999). The state of boosting. Comp. Sci. Stat. 31:172-181

Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis and Torsten Hothorn (2007). Bias in Random Forest variable Importance Measures: Illustrations, Sources and a Solution. BMC Bioinformatics, 8, 25. http://www.biomedcentral.co,/1471-2105/8/25

Carolin Strobl, James Malley and Gerhard Tutz (2009). An Introduction to Recursive Partitioning: Rationale, Application, and Characteristics of Classification and Regression Trees, Bagging, and Random forests. Phsycological Methods, 14(4), 323-348.

Torsten Hothorn, Berthold Lausen, Axel Benner and Martin Radespiel-Troeger (2004). Bagging Survival Trees. Statistics in Medicine, 23(1), 77-91.

Torsten Hothorn, Peter Buhlmann, Sandrine Dudoit, Annette Molinaro and Mark J. ven der Laan (2006a). Survival Ensembles. Biostatistics, 7(3), 355-373.

Torston Hothorn, Kurt Hornik and Achim Zeileis (2006b). Unbiased Recursive Partitioning: A Conditional Inference Framework. JOurnal of Computational and Graphical Statistics, 15(3), 651-674. Preprint available from http://statmath.wu-wein.ac.at/~zeileis/papers/Hothorn+Hornik+Zeileis-2006.pdf

Internal ModelMap Functions

Description

Internal Subfunctions

Details

These are not to be called by the user.

Author(s)

Elizabeth Freeman

Build a raster Look-UP-Table for training dataset

Description

GUI prompts will help the user build a Look-Up-Table to associated predictor variable with their corresponding spatial rasters.

Usage

build.rastLUT(imageList=NULL,predList=NULL,qdata.trainfn=NULL,rastLUTfn=NULL,folder=NULL)

Arguments

Details

This function helps the user create a raster Look-Up-Table to be used later bymodel.mapmake(). Currently this function only works in a Windows environment.

First, if"folder" is not given, the user selects the output folder for the Look-UP-Table.

Second, if"predList" or"qdatatrainfn" are not given, the user selects the file containing the training data. The header of the file is used to generate a selection list of possible predictor variables.

Third, if"imageList" is not provided, the user selects the rasters.

Finally, the function steps through each band of each raster, and the user selects the appropriate predictor.

Value

Returns a data frame containing the raster Look-Up-Table. Also Writes a.csv file containing the raster Look-Up-Table.

Author(s)

Elizabeth Freeman

Examples

## Not run: folder<-system.file("extdata", "helpexamples", package = "ModelMap")qdata.trainfn = paste(folder,"/DATATRAIN.csv",sep="")#build.rastLUT(qdata.trainfn=qdata.trainfn,#folder=folder)## End(Not run) # end dontrun

colors to transparent colors

Description

transform color names to transparent versions of rgb color codes

Usage

col2trans(col.names,alpha=0.5)

Arguments

Details

Translates a vector of color names to a vector of transparent rgb color codes. Color names must be from names given bycolors.

Value

Outputs a vector of transparent color codes.

Author(s)

Elizabeth Freeman

Examples

col.names=c("blue","violetred4","thistle3","yellowgreen")col2trans(col.names,alpha=.2)###to see effect of alpha###alpha<-(0:10)/10colmat<-matrix(1:(length(alpha)*length(col.names)),nrow=length(alpha),ncol=length(col.names),byrow=TRUE)color.codes<-vector("character",0)for(i in 1:length(alpha)){color.codes<-c(color.codes,col2trans(col.names,alpha=alpha[i]))}#make plot#plot(c(0,1),c(0,1),type="n",xlab="alpha",ylab="color name",yaxt="n",xaxs="i",yaxs="i")abline(h=(0:100)/100)image(z=colmat,x=(0:length(alpha))/length(alpha),y=(0:length(col.names))/length(col.names),col=color.codes,add=TRUE)op<-par(xpd=TRUE)text(col.names,x=-.08,y=(1:length(col.names)-.5)/length(col.names),srt=90)par(op)

Randomly Divide Data into Training and Test Sets

Description

Uses random selection to split a dataset into training and test data sets

Usage

get.test(proportion.test, qdatafn = NULL, seed = NULL, folder=NULL, qdata.trainfn = paste(strsplit(qdatafn, split = ".csv")[[1]], "_train.csv", sep = ""), qdata.testfn = paste(strsplit(qdatafn, split = ".csv")[[1]], "_test.csv", sep = ""))

Arguments

Details

This function should be run once, before starting analysis to create training and test sets. If the cross validation option is to be used with RF or SGB models, or if the OOB option is to be used for RF models, then this step is unnecessary.

Value

Outputs a training data file and test data file. Unlessqdata.trainfn orqdata.testfn are specified, the output will be located infolder. The output will have the same rows and columns as the original data.

Author(s)

Elizabeth Freeman

Examples

## Not run: qdatafn<-system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")qdata<-read.table(file=qdatafn,sep=",",header=TRUE,check.names=FALSE)get.test(proportion.test=0.2,qdatafn=qdatafn,seed=42,folder=getwd(),qdata.trainfn="example.train.csv",qdata.testfn="example.test.csv")## End(Not run) # end dontrun

Model Building

Description

Create sophisticated models using Random Forest, Quantile Regression Forests, Conditional Forests, or Stochastic Gradient Boosting from training data

Usage

 model.build(model.type = NULL, qdata.trainfn = NULL, folder = NULL,MODELfn = NULL, predList = NULL, predFactor = FALSE, response.name = NULL,response.type = NULL, unique.rowname = NULL, seed = NULL, na.action = NULL,keep.data = TRUE, ntree = switch(model.type,RF=500,QRF=1000,CF=500,500),mtry = switch(model.type,RF=NULL,QRF=ceiling(length(predList)/3),CF = min(5,length(predList)-1),NULL), replace = TRUE, strata = NULL,sampsize = NULL, proximity = FALSE, importance=FALSE, quantiles=c(0.1,0.5,0.9), subset = NULL, weights = NULL,controls = NULL, xtrafo = NULL, ytrafo = NULL, scores = NULL)

Arguments

Details

This package provides a push button approach to complex model building and production mapping. It contains three main functions:model.build,model.diagnostics, andmodel.mapmake.

In addition it contains a simple functionget.test that can be used to randomly divide a training dataset into training and test/validation sets;build.rastLUT that uses GUI prompts to walk a user through the process of setting up a Raster look up table to link predictors from the training data with the rasters used for map contruction;model.explore, for preliminary data exploration; and,model.importance.plot andmodel.interaction.plot for interpreting the effects of individual model predictors.

These functions can be run in a traditional R command mode, where all arguments are specified in the function call. However they can also be used in a full push button mode, where you type in, for example, the simple commandmodel.build, and GUI pop up windows will ask questions about the type of model, the file locations of the data, etc...

When running theModelMap package on non-Windows platforms, file names and folders need to be specified in the argument list, but other pushbutton selections are handled by theselect.list() function, which is platform independent.

Binary, categorical, and continuous response models are supported for Random Forest and Conditional Forest. Quantile Random Forest is appropriate for only continuous response models.

Random Forest is implemented through therandomForest package withinR. Random Forest is more user friendly than Stochastic Gradient Boosting, as it has fewer parameters to be set by the user, and is less sensitive to tuning of these parameters. A Random Forest model consists of multiple trees that vote on predictions. For each tree a random subset of the training data is used to construct the tree, with the remaining data points used to construct out-of-bag (OOB) error estimates. At each node of the tree a random selection of predictors is chosen to determine the split. The number of predictors used to select the splits (argumentmtry) is the primary user specified parameter that can affect model performance.

By defaultmtry will be automatically optimized using therandomForest packagetuneRF() function. Note that this is a stochastic process. If there is a chance that models may be combined later with therandomForest packagecombine function then for consistency it is important to provide themtry argument rather that using the default optimization process.

Random Forest will not over fit data, therefore the only penalty of increasing the number of trees is computation time. Random Forest can compute variable importance, an advantage over some "black box" modeling techniques if it is important to understand the ecological relationships underlying a model (Brieman, 2001).

Quantile Regression Forests is implemented through thequantregForest package.

Conditional Forests is implemented with thecforest() function in theparty package. As stated in theparty package, ensembles of conditional inference trees have not yet been extensively tested, so this routine is meant for the expert user only and its current state is rather experimental.

For CF models,ModelMap currently only supports binary, categorical and continuous response models. Also, for some CF model parameters (subset,weights, andscores)ModelMap only provides OOB and independent test set diagnostics, and does not support cross validation diagnostics.

Stochastic gradient boosting is not currently supported byModelMap.

Value

The function will return the model object. Additionally, it will write a text file to disk, in the folder specified byfolder. This file lists the values of each argument as chosen from GUI prompts used for the function call.

Author(s)

Elizabeth Freeman and Tracey Frescino

References

Breiman, L. (2001) Random Forests. Machine Learning, 45:5-32.

Elith, J., Leathwick, J. R. and Hastie, T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology. 77:802-813.

Liaw, A. and Wiener, M. (2002). Classification and Regression by randomForest. R News 2(3), 18–22.

N. Meinshausen (2006) "Quantile Regression Forests", Journal of Machine Learning Research 7, 983-999 http://jmlr.csail.mit.edu/papers/v7/

Ridgeway, G., (1999). The state of boosting. Comp. Sci. Stat. 31:172-181

Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis and Torsten Hothorn (2007). Bias in Random Forest variable Importance Measures: Illustrations, Sources and a Solution. BMC Bioinformatics, 8, 25. http://www.biomedcentral.co,/1471-2105/8/25

Carolin Strobl, James Malley and Gerhard Tutz (2009). An Introduction to Recursive Partitioning: Rationale, Application, and Characteristics of Classification and Regression Trees, Bagging, and Random forests. Phsycological Methods, 14(4), 323-348.

Torsten Hothorn, Berthold Lausen, Axel Benner and Martin Radespiel-Troeger (2004). Bagging Survival Trees. Statistics in Medicine, 23(1), 77-91.

Torsten Hothorn, Peter Buhlmann, Sandrine Dudoit, Annette Molinaro and Mark J. ven der Laan (2006a). Survival Ensembles. Biostatistics, 7(3), 355-373.

Torston Hothorn, Kurt Hornik and Achim Zeileis (2006b). Unbiased Recursive Partitioning: A Conditional Inference Framework. Journal of Computational and Graphical Statistics, 15(3), 651-674. Preprint available from http://statmath.wu-wein.ac.at/~zeileis/papers/Hothorn+Hornik+Zeileis-2006.pdf

See Also

get.test,model.diagnostics,model.mapmake

Examples

## Not run: ######################################################################################################## Run this set up code: ################################################################################################### set seed:seed=38# Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")# Define folder for all output:folder=getwd()#identifier for individual training and test data pointsunique.rowname="ID"######################################################################################### Pick one of the following sets of definitions: ################################################################################################## Continuous Response, Continuous Predictors #############file name:MODELfn="RF_Bio_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="BIO"response.type="continuous"########## binary Response, Continuous Predictors #############file name to store model:MODELfn="RF_CONIFTYP_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="CONIFTYP"# This variable is 1 if a conifer or mixed conifer type is present, # otherwise 0.response.type="binary"########## Continuous Response, Categorical Predictors ############# In this example, NLCD is a categorical predictor.## You must decide what you want to happen if there are categories# present in the data to be predicted (either the validation/test set# or in the image file) that were not present in the original training data.# Choices:#       na.action = "na.omit"#                    Any validation datapoint or image pixel with a value for any#                    categorical predictor not found in the training data will be#                    returned as NA.#       na.action = "na.roughfix"#                    Any validation datapoint or image pixel with a value for any#                    categorical predictor not found in the training data will have#                    the most common category for that predictor substituted,#                    and the a prediction will be made.# You must also let R know which of the predictors are categorical, in other# words, which ones R needs to treat as factors.# This vector must be a subset of the predictors given in predList#file name to store model:MODELfn="RF_BIO_TCandNLCD"#predictors:predList=c("TCB","TCG","TCW","NLCD")#define which predictors are categorical:predFactor=c("NLCD")# Response name and type:response.name="BIO"response.type="continuous"###################################################################################################### build model: ################################################################################################################ create model before batching (only run this code once ever!) ###model.obj = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed,                       na.action="na.roughfix")## End(Not run) # end dontrun

Model Predictions and Diagnostics

Description

Takes model object and makes predictions, runs model diagnostics, and creates graphs and tables of the results.

Usage

model.diagnostics(model.obj = NULL, qdata.trainfn = NULL, qdata.testfn = NULL, folder = NULL, MODELfn = NULL, response.name = NULL, unique.rowname = NULL, diagnostic.flag=NULL, seed = NULL, prediction.type=NULL, MODELpredfn = NULL,  na.action = NULL, v.fold = 10, device.type = NULL, DIAGNOSTICfn = NULL,  res=NULL, jpeg.res = 72, device.width = 7,  device.height = 7, units="in",  pointsize=12, cex=par()$cex, req.sens, req.spec, FPC, FNC, quantiles=NULL,  all=TRUE, subset = NULL, weights = NULL, mtry = NULL, controls = NULL,  xtrafo = NULL, ytrafo = NULL, scores = NULL)

Arguments

Details

model.diagnostics()takes model object and makes predictions, runs model diagnostics, and creates graphs and tables of the results.

model.diagnostics() can be run in a traditional R command mode, where all arguments are specified in the function call. However it can also be used in a full push button mode, where you type in the simple commandmodel.map(), and GUI pop up windows will ask questions about the type of model, the file locations of the data, etc...

When runningmodel.map() on non-Windows platforms, file names and folders need to be specified in the argument list, but other pushbutton selections are handled by theselect.list() function, which is platform independent.

Diagnostic predictions are made my one of four methods, and a text file is generated consisting of three columns: Observation ID, observed values and predicted values. Ifpredition.type = "CV") an additional column indicates which cross-fold each observation fell into. If the models response type is categorical then in addition a column giving the category predicted by majority vote, there are also categories for each possible response category giving the proportion of trees that predicted that category.

A variable importance graph is made. Ifresponse.type = "categorical", category specific graphs are generated for variable importance. These show how much the model accuracy for each category is affected when the values of each predictor variable is randomly permuted.

The packagecorrplot is used to generate a plot of correlation between predictor variables. If there are highly correlated predictor variables, then the variable importances of"RF" and"QRF" models need to be interpreted with care, and users may want to consider looking at the conditional variable importances available for"CF" models produced by theparty package.

Ifmodel.type = "RF", the OOB error is plotted as a function of number of trees in the model. Ifresponse.type = "binary" or Ifresponse.type = "categorical" category specific graphs are generated for OOB error as a function of number of trees.

Ifresponse.type = "binary", a summary graph is made using thePresenceAbsence package and a*.csv spreadsheets are created of optimized thresholds by several methods with their associated error statistics, and predicted prevalence.

Ifresponse.type = "continuous" a scatterplot of observed vs. predicted is created with a simple linear regression line. The graph is labeled with slope and intercept of this line as well as Pearson's and Spearman's correlation coefficients.

Ifresponse.type = "categorical", a confusion matrix is generated, that includes erros of ommission and comission, as well as Kappa, Percent Correctly Classified (PCC) and the Multicategorical Area Under the Curve (MAUC) as defined by Hand and Till (2001) and calculated by the packageHandTill2001.

Value

The function will return a dataframe of the row ID, and the Observed and predicted values.

For Binary response models the predicted probability of presence is returned.

For Categorical Response models the predicted category (by majority vote) is returned as well as a column for each category giving the probability of that category. If necessary,make.names is applied to the categories to create valid column names.

For Continuous response models the predicted value is returned.

Ifprediction.type = "CV" the dataframe also includes a column indicating which cross-validation fold each datapoint was in.

Note

Importance currently unavailable for QRF models.

If you are running cross validation diagnostics on a CF model, the model parameters will NOT automatically be passed tomodel.diagnostics(). For cross validation, it is the users responsibility to be certain that the CF arguments are the same inmodel.build() andmodel.diagnostics().

Also, for some CF model parameters (subset,weights, andscores)ModelMap only provides OOB and independent test set diagnostics, and does not support cross validation diagnostics.

Author(s)

Elizabeth Freeman and Tracey Frescino

References

Breiman, L. (2001) Random Forests. Machine Learning, 45:5-32.

Elith, J., Leathwick, J. R. and Hastie, T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology. 77:802-813.

Hand, D. J., & Till, R. J. (2001). A simple generalisation of the area under the ROC curve for multiple class classification problems. Machine Learning, 45(2), 171-186.

Liaw, A. and Wiener, M. (2002). Classification and Regression by randomForest. R News 2(3), 18–22.

Ridgeway, G., (1999). The state of boosting. Comp. Sci. Stat. 31:172-181

See Also

get.test,model.build,model.mapmake

Examples

## Not run: ######################################################################################################## Run this set up code: ################################################################################################### set seed:seed=38# Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")qdata.testfn = system.file("extdata", "helpexamples","DATATEST.csv", package = "ModelMap")# Define folder for all output:folder=getwd()#identifier for individual training and test data pointsunique.rowname="ID"######################################################################################### Pick one of the following sets of definitions: ################################################################################################## Continuous Response, Continuous Predictors #############file name to store model:MODELfn="RF_Bio_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="BIO"response.type="continuous"########## binary Response, Continuous Predictors #############file name to store model:MODELfn="RF_CONIFTYP_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="CONIFTYP"# This variable is 1 if a conifer or mixed conifer type is present, # otherwise 0.response.type="binary"########## Continuous Response, Categorical Predictors ############# In this example, NLCD is a categorical predictor.## You must decide what you want to happen if there are categories# present in the data to be predicted (either the validation/test set# or in the image file) that were not present in the original training data.# Choices:#       na.action =  "na.omit"#                    Any validation datapoint or image pixel with a value for any#                    categorical predictor not found in the training data will be#                    returned as NA.#       na.action =  "na.roughfix"#                    Any validation datapoint or image pixel with a value for any#                    categorical predictor not found in the training data will have#                    the most common category for that predictor substituted,#                    and the a prediction will be made.# You must also let R know which of the predictors are categorical, in other# words, which ones R needs to treat as factors.# This vector must be a subset of the predictors given in predList#file name to store model:MODELfn="RF_BIO_TCandNLCD"#predictors:predList=c("TCB","TCG","TCW","NLCD")#define which predictors are categorical:predFactor=c("NLCD")# Response name and type:response.name="BIO"response.type="continuous"###################################################################################################### build model: ################################################################################################################ create model ###model.obj = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed,                       na.action="na.roughfix")############################################################################### Then Run this code make validation predictions and diagnostics: #################################################################################### for Out-of-Bag predictions ###MODELpredfn<-paste(MODELfn,"_OOB",sep="")PRED.OOB<-model.diagnostics( model.obj=model.obj,qdata.trainfn=qdata.trainfn,                   folder=folder,                   unique.rowname=unique.rowname,                # Model Validation Arguments                   prediction.type="OOB",                   MODELpredfn=MODELpredfn,                   device.type=c("default","jpeg","pdf"),                   na.action="na.roughfix")PRED.OOB### for Cross-Validation predictions ####MODELpredfn<-paste(MODELfn,"_CV",sep="")#PRED.CV<-model.diagnostics( model.obj=model.obj,#                   qdata.trainfn=qdata.trainfn,#                   folder=folder,#                   unique.rowname=unique.rowname,#                   seed=seed,#                # Model Validation Arguments#                   prediction.type="CV",#                   MODELpredfn=MODELpredfn,#                   device.type=c("default","jpeg","pdf"),#                   v.fold=10,#                   na.action="na.roughfix"#)#PRED.CV### for Independent Test Set predictions ####MODELpredfn<-paste(MODELfn,"_TEST",sep="")#PRED.TEST<-model.diagnostics( model.obj=model.obj,#                   qdata.testfn=qdata.testfn,#                   folder=folder,#                   unique.rowname=unique.rowname,#                # Model Validation Arguments#                   prediction.type="TEST",#                   MODELpredfn=MODELpredfn,#                   device.type=c("default","jpeg","pdf"),#                   na.action="na.roughfix"#)#PRED.TEST)## End(Not run) # end dontrun

Exploratory data analysis

Description

Graphically explores the relationships between the training data and the predictor rasters.

Usage

model.explore(qdata.trainfn = NULL, folder = NULL, predList = NULL, predFactor = FALSE, response.name = NULL, response.type = NULL, response.colors = NULL, unique.rowname = NULL, OUTPUTfn = NULL, device.type = NULL, allow.default.graphics=FALSE, res=NULL, jpeg.res = 72, MAXCELL=100000, device.width = NULL, device.height = NULL, units="in", pointsize=12, cex=1, rastLUTfn = NULL, create.extrapolation.masks = FALSE, na.value = -9999, col.ramp = rainbow(101, start = 0, end = 0.5), col.cat = palette()[-1])

Arguments

Details

Themodel.explore function is intended to aid with preliminary data exploration before model building. It includes graphical tools to explore the relationships between the training data (both predictors and responses) as well as the predictor rasters. It uses thecorrplot package to create a correlation plot of the continuous predictor. This can aid in interpreting themodel.importance.plot output from the models, as Random Forest models divide importance between correlated predictors, while Stochastic Gradient Boosting models assing the majority of the importance to the correlated predictor that is used earlies in the model.

Themodel.explore function also can aid in identifying if additional training data is needed. For example, the maps of the extrapolation masks for the predictor rasters help spot areas of the map where the pixels lie outside the range of the training data, and therefore any model predictions will be extrapolations, and possibly unreliable. The user can decide to either collect additional training data, or mask out these areas of the final prediction output ofmodel.mapmake.

To increase speed, the default behavior for large predictor rasters is to create the graphics from subsampled rasters. (Note: for categorical predictors, the full raster is always used to identify all categories found in the map area.) Ifcreate.extrapolation.masks=TRUE, then the full rasters are used for the extrapolation masks, regardless of size of the reasters. This option runs much slower, as large rasters need to be read into R a block at a time.

Value

Function does not return a value, but does create files.

Graphical files are created for each predictor variable, with file type determined bydevice.type. In addition, ifcreate.extrapolation.masks, an extrapolation mask raster is created for each predictor as well as an overall extrapolation mask, with the value1 for pixels with predictor values within the range of the training data, or categories found in the training data, and the value0 for pixels outside the range of the training data, categories not found in the training data, or NA value. The overall extrapolation mask has0 if any of the predictors for that pixel are extrapolated. Note that this option is much slower to run.

Note

The default graphics device is disabled unlessallow.default.graphics is set toTRUE. These graphics can be slow to produce, and if the on screen graphics device is moved or closed while the graphic is in progress, it can crash R. It is recomended that graphics be written to a file by using jpeg, pdf, etc...device.type.

Author(s)

Elizabeth Freeman

Examples

## Not run: ######################################################################################################## Run this set up code: #####################################################################################################Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")###Define folder for all output:folder=getwd()###identifier for individual training and test data pointsunique.rowname="ID"###predictors:predList=c("TCB","TCG","TCW","NLCD")###define which predictors are categorical:predFactor=c("NLCD")###Create a the filename (including path) for the rast Look up Tables ###rastLUTfn.2001 <- system.file("extdata","helpexamples","LUT_2001.csv",package="ModelMap")###Load rast LUT table, and add path to the predictor raster filenames in column 1 ###rastLUT.2001 <- read.table(rastLUTfn.2001,header=FALSE,sep=",",stringsAsFactors=FALSE)for(i in 1:nrow(rastLUT.2001)){rastLUT.2001[i,1] <- system.file("extdata","helpexamples",rastLUT.2001[i,1],package="ModelMap")}                                 #################Continuous Response######################Response name and type:response.name="BIO"response.type="continuous"###file name to store model:OUTPUTfn="BIO_TCandNLCD.img"###run model.exploremodel.explore(qdata.trainfn=qdata.trainfn,folder=folder,predList=predList,predFactor=predFactor,response.name=response.name,response.type=response.type,unique.rowname=unique.rowname,OUTPUTfn=OUTPUTfn,device.type="jpeg",jpeg.res=144,# Raster argumentsrastLUTfn=rastLUT.2001,na.value=-9999,# colors for continuous predictorscol.ramp=rainbow(101,start=0,end=.5), # colors for categorical predictorscol.cat=c("wheat1","springgreen2","darkolivegreen4",  "darkolivegreen2","yellow","thistle2",  "brown2","brown4"))## End(Not run) # end dontrun

Compares the variable importance of two models with a back to back barchart.

Description

Takes two models and produces a back to back bar chart to compare the importance of the predictor variables. Models can be any combination of Random Forest or Stochastic Gradient Boosting, as long as both models have the same predictor variables.

Usage

model.importance.plot(model.obj.1 = NULL, model.obj.2 = NULL, model.name.1 = "Model 1", model.name.2 = "Model 2", imp.type.1 = NULL, imp.type.2 = NULL, type.label=TRUE, class.1 = NULL, class.2 = NULL, quantile.1=NULL, quantile.2=NULL,col.1="grey", col.2="black", scale.by = "sum", sort.by = "model.obj.1", cf.mincriterion.1 = 0, cf.conditional.1 = FALSE, cf.threshold.1 = 0.2, cf.nperm.1 = 1, cf.mincriterion.2 = 0, cf.conditional.2 = FALSE, cf.threshold.2 = 0.2, cf.nperm.2 = 1, predList = NULL, folder = NULL, PLOTfn = NULL, device.type = NULL, res=NULL, jpeg.res = 72, device.width = 7,  device.height = 7, units="in", pointsize=12, cex=par()$cex,...)

Arguments

Details

The importance measures used in this plot depend on the model type (RF verses SGB) and the response type (continuous, categorical, or binary).

Importance type 1 is permutation based, as described in Breiman (2001). Importance is calculated by randomly permuting each predictor variable and computing the associated reduction in predictive performance using Out Of Bag error for RF models and training error for SGB models. Note that for SGB models permutation based importance measures are still considered experimental. Importance type 2 is model based. For RF models, importance type 2 is calculated by the decrease in node impurities attributable to each predictor variable. For SGB models, importance type 2 is the reduction attributable to each variable in predicting the gradient on each iteration as described in described in Friedman (2001).

For RF models:

For Random Forest models, ifimp.type not specified, importance type defaults toimp.type of1 - permutation importance. For SGB models, permutation importance is considered experimental so importance defaults toimp.type of2 - reduction of gradient of the loss function.

Also, for binary and categorical Random Forest models, class specific importance plots can be generated by the use of theclass argument. Note that class specific importance is only available for Random Forest models with importance type 1.

For CF models:

For binary CF models, ifimportance.type = 2, function uses AUC-based variables importances as described by Janitza et al. (2012). Here, the area under the curve instead of the accuracy is used to calculate the importance of each variable. This AUC-based variable importance measure is more robust towards class imbalance.

Also, for CF models, ifcf.conditional = TRUE, the importance of each variable is computed by permuting within a grid defined by the covariates that are associated (with 1 - p-value greater than threshold) to the variable of interest. The resulting variable importance score is conditional in the sense of beta coefficients in regression models, but represents the effect of a variable in both main effects and interactions. See Strobl et al. (2008) for details. Conditional improtance can be slow for large datasets.

Value

The function returns a two element list:IMP1 is the variable importance formodel.obj.1; and,IMP2 is the variable importance formodel.obj.2. This is mostly intended for CF models, where calculating the conditional importance can represent a considerable time investment. For other model types it would be just as easy to recalcuate importances on the fly as needed.

Note

Importance currently unavailable for QRF models.

Author(s)

Elizabeth Freeman

References

Breiman, L. (2001) Random Forests. Machine Learning, 45:5-32.

Alexander Hapfelmeier, Torsten Hothorn, Kurt Ulm, and Carolin Strobl (2012). A New Variable Importance Measure for Random Forests with Missing Data. Statistics and Computing, http://dx.doi.org/10.1007/s11222-012-9349-1

Torsten Hothorn, Kurt Hornik, and Achim Zeileis (2006b). Unbiased Recursive Partitioning: A Conditional Inference Framework. Journal of Computational and Graphical Statistics, 15 (3), 651-674. Preprint available from http://statmath.wu-wien.ac.at/~zeileis/papers/Hothorn+Hornik+Zeileis-2006.pdf

Silke Janitza, Carolin Strobl and Anne-Laure Boulesteix (2013). An AUC-based Permutation Variable Importance Measure for Random Forests. BMC Bioinformatics.2013, 14 119. http://www.biomedcentral.com/1471-2105/14/119

Carolin Strobl, Anne-Laure Boulesteix, Thomas Kneib, Thomas Augustin, and Achim Zeileis (2008). Conditional Variable Importance for Random Forests. BMC Bioinformatics, 9, 307. http://www.biomedcentral.com/1471-2105/9/307

See Also

model.build

Examples

## Not run: ######################################################################################################## Run this set up code: ################################################################################################### set seed:seed=38# Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")# Define folder for all output:folder=getwd()#identifier for individual training and test data pointsunique.rowname="ID"############################################################################ Continuous Response, Continuous Predictors ###############################################################################file names:MODELfn.RF="RF_Bio_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="BIO"response.type="continuous"########## Build Models #################################model.obj.RF = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn.RF,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed)########## Make Imortance Plot - RF Importance type 1 vs 2 #######model.importance.plot(model.obj.1=model.obj.RF, model.obj.2=model.obj.RF, model.name.1="PercentIncMSE", model.name.2="IncNodePurity",imp.type.1=1,imp.type.2=2,scale.by="sum",sort.by="predList", predList=predList,main="Imp type 1 vs Imp type 2",device.type="default")############################################################################ Categorical Response, Continuous Predictors #############################################################################file name:MODELfn="RF_NLCD_TC"predictors:predList=c("TCB","TCG","TCW")define which predictors are categorical:predFactor=FALSE Response name and type:response.name="NLCD"response.type="categorical"########## Build Model #################################model.obj.NLCD = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed)############## Make Imortance Plot ###################model.importance.plot(model.obj.1=model.obj.NLCD, model.obj.2=model.obj.NLCD, model.name.1="NLCD=41", model.name.2="NLCD=42",class.1="41",class.2="42",scale.by="sum",sort.by="predList", predList=predList,main="Class 41 vs. Class 42",device.type="default")################################################################################ Conditonal inference forest models ###################################################################################predictors:predList=c("TCB","TCG","TCW","NLCD")#define which predictors are categorical:predFactor=c("NLCD")#binary responseresponse.name="CONIFTYP"response.type="binary"MODELfn.CF="CF_CONIFTYP_TCandNLCD"####################### Build Model ##############################model.obj.CF = model.build( model.type="CF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn.CF,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed)################## Make Imortance Plot ###########################Conditional vs. Unconditional importance#model.importance.plot(model.obj.1=model.obj.CF, model.obj.2=model.obj.CF, model.name.1="conditional",model.name.2="unconditional",imp.type.1=1,imp.type.2=1, cf.conditional.1=TRUE,cf.conditional.2=FALSE,scale.by="sum",sort.by="predList", predList=predList,main="Conditional verses Unconditional",device.type="default")## End(Not run) # end dontrun

plot of two-way model interactions

Description

Image or Perspective plot of two-way model interactions. Ranges of two specified predictor variables are plotted on X and Y axis, and fitted model values are plotted on the Z axis. The remaining predictor variables are fixed at their mean (for continuous predictors) or their most common value (for categorical predictors).

Usage

model.interaction.plot(model.obj = NULL, x = NULL, y = NULL, response.category=NULL, quantiles=NULL, all=FALSE, obs=1, qdata.trainfn = NULL, folder = NULL, MODELfn = NULL,  PLOTfn = NULL, pred.means = NULL, xlab = NULL, ylab = NULL, x.range = NULL, y.range = NULL, z.range = NULL, ticktype = "detailed", theta = 55, phi = 40, smooth = "none", plot.type = NULL, device.type = NULL, res=NULL, jpeg.res = 72, device.width = 7,  device.height = 7, units="in", pointsize=12, cex=par()$cex, col = NULL, xlim = NULL, ylim = NULL, zlim = NULL, ...)

Arguments

Details

This 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 of predictor values over the range of both variables. The remaining predictor variables from the model 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.

This function works with both continuous and categorical predictors, though the perspective plot should be interpreted with care for categorical predictors. In particular, thesmooth option is not appropriate if either of the two selected predictor variables is categorical.

For categorical response models, a particular value must be specified for the response using theresponse.category argument.

Author(s)

Elizabeth Freeman

References

This function is adapted fromgbm.perspec version 2.9 April 2007, J Leathwick/J Elith. See appendix S3 from:

Elith, J., Leathwick, J. R. and Hastie, T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology. 77:802-813.

Examples

## Not run: ######################################################################################################## Run this set up code: ################################################################################################### set seed:seed=38# Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")qdata.testfn = system.file("extdata", "helpexamples","DATATEST.csv", package = "ModelMap")# Define folder for all output:folder=getwd()########## Continuous Response, Categorical Predictors #############file name to store model:MODELfn="RF_BIO_TCandNLCD"#predictors:predList=c("TCB","TCG","TCW","NLCD")#define which predictors are categorical:predFactor=c("NLCD")# Response name and type:response.name="BIO"response.type="continuous"#identifier for individual training and test data pointsunique.rowname="ID"###################################################################################################### build model: ################################################################################################################ create model ###model.obj = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed,                       na.action=na.roughfix)################################################################################################# make interaction plots: ################################################################################################################################### Perspective Plots ############################### specify first and third  predictors in 'predList (both continuous) ###model.interaction.plot(model.obj,x=1,y=3, main=response.name, plot.type="persp", device.type="default") ### specify predictors in 'predList' by name (one continuous one factored) ###model.interaction.plot(model.obj,x="TCB", y="NLCD", main=response.name, plot.type="persp", device.type="default") ###################### Image Plots ######################### same as previous example, but image plot ###l <- seq(100,0,length.out=101)c <- seq(0,100,length.out=101)col.ramp <- hcl(h = 120, c = c, l = l) model.interaction.plot(model.obj,x="TCB", y="NLCD",main=response.name,plot.type="image", device.type="default",col = col.ramp) ############################ 3-way Interaction ############################### use 'pred.means' argument to fix values of additional predictors ###### factored 3rd predictor ###interaction between TCG and TCW for 3 most common values of NLCDnlcd<-levels(model.obj$predictor.data$NLCD)nlcd.counts<-table(model.obj$predictor.data$NLCD)nlcd.ordered<-nlcd[order(nlcd.counts,decreasing=TRUE)]for(i in nlcd.ordered[1:3]){pred.means=list(NLCD=i)model.interaction.plot(model.obj,x="TCG", y="TCW",  main=paste("NLCD=",i," (",nlcd.counts[i]," plots)", sep=""),pred.means=pred.means, z.range=c(0,110),theta=290,plot.type="persp", device.type="default") }### continuos 3rd predictor ###tcb<-seq(min(model.obj$predictor.data$TCB),max(model.obj$predictor.data$TCB),length=3)tcb<-signif(tcb,2)for(i in tcb){pred.means=list(TCB=i)model.interaction.plot(model.obj,x="TCG", y="TCW",  main=paste("TCB =",i),pred.means=pred.means, z.range=c(0,120),theta=290,plot.type="persp", device.type="default") }### 4-way Interesting combos ###tcb=c(1300,2900,3400)nlcd=c(11,90,95)for(i in 1:3){pred.means=list(TCB=tcb[i],NLCD=nlcd[i])model.interaction.plot(model.obj,x="TCG", y="TCW",  main=paste("TCB =",tcb[i],"        NLCD =",nlcd[i]),pred.means=pred.means, z.range=c(0,120),theta=290,plot.type="persp", device.type="default") }## End(Not run) #end dontrun

Map Making

Description

Applies models to either ERDAS Imagine image (.img) files or ESRI Grids of predictors to create detailed prediction surfaces. It will handle large predictor files for map making, by reading in the.img files in rows, and output to the.img file the prediction for each data row, before reading the next row of data.

Usage

model.mapmake(model.obj= NULL, folder = NULL, MODELfn = NULL, rastLUTfn = NULL, na.action = NULL, na.value=-9999, keep.predictor.brick=FALSE, map.sd = FALSE, OUTPUTfn = NULL, quantiles=NULL)

Arguments

Details

model.mapmake() can be run in a traditional R command mode, where all arguments are specified in the function call. However it can also be used in a full push button mode, where you type in the simple commandmodel.mapmake(), and GUI pop up windows will ask questions about the type of model, the file locations of the data, etc...

When runningmodel.mapmake() on non-Windows platforms, file names and folders need to be specified in the argument list, but other pushbutton selections are handled by theselect.list() function, which is platform independent.

The R packageraster is used to read spatial rasters into R. The data for production mapping should be in the form of pixel-based raster layers representing the predictors in the model. If there is more than one predictor or raster layer, the layers must all have the same number of columns and rows. The layers must also have the same extent, projection, and pixel size, for effective model development and accuracy. Theraster package functioncompareRaster() is used to check predictor layers for consistency.

The layers must also be in (single or multi-band) raster data formats that can be read by packageraster, for example ESRI Grid or ERDAS Imagine image files. The predictor layers must have continuous or categorical data values. SeewriteRaster for a list of available formats.

To improve processing speed, theraster package is used to create a raster brick object with a layer for each predictor in the model. By default, this brick is a temporary file that is automatically deleated as soon as the map is completed. Ifkeep.predictor.brick=TRUE, the predictor brick with be saved as a nativeraster package file, with a file name created by appending'_brick' to theOUTPUTfn. Warning: these bricks can be quite large, as they contain all the predictor data for every pixel in the map.

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 in the predictions, thus saving computing time.

The functionmodel.mapmake() outputs an rater file of map information suitable to be imported into a GIS. Maps can also be imported back into R using the functionraster() from theraster package. The file extension ofOUTPUTfn determines the write format. IfOUTPUTfn does not include a file extension, output will default to an ERDAS Imagine image file with extension".img"

For Binary response models the output is in the form of predicted probability of presence for each pixel. For Continuous response models the output is the predicted value for each pixel. For Categorical response models the map output depends on the category labels. If the categorical response variable is numeric, the map output will use the original numeric categories. If the categories are non-numeric (for example, character strings), map output is in the form of integer class codes for each pixel, coded for each level of the factored response, and a CSV file containing a look up table is also generated to associate the integer codes with the original values of the response categories.

The first predictor frompredList is used to determine projection of output Imagine Image file.

Value

Themodel.mapmake() function does not return a value, instead it writes a raster file of map information (suitable for importing into a GIS) to the specified folder. The output raster is saved in the format specifed by the file extension ofOUTPUTfn

Themodel.mapmake() function also writes a text file listing the projections of all predictor rasters.

For categorical response models, a csv file map key is written giving the integer code associated with each response category.

Ifkeep.predictor.brick = TRUE then a raster brick of all the predictor rasters from the model is also saved to the specified folder. Ifkeep.predictor.brick = FALSE (the default) then the predictor brick is written to a temprary file, and deleted. Warning: the predictor bricks can be quite large, and saving them can require quite a bit of memory.

Note

Ifmodel.mapmake() is interupted it may leave orphan.gri and.grd files in your temporary directory. Theraster package functionsshowTmpFiles andremoveTmpFiles can be used to locate and remove these files, or they can be deleated manually from the temporary directory.

Author(s)

Elizabeth Freeman and Tracey Frescino

References

Breiman, L. (2001) Random Forests. Machine Learning, 45:5-32.

Liaw, A. and Wiener, M. (2002). Classification and Regression by randomForest. R News 2(3), 18–22.

Ridgeway, G., (1999). The state of boosting. Comp. Sci. Stat. 31:172-181

Simpson, E. H. (1949). Measurement of diversity. Nature.

See Also

get.test,model.build,model.diagnostics,compareRaster,writeRaster

Examples

## Not run: ######################################################################################################## Run this set up code: ################################################################################################### set seed:seed=38# Define training and test files:qdata.trainfn = system.file("extdata", "helpexamples","DATATRAIN.csv", package = "ModelMap")# Define folder for all output:folder=getwd()#identifier for individual training and test data pointsunique.rowname="ID"################################################################################################### Define the model: ##################################################################################################################### Continuous Response, Continuous Predictors #############file name to store model:MODELfn="RF_Bio_TC"#predictors:predList=c("TCB","TCG","TCW")#define which predictors are categorical:predFactor=FALSE# Response name and type:response.name="BIO"response.type="continuous"###################################################################################################### build model: ################################################################################################################ create model ###model.obj = model.build( model.type="RF",                       qdata.trainfn=qdata.trainfn,                       folder=folder,                       unique.rowname=unique.rowname,                       MODELfn=MODELfn,                       predList=predList,                       predFactor=predFactor,                       response.name=response.name,                       response.type=response.type,                       seed=seed,                       na.action="na.roughfix")####################################################################################### Then Run this code to predict map pixels ################################################################################################### Create a the filename (including path) for the rast Look up Tables ###rastLUTfn.2001 <- system.file("extdata","helpexamples","LUT_2001.csv",package="ModelMap")                                          ### Load rast LUT table, and add path to the predictor raster filenames in column 1 ###rastLUT.2001 <- read.table(rastLUTfn.2001,header=FALSE,sep=",",stringsAsFactors=FALSE)for(i in 1:nrow(rastLUT.2001)){rastLUT.2001[i,1] <- system.file("extdata","helpexamples",rastLUT.2001[i,1],package="ModelMap")}                                 ### Define filename for map  output ###OUTPUTfn.2001 <- "RF_BIO_TCandNLCD_01.img"OUTPUTfn.2001 <- paste(folder,OUTPUTfn.2001,sep="/")### Create image files of predicted map data ###model.mapmake( model.obj=model.obj,               folder=folder,               rastLUTfn=rastLUT.2001,           # Mapping arguments               OUTPUTfn=OUTPUTfn.2001               )########################################################################################### run this code to create maps in R ###################################################################################################### Define Color Ramp ###l <- seq(100,0,length.out=101)c <- seq(0,100,length.out=101)col.ramp <- hcl(h = 120, c = c, l = l)### read in map data ###mapgrid.2001 <- raster(OUTPUTfn.2001)#mapgrid.2001 <- setMinMax(mapgrid.2001)### create map ###dev.new(width = 5, height = 5)opar <- par(mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)zlim <- c(0,max(maxValue(mapgrid.2001)))legend.label<-rev(pretty(zlim,n=5))legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]image(  mapgrid.2001, col = col.ramp, zlim=zlim, asp=1, bty="n",         xaxt="n", yaxt="n", main="", xlab="", ylab="")mtext("2001 Imagery",side=3,line=1,cex=1.2)legend(x=xmax(mapgrid.2001),y=ymax(mapgrid.2001),legend=legend.label,fill=legend.colors,bty="n",cex=1.2)mtext("Predictions",side=3,line=1,cex=1.5,outer=TRUE)par(opar)## End(Not run) # end dontrun