---
title: "ADA2: Homework 10, Ch 11, Logistic Regression"
author: "Name Here"
date: "mm/dd/yyyy"
output:
pdf_document:
number_sections: yes
toc: yes
html_document:
toc: true
number_sections: true
code_folding: show
---
```{R, echo=FALSE}
# I set some GLOBAL R chunk options here.
# (to hide this message add "echo=FALSE" to the code chunk options)
knitr::opts_chunk$set(comment = NA, message = FALSE, warning = FALSE, width = 100)
knitr::opts_chunk$set(fig.align = "center", fig.height = 4, fig.width = 6)
knitr::opts_chunk$set(cache = TRUE, autodep=TRUE) #$
```
# [Galapagos Island Species Data](http://www.statsci.org/data/general/galapagos.html)
The Galapagos Islands about 600 miles off the coast of Ecuador provide an
excellent laboratory for studying the factors that influence the development
and survival of different life species. They were the site of much of Charles
Darwin's original research leading later to publication of his "Origin of
Species". Descending from a few stranded ancestors and cut off from the rest of
the world, the Galapagos animals offer much more obvious proofs of the fact of
evolution than can be seen in the more intricate complexities of life in most
environments. Darwin wrote:
_The natural history of these islands is eminently curious, and well deserves
attention. Most of the organic productions are aboriginal creations, found
nowhere else; there is even a difference between the inhabitants of the
different islands; yet all show a marked relationship with those of America,
though separated from that continent by an open space of ocean, between 500 and
600 miles in width. The archipelago is a little world in itself, or rather a
satellite attached to America, whence it has derived a few stray colonists and
has received the general character of its indigenous productions. Considering
the small size of the islands, we feel the more astonished at the number of
their aboriginal beings, and at their confined range. Seeing every height
crowned with its crater, and the boundaries of most of the lava-streams still
distinct, we are led to believe that within a period geologically recent the
unbroken ocean was here spread out. Hence, both in space and time, we seem to
be brought somewhere near to that great fact---that mystery of mysteries---the
first appearance of new beings on earth._
And from elsewhere in Darwin's diary:
_I never dreamed that islands 50 or 60 miles apart, and most of them in sight of
each other, formed of precisely the same rocks, placed under a quite similar
climate, rising to a nearly equal height, would have been differently tenanted...
It is the circumstance that several of the islands possess their own species
of the tortoise, mocking-thrush, finches and numerous plants, these species
having the same general habits, occupying analogous situations, and obviously
filling the same place in the natural economy of the archipelago, that strikes
me with wonder._
M.P. Johnson and P.H. Raven, "Species number and endemism: The Galapagos
Archipelago revisited", Science, 179, 893-895 (1973), have presented data
giving the number of plant species and related variables for 29 different
islands. __Counts are given for both the total number of species and the number
of species that occur only in the Galapagos (the endemics).__
Elevations for Baltra and Seymour obtained from web searches. Elevations for
four other small islands obtained from large-scale maps.
```
Variable Description
Island Name of Island
Plants Number of plant species
PlantEnd Number of endemic plant species
Finches Number of finch species
FinchEnd Number of endemic finch species
FinchGenera Number of finch genera
Area Area (km^2)
Elevation Maximum elevation (m)
Nearest Distance from to nearest island (km)
StCruz Distance to Santa Cruz Island (km)
Adjacent Area of adjacent island (km^2)
```
__Goal:__
To build a model to predict the proportion of endemic plants on an island based
on the island characteristics.
```{R}
fn.data <- "http://statacumen.com/teach/ADA2/homework/ADA2_HW_10_Galapagos.txt"
gal <- read.table(fn.data, skip=28, header=TRUE)
# to label plots in plots
gal$id <- 1:nrow(gal)
```
Artificially remove the responses for three islands to predict later.
```{R}
# list of islands to predict
island.pred.list <- c("Gardner2", "Santa.Fe", "Wolf")
## capture the observed probabilities
gal.pred.true <- subset(gal, Island %in% island.pred.list)
# Set these islands with missing response variables
gal[(gal$Island %in% island.pred.list), c("Plants", "PlantEnd")] <- NA
```
Compute the observed proportion and empirical logits of endemic plants on each island.
```{R}
# observed proportions
gal$p.hat <- gal$PlantEnd / gal$Plants
# emperical logits
gal$emp.logit <- log((gal$p.hat + 0.5/gal$Plants) /
(1 - gal$p.hat + 0.5/gal$Plants))
```
__Data modifications__
```{R}
## RETURN HERE TO SUBSET AND TRANSFORM THE DATA
# gal <- subset(gal, !(X > Z)) # exclude observations of this type
# gal <- gal[!(gal$X > Z),] # exclude observations of this type
### SOLUTION
```
## __(2 p)__ Interpret plot of observed proportions against predictors
```{R, fig.height = 3, fig.width = 8}
# Create plots for proportion endemic for each variable
library(reshape2)
gal.long <- melt(subset(gal, select = c(Island, id, p.hat, emp.logit, Area, Elevation, Nearest, StCruz, Adjacent)), id.vars = c("Island", "id", "p.hat", "emp.logit"))
# Plot the data using ggplot
library(ggplot2)
p <- ggplot(gal.long, aes(x = value, y = p.hat, label = id))
p <- p + geom_hline(yintercept = c(0,1), alpha = 0.25)
p <- p + geom_text(hjust = 0.5, vjust = -0.5, alpha = 0.25, colour = 2)
p <- p + geom_point()
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + facet_wrap( ~ variable, scales = "free_x", nrow = 1)
print(p)
# Plot the data using ggplot
library(ggplot2)
p <- ggplot(gal.long, aes(x = value, y = emp.logit, label = id))
p <- p + geom_text(hjust = 0.5, vjust = -0.5, alpha = 0.25, colour = 2)
p <- p + geom_point()
p <- p + facet_wrap( ~ variable, scales = "free_x", nrow = 1)
print(p)
```
__Comment__ on how the proportion of endemic plants depends on each variable (in terms of increase/decrease).
Also, __interpret__ the plot regarding whether the empirical logits appear linear (or any trends).
Note that the marginal empirical logit plots _do not_ have to be linear, but the model in 6-dimensional space should be roughly "linear".
Indicate whether any observations are gross outliers and should be dropped, and
whether variables are obvious candidates for transformation.
Then, drop outliers, transform and update your comments, and repeat until satisfied.
### Solution
__Proportion plots: __
[answer]
__Empirical logit plots:__
[answer]
## __(0 p)__ Predictor scatterplot matrix
For further information, the relationship between predictors is plotted.
```{R, fig.height = 6, fig.width = 6}
# relationships between predictors
library(ggplot2)
library(GGally)
p <- ggpairs(subset(gal, select = c(Area, Elevation, Nearest, StCruz, Adjacent)))
print(p)
```
## __(1 p)__ Fit a logistic regression model, interpret deviance lack-of-fit
Fit a logistic model relating the probability of endemic plants to the predictors.
Decide which predictors to use.
```{R}
### SOLUTION
# Don't include both Area and Elevation, since highly correlated
glm.g.aensj <- glm(cbind(PlantEnd, Plants - PlantEnd) ~ Area + Nearest + StCruz + Adjacent, family = binomial, gal)
# Test residual deviance for lack-of-fit (if > 0.10, little-to-no lack-of-fit)
dev.p.val.full <- 1 - pchisq(glm.g.aensj$deviance, glm.g.aensj$df.residual)
dev.p.val.full
## Stepwise selection
# option: trace = 0 doesn't show each step of the automated selection
glm.g.aensj.red.AIC <- step(glm.g.aensj, direction="both", trace = 0)
# the anova object provides a summary of the selection steps in order
glm.g.aensj.red.AIC$anova
coef(glm.g.aensj.red.AIC)
# Test residual deviance for lack-of-fit (if > 0.10, little-to-no lack-of-fit)
dev.p.val.red <- 1 - pchisq(glm.g.aensj.red.AIC$deviance, glm.g.aensj.red.AIC$df.residual)
dev.p.val.red
```
Look at the residual deviance lack-of-fit statistic for the __full model__.
__Is there__ evidence of any gross deficiencies with the model?
How about for the __reduced model__?
__(Regardless of lack of fit result, continue with the assignment.
This is a realistic example and not everything may work out nicely.)__
### Solution
[answer]
Full:
Reduced:
## __(2 p)__ Interpret logistic regression coefficients
Which variables appear to be a useful predictor of the probability of endemic plants?
__Interpret__ the hypothesis test(s).
```{R}
summary(glm.g.aensj.red.AIC)
```
### Solution
[answer]
`XXX` is significant and it's negative coefficient of
`r signif(glm.g.aensj.red.AIC$coefficients[2], 3)` suggests that as `XXX`
increases, the proportion of endemic plants decreases.
etc.
## __(1 p)__ Write model equation
__Provide__ an equation relating the fitted probability of endemic plants to
the selected predictor variables on the probability scale.
### Solution
[answer]
The logistic equation is
$$
\tilde{p}
=
\frac{ \exp( XXX )}
{ 1 + \exp( XXX ) }
.
$$
## __(2 p)__ Plot the fitted probabilities as a function of the selected predictor variables.
Note that if there are more than one predictor,
these plots may be jagged.
That is because the predictions we're getting for each observation is
conditional on the _other_ variables in the model.
This is not an ideal way of plotting the data and model,
but it will give some sense of whether the predictions are close to the observed proportions.
```{R}
# put the fitted values in the data.frame
#gal$fitted.values <- glm.g.aensj.red.AIC$fitted.values
# predict() uses all the Load values in dataset, including appended values
fit <- predict(glm.g.aensj.red.AIC, subset(gal, select = c(Area, Elevation, Nearest, StCruz, Adjacent)), type = "link", se.fit = TRUE)
gal$fit <- fit$fit
gal$se.fit <- fit$se.fit
# CI for fitted values
gal <- within(gal, {
# added "fitted" to make predictions at appended Load values
fitted = exp(fit) / (1 + exp(fit))
fit.lower = exp(fit - 1.96 * se.fit) / (1 + exp(fit - 1.96 * se.fit))
fit.upper = exp(fit + 1.96 * se.fit) / (1 + exp(fit + 1.96 * se.fit))
})
```
### Solution
[answer]
```
{R, fig.height = 4, fig.width = 10}
library(ggplot2)
p1 <- ggplot(gal, aes(x = XXX, y = p.hat))
# predicted curve and point-wise 95% CI
p1 <- p1 + geom_ribbon(aes(x = XXX, ymin = fit.lower, ymax = fit.upper), alpha = 0.2)
p1 <- p1 + geom_line(aes(x = XXX, y = fitted), colour="red")
# fitted values
p1 <- p1 + geom_point(aes(y = fitted), size=2, colour="red")
# observed values
p1 <- p1 + geom_point(size = 2)
p1 <- p1 + scale_y_continuous(limits = c(0,1))
p1 <- p1 + ylab("Probability")
#p1 <- p1 + labs(title = "Observed and predicted probability of endemic")
#print(p1)
library(ggplot2)
p2 <- ggplot(gal, aes(x = XXX, y = p.hat))
# predicted curve and point-wise 95% CI
p2 <- p2 + geom_ribbon(aes(x = XXX, ymin = fit.lower, ymax = fit.upper), alpha = 0.2)
p2 <- p2 + geom_line(aes(x = XXX, y = fitted), colour="red")
# fitted values
p2 <- p2 + geom_point(aes(y = fitted), size=2, colour="red")
# observed values
p2 <- p2 + geom_point(size = 2)
p2 <- p2 + scale_y_continuous(limits = c(0,1))
p2 <- p2 + ylab("Probability")
#p2 <- p2 + labs(title = "Observed and predicted probability of endemic")
#print(p2)
library(gridExtra)
grid.arrange(grobs = list(p1, p2), nrow=1, top = "Observed and predicted probability of endemic")
```
## __(1 p)__ Interpret the prediction with 95% CI at the three islands we didn't use to build the model
Compute the estimated probability of endemic for these islands:
```{R}
island.pred.list
```
__Provide and interpret__ the 95% CIs for this probability.
We have already augmented the data set with the values to predict, so the
`predict()` function above has already done the calculations for us.
### Solution
[answer]
## __(1 p)__ Caveats
What limitations may exist in this analysis?
Do you have reason to expect that model predictions may not be accurate?
### Solution
[answer]