#########################################################
# Script to use GLMs to investigate the relationship between
# douglas fir presence and envrionmental variables
# Author: Jim Graham
# Date: 4/24/2012, updated 2/21/2017
#########################################################

#########################################################
# Read and clearn up the data

TheData <- read.csv("f:/ProjectsModeling/Cali_5Minute/DougFir_CA_Predict 4.csv")

# Remove any NA (null) values, the data should already be loaded
TheData=na.omit(TheData)

boxplot(TheData$AnnualPrecip~TheData$Present, 
        main="Annual Preciptation vs. Presence", 
        xlab="Presence", ylab="Height")

#########################################################
# Explore the data with histograms and plots

# Histogram the height values
hist(TheData$Height,probability=TRUE,main="Tree Height",breaks=400) 

# Histogram the diameter values
hist(TheData$MaxTemp,probability=TRUE,main="Max Temp",breaks=400) 

# Histogram the precipitation values
hist(TheData$AnnualPrecip,probability=TRUE,main="Annual Precip",breaks=400) 

# Histogram the mean temperature values
hist(TheData$MeanTemp,probability=TRUE,main="Mean Temp",breaks=400) 

# Plot the height values against the diameter values
plot(TheData$Height~TheData$MeanTemp,main="Height vs. Mean Temp")

# Plot the height values against the diameter values
plot(TheData$Height~TheData$DEM,main="Height vs. DEM")

# Plot the height values against the diameter values
plot(TheData$Height~TheData$AnnualPrecip,main="Height vs. AnnualPrecip")

##########################################################################
# Linear regressoin

# Add a linear regression line to the plot
RegressionModel<-lm(TheData$Height~TheData$AnnualPrecip)
abline(RegressionModel,col="red")
summary(RegressionModel)

# Create a response vector based on the model
ResponsePrediction=predict(RegressionModel)

# Add the prediction back into the original dataset
TheData2=TheData
TheData2[10]=ResponsePrediction

# Add the response to the plot
points(TheData$AnnualPrecip,ResponsePrediction,  type="p", col="red", lwd=2)

# write out the final results
write.csv(TheData2,"C:\\ProjectsModeling\\Cali_5Minute\\Final.csv")

########################################
# use GLMs with the entire dataset, one predictor

# Plot the data vs. predictors
plot(TheData$AnnualPrecip,TheData$Present,main="Present vs. AnnualPrecip")

# Create the GLM
TheGLM <- glm(TheData$Present~TheData$AnnualPrecip,family=binomial)

# Display the summary statistics
summary(TheGLM)

# Show the plots
plot(TheGLM)

# Get a look at the analysis of deviance
anova(TheGLM,test="Chi")

# create a predicted response
ResponsePrediction=predict(TheGLM)

# plot the result
plot(TheData$AnnualPrecip,TheData$Present,main="Height vs. AnnualPrecip")
points(TheData$AnnualPrecip,fitted(TheGLM),col=2,pch="+")

########################################
# use GLMs with the entire dataset, two predictors

# Create the GLM
TheGLM <- glm(TheData$Present~TheData$AnnualPrecip+TheData$MeanTemp,family=binomial)
summary(TheGLM)
anova(TheGLM,test="Chi")
plot(TheGLM)

# create a predicted response
ResponsePrediction=predict(TheGLM)

# Plot the data vs. annual precip
plot(TheData$Present~TheData$MeanTemp,main="Height vs. MeanTemp")
points(TheData$MeanTemp,fitted(TheGLM),col=2,pch="+")

# Plot the data vs. annual precip
plot(TheData$Present~TheData$AnnualPrecip,main="Height vs. AnnualPrecip")
points(TheData$AnnualPrecip,fitted(TheGLM),col=3,pch="+")