rm(list = ls())
chap <- read.table(url(
"http://stat.unm.edu/~fletcher/LLM/DATA/CHAPMAN.DAT"),
#"C:\\E-drive\\Books\\ANREG2\\newdata\\chapman.dat",
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\chapman.dat",
sep="",col.names=
c("Case","Ag","S","D","Ch","H","W","y"))
attach(chap)
chap
summary(chap)

#Summary tables
ch = glm(y ~ Ag+S+D+Ch+H+W,family=binomial)
chp=summary(ch)
chp
#anova(ch)

rwt=ch$fit*(1-ch$fit)
yy=log(ch$fit/(1-ch$fit))+(y-ch$fit)/rwt
# If Bin(n_i,p_i)s have n_i different from 1,
# multiply rwt and second term in yy by by n_i

ch1 <- lm(yy ~ Ag+S+D+Ch+H+W,weights=rwt)
ch1p=summary(ch1)
ch1p
anova(ch1)
# Note the agreement between the glm and lm fits!!!

# assign number of best models and number of
# predictor variables.

#install.packages("leaps")
library(leaps)
x <- model.matrix(ch1)[,-1]
nb=3
xp=ch1p$df[1]-1
dfe=length(y)- 1- c(rep(1:(xp-1),each=nb),xp)
g <- regsubsets(x,yy,nbest=nb,weights=rwt)
gg = summary(g)
tt=c(gg$rsq,gg$adjr2,gg$cp,sqrt(gg$rss/dfe))
tt1=matrix(tt,nb*(xp-1)+1,4,
	dimnames = list(NULL,c("R2","AdjR2","Cp","RootMSE")))
tab1=data.frame(tt1,gg$outmat)
tab1