R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Maps Predicted Values and Clusters for modgam Objects

colormap

R Documentation

Maps Predicted Values and Clusters for modgam Objects

Description

Displays a color image map, including a legend, scale bar, and optional North arrow, showing crude or adjusted odds ratios (or linear predictors) for a grid of points. Irregular grids are allowed. Also draws contour lines for regions of signficantly increased or decreased values of the outcome variable ("clusters"), if permutation ranks are provided. Designed to display modgam objects but can be used with other model results if the modgamobj list contains the correct elements.

Usage

colormap(modgamobj, map = NULL, add=F, contours="none", mapmin = NULL, mapmax = NULL, 
         arrow = T, axes = F, ptsize = 0.9, alpha = 0.05)

Arguments

`modgamobj`	(Required) A list containing at least these two elements: "grid" (a 2 column data frame with X and Y coordinates) and "fit" (a vector of linear predictors, with length equal to the number of grid points). If the list contains an element named "OR" (a vector of odds ratios, with length equal to the number of grid points) then "OR" will be used instead of "fit". If the list contains the element "pointwise" (a vector of percentile ranks generated by permutation test) then those values will be used to generate contours. The correct list format is provided as output from the modgam function (see examples).
`map`	Can be used to map predicted values on a base map from the `map` function in the maps package, or on a base map produced by the `readShapePoly` function in maptools package. ReadShapePoly reads maps from external files in the ESRI shapefile format. `map=NULL` produces a color image without any base map.
`add`	Use `add=T` to add the color map to an existing plot. This will often result in loss of the legend and scale, which are added outside of the normal map boundaries. `add` is ignored when a map is provided using the `map` argument.
`contours`	Use `contours="response"` to add contour lines for the predicted response, for example to draw isoboles for mixtures of exposures. Use `contours="permrank"` to add contour lines for pointwise p-values computed from the permutation ranks, at alpha/2 and (1-alpha)/2. The default is "none" which produces no contour lines.
`mapmin`	The minimum value for the color scale legend
`mapmax`	The maximum value for the color scale legend
`arrow`	Use arrow=T to add a North arrow to the map.
`axes`	Use axes=T to add axes to the map (useful for chemical mixture isoboles).
`ptsize`	The size of the points used to fill in colors on the map. Increase to remove white space inside the map or decrease to remove color outside the map boundaries. NOTE: white space can also be eliminated by increasing the grid size in `predgrid`, which is often preferable as it results in a higher resolution map.
`alpha`	The nominal pointwise type I error rate; only used when `contours="permrank"`.

Value

Produces an image map showing crude or adjusted linear predictors (or odds ratios). If the base map is in readShapePoly format a scale bar is included. The scale bar assumes that the X and Y coordinates are provided in meters.

Warning

Note that the contour lines use a pointwise nominal type I error rate of alpha; the chance of a type I error occurring for at least one point on the map is much higher, typically approaching 100% at alpha=0.05, because a spatial prediction grid generally contains many points.

Author(s)

Veronica Vieira, Scott Bartell, and Robin Bliss

Send bug reports to sbartell@uci.edu.

Examples


### Load readShapePoly base map and data 
data(MAmap)		
data(MAdata) 
# Create a grid on the base map (PBSmapping package recommended)
if(require(PBSmapping)) MAgrid <- predgrid(MAdata, MAmap) else 
	MAgrid <- predgrid(MAdata)		
# fit crude GAM model to the MA data using span size of 50%
fit1 <- modgam(MAdata, MAgrid, m="crude", sp=0.5)  
# Plot a map showing crude odds ratios
colormap(fit1, MAmap)					

