setwd("D:\\study\\572")

library(arm)

#hw question
samara=read.table("samara.txt",header=T)
str(samara)

lm1=lm(velocity~1,data=samara)
lm2=lm(velocity~sqrt(loading),data=samara)
lm3=lm(velocity~sqrt(loading)+symmetry,data=samara)
lm4=lm(velocity~sqrt(loading)+genus,data=samara)
lm5=lm(velocity~sqrt(loading)+species,data=samara)

anova(lm1,lm2,lm3,lm4,lm5)

#GLM example
glove=read.table("glove.txt",header=T)
str(glove)

library(lattice)
xyplot(Per~Experience|Period,glove.b)


glove.b=glove
glove.b$Per=with(glove,Gloves/Observed)
str(glove.b)

glm1=glm(Per~as.factor(Period)*Experience,data=glove.b,weights=Observed,family=binomial)
display(glm1,digits=3)

glm2=glm(Per~as.factor(Period)+Experience,data=glove.b,weights=Observed,family=binomial)
display(glm2,digits=3)

sim1=sim(glm2,1000)

#(a)
apply(sim1$beta,2,mean)
apply(sim1$beta,2,sd)

mean(sim1$sigma)

#prediction
invlogit(predict(glm2,data.frame(Period=1,Experience=5)))
invlogit(predict(glm2,data.frame(Period=2,Experience=5)))
invlogit(predict(glm2,data.frame(Period=1,Experience=20)))
invlogit(predict(glm2,data.frame(Period=2,Experience=20)))

coef=apply(sim1$beta,2,mean)
x11.new=c(1,0,0,5)
invlogit(coef%*%x11.new)
x12.new=c(1,1,0,5)
invlogit(coef%*%x12.new)

x21.new=c(1,0,0,20)
invlogit(coef%*%x21.new)
x22.new=c(1,1,0,20)
invlogit(coef%*%x22.new)

#interval
p11=invlogit(sim1$beta%*%x11.new)
range(p11)
quantile(p11,c(0.025,0.975))


p21=invlogit(sim1$beta%*%x21.new)
range(p21)
quantile(p21,c(0.025,0.975))

#confidence interval (this part is incorrect)
pre=invlogit(sim1$beta%*%as.matrix(c(1,0,0,5))+rnorm(1000,0,sim1$sigma))
print(quantile(pre,c(.025,.975)),digit=3)
print(predict(glm2,data.frame(Period=1,Experience=5)),digit=3)
display(glm2)
h=c(-0.577-2*63.8,-0.577+2*63.8)
print(invlogit(h),digit=3)