donner=read.table("donner.txt")
survive=donner[,3]
age=donner[,1]
sex=donner[,2]

### Table sex*survive
table=as.matrix(table(sex,survive))
contingency_table=list(Frequency=table,Expected=chisq.test(table)$expected)
contingency_table
chisq.test(table, correct=FALSE) ## chisq. test of independence 
LM=lm(survive~age) ## linear regression
summary(LM)
LMANOVA=anova(LM) ## anova
LMANOVA

### Plot Survive*Age

plot(age,survive,xlim=c(15,70),ylim=c(-1.5,2.0),main="survive v.s. age")
abline(LM,col="red")
abline(confint(LM)[1],confint(LM)[2],col="green")
abline(confint(LM)[3],confint(LM)[4],col="purple")
### Plot Predicted*Age
plot(age,fitted(LM),main="Predicted v.s. Age")

### Q-Q plot
qqnorm(residuals(LM),main="Q-Q Plot")

### Studentized residuals v.s Observation
plot(rstudent(LM),main="Studentized residual v.s. observation")
abline(h=0)

###----------fitting logistic regression survive~age
result=glm(survive~age,family=binomial("logit"))
summary(result)

## to get the specific coefficient, this command will produce a vector
coefficients(result)
## to access the estimated slope and turn it into the odds-ratio
exp(coefficients(result)[2])

confint(result) ## confidence interval for parameters
exp(confint(result)) ## exponentiate to get on the odds-scale

### Diagnostics Measures

lm.influence(result)

### Pearson Residuals v.s. observation 
plot(residuals(result,type="pearson"),main="pearson residual plot")

### Deviance Residuals v.s. observation
plot(residuals(result,type="deviance"),main="deviance residual plot")

### Hat Diagonal Plot
plot(hatvalues(result),ylab="H",xlab="Case Number Index")

### Intercept DfBeta Plot
plot(dfbetas(result)[,1],ylab="DFBETA0",xlab="Case Number Index")

### Intercept DfBeta Plot
plot(dfbetas(result)[,2],ylab="DFBETA1",xlab="Case Number Index")

### Table age*survive
table=as.matrix(table(age,survive))
contigency_table=list(Frequency=table,Expected=chisq.test(table)$expected,Percent=prop.table(table),RowPct=prop.table(table,1),ColPct=prop.table(table,2))
contigency_table
chisq.test(table)

###----------fitting logistic regression survive~age+sex
result=glm(survive~age+sex,family=binomial("logit"))
summary(result)
confint(result) ## confidence interval for the parameters 

### Diagnostics Measures

lm.influence(result)

### Pearson Residuals v.s. observation 
plot(residuals(result,type="pearson"),main="pearson residual plot")

### Deviance Residuals v.s. observation
plot(residuals(result,type="deviance"),main="deviance residual plot")

### Hat Diagonal Plot
plot(hatvalues(result),ylab="H",xlab="Case Number Index")

### Intercept DfBeta Plot
plot(dfbetas(result)[,1],ylab="DFBETA0",xlab="Case Number Index")

### Intercept DfBeta Plot
plot(dfbetas(result)[,2],ylab="DFBETA1",xlab="Case Number Index")

###----------fitting logistic regression survive~age+sex+age*sex
donner=as.data.frame(donner)
sort(donner,c(donner$V1,donner$V2))
result=glm(survive~age+sex+age*sex,family=binomial("logit"))
out=data.frame(survive,age,sex,pi=result$fitted)
out

### Sort by sex age and Plot by sex group
group1=(out[which(sex==0),])[sort.int(age[which(sex==0)],index=TRUE)$ix,]
group2=(out[which(sex==1),])[sort.int(age[which(sex==1)],index=TRUE)$ix,]
plot(group1$age,group1$pi,col="red",type="l",xlim=c(10,70),ylim=c(0,1),ylab="Estimated Probability",xlab="age")
lines(group2$age,group2$pi,col="blue",type="c")

##what if we have two categorical predictors and want to add an interaction term?