#Supplement 1

# Fricker, G.A., Wolf, J.A., Saatchi, S.S., Gillespie, T.W. 
#Predicting spatial variations of tree species richness in tropical forests from high resolution remote sensing. Ecological Applications.
# File name: 2015_02_03_ols_gls_optimization_Fig5.r
#is the r-code for performing the Generalize Least Squares (GLS) spatial model optimization.  
# The purpose of this code is to demonstrate how GLS spatial models are optimized and also to create Figure 5.

######GLS Models
#Selected variables from the OLS regression are tested as a GLS model with utm coordinates
#   *note: this is only the variables used to predict all stems (>1cm dbh), this process was performed for each discrete tree size class
    #Some remote sensing had contradicting GLS and OLS model fits and were excluded (from the prediction process and this script)

#This code shows how the GLS spatial model is optimized for each remote sensing variable:

#GLS optimization
#three varigrom models are fit: Spherical, Gaussian, Exponential
#Lowest AIC score is the optimized model.  
#If the slope of the OLS (grey) model fit and the GLS (black) model fit are the same sign (postive or negative), 
#the GLS model does not show strong spatial autocorrelation and can be considered valid for a non-spatial OLS model

#INPUT TABLE - all stems > 10 mm dbh (original analysis performed in mm)
d10 = read.csv(file = 'D:\\bci_publication\\r\\in\\gls\\10.csv', header=T)

#In this script, the variables are:
#twi_sd (Topographic Wetness Index -standard deviation)
#curv (mean terrain curvature)
#mch_sd (Mean Canopy Height - standard deviation)
#evi_sd (Enhanced Vegetation Index - standard deviation)

#OUTPUT DIRECTORY (to print pdf)
setwd("D:\\bci_publication\\r\\out")

#nlme: Linear and Nonlinear Mixed Effects Models
#http://cran.r-project.org/web/packages/nlme/index.html
library(nlme)

# GLS models
#Species Richness
#All Stems (greater than 10 mm DBH) 

#Figure 5:
# open pdf to write to - does not show up on screen. saved directly to file
pdf("2014_OLS_GLS_d10_EA.pdf", width=8, height=8)
par(mfrow=c(2,2))

