# R code for class 2/7/2019

library(faraway)
data(orings)

plot(damage/6~temp,orings,xlim=c(25,85),ylim=c(0,1),xlab='Temperature',
ylab='Proportion of Damaged Rings')

lmod=glm(cbind(damage,6-damage)~temp,family=binomial,orings)
sumary(lmod)

x=seq(25,85,1)
lines(x,ilogit(11.6630-0.2162*x))

# alternative model including squared temperature
temp2=with(orings,temp^2)
lmod2=glm(cbind(damage,6-damage)~temp+temp2,family=binomial,orings)
sumary(lmod2)

x=seq(25,85,1)
lines(x,ilogit(53.091126-1.531653*x+0.010293*x*x),col='red')

ilogit(11.6630-0.2162*31)
ilogit(53.091126-1.531653*31+0.010293*31*31)

# formal test of model 2 against model 1
anova(lmod,lmod2)
pchisq(2.3864,1,lower=F)

# deviance of 2.3864, df=1, pvalue=0.12
# Addition of quadratic term does NOT significantly improve the model



pchisq(deviance(lmod),df.residual(lmod),lower=F)

pchisq(38.9,22,lower=F)

deviance(lmod)
# pearson X^2 statistic
sum(residuals(lmod,type='pearson')^2)

pchisq(sum(residuals(lmod,type='pearson')^2),df.residual(lmod),lower=F)

# estimate overdispersion
sum(residuals(lmod,type='pearson')^2)/df.residual(lmod)

# alternatively, do this:
lmodod=glm(cbind(damage,6-damage)~temp,family=quasibinomial,orings)

summary(lmodod)

# end of class 2/7/2019


data(mammalsleep)
mammalsleep$pdr=with(mammalsleep,dream/sleep)
summary(mammalsleep$pdr)

modl=glm(pdr~log(body)+log(brain)+log(lifespan)+log(gestation)+predation+exposure+danger,
family=quasibinomial,mammalsleep)
summary(modl)

# do backward selection dropping one variable at a time

modl=glm(pdr~log(body)+log(brain)+log(lifespan)+log(gestation)+exposure+danger,
family=quasibinomial,mammalsleep)
summary(modl)

modl=glm(pdr~log(body)+log(brain)+log(lifespan)+log(gestation)+danger,
family=quasibinomial,mammalsleep)
summary(modl)

modl=glm(pdr~log(body)+log(lifespan)+log(gestation)+danger,family=quasibinomial,mammalsleep)
summary(modl)

modl=glm(pdr~log(body)+log(lifespan)+danger,family=quasibinomial,mammalsleep)
summary(modl)

# diagnostic plots page 63

ll=row.names(na.omit(mammalsleep[,c(1,6,10,11)]))
halfnorm(cooks.distance(modl),labs=ll)

plot(predict(modl),residuals(modl,type='pearson'),xlab='Linear Predictor',
ylab='Pearson Residuals',pch=20)


# omit Asian elephant and refit model
mammalsleep[5,11]=NA
mammalsleep

modl2=glm(pdr~log(body)+log(lifespan)+danger,family=quasibinomial,mammalsleep)
summary(modl2)

# restore Asian elephant and do beta model fit of section 3.6

data(mammalsleep)
mammalsleep$pdr=with(mammalsleep,dream/sleep)
library(mgcv)
modb=gam(pdr~log(body)+log(lifespan),family=betar(),mammalsleep)
summary(modb)


plot(modl$fitted,modb$fitted,pch=20)
lines(c(-1000,1000),c(-1000,1000))