# R script for lecture 6 in STK4900/9900

# Slide 12: Kittiwake example --------------------------------------------------

kittiwake <- matrix(c(100,175,60,434), nrow=2)

rownames(kittiwake) <- c("unsuccessful","successful")
colnames(kittiwake) <- c("divorced","not_divorced")
kittiwake

chisq.res <- chisq.test(kittiwake, correct=FALSE)
chisq.res
sum(((chisq.res$observed-chisq.res$expected)^2)/chisq.res$expected)

prop.res <- prop.test(kittiwake, correct=FALSE)
prop.res

# Slide 15: Contingency table --------------------------------------------------

bloodpr <- matrix(c(14,11,6,11,11,10,8,9,12), nrow=3)
rownames(bloodpr) <- c("F.low","F.middle","F.upper")
colnames(bloodpr) <- c("C.low","C.middle","C.upper")
bloodpr

chisq.res <- chisq.test(bloodpr, correct=FALSE)
chisq.res
chisq.res$observed
chisq.res$expected


# Slide 27-28: Logistic regression - CHD data ----------------------------------

wcgs <- read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/data/wcgs.txt", sep="\t", header=TRUE, na.strings=".")

# Group data by mean age in each age category
wcgs$agem <- 39.5*(wcgs$agec==0)+42.9*(wcgs$agec==1)+47.9*(wcgs$agec==2)+52.8*(wcgs$agec==3)+57.3*(wcgs$agec==4)
table(wcgs$agem,wcgs$chd69)

fit.binary <- glm(chd69~agem, data=wcgs, family=binomial)
summary(fit.binary)
predict(fit.binary, type="response", data.frame(agem=50))

tmp <- table(wcgs$agem,wcgs$chd69)
chd.grouped <- data.frame(chd=tmp[,2], no.chd=tmp[,1], agem=as.numeric(rownames(tmp)))
row.names(chd.grouped) <- 1:5

fit.grouped <- glm(cbind(chd,no.chd)~agem, data=chd.grouped, family=binomial)
summary(fit.grouped)