#twi_sd  
##GLS Optimization
SML_gls_twi_sd_spher = gls(sr_all ~ twi_sd,data = d10,correlation = corSpher(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_twi_sd_spher)#R. AIC= 363.3605
SML_gls_twi_sd_gau = gls(sr_all ~ twi_sd,data = d10,correlation = corGaus(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_twi_sd_gau)#R. AIC=362.0129
SML_gls_twi_sd_exp = gls(sr_all ~ twi_sd,data = d10,correlation = corExp(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_twi_sd_exp) #R. AIC=362.3698

olsfit = lm(sr_all ~ twi_sd,data = d10)
summary(olsfit)

# Lowest AIC GLS model
palette(c("yellow3","chartreuse","darkgreen"))
plot(d10$sr_all~d10$twi_sd, xlab='Topographic Wetness Index (standard deviation)', ylab='species/ha', pch=20,cex=1.5,col=as.factor(d10$forest_age), main = "1.0 ha")
fit100= lm(d10$sr_all~d10$twi_sd)
abline(fit100,lwd=2,col = "grey") #OLS FIT)
glsfit = SML_gls_twi_sd_gau
abline(glsfit,lwd=2,col = "black") #GLS FIT)

r2 = round(summary(fit100)$r.squared, 2)
my.p = summary(fit100)$coefficients[2,4]
mylabel = bquote(italic(R)^2 == .(format(r2, digits =2)))
rp = vector('expression',2)
rp[1] = substitute(expression(italic(R)^2 == r2), 
                   list(r2 = format(r2,dig=2)))[2]
rp[2] = substitute(expression(italic(p) == my.p), 
                   list(my.p = format(my.p, digits = 2)))[2]
legend('bottomright', legend = rp, bty = 'n')

#curv
##GLS Optimization
SML_gls_curv_spher = gls(sr_all ~ curv,data = d10,correlation = corSpher(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_curv_spher)#R. AIC= 365.6982
SML_gls_curv_gau = gls(sr_all ~ curv,data = d10,correlation = corGaus(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_curv_gau)#R. AIC= 363.9268
SML_gls_curv_exp = gls(sr_all ~ curv,data = d10,correlation = corExp(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_curv_exp) #R. AIC=364.394

olsfit = lm(sr_all ~ curv,data = d10)
summary(olsfit)

# Lowest AIC GLS model
palette(c("yellow3","chartreuse","darkgreen"))
plot(d10$sr_all~d10$curv, xlab='Terrain curvature', ylab='species/ha', pch=20,cex=1.5,col=as.factor(d10$forest_age), main = "1.0 ha")
fit100= lm(d10$sr_all~d10$curv)
abline(fit100,lwd=2,col = "grey") #OLS FIT)
glsfit = SML_gls_curv_gau
abline(glsfit,lwd=2,col = "black") #GLS FIT
r2 = round(summary(fit100)$r.squared, 2)
my.p = summary(fit100)$coefficients[2,4]
mylabel = bquote(italic(R)^2 == .(format(r2, digits =2)))
rp = vector('expression',2)
rp[1] = substitute(expression(italic(R)^2 == r2), 
                   list(r2 = format(r2,dig=2)))[2]
rp[2] = substitute(expression(italic(p) == my.p), 
                   list(my.p = format(my.p, digits = 2)))[2]
legend('bottomright', legend = rp, bty = 'n')

#mch_sd
##GLS Optimization
SML_gls_mch_sd_spher = gls(sr_all ~ mch_sd,data = d10,correlation = corSpher(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_mch_sd_spher)#R. AIC= 364.913
SML_gls_mch_sd_gau = gls(sr_all ~ mch_sd,data = d10,correlation = corGaus(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_mch_sd_gau)#R. AIC=365.2666
SML_gls_mch_sd_exp = gls(sr_all ~ mch_sd,data = d10,correlation = corExp(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_mch_sd_exp) #R. AIC=365.9774

olsfit = lm(sr_all ~ mch_sd,data = d10)
summary(olsfit)

# Lowest AIC GLS model
palette(c("yellow3","chartreuse","darkgreen"))
plot(d10$sr_all~d10$mch_sd, xlab='Mean Canopy Height (standard deviation)', ylab='species/ha', pch=20,cex=1.5,col=as.factor(d10$forest_age), main = "1.0 ha")
fit100= lm(d10$sr_all~d10$mch_sd)
abline(fit100,lwd=2,col = "grey") #OLS FIT)
glsfit = SML_gls_mch_sd_spher
abline(glsfit,lwd=2,col = "black") #GLS FITr2 = round(summary(fit100)$r.squared, 2)
my.p = summary(fit100)$coefficients[2,4]
mylabel = bquote(italic(R)^2 == .(format(r2, digits =2)))
rp = vector('expression',2)
#rp[1] = substitute(expression(italic(R)^2 == r2),
rp[1] = substitute(expression(italic(R)^2 == 0.15),
                   list(r2 = format(r2,dig=2)))[2]
rp[2] = substitute(expression(italic(p) == my.p), 
                   list(my.p = format(my.p, digits = 2)))[2]
legend('bottomright', legend = rp, bty = 'n')

#evi_sd
##GLS Optimization
SML_gls_evi_sd_spher = gls(sr_all ~ evi_sd,data = d10,correlation = corSpher(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_evi_sd_spher)#R. AIC= 375.5839
SML_gls_evi_sd_gau = gls(sr_all ~ evi_sd,data = d10,correlation = corGaus(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_evi_sd_gau)#R. AIC=374.8626
SML_gls_evi_sd_exp = gls(sr_all ~ evi_sd,data = d10,correlation = corExp(form = ~utm_e+utm_n,nugget=T))
summary(SML_gls_evi_sd_exp) #R. AIC=375.2556

olsfit = lm(sr_all ~ evi_sd,data = d10)
summary(olsfit)

# Lowest AIC GLS model
palette(c("yellow3","chartreuse","darkgreen"))
plot(d10$sr_all~d10$evi_sd, xlab='Enhanced Vegetation Index (standard deviation)', ylab='species/ha', pch=20,cex=1.5,col=as.factor(d10$forest_age), main = "1.0 ha")
fit100= lm(d10$sr_all~d10$evi_sd)
abline(fit100,lwd=2,col = "grey") #OLS FIT)
glsfit = SML_gls_evi_sd_gau
abline(glsfit,lwd=2,col = "black") #GLS FIT
r2 = round(summary(fit100)$r.squared, 2)
my.p = summary(fit100)$coefficients[2,4]
mylabel = bquote(italic(R)^2 == .(format(r2, digits =2)))
rp = vector('expression',2)
rp[1] = substitute(expression(italic(R)^2 == r2), 
                   list(r2 = format(r2,dig=2)))[2]
rp[2] = substitute(expression(italic(p) == my.p), 
                   list(my.p = format(my.p, digits = 2)))[2]
legend('bottomright', legend = rp, bty = 'n')

# turn output device off
dev.off()