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, the smooth 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 the response.category argument.

Author(s)

Elizabeth Freeman

References

This function is adapted from gbm.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

###########################################################################
############################# 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 points

unique.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") 

## Not run: 
### 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 NLCD

nlcd<-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

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(ModelMap)
Loading required package: randomForest
randomForest 4.6-12
Type rfNews() to see new features/changes/bug fixes.
Loading required package: raster
Loading required package: sp
Loading required package: rgdal
rgdal: version: 1.1-10, (SVN revision 622)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.3, released 2015/09/16
 Path to GDAL shared files: /usr/share/gdal/1.11
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-3 
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ModelMap/model.interaction.plot.Rd_%03d_medium.png", width=480, height=480)
> ### Name: model.interaction.plot
> ### Title: plot of two-way model interactions
> ### Aliases: model.interaction.plot
> ### Keywords: models
> 
> ### ** Examples
> 
> ###########################################################################
> ############################# 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 points
> 
> unique.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
+ )
mtry = 1  OOB error = 1819.048 
Searching left ...
Searching right ...
mtry = 2 	OOB error = 1907.638 
-0.04870136 0.05 
> 
> ###########################################################################
> ###################### 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") 
maximum value =  116.61 
Error in dev.new(width = device.width, height = device.height, pointsize = pointsize) : 
  no suitable unused file name for pdf()
Calls: model.interaction.plot -> initialize.device -> dev.new
Execution halted