Lib_HW2

Load Packages

library(sm)           # for sm.density in diagnostics

## Package 'sm', version 2.2-5.5: type help(sm) for summary information

library(leaps)        # to provide combinations
library(MPV)          # for calculating PRESS
library(akima)        # for interp to smooth for plotting
library(fields)       # for surface to plot surface plot

## Loading required package: spam

## Loading required package: dotCall64

## Loading required package: grid

## Spam version 2.2-0 (2018-06-19) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: maps

## See www.image.ucar.edu/~nychka/Fields for
##  a vignette and other supplements.

library(locfit)       # for fitting local polynomial

## locfit 1.5-9.1    2013-03-22

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:MPV':
## 
##     cement

## The following object is masked from 'package:sm':
## 
##     muscle

library(VGAM)         # useful for multynomial fit

## Loading required package: stats4

## 
## Attaching package: 'stats4'

## The following object is masked from 'package:spam':
## 
##     mle

## Loading required package: splines

library(nnet)         # to fit multinomial regression
library(verification)

## Loading required package: boot

## 
## Attaching package: 'boot'

## The following objects are masked from 'package:VGAM':
## 
##     logit, simplex

## The following object is masked from 'package:MPV':
## 
##     motor

## The following object is masked from 'package:sm':
## 
##     dogs

## Loading required package: CircStats

## 
## Attaching package: 'CircStats'

## The following objects are masked from 'package:VGAM':
## 
##     dcard, rcard

## Loading required package: dtw

## Loading required package: proxy

## 
## Attaching package: 'proxy'

## The following object is masked from 'package:spam':
## 
##     as.matrix

## The following objects are masked from 'package:stats':
## 
##     as.dist, dist

## The following object is masked from 'package:base':
## 
##     as.matrix

## Loaded dtw v1.20-1. See ?dtw for help, citation("dtw") for use in publication.

## 
## Attaching package: 'verification'

## The following object is masked from 'package:VGAM':
## 
##     exponential

library(latex2exp)    #  LaTeX math formulas to R's
library(maps)         # useful for graphs function map
library(rgdal)        # for converting projection of coordinates

## Loading required package: sp

## rgdal: version: 1.3-4, (SVN revision 766)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/DELL/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/DELL/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1

library(ggplot2)      # for plotting spatial map
library(reshape2)     # for using the melt function (Convert an object into a molten data frame)
library(scales)      # for Visualization spatial map
library(spBayes)     # for bayesian kriging

## Loading required package: coda

## 
## Attaching package: 'coda'

## The following object is masked from 'package:VGAM':
## 
##     nvar

## Loading required package: magic

## Loading required package: abind

## Loading required package: Formula

library(geoR)

## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------

library(mclust)

## Package 'mclust' version 5.4.1
## Type 'citation("mclust")' for citing this R package in publications.

## 
## Attaching package: 'mclust'

## The following object is masked from 'package:maps':
## 
##     map

library(tree)        # classification and regression trees
library(randomForest)  #Random Forest

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

library(RColorBrewer) # Color schemes for maps
library(Hmisc)        # Add minor tick marks

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:boot':
## 
##     melanoma

## Loading required package: survival

## 
## Attaching package: 'survival'

## The following object is masked from 'package:boot':
## 
##     aml

## 
## Attaching package: 'Hmisc'

## The following object is masked from 'package:fields':
## 
##     describe

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(wesanderson)
library(cluster)      # To perform clustering of the Extremes

## 
## Attaching package: 'cluster'

## The following object is masked from 'package:maps':
## 
##     votes.repub

library(kohonen) #required to perform som

## 
## Attaching package: 'kohonen'

## The following object is masked from 'package:mclust':
## 
##     map

## The following object is masked from 'package:maps':
## 
##     map

Quilt plotting function

Quilt_plotting=function(lon, lat, ypred, lon1, lat1, yob,type=""){
    if (type=="binary"){
      par(mfrow=c(3,1))
      quilt.plot(lon, lat,ypred$fit,xlab="Longitude (m)",ylab="Latitude (m)",main='Posterior Binary Precipitation on DEM Grid',zlim=range(ypred$fit,yob))
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
      
      quilt.plot(lon1, lat1,yob,xlab="Longitude (m)",ylab="Latitude (m)",main='Observed Binary Precipitation',zlim=range(ypred$fit,yob))
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
      
      quilt.plot(lon, lat,ypred$se,xlab="Longitude (m)",ylab="Latitude (m)",main='Standard Error on DEM Grid')
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
    } else{
      par(mfrow=c(3,1))
      quilt.plot(lon, lat,ypred$fit,xlab="Longitude (m)",ylab="Latitude (m)",main='Posterior Precipitation on DEM Grid (mm)',zlim=range(ypred$fit,yob))
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
      
      quilt.plot(lon1, lat1,yob,xlab="Longitude (m)",ylab="Latitude (m)",main='Observed Precipitation (mm)',zlim=range(ypred$fit,yob))
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
      
      quilt.plot(lon, lat,ypred$se,xlab="Longitude (m)",ylab="Latitude (m)",main='Standard Error on DEM Grid (mm)')
      grid(col="gray70",lty=2)
      US(add=TRUE, col="gray50", lwd=2,xlim=range(-125,-100))
      box()
    }
  }

Drop-10% CV function

Drop_10_pred = function(X,bestTree,type) {
  mod_data = X
  N = length(X$prec)
  drop = 10
  nsample = 500
  i_full = 1:N
  # initialize skill score vectors
  skill_rmse = vector(mode="numeric", length=nsample)
  skill_cor = vector(mode="numeric", length=nsample)
  for (i in 1:nsample){
    if (type=="tree"){
      i_drop = sample(i_full,N*drop/100)            # can add argument replace=TRUE
      drop_dt = mod_data[-i_drop,] # drop 10% precip value
      myTree =tree(prec ~ PC1+PC2+PC3, data = drop_dt, model = T)
      drop_mod= prune.tree(myTree, best = bestTree)
      drop_pred=predict(drop_mod,newdata=mod_data[i_drop,])
      drop_actual = mod_data[i_drop,1]
      skill_rmse[i] = sqrt(mean((drop_actual - drop_pred)^2))
      skill_cor[i] = cor(drop_actual,drop_pred)
    }else{
      i_drop = sample(i_full,N*drop/100)            # can add argument replace=TRUE
      drop_dt = mod_data[-i_drop,] # drop 10% precip value
      drop_mod=randomForest(prec ~ PC1+PC2+PC3, data = drop_dt)
      drop_pred=predict(drop_mod,newdata=mod_data[i_drop,])
      drop_actual = mod_data[i_drop,1]
      skill_rmse[i] = sqrt(mean((drop_actual - drop_pred)^2))
      skill_cor[i] = cor(drop_actual,drop_pred)
    }
    
  }
  CV=as.data.frame(cbind(skill_rmse,skill_cor))
  names(CV)=c("rmse","cor")
  return(CV)
}

logit2prob = function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}

Function to prune the tree

treeFun = function(myTree, toPlot = F, title = ""){

  #Perform CV on tree object
  cvTree <- cv.tree(myTree)
  optTree <- which.min(cvTree$dev)
  bestTree <- cvTree$size[optTree]
  
  #prune Tree based on CV results
  pruneTree <- prune.tree(myTree, best = bestTree)
  
  #If plotting is selected
  if(toPlot){
    
    #Plot unpruned Tree
    plot(myTree)
    text(myTree, cex = .75)
    title(main = paste("Unpruned Tree for", title))
    
    #Plot CV
    plot(cvTree$size, cvTree$dev, type = "b", 
         main = paste("Cross Validation for", title))
    
    #Plot Prunned Tree
    plot(pruneTree)
    text(pruneTree, cex = .75)
    title(main = paste("Pruned Tree for", title))
  } 
  pruneTree$besTree=bestTree
  return(pruneTree)
}

Extreme clustering function

pam_fmado_ll = function (x, k, ll) {
  
  # x - Design matrix, typically colums are stations, rows are time 
  # k - Number of clusters
  # ll - two column matrix of lat long points (or preferably projected) with N rows
  
  N = ncol(x) # number of stations 
  T = nrow(x) # number of time points
  
  # compute the F-madogram distance
  V = array(NaN, dim = c(T, N))
  for (p in 1:N) {
    x.vec = as.vector(x[, p])
    Femp = ecdf(x.vec)(x.vec)
    V[, p] = Femp
  }
  DD = dist(t(V), method = "manhattan", diag = TRUE, upper = TRUE)/(2 * T)
  
  # weight by physical distance
  DDll = dist(ll,method='manhattan')
  DDw = as.matrix(DD) + t(t(as.matrix(DDll))/apply(as.matrix(DDll),2,max))*max(as.matrix(DD))
  
  # do the clustering
  output = pam(DDw, k, diss = TRUE, medoids = NULL)
  return(output)
}

multiplot function

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

Lib_HW2

Alvaro Ossandon

November 24, 2018

Load Packages

Quilt plotting function

Drop-10% CV function

Function to prune the tree

Extreme clustering function

multiplot function