#  Cochran-Armitage trend test for 2 x I tables where I have scores
# Edps 689
# Fall 2018
# C.J.Anderson
#
# Framingham Heart Study  (Agresti, 1990, p 93);
#
library(DescTools)            # needed for Cochran-Armitage trend test
library(vcd)

setwd("D:\\Dropbox\\edps 589\\2 Chi-square")

hs <- read.table("framingham_heart_data.txt",header=TRUE)
hs

hs.tab <- xtabs(count ~ bp + heart,data=hs)

CochranArmitageTest(hs.tab, alternative = c("two.sided", "increasing", "decreasing"))


hs.mat <-as.matrix(hs.tab)
n.row <- hs.mat[,1] + hs.mat[,2]
new <- cbind(hs.mat,n.row)
new.data <- as.data.frame(new)
new.data$p.yes<- new.data$yes/new.data$n.row
new.data$p.no <- new.data$no/new.data$n.row
new.data$scores <- c(1,10,2,3,4,5,5.5,8)

# Figure of data, and 2 smoothing curves
scatter.smooth(new.data$scores,new.data$p.yes,
     main="Framingham Heart Study",
     xlab="Blood Pressure Scores",
	 ylab="Proportion Who Had Heart Attack",
	 pch=20,
	 cex=1.5)
fitted.lm <- fitted(lm(new.data$p.yes ~ new.data$scores))
lines(new.data$scores,fitted.lm,col="blue")
text(1.45,.175,"smooth",col="black")
text(2,.185,"linear regression",col="blue")

	
# Proportation heart attack:
prop.trend.test(new.data$yes, n.row, score=c(1,10,2,3,4,5,5.5,8))
# which is the same as not heart attack:
prop.trend.test(new.data$no, n.row, score=c(1,10,2,3,4,5,5.5,8))