#Program to analyze the data in Table 10.1 for Example 10.1
yd=read.table("C:\\Research\\Book\\Data\\table101.txt")
dose = yd[,1]
died = yd[,2]
count =yd[,3]
alive = count-died
l.dose=log(dose)
yy=cbind(died,alive)
prop=died/count
md1 = yy~l.dose
md2 = prop~l.dose

#Linear regression analysis
lreg = lm(md2)
resid=residuals(lreg); pred=fitted(lreg)
summary(lreg)
anova(lreg)
cbind(dose,died,count,l.dose,prop,pred,resid)

#Logistic regression analysis
gl1 = glm(md1, family=binomial("logit"))
resid1=residuals(gl1); phat1=fitted.values(gl1)
yhat1=count*phat1
anova(gl1)
summary(gl1)
cbind(died,count,yhat1,phat1,resid1)

gl2 = glm(md2, family=binomial("logit"),weights=count)
resid2=residuals(gl2); phat2=fitted.values(gl2)
yhat2=count*phat2
anova(gl2)
summary(gl2)
cbind(died,count,yhat2,phat2,resid2)

###########################################################################

#Program to analyze the data in Table 10.2 as part of Example 10.1
yd=read.table("C:\\Research\\Book\\Data\\table102.txt")
dose = yd[,1]
count= yd[,2]
stat0= yd[,3]
ldose=log10(dose)
status=rep(0,length(stat0))
for (k in 1:length(stat0))
{if (stat0[k]=="died") status[k]=1}
#resp=as.factor(resp)

cbind(dose,ldose,count,status)

#Logistic regression analysis
md = status~ldose
gl1 = glm(md, family=binomial("logit"), weights=count)
anova(gl1, test="Chisq")
summary(gl1)

###########################################################################

#Program to analyze the data in Table 10.3 for Example 10.2
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table103.txt")
grade= yd[,1]
resp = yd[,2]
resp=as.factor(resp)

#Logistic regression analysis
md = resp~grade
gl1 = glm(md, family=binomial("logit"))
anova(gl1, test="Chisq")
summary(gl1)
deviance(gl1)
predicted.prob=predict(gl1,type="response")
cbind(grade,resp,predicted.prob)

#Multiple logistic regression analysis
grade2 = grade^2
md2 = resp~grade+grade2
gl2 = glm(md2, family=binomial("logit"))
anova(gl2, test="Chisq")
summary(gl2)
deviance(gl2)
linearHypothesis(gl2, "grade2=0")

###########################################################################

#Program to analyze the data for Example 10.3
library(car)
yd=read.table("C:\\Research\\Book\\Data\\example103.txt")
logdose=yd[,1]
count= yd[,2]
died = yd[,3]
dose2=logdose^2
alive = count-died
yy=cbind(died,alive)

#Logistic regression analysis
md = yy~logdose
gl1 = glm(md, family=binomial("logit"))
anova(gl1, test="Chisq")
summary(gl1)

#Multiple logistic regression analysis
md2 = yy~logdose+dose2
gl2 = glm(md2, family=binomial("logit"))
anova(gl2, test="Chisq")
summary(gl2)
linearHypothesis(gl2, "dose2=0")

#Logistic regression analysis
md3 = yy~logdose
gl3 = glm(md3, family=binomial("cloglog"))
anova(gl3, test="Chisq")
summary(gl3)

###########################################################################

#Program to analyze the data in Table 10.4 for Example 10.4
library(car)
yd=read.table("C:\\Research\\Book\\Data\\table104.txt")
idno = yd[,1]
cons = yd[,2]
volu = yd[,3]
rate = yd[,4]
inte = volu*rate

#Logistic regression analysis
md1 = cons ~ volu+rate
gl1 = glm(md1, family=binomial("logit"))
anova(gl1, test="Chisq")
summary(gl1)

#Multiple logistic regression analysis
md2 = cons ~ volu+rate+inte
gl2 = glm(md2, family=binomial("logit"))
anova(gl2, test="Chisq")
summary(gl2)
linearHypothesis(gl2, "inte=0")

###########################################################################

#Program to analyze the data in Table 10.5 for Example 10.5
yd=read.table("C:\\Research\\Book\\Data\\table105.txt")
surv=yd[,1]
count=yd[,2]
yr=seq(1938,1947,1); 
year = rep(yr, each=2)
gender=rep(1, length(yr)*2)
x1w=rep(0, length(yr)*2)
sex=rbind(cbind(gender),cbind(x1w))
x2as=cbind(1,0); fac=rep(x2as, length(yr)*2)
y=cbind(year,sex,fac,surv,count)
yv=rbind(y[1:16,],y[18:length(count),])
gender=as.factor(yv[,2]); fac=as.factor(yv[,3])
surv=yv[,4]; count=yv[,5]; year=yv[,1]

#Logistic regression analysis for year with gender and faculty as factors
died = count-surv
yy=cbind(surv,died)
md1 = yy~year+gender+fac
gl1 = glm(md1, family=binomial("logit"))
anova(gl1)
summary(gl1)

#Logistic regression analysis with gender and faculty as factors
md2 = yy~gender+fac
gl2 = glm(md2, family=binomial("logit"))
anova(gl2)
summary(gl2)

#Logistic regression analysis with gender, faculty, and year as factors
#Note that 1938 is taken as the reference year in this analysis
year.n=as.factor(year)
md3 = yy~gender+fac+year.n
#Logistic regression analysis
gl3 = glm(md3, family=binomial("logit"))
anova(gl3)
summary(gl3)

###########################################################################

#Program to analyze the data used in section 10.5 (Arthritis data)
yd=read.table("C:\\Research\\Book\\Data\\section105.txt")
idno =yd[,1]
treat=yd[,2]
gender=yd[,3]
age = yd[,4]
imp = yd[,5]

no=length(imp)
better=rep(0,no); trt=rep(0,no); sex=rep(0,no) 
for (k in 1:no)
{if (imp[k] > 0) better[k]=1
if (treat[k]=="Treat") trt[k]=1
if (gender[k]=="F") sex[k]=1
}

#Multiple logistic regression analysis
md1 = better ~ trt+sex+age
gl1 = glm(md1, family=binomial("logit"))
anova(gl1, test="Chisq")
summary(gl1)