---
title: "logistic regression"
output: 
  html_notebook:
    self_contained: yes
  html_document: default
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
library(foreign)
library(MASS)
donner <- read.csv("~/Dropbox/WorkDocuments/Courses/Stat 224/Data/donner.txt", sep="")
attach(donner)
colors=c("red","blue")
plot(Survived~Age, pch = 16,col=colors[Sex+1])
```
```{r}
fSex=factor(Sex,labels=c("Female","Male"))
ggplot(donner, aes(Age, as.numeric(Survived), color=fSex)) +
  stat_smooth(method="loess", formula=y~x,
              alpha=0.2, size=2, aes(fill=fSex)) +
  geom_point(position=position_jitter(height=0.03, width=0)) +
  xlab("Age") + ylab("Pr (survived)")
```

```{r}
fit<-glm(Survived~Age,family=binomial())
summary(fit)
plot(fit)
```

```{r}
range(Age)
newage=seq(15,65,.1)
predicted_probs=predict(fit,list(Age=newage),type="response")
plot(Survived~Age, pch = 16,col=colors[Sex+1])
lines(newage,predicted_probs)
legend("right",y=c('female','male'),pch = 16,col=colors[1:2])
```
```{r}
ggplot(donner[c("Survived","Age")], aes(Age, as.numeric(Survived))) +
  stat_smooth(method="glm", family=binomial, formula=y~x,
              alpha=0.2, size=2) +
  geom_point(position=position_jitter(height=0.03, width=0)) +
  xlab("Age") + ylab("Pr (survived)")
```
```




```{r}
fit2<-glm(Survived~Age+Sex,family=binomial())
summary(fit2)
female_frame=data.frame(Age=newage,Sex=rep(0,length(newage)))
male_frame=data.frame(Age=newage,Sex=rep(1,length(newage)))
predicted_probs_F=predict(fit2,female_frame,type="response")
predicted_probs_M=predict(fit2,male_frame,type="response")
plot(Survived~Age, pch = 16,col=colors[Sex+1])
lines(newage,predicted_probs_F,col=colors[1])
lines(newage,predicted_probs_M,col=colors[2])
legend("right",y=c('female','male'),pch = 16,col=colors[1:2])
```
```{r}
ggplot(donner, aes(Age, as.numeric(Survived), color=fSex)) +
  stat_smooth(method="glm", family=binomial, formula=y~x,
              alpha=0.2, size=2, aes(fill=fSex)) +
  geom_point(position=position_jitter(height=0.03, width=0)) +
  xlab("Age") + ylab("Pr (survived)")
```

```{r}
data(ships)
head(ships)
?ships
```
`type
type: "A" to "E".

year
year of construction: 1960–64, 65–69, 70–74, 75–79 (coded as "60", "65", "70", "75").

period
period of operation : 1960–74, 75–79.

service
aggregate months of service.

incidents
number of damage incidents.`

```{r}
# exclude ships with 0 months of service
ships2 <- subset(ships, service > 0)
# convert the period and year variables to factors
ships2$year <- as.factor(ships2$year)
ships2$period <- as.factor(ships2$period)

# consider a log-linear model including all the variables
# "offset" means a term with a fixed coefficient of 1
glm1 <- glm(incidents ~ type + year + period, 
    family = poisson(link = "log"), data = ships2, offset = log(service))
summary(glm1)
par(mfrow=c(2,2))
plot(glm1)
```
```{r}
ggplotRegression <- function (fit) {

require(ggplot2)

ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
                     "Intercept =",signif(fit$coef[[1]],5 ),
                     " Slope =",signif(fit$coef[[2]], 5),
                     " P =",signif(summary(fit)$coef[2,4], 5)))
}
```

```{r}
ggplotRegression(fit)
```