#### A detailed example including map projections and data trimming
# NOTE: this example requires the maps, mapproj, and PBSmapping packages
# Convert base map and beer tweet data locations to US Albers projection 
# projected coords yield better distance estimates than (lat,long)
if(require(maps) & require(mapproj) & require(PBSmapping)) {
	USmap <- map("state",projection="albers",parameters=c(29.5,45.5),
				plot=FALSE,fill=TRUE,col="transparent")
	data(beertweets)
	case <- beertweets$beer
    # Reuse last map projection to convert data coordinates	
	XY <- mapproject(beertweets$longitude,beertweets$latitude)[1:2]  
	beerproj <- data.frame(case, XY[1], XY[2])			 
	# Generate grid on the US map, trimmed to range of beer data
	USgrid <- predgrid(beerproj, USmap)						
    # Fit unadjusted model--geolocation only
	fit2 <- modgam(beerproj, USgrid, m="unadjusted", sp=0.2)	
	dev.new(width=7,height=5)
    colormap(fit2, USmap)
	title(main="Beer Tweet Odds Ratios")	
}

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(MapGAM)
Loading required package: sp
Loading required package: gam
Loading required package: splines
Loading required package: foreach
Loaded gam 1.12

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MapGAM/colormap.Rd_%03d_medium.png", width=480, height=480)
> ### Name: colormap
> ### Title: Maps Predicted Values and Clusters for modgam Objects
> ### Aliases: colormap
> ### Keywords: hplot misc smooth
> 
> ### ** Examples
> 
> ## No test: 
> ### Load readShapePoly base map and data 
> data(MAmap)		
> data(MAdata) 
> # Create a grid on the base map (PBSmapping package recommended)
> if(require(PBSmapping)) MAgrid <- predgrid(MAdata, MAmap) else 
+ 	MAgrid <- predgrid(MAdata)		
Loading required package: PBSmapping

-----------------------------------------------------------
PBS Mapping 2.69.76 -- Copyright (C) 2003-2016 Fisheries and Oceans Canada

PBS Mapping comes with ABSOLUTELY NO WARRANTY;
for details see the file COPYING.
This is free software, and you are welcome to redistribute
it under certain conditions, as outlined in the above file.

A complete user guide 'PBSmapping-UG.pdf' is located at 
/home/ddbj/local/lib64/R/library/PBSmapping/doc/PBSmapping-UG.pdf

Packaged on 2015-04-23
Pacific Biological Station, Nanaimo

All available PBS packages can be found at
http://code.google.com/p/pbs-software/

To see demos, type '.PBSfigs()'.
-----------------------------------------------------------


> # fit crude GAM model to the MA data using span size of 50%
> fit1 <- modgam(MAdata, MAgrid, m="crude", sp=0.5)  
The unadjusted model is: 
Case ~ lo(Xcoord, Ycoord, span = 0.5)
Family: binomial Link: logit
> # Plot a map showing crude odds ratios
> colormap(fit1, MAmap)					
> 
> #### A detailed example including map projections and data trimming
> # NOTE: this example requires the maps, mapproj, and PBSmapping packages
> # Convert base map and beer tweet data locations to US Albers projection 
> # projected coords yield better distance estimates than (lat,long)
> if(require(maps) & require(mapproj) & require(PBSmapping)) {
+ 	USmap <- map("state",projection="albers",parameters=c(29.5,45.5),
+ 				plot=FALSE,fill=TRUE,col="transparent")
+ 	data(beertweets)
+ 	case <- beertweets$beer
+     # Reuse last map projection to convert data coordinates	
+ 	XY <- mapproject(beertweets$longitude,beertweets$latitude)[1:2]  
+ 	beerproj <- data.frame(case, XY[1], XY[2])			 
+ 	# Generate grid on the US map, trimmed to range of beer data
+ 	USgrid <- predgrid(beerproj, USmap)						
+     # Fit unadjusted model--geolocation only
+ 	fit2 <- modgam(beerproj, USgrid, m="unadjusted", sp=0.2)	
+ 	dev.new(width=7,height=5)
+     colormap(fit2, USmap)
+ 	title(main="Beer Tweet Odds Ratios")	
+ }
Loading required package: maps

 # maps v3.1: updated 'world': all lakes moved to separate new #
 # 'lakes' database. Type '?world' or 'news(package="maps")'.  #


Loading required package: mapproj
The unadjusted model is: 
case ~ lo(x, y, span = 0.2)
Family: binomial Link: logit
Error in dev.new(width = 7, height = 5) : 
  no suitable unused file name for pdf()
Execution halted