
```{r }
library(DAAG)
data(frogs)
plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1],
     xlab="Meters east of reference point", ylab="Meters north")
library(DAAG)
data(spam7)
head(spam7)
spam = spam7
m1 = glm(pres.abs ~ altitude + distance + NoOfPools + NoOfSites + avrain + meanmin + meanmax, family = binomial, data = frogs)
summary(m1)
x0 = c(1, 1700, 400, 30, 8, 150, 4, 16)
sum(x0 * m1$coefficients)
exp(sum(x0 * m1$coefficients))
exp(sum(x0 * m1$coefficients))/(1 + exp(sum(x0 * m1$coefficients)))
m1$fitted.values
i = 45
rrg = c(1, frogs$altitude[i], frogs$distance[i], frogs$NoOfPools[i], frogs$NoOfSites[i], frogs$avrain[i], frogs$meanmin[i], frogs$meanmax[i])
#or use model.matrix(m1)[i,]
eta = sum(rrg*m1$coefficients)
prr = exp(eta)/(1 + exp(eta))
#This agrees with m1$fitted.values[i]
c(prr, m1$fitted.values[i])
plot(frogs$pres.abs, m1$fitted.values)
```

Checking the unusual points:

```{r }
sel <- frogs$pres.abs==0 & m1$fitted >.7
plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1], col=sel+1, xlab="Meters east of reference point", ylab="Meters north") 
```



```{r }
sel <- frogs$pres.abs==1 & m1$fitted <.2
plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1], col=sel+1, xlab="Meters east of reference point", ylab="Meters north") 
thr = 0.5
y = frogs$pres.abs
yhat = m1$fitted.values
a <- sum((!y) & (yhat<=thr))
b <- sum((!y) & (yhat>thr))
c <- sum((y) & (yhat<=thr)) 
d <- sum((y) & (yhat>thr))
c(a, b, c, d)
thr = 0.3
y = frogs$pres.abs
yhat = m1$fitted.values
a <- sum((!y) & (yhat<=thr))
b <- sum((!y) & (yhat>thr))
c <- sum((y) & (yhat<=thr)) 
d <- sum((y) & (yhat>thr))
c(a, b, c, d)
conf <- matrix(0,nrow=21,ncol=5)
colnames(conf) <- c("thr","a","b","c","d")
conf[,1] <- seq(0,1,by=.05)
y <- frogs$pres.abs
yhat <- m1$fitted.values
for ( i in 1:21)
  {
    a <- sum((!y) & (yhat<=conf[i,1]))
    b <- sum((!y) & (yhat>conf[i,1]))
    c <- sum((y) & (yhat<=conf[i,1])) 
    d <- sum((y) & (yhat>conf[i,1]))
    conf[i,2:5] <- c(a,b,c,d)
  }
conf
plot(conf[,1],conf[,3]+conf[,4], xlab="threshold", ylab="b+c")
library(DAAG)
data(spam7)
head(spam7)
summary(spam)
par(mfrow=c(3,2))
for (i in 1:5)
  hist(spam[,i],main="",xlab=names(spam)[i],breaks=10000)
par(mfrow=c(1,1))
pairs(spam, cex=0.5)
s = 0.001 #correction for zero as log(0) is infinity. 
pairs(~log(crl.tot)+log(dollar+s) +log(bang+s)+log(money+s)+log(n000+s)+log(make+s) + yesno,       data=spam, cex=.5)
spam.glm <- glm(yesno~ log(crl.tot) + log(dollar+s) + log(bang+s)
                +log(money+s) + log(n000+s) + log(make+s),
                family=binomial, data=spam)
summary(spam.glm)
spam.lm <-  lm(as.numeric(yesno=="y")~ log(crl.tot) + log(dollar+s)
               + log(bang+s) +log(money+s) + log(n000+s) + log(make+s)
               ,data=spam)
summary(spam.lm)
plot(spam.lm$fitted.values, spam.glm$fitted.values)
abline(c(0, 1), col = "red")
confusion <- function (y, yhat, thres)
  {
    n <- length(thres)
    conf <- matrix(0,length(thres),ncol=4)
    colnames(conf) <- c("a","b","c","d")
    for ( i in 1:n)
      {
        a <- sum((!y) & (yhat<=thres[i]))
        b <- sum((!y) & (yhat>thres[i]))
        c <- sum((y) & (yhat<=thres[i])) 
        d <- sum((y) & (yhat>thres[i]))
        conf[i,] <- c(a,b,c,d)
      }
    return(conf)
  }
v <- seq(0.05,0.95,by=0.05)
y <- as.numeric(spam$yesno=="y")
glm.conf <- confusion(y,spam.glm$fitted,v)
glm.conf
lm.conf <- confusion(y,spam.lm$fitted,v)
lm.conf
matplot(v,cbind((glm.conf[,"b"]+glm.conf[,"c"])/4601, (lm.conf[,"b"]+lm.conf[,"c"])/4601), xlab="threshold", ylab="b+c",type="l")
legend(.8, .4, lty=1:2, col=1:2,c("glm","lm"))
spam.lm1 = lm(as.numeric(yesno=="y")~ crl.tot + dollar
               + bang + money + n000 + make, data=spam)
summary(spam.lm1)
summary(spam.lm)
```

Precision and recall for the GLM:

```{r }
cbind(glm.conf[,"d"]/(glm.conf[,"b"]+glm.conf[,"d"]),glm.conf[,"d"]/(glm.conf[,"c"]+glm.conf[,"d"]) )
plot(glm.conf[,"d"]/(glm.conf[,"b"]+glm.conf[,"d"]),glm.conf[,"d"]/(glm.conf[,"c"]+glm.conf[,"d"]) )
matplot(x = cbind(glm.conf[,"d"]/(glm.conf[,"b"]+glm.conf[,"d"]),
              lm.conf[,"d"]/(lm.conf[,"b"]+lm.conf[,"d"])),
        y = cbind(glm.conf[,"d"]/(glm.conf[,"c"]+glm.conf[,"d"]),
              lm.conf[,"d"]/(lm.conf[,"c"]+lm.conf[,"d"])),
        xlab="Precision", ylab="Recall",type="l")
legend(.8, .4, lty=1:2, col=1:2,c("glm","lm"))
```



```{r }
library(DAAG)
data(frogs)
m1 = glm(pres.abs ~ altitude + distance + NoOfPools + NoOfSites + avrain + meanmin + meanmax, family = binomial, data = frogs)
summary(m1)
```



```{r }
deviance(m1)
rd = sum(frogs$pres.abs * (-2) * log(m1$fitted.values)) + sum((1 - frogs$pres.abs)*(-2) * log(1 - m1$fitted.values))  
rd
m2 = glm(pres.abs ~ altitude + distance + NoOfSites + avrain + meanmin + meanmax, family = binomial, data = frogs) 
summary(m2)
```

AIC

```{r }
m1$aic
m1$deviance + 2*(7+1)
```



```{r }
step(m1)
step(m1, direction = "both")
```

