setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
hemodata<-read.table(file="Hemophilia-dat.txt")
colnames(hemodata)<-c("group", "AHF activity", "AHF-like antigen")
normalgroup<-hemodata[hemodata[,1]==1,]
carriergroup<-hemodata[hemodata[,1]==2,]
plot(normalgroup[,2],normalgroup[,3],xlim=c(-0.7,0.2),ylim=c(-0.4,0.4),xlab="log(AHF activity)", ylab="log(AHF-like antigen)")
points(carriergroup[,2],carriergroup[,3],pch=3)

## Performing the linear classification

library(MASS)
train<-sample(1:dim(hemodata)[1],60)
z<-lda(group~ ., hemodata,prior=c(0.5,0.5),subset=train)
predgroups<-predict(z,hemodata[-train,])$class
truth<-hemodata[-train,1]
predgroups==truth
table(truth,predgroups)

## Linear classifier by “hand”

traindata<-hemodata[train,]
groupmeans<-aggregate(traindata[,2:3], list(traindata$group), mean)
group1<-(traindata$group==1)
group2<-(traindata$group==2)
S1<-cov(traindata[group1,2:3])
S2<-cov(traindata[group2,2:3])
n1<-length(group1)
n2<-length(group2)
Spooled<-((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)
invSn<-solve(Spooled)
ahat<-as.numeric(groupmeans[1,2:3]-groupmeans[2,2:3])%*%invSn
mhat<-ahat%*%as.numeric(groupmeans[1,2:3]+groupmeans[2,2:3])/2
testdata<-sapply(hemodata[-train,2:3],as.numeric)
linearrule<-as.vector(testdata%*%t(ahat))-mhat
predictions<-1+(linearrule<0)*1

normalgroup0<-traindata[traindata[,1]==1,]
carriergroup0<-traindata[traindata[,1]==2,]
plot(normalgroup0[,2],normalgroup0[,3],xlim=c(-0.7,0.2),ylim=c(-0.4,0.4),xlab="log(AHF activity)", ylab="log(AHF-like antigen)")
points(carriergroup0[,2],carriergroup0[,3],pch=3)
abline(a=mhat/ahat[2],b=-ahat[1]/ahat[2],col=2)