#LESSON 3 (DESCRIPTIVE STATISTICS & CONTINGENCY TABLES) (12/04/2021)
library(ISwR)
attach(thuesen)
attach(red.cell.folate)
red.cell.folate

tapply(red.cell.folate)
tapply(folate, ventilation, mean)
tapply(folate, ventilation, sd)
tapply(folate, ventilation, length)

#Try to something nicer
xbar<-tapply(folate, ventilation, mean)
s<-tapply(folate, ventilation, sd)
n<-tapply(folate, ventilation, length)
cbind(mean=xbar, std.dev=s, n=n)

data(juul)
attach(juul)
tapply(igf1, tanner, mean)
tapply(igf1, tanner, mean, na.rm=T)

caff.marital<-matrix(c(652, 1537, 598, 242, 36, 46, 38, 21, 218, 327, 106, 67), nrow=3, byrow=T)
caff.marital

colnames(caff.marital)<-c("0", "1-150", "151-300", ">300")
rownames(caff.marital)<-c("Married", "Prev.married", "Single")
caff.marital

#Tabulation
table(sex)
table(sex,menarche)
table(menarche, tanner)
table(caff.marital)
t(caff.marital)

table(tanner, sex)
tanner.sex<-table(tanner, sex)
tanner.sex

margin.table<-table(tanner.sex,1)

juul<-transform(juul, sex=factor(sex, labels=c("M", "F")), menarche=factor(menarche, labels=c("No", "Yes")), tanner=factor(tanner, labels=c("I", "II", "III", "IV", "V")))
margin.table<-table(tanner.sex,1)
detach(juul)
juul$sex<-factor(juul$sex, labels=c("M", "F"))
juul$menarche<-factor(juul$menarche, labels=c("No", "Yes"))
juul$tanner<-factor(juul$tanner, labels=c("I", "II", "III", "IV", "V"))
attach(juul)

#COUNT
margin.table(tanner.sex, 1)
margin.table(tanner.sex, 2)

#PROPORTION
prop.table(tanner.sex, 1)
prop.table(tanner.sex, 2)

total.caff<-margin.table(caff.marital, 2)
total.caff

barplot(total.caff, col="white")


par(mfrow=c(2,2))
barplot(caff.marital, col="white")
barplot(t(caff.marital), col="white")
barplot(t(caff.marital), col="white", beside=T)
barplot(prop.table(t(caff.marital), 2), col="white", beside=T)
 

caff.marital
prop.table(t(caff.marital))

par(mfrow=c(1,1))
barplot(prop.table(t(caff.marital), 2), beside=T, legend.text=colnames(caff.marital), col=c("white", "grey80", "grey50", "black"))


data(energy)
attach(energy)
expend.lean<-expend[stature=="lean"]
expend.obese<-expend[stature=="obese"]

par(mfrow=c(2,1))
hist(expend.lean, breaks=10, xlim=c(5,13), ylim=c(0,4), col="white")
hist(expend.obese, breaks=10, xlim=c(5,13), ylim=c(0,4), col="grey")

boxplot(expend.lean, expend.obese)

#These commands show only (P VALUE)
chisq.test(caff.marital) 
chisq.test(tanner.sex)
 
#CALCULATION
E <- chisq.test(caff.marital)$expected
O <- chisq.test(caff.marital)$observed
(O-E)^2/E


data(juul)
attach(juul)

chisq.test(tanner,sex)
fisher.test(caff.marital)