# NRES 798 201501 Lab 6
# Objectives of today's labs
# Implement a  2-way ANCOVA and interpret the results
# Develop a binomial GLM: plot and interpret the results
# Develop a Poisson GLM: plot and interpret the results
# Run a GAM and use the output to evaluate the functional relationship
# Compare the GAM fitted model to GLM and LM model fits

#########################################
# ANCOVA
#########################################
# Number of eggs produced, mass of eggs produced, length of egg chaimber, number of mates for mountain pine beetle
# Two treatments 
# condition: beetles maintained on low, medium, and high food treatments
# dispersal flight: beetles subjected to a dispersal flight prior to reproduction (0 no, 1 yes)
# covariates: body mass, elytra length

# Hypothesis to test:
# 1: Female beetle condition (food prior to reproduction) influences the length of the egg chamber, number of eggs produced, and the total mass of eggs produced
# 2: Females that undergo a dispersal flight produce fewer egg, smaller egg mass, and shorter egg chambers
# 3: The condition of females prior to dispersal influences how dispersal affects their reproductive investment

# TODO: Read in the data set "beetle.txt"
# TODO: Make sure "condit" and "disper" are recognized as factors (i.e. categorical variables, not continuous)
# TODO: Run simple linear models for "chamber","eggnum","eggmass" using "condit","disper","elytra", and "mass" and evaluate the residuals
# Are there any problems with the residuals?

# TODO: Run a two-way ANOVA (no-interaction) and evaluate how beetle condition and dispersal influence each of the reproduction metrics
# TODO: Include the interaction term and determine if the fit of the model is better for each of the reproduction metrics

# Note: the function anova() can be used to evaluate whether or not the inclusion of additional terms improves model fit
# e.g. m1 <- lm(y~x); m2 <- lm(y~x+z)
# anova(m1,m2)  # The test evaluates the two models deviance values

# TODO: Interpret any significant interaction terms
# hint: interaction.plot()

# "elytra" and "mass" (the size and mass of each female beetle) can potentially influence beetle reproductive effort
# TODO: Run a two-way ANCOVA using "elytra" and "mass" and evaluate whether or not the model is improved for each of the three response variables

# TODO: Create a model that includes all two-way interactions (between categorical and continuous variables also)
# hint: what does lm(y~(x+y+z)^3) do?

# TODO: evaluate each term in the above specified full model and reduce the model to the optimal size/complexity
# TODO: Using this optimal model evaluate the three state

#########################################
# Generalized Linear Models (GLM): Binary
#########################################
# Install package "labdsv" that included data to be used in GLM and GAM analysis
install.packages("labdsv")
library(labdsv)

# Dave Roberts Bryce Data
data(bryceveg)
data(brycesite)

attach(bryceveg)
attach(brycesite)

# This data set is from Bryce Canyon and includes two components:
# "brycesite" contain plot information from 160 sites
# "bryceveg" is a presence absence data set for a wide range of plant species
# We will be focusing on the data for Berberis repens ("berrep"), a common shrub, and how its occurrence if influenced by elevation and aspect

# TODO: Inspect the data for any inconsistencies or problems (focus on "berrep")
# TODO: how should the inconsistencies be dealt with.

# TODO: run a linear model on "berrep" with elevation as the independent variable
# TODO: Inspect the residuals.  Does this pattern make sense? What distribution defines the data?

# TODO: run a glm on the berrep occurrence data using elevation as the predictor variable
# Note: the general form for a glm is <- glm (y > 0 ~ x, family = binomial)

# TODO: plot the fitted values from the model
# e.g. plot(elev,fitted(model),xlab='Elevation', ylab='Probability of Occurrence')

# TODO: evaluate the fitted glm model (summary(model)
# TODO: What do the residual deviance and df values tell us?

# TODO: run a glm model that includes both elevation "elev" and aspect "av", and evaluate if the fit is improved
# TODO: use anova(model,test="Chi") to evaluate how each term in the model reduces deviance

#########################################
# Generalized Linear Models (GLM): Poisson
#########################################

# In the "beetle' data set the number of mates each female attracted was also recorded.
# TODO: What distribution likely best represents the number of mates?

# TODO: run linear first order models on beetle$ummates using first only the categorical variables, and then the categorical and continuous variables

# TODO: evaluate the overall fits of these models

# TODO: specify glm models that have the same structure as the liner model used above, using an appropriate linker function

# TODO: evaluate the fit of the glm models compared to the linear models


#########################################
# Generalized Additive Models (GAM)
#########################################
# The occurrence of Berberis repens in the Bryce Canyon data set can also be analysed using a GAM

# TODO: run a gam on "berrep" that just includes elevation as a smoothed predictor variable
# Note: the general form for a gam is <- gam(y > 0 ~ s(x), family=binomial)

# TODO: use summary(model) to examine the gam model fit
# TODO: plot the fitted gam model to visualize the smoothing function
# Note: plot(gam.model,se=TRUE) is the easiest way to do this.
# TODO: how does the gam fitted function visually compare to the glm fitted function?

# TODO: runa a GAM that includes both elevation and aspect 
# TODO: compare the elevation and the elevation and aspect models using anova
# Note: the structure is anova(mod1,mod2,test='Chi')