CoDA week4 PA.R

# wk4 PA.R  R script for Prog Assign - week 4

# load the data set and create contigency table c
# for grouped petal width (col 4) and species col 5

fname <- 'wk4 dataset_390_1.txt'
colnames <- c('SL','SW','PL','PW','Species')
ds <- read.csv( fname, header = FALSE, sep = "\t", col.names=colnames)
breaks <- c(0.0, 0.5, 1.5, 2.0, 99.0)   # right border included
ds$PWc <- cut(ds$PW, breaks=breaks, right = TRUE)
c <- table(ds$Species, ds$PWc)

cat("\nContingency Frequency Table\n")
print(addmargins(c))

cat("\nContingency Probability Table\n")
# = relative frequencies totaling to 1
p <- prop.table(c)
print(addmargins(p))

cat("\nEpected Values Table (as probabilities)\n")
Xsq <- chisq.test(p)
ep <- Xsq$expected
print(Xsq$expected)

cat("\nQuetelet Index Table\n")
# quetelet: q[r,c] = (p[r,c] / (p[r]*p[c])) - 1
q <-  p / Xsq$expected - 1
print(q)


cat("\nPearson Index Table = Pearson Residuals Xsq$residuals\n")
# pearson: x[r,c] = (p[r,c] - (p[r]*p[c])) / sqrt((p[r]*p[c]))
#x <-  (p - Xsq$expected) / sqrt(Xsq$expected)
print(Xsq$residuals)


cat("\noutput for grading: Quetelet || Pearson\n\n")
cat(sprintf("%.3f %.3f %.3f %.3f\n",q[1,],q[2,],q[3,],q[4]),sep='')
cat(sprintf("%.3f %.3f %.3f %.3f\n",Xsq$residuals[1,],Xsq$residuals[2,],Xsq$residuals[3,],Xsq$residuals[4]),sep='')

updated for proper printing at end