# Rscript

############################################################
#################### 60 Characters Wide ####################
# 06-05-10

x <- read.table('~/Desktop/LittleCont/data.txt', header = TRUE)
plot(x$PERENT, x$TEMP)

############################################################
# Least squares linear model

abline(lm(x$TEMP ~ x$PERENT)$coef, col = 'black')
text(0.6, 12, labels = 'least squares', col = 'black')

############################################################
# Lowess model

lines(lowess(x$PERENT, x$TEMP, f = 2/3), col = 'red')
text(0.6, 8, labels = 'lowess f = 2/3', col = 'red')

lines(lowess(x$PERENT, x$TEMP, f = 1/3), col = 'green')
text(0.6, 6, labels = 'lowess f = 1/3', col = 'green')

############################################################
# Fitted cosine

x.nls <- nls(TEMP ~ A * cos(B * (PERENT + C)) + D,
  data = x, start = list(A = 1, B = 1, C = 1, D = 1),
  trace = TRUE, alg = 'port')
plot(x, xlim = c(-1.2, 1.2), ylim = c(-5, 30))
curve(summary(x.nls)$parameters[1] *
      cos(summary(x.nls)$parameters[2] * (x +
      summary(x.nls)$parameters[3])) +
      summary(x.nls)$parameters[4],
      from = -10, to = 10, col = 'blue', add = TRUE, n = 10001)
points(x$PERENT,fitted(x.nls), col = 'red', pch = '+')

############################################################
# Fitted cubic

x.nls <- nls(TEMP ~ A * PERENT^3 + B * PERENT^2 + C * PERENT + D,
  data = x, start = list(A = 0, B = 0, C = 0, D = 0),
  trace = TRUE, alg = 'port')

curve(summary(x.nls)$parameters[1] * x^3 +
      summary(x.nls)$parameters[2] * x^2 +
      summary(x.nls)$parameters[3] * x +
      summary(x.nls)$parameters[4],
      from = -10, to = 10, col = 'green', add = TRUE, n = 10001)
points(x$PERENT,fitted(x.nls), pch = '.')

############################################################
# Why isn't this curve a better fit?

curve(120 * (x - 0.7)^3 + 3 * x^2 + 21, from = -10, to = 10, col = 'yellow', add = TRUE, n = 100001)

curve(120 * x^3 - 250 * x^2 + 175 * x - 19, from = -10, to = 10, add = TRUE, n = 100001, col = 'pink')

############################################################
# Still doesn't work, even changing the starting values of
#  the parameters

x.nls <- nls(TEMP ~ A * PERENT^3 + B * PERENT^2 + C * PERENT + D,
  data = x, start = list(A = 120, B = -250, C = 175, D = -19),
  trace = TRUE)

curve(summary(x.nls)$parameters[1] * x^3 +
      summary(x.nls)$parameters[2] * x^2 +
      summary(x.nls)$parameters[3] * x +
      summary(x.nls)$parameters[4],
      from = -10, to = 10, col = 'grey',
      lwd = 5, add = TRUE, n = 10001)

############################################################
# A synthetic test?

x <- -49:50
y <- 5*(x+5)^3 + 100000 * rnorm(100)
plot(x,y)

x.nls <- nls(y ~ A * x^3 + B * x^2 + C * x + D,
  start = list(A = 0, B = 0, C = 0, D = 0),
  trace = TRUE)

curve(summary(x.nls)$parameters[1] * x^3 +
      summary(x.nls)$parameters[2] * x^2 +
      summary(x.nls)$parameters[3] * x +
      summary(x.nls)$parameters[4],
      from = -10, to = 10, col = 'grey',
      lwd = 5, add = TRUE, n = 10001)

# This works fine....so the problem with the real data
#  must be that it is not uniformly distributed in x so
#  so the least squares algorithm is willing to ignore
#  a poor fit in the sparse (right hand) part of the graph
#  in order to get a small improvement in the densely
#  occupied part of the space.

############################################################
# Kernel ridge regression

x <- read.table('~/Desktop/LittleCont/data.txt', header = TRUE)
plot(x$PERENT, x$TEMP)

n <- nrow(x)
X <- matrix(x$PERENT,ncol=1)
y <- matrix(x$TEMP,ncol=1)
plot(X,y)

#K <- function(u,v) exp(-(u-v)^2)
K <- function(u,v) (u-v)^3

Km <- matrix(nrow=n, ncol=n)
for (i in 1:n) for (j in (i:n)) Km[i,j] <- Km[j,i] <- K(X[i,], X[j,])

gamma <- 0.001
alphaHat <- solve( Km + gamma*diag(n)) %*% y

m<-100
Xnew <- matrix( seq(0,1, length=m), ncol=1)

fitted.values <- rep(0,m)
for (i in 1:m) for (j in 1:n) fitted.values[i] <-
        fitted.values[i] + K(Xnew[i,],X[j,]) * alphaHat[j]

points(Xnew,fitted.values, type="l", col="red")

###########

K <- function(u,v) exp(-(u-v)^2)

Km <- matrix(nrow=nrow(x), ncol=nrow(x))
for (i in 1:nrow(x))
  for (j in (1:i))
    Km[i,j] <- Km[j,i] <- K(x[i,], x[j,])

gamma <- 0.001
alphaHat <- solve(Km + gamma*diag(nrow(x))) %*% x$TEMP

m <- 100
Xnew <- matrix(seq(0,1, length = m), ncol = 1)

fitted.values <- rep(0,m)
for (i in 1:m) for (j in 1:n) fitted.values[i] <-
        fitted.values[i] + K(Xnew[i,],x[j,]) * alphaHat[j]

points(Xnew,fitted.values, type="l", col="red")

############################################################
#################### 60 Characters Wide ####################
# 06-05-15

x <- read.table('~/Desktop/LittleCont/data7studies.csv',
   header = TRUE, sep = '\t')
plot(x$P, x$MATdegC, type = 'n',
     main = 'Comparisons of Leaf Margin Analysis Training Data',
     xlab = 'Percentage of Entire Margined Leaves',
     ylab = 'Mean Annual Temperature in degrees C')
for(i in c(5, 4, 7, 6, 3, 2, 1)){
  points(x$P[x$STUDY == levels(x$STUDY)[i]],
     x$MATdegC[x$STUDY == levels(x$STUDY)[i]],
     col = as.numeric(x$STUDY[x$STUDY == levels(x$STUDY)[i]]),
     pch = 19, cex = 1.2)
  text(100, c(-2, 0, 2, 4, 6, 8, 10)[i], labels = levels(x$STUDY)[i],
     col = (1:length(levels(x$STUDY)))[i], adj = 1)
  lines(lowess(x$P[x$STUDY == levels(x$STUDY)[i]],
     x$MATdegC[x$STUDY == levels(x$STUDY)[i]], f = 2/3),
     col = (1:length(levels(x$STUDY)))[i], lwd = 3)
}
text(0, 25, 'curves are loess fits with bandwidth = 2/3', adj = 0)

############################################################
#################### 60 Characters Wide ####################
# 06-05-31

y <- read.table('~/Desktop/LittleCont/all.data.csv',
   header = TRUE, sep = '\t')
par(mfrow = c(1,2))
plot(y$PERENT, y$TEMP, col = hsv(y$PERENT, 1, 1), cex = 1)
plot(y$X_COORD, y$Y_COORD, pch = 20, col = hsv(y$PERENT, 1, 1), cex = 1)

############################################################
#################### 60 Characters Wide ####################
# 06-06-08

z <- read.table('~/Desktop/LittleCont/result_walton_wholecontinent.csv',
   header = TRUE, sep = '\t')

#normalizing to [0,1]
x.norm <- (z$X_COORD - min(z$X_COORD, na.rm = TRUE)) /
   (max(z$X_COORD, na.rm = TRUE) - min(z$X_COORD, na.rm = TRUE))
y.norm <- (z$Y_COORD - min(z$Y_COORD, na.rm = TRUE)) /
   (max(z$Y_COORD, na.rm = TRUE) - min(z$Y_COORD, na.rm = TRUE))
T.norm <- (z$TEMP - min(z$TEMP, na.rm = TRUE)) /
   (max(z$TEMP, na.rm = TRUE) - min(z$TEMP, na.rm = TRUE))

par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(T.norm, 1, 1),
   cex = 1, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by temperature')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(T.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by temperature')

par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(x.norm, 1, 1),
   cex = 1, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by latitude')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(x.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by latitude')

par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(y.norm, 1, 1),
   cex = 1, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by longitude')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(y.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by longitude')

par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(z$Percent.of.entire, 1, 1),
   cex = 1, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by P')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(z$Percent.of.entire, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by P')

w <- read.table('~/Desktop/LittleCont/result_walton_east.csv',
   header = TRUE, sep = '\t')

T.norm <- (w$TEMP - min(w$TEMP, na.rm = TRUE)) /
   (max(w$TEMP, na.rm = TRUE) - min(w$TEMP, na.rm = TRUE))

par(mfrow = c(1,2))
plot(w$Percent_entire, w$TEMP, col = hsv(w$Percent_entire, 1, 1),
   cex = 1, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by P')
plot(w$X_COORD, w$Y_COORD, pch = 20, col = hsv(w$Percent_entire, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by P')

v <- read.table('~/Desktop/LittleCont/eastover20less300m.csv',
   header = TRUE, sep = '\t')

wa <- read.table('~/Desktop/LittleCont/westamerica.csv',
   header = TRUE, sep = '\t')
wa1000 <- read.table('~/Desktop/LittleCont/westover1000rain.csv',
   header = TRUE, sep = '\t')

plot(wa$PERENTIRE, wa$TEMP, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (degrees C)', main = 'Western North America')
abline(lm(wa$TEMP ~ wa$PERENTIRE)$coef)

summary(lm(wa$TEMP ~ wa$PERENTIRE))

pv <- format(pf(summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[1],
         summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[2],
         summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[3], lower.tail = FALSE),
         digits = 2)

rs <- format(summary(lm(wa$TEMP ~ wa$PERENTIRE))$r.squared, digits = 2)
text(0.5, 3, paste('n = ', nrow(wa)), adj = 0)
text(0.5, 1.5, paste('r^2 = ', rs), adj = 0)
text(0.5, 0, paste('p = ', pv), adj = 0)

plot(wa1000$PERENTIRE, wa1000$TEMP, xlab = 'Percent Entire Leaves', ylab = 'Mean Annual Temperature (degrees C)', main = 'Western North America with > 1 m annual rainfall')
abline(lm(wa1000$TEMP ~ wa1000$PERENTIRE)$coef)

summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))

pv <- format(pf(summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[1],
         summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[2],
         summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[3], lower.tail = FALSE),
         digits = 2)

rs <- format(summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$r.squared, digits = 2)
text(0.15, 15, paste('n = ', nrow(wa1000)), adj = 0)
text(0.15, 13.5, paste('r^2 = ', rs), adj = 0)
text(0.15, 12, paste('p = ', pv), adj = 0)

############################################################
#################### 60 Characters Wide ####################
# 06-10-19

# This is to create all the figures for the revised MS

PATH <- '/Users/wag/Desktop/LittleCont/AdamsCollab/JBiogeogSubmission/revision/data/'

# Using loess()
y <- read.table(paste(PATH,'all.data.csv', sep = ''),
   header = TRUE, sep = '\t')
par(mfrow = c(1,2))
plot(y$PERENT, y$TEMP, col = hsv(y$PERENT, 1, 1), cex = 1)
y.l <- loess(y$TEMP ~ y$PERENT, span = 2/3, se = TRUE)
print(y.l)
pred <- predict(y.l, data.frame(PERENT = seq(0, 1, 0.01)))
lines(seq(0, 1, 0.01), pred)
text(0.6, 8, labels = 'lowess bandwidth = 2/3')
plot(y$X_COORD, y$Y_COORD, pch = 20, col = hsv(y$PERENT, 1, 1), cex = 1)

# Using lowess()
y <- read.table(paste(PATH, 'all.data.csv', sep = ''),
   header = TRUE, sep = '\t')
par(mfrow = c(1,2))
plot(y$PERENT, y$TEMP, col = hsv(y$PERENT, 1, 1), cex = 1)
lines(lowess(y$PERENT, y$TEMP, f = 2/3))
text(0.6, 8, labels = 'lowess bandwidth = 2/3')
plot(y$X_COORD, y$Y_COORD, pch = 20, col = hsv(y$PERENT, 1, 1), cex = 1)

# Why so different?

z <- read.table(paste(PATH, 'result_walton_wholecontinent.csv', sep = ''),
   header = TRUE, sep = '\t')

#normalizing to [0,1]
x.norm <- (z$X_COORD - min(z$X_COORD, na.rm = TRUE)) /
   (max(z$X_COORD, na.rm = TRUE) - min(z$X_COORD, na.rm = TRUE))
y.norm <- (z$Y_COORD - min(z$Y_COORD, na.rm = TRUE)) /
   (max(z$Y_COORD, na.rm = TRUE) - min(z$Y_COORD, na.rm = TRUE))
T.norm <- (z$TEMP - min(z$TEMP, na.rm = TRUE)) /
   (max(z$TEMP, na.rm = TRUE) - min(z$TEMP, na.rm = TRUE))

#################### 60 Characters Wide ####################
# FIGURE 1

pdf('fig1.pdf', width = 7.5, height = 3.75, family = 'Times', pointsize = 8)
par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(T.norm, 1, 1),
   cex = 1, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by temperature')
lines(lowess(z$Percent.of.entire, z$TEMP, f = 2/3))
text(0.7, 8, labels = 'lowess bandwidth = 2/3')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(T.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by temperature')
dev.off()
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# EXTRA FIG
par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(x.norm, 1, 1),
   cex = 1, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by longitude')
lines(lowess(z$Percent.of.entire, z$TEMP, f = 2/3))
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(x.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by latitude')
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# EXTRA FIG
par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(y.norm, 1, 1),
   cex = 1, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by latitude')
lines(lowess(z$Percent.of.entire, z$TEMP, f = 2/3))
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(y.norm, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by longitude')
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# FIG 2
pdf('fig2.pdf', width = 7.5, height = 3.75, family = 'Times', pointsize = 8)
par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, col = hsv(z$Percent.of.entire, 1, 1),
   cex = 1, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by P')
lines(lowess(z$Percent.of.entire, z$TEMP, f = 2/3))
text(0.7, 8, labels = 'lowess bandwidth = 2/3')
plot(z$X_COORD, z$Y_COORD, pch = 20, col = hsv(z$Percent.of.entire, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by P')
dev.off()
#################### 60 Characters Wide ####################

w <- read.table(paste(PATH, 'result_walton_east.csv', sep = ''),
   header = TRUE, sep = '\t')

T.norm <- (w$TEMP - min(w$TEMP, na.rm = TRUE)) /
   (max(w$TEMP, na.rm = TRUE) - min(w$TEMP, na.rm = TRUE))

#################### 60 Characters Wide ####################
# FIG 3
pdf('fig3.pdf', width = 7.5, height = 3.75, family = 'Times', pointsize = 8)
par(mfrow = c(1,2))
plot(w$Percent_entire, w$TEMP, col = hsv(w$Percent_entire, 1, 1),
   cex = 1, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by P')
lines(lowess(w$Percent_entire, w$TEMP, f = 2/3))
plot(w$X_COORD, w$Y_COORD, pch = 20, col = hsv(w$Percent_entire, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by P')
dev.off()
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# FIG 6
v <- read.table(paste(PATH, 'eastover20less300m.csv', sep = ''),
   header = TRUE, sep = '\t')

wa <- read.table(paste(PATH, 'westamerica.csv', sep = ''),
   header = TRUE, sep = '\t')
wa1000 <- read.table(paste(PATH, 'westover1000rain.csv', sep = ''),
   header = TRUE, sep = '\t')

pdf('fig6.pdf', width = 3.75, height = 3.75, family = 'Times', pointsize = 8)
plot(wa$PERENTIRE, wa$TEMP, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (degrees C)', main = 'Western North America')
abline(lm(wa$TEMP ~ wa$PERENTIRE)$coef)

summary(lm(wa$TEMP ~ wa$PERENTIRE))

pv <- format(pf(summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[1],
         summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[2],
         summary(lm(wa$TEMP ~ wa$PERENTIRE))$fstatistic[3], lower.tail = FALSE),
         digits = 2)

rs <- format(summary(lm(wa$TEMP ~ wa$PERENTIRE))$r.squared, digits = 2)
text(0.5, 3, paste('n = ', nrow(wa)), adj = 0)
text(0.5, 1.5, paste('r^2 = ', rs), adj = 0)
text(0.5, 0, paste('p = ', pv), adj = 0)
dev.off()
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# FIG 5
pdf('fig5.pdf', width = 3.75, height = 3.75, family = 'Times', pointsize = 8)
plot(wa1000$PERENTIRE, wa1000$TEMP, xlab = 'Proportion Entire Leaves', ylab = 'Mean Annual Temperature (degrees C)', main = 'Western North America with > 1 m annual rainfall')
abline(lm(wa1000$TEMP ~ wa1000$PERENTIRE)$coef)

summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))

pv <- format(pf(summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[1],
         summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[2],
         summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$fstatistic[3], lower.tail = FALSE),
         digits = 2)

rs <- format(summary(lm(wa1000$TEMP ~ wa1000$PERENTIRE))$r.squared, digits = 2)
text(0.15, 15, paste('n = ', nrow(wa1000)), adj = 0)
text(0.15, 13.5, paste('r^2 = ', rs), adj = 0)
text(0.15, 12, paste('p = ', pv), adj = 0)
dev.off()
#################### 60 Characters Wide ####################

#################### 60 Characters Wide ####################
# FIG 4
pdf('fig4.pdf', width = 7.5, height = 5, family = 'Times', pointsize = 8)
x <- read.table(paste(PATH, 'data7studies.csv', sep = ''),
   header = TRUE, sep = '\t')
plot(x$P, x$MATdegC, type = 'n',
     main = 'Comparisons of Leaf Margin Analysis Training Data',
     xlab = 'Percent of Entire Margined Leaves',
     ylab = 'Mean Annual Temperature in degrees C')
for(i in c(5, 4, 7, 6, 3, 2, 1)){
  points(x$P[x$STUDY == levels(x$STUDY)[i]],
     x$MATdegC[x$STUDY == levels(x$STUDY)[i]],
     col = as.numeric(x$STUDY[x$STUDY == levels(x$STUDY)[i]]),
     pch = 19, cex = 1.2)
  text(100, c(-2, 0, 2, 4, 6, 8, 10)[i], labels = levels(x$STUDY)[i],
     col = (1:length(levels(x$STUDY)))[i], adj = 1)
  lines(lowess(x$P[x$STUDY == levels(x$STUDY)[i]],
     x$MATdegC[x$STUDY == levels(x$STUDY)[i]], f = 2/3),
     col = (1:length(levels(x$STUDY)))[i], lwd = 3)
}
text(0, 25, 'curves are lowess fits with bandwidth = 2/3', adj = 0)
dev.off()
#################### 60 Characters Wide ####################

g <- read.table(paste(PATH, 'walton1017.csv', sep = ''),
                header = TRUE, sep = '\t')

decid <- g$EVER
decid[is.na(decid)] <- 0

#################### 60 Characters Wide ####################
# FIG 7
pdf('fig7.pdf', width = 7.5, height = 7.5, family = 'Times', pointsize = 8)
par(mfrow = c(2,2))
plot(g$EVER, g$TEMP, col = hsv(decid, 1, 1), ylim = c(-5, 25),
   cex = 1, xlab = 'Proportion Broadleaf Deciduous Species', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by proportion broadleaf deciduous')

plot(g$X_COORD, g$Y_COORD, pch = 20, col = hsv(decid, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by proportion broadleaf deciduous', xlim = c(-130, -60), ylim = c(20, 55))

plot(g$EVER, g$PERENT, col = hsv(decid, 1, 1), ylim = c(0,1),
   cex = 1, xlab = 'Proportion Broadleaf Deciduous Species', ylab = 'Proportion Entire Leaves', main = 'Scatter plot, colored by proportion broadleaf deciduous')

plot(g$EVER, g$Y_COORD, col = hsv(decid, 1, 1), ylim = c(20, 55),
   cex = 1, xlab = 'Proportion Broadleaf Deciduous Species', ylab = 'North Latitude', main = 'Scatter plot, colored by proportion broadleaf deciduous')
dev.off()

h <- g
h <- h[h$TOTSPECIES > 20,]
hdecid <- h$EVER
#hdecid[is.na(decid)] <- 0

pdf('fig7a.pdf', width = 7.5, height = 7.5, family = 'Times', pointsize = 8)
par(mfrow = c(2,2))
plot(h$EVER, h$TEMP, col = hsv(hdecid, 1, 1), ylim = c(-5, 25),
   cex = 1, xlab = 'Proportion Broadleaf Evergreen Species', ylab = 'Mean Annual Temperature (deg C)', main = 'Scatter plot, colored by proportion broadleaf evergreen')

plot(h$X_COORD, h$Y_COORD, pch = 20, col = hsv(hdecid, 1, 1), cex = 1, xlab = 'West Longitude', ylab = 'North Latitude', main = 'Map, colored by proportion broadleaf evergreen', xlim = c(-130, -60), ylim = c(25, 55))

plot(h$EVER, h$PERENT, col = hsv(hdecid, 1, 1), ylim = c(0,1),
   cex = 1, xlab = 'Proportion Broadleaf Evergreen Species', ylab = 'Proportion Entire Leaves', main = 'Scatter plot, colored by proportion broadleaf evergreen')

plot(h$EVER, h$Y_COORD, col = hsv(hdecid, 1, 1), ylim = c(20, 55),
   cex = 1, xlab = 'Proportion Broadleaf Evergreen Species', ylab = 'North Latitude', main = 'Scatter plot, colored by proportion broadleaf evergreen')
dev.off()

############################################################
# 07-06-14 NEW B&W FIGURES?

pdf('fig1bw.pdf', width = 7.5, height = 3.75, family = 'Times', pointsize = 8)
par(mfrow = c(1,2))
plot(z$Percent.of.entire, z$TEMP, pch = 21,
     cex = T.norm[length(T.norm):1], xlab = 'Proportion Entire Leaves',
     ylab = 'Mean Annual Temperature (deg C)',
     main = 'Scatter plot, plotting symbol size inversely proportional to temperature')
lines(lowess(z$Percent.of.entire, z$TEMP, f = 2/3))
text(0.7, 8, labels = 'lowess bandwidth = 2/3')
plot(z$X_COORD, z$Y_COORD, pch = 21, cex = T.norm[length(T.norm):1],
     xlab = 'West Longitude', ylab = 'North Latitude',
     main = 'Map, plotting symbol size inversely proportional to temperature')
dev.off()