### Data from Zhao and  Boerwinkle (http://genome.cshlp.org/content/12/11/1679.full)
AG = 844427
CT = 845441
AC = 231506
GT = 233387
AT = 192285
CG = 229857

snp = c(AG,CT,AC,GT,AT,CG)
names(snp) = c("A/G","C/T","A/C","G/T","A/T","C/G")
barplot(snp,col=c("blue","blue",rep("red",4)),ylab="Count")
barplot(snp/sum(snp),col=c("blue","blue",rep("red",4)),ylab="Freq")
pie(snp)



### Iris data
### library(MASS)
data(iris)
par(mfrow=c(2,2))
hist(iris[,1],main="Histogram of Sepal.Length",xlab="Sep.Length")
hist(iris[,2],main="Histogram of Sepal.Width",xlab="Sepal.Width")
hist(iris[,3],main="Histogram of Petal.Length",xlab="Petal.Length")
hist(iris[,4],main="Histogram of Petal.Width",xlab="Petal.Width")

plot(iris[,c(3,4)])
plot(iris[,-5])


plot(iris[,3],iris[,4],pch=19,xlab="Petal.Length",ylab="Petal.Width",col=c("red","green3","blue")[unclass(iris$Species)])
legend("bottomright",legend=c("setosa","versicolor","virginica"),col=c("red","green3","blue"),pch=15)



pairs(iris[1:4],main="Iris Data (red=setosa,green=versicolor,blue=virginica)", pch=21,
  bg=c("red","green3","blue")[unclass(iris$Species)])

par(mfrow=c(1,2))
boxplot(Petal.Length~ Species,data=iris,col=c("red","green3","blue"),main="Boxplot of Petal.Length")
boxplot(Petal.Width~ Species,data=iris,col=c("red","green3","blue"),main="Boxplot of Petal.Width")

library("beeswarm")

beeswarm(Petal.Length~ Species,data=iris,col=c("red","green3","blue"),pch=20,method="square")

iris.tmp = iris
for(i in 1:4) iris.tmp[,i] = iris.tmp[,i] + rnorm(nrow(iris),mean=0,sd=0.01)
boxplot(Petal.Length~ Species,data=iris.tmp,main="Boxplot of Petal.Length")
beeswarm(Petal.Length~ Species,data=iris.tmp,col=c("red","green3","blue"),add=TRUE,pch=20,method="square")

boxplot(Petal.Width~ Species,data=iris.tmp,main="Boxplot of Petal.Width")
beeswarm(Petal.Width~ Species,data=iris.tmp,col=c("red","green3","blue"),add=TRUE,pch=20,method="square")


#m1 = mean(iris[,2])
#sd = sd(iris[,2])
#x1 = seq(0,6,by=0.01)
#y = dnorm(x1,mean=m1,sd=sd)
#hist(iris[,2],main="Histogram of Sepal.Width",xlab="Sepal.Width",freq=F)
#lines(x1,y,col="red",lwd=2)


###Aspirin data
x = matrix(c(104,10933,189,10845),nrow=2)
colnames(x) = c("aspirin","placebo")
rownames(x) = c("HeartAttack","NoHeartAttach")
fisher.test(x,alternative="less")


#### normal fitting
library(MASS)
data(Pima.tr)
m1 = mean(Pima.tr$bmi)
sd = sd(Pima.tr$bmi)
x1 = seq(0,60,by=0.1)
y = dnorm(x1,mean=m1,sd=sd)
hist(Pima.tr$bmi,main="Histogram of BMI",xlab="Sepal.Width",freq=F)
lines(x1,y,col="red",lwd=2)
legend("topleft",paste(c("mu=","sd="),format(c(m1,sd),digits=2)),bty="n")

### quantile-quantile plot (QQ-plot)
### BMI
m1 = mean(Pima.tr$bmi)
sd = sd(Pima.tr$bmi)
prb = (1:nrow(Pima.tr))/nrow(Pima.tr)
q.n = qnorm(prb,mean=m1,sd=sd)
plot(sort(Pima.tr$bmi),q.n,xlab="Observed Quantile",ylab="Theoretical Quantile",main="QQ-plot of BMI in Pima.tr")
abline(a=0,b=1,col="red",lwd=2)
###Petal.Length
m1 = mean(iris$Petal.Length)
sd = sd(iris$Petal.Length)
prb = (1:nrow(iris))/nrow(iris)
q.n = qnorm(prb,mean=m1,sd=sd)
plot(sort(iris$Petal.Length),q.n,xlab="Observed Quantile",ylab="Theoretical Quantile",main="QQ-plot of BMI in Pima.tr")
abline(a=0,b=1,col="red",lwd=2)