pdf("ex0509a.pdf")

# Read in data from ASCII formatted (Fixed With File)

ex0509dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0509.dat", width=c(8,8), 
        col.names=c("dose","exhaustion"))


# Create a dataset (frame) from input data

ex0509 <- data.frame(ex0509dat)

# Attach the dataset for analysis 

attach(ex0509)

numrows <- max(dose)

numcols <- max(exhaustion)

# Obtain crosstabulation of dose and exhaustion

freqde <- table(dose,exhaustion)

freqde

# Obtain Counts Concordant and Discordant Pairs for each cell

concord <- matrix(rep(0,numrows*numcols),nrow=numrows,byrow=T)
discord <- matrix(rep(0,numrows*numcols),nrow=numrows,byrow=T)


for (i1 in 1:numrows) {
 for (i2 in 1:numcols) {
  for (i3 in 1:numrows)   {
   for (i4 in 1:numcols)  {
    if ((i3>i1) & (i4>i2)) concord[i1,i2] <- concord[i1,i2] + freqde[i3,i4]
    if ((i3>i1) & (i4<i2)) discord[i1,i2] <- discord[i1,i2] + freqde[i3,i4]
    if ((i3<i1) & (i4<i2)) concord[i1,i2] <- concord[i1,i2] + freqde[i3,i4]
    if ((i3<i1) & (i4>i2)) discord[i1,i2] <- discord[i1,i2] + freqde[i3,i4]
 } } } }

concord

discord

# Obtain P and Q (Twice C and D)

pconc <- 0.0
qdisc <- 0.0

for (i1 in 1:numrows) {
  for (i2 in 1:numcols) {
   pconc <- pconc + (freqde[i1,i2]*concord[i1,i2])
   qdisc <- qdisc + (freqde[i1,i2]*discord[i1,i2])
   }  }

pconc

qdisc


# Obtain row and column totals and sums of squares for Kendall's tau 

rowtot <- c(rep(0,numrows))

coltot <- c(rep(0,numcols))

for (i1 in 1:numrows) {
 for (i2 in 1:numcols) {
  rowtot[i1] <- rowtot[i1]+freqde[i1,i2]
  } }

for (i1 in 1:numcols) {
 for (i2 in 1:numrows) {
  coltot[i1] <- coltot[i1]+freqde[i2,i1]
  } }

alltot <- sum(freqde)

sumrowtotsq <- 0
sumcoltotsq <- 0

for (i1 in 1:numrows) {
 sumrowtotsq <- sumrowtotsq + (rowtot[i1]^2)
}


for (i1 in 1:numcols) {
 sumcoltotsq <- sumcoltotsq + (coltot[i1]^2)
}

# Denominator for Kendall's Tau

kentaudenr <- alltot^2-sumrowtotsq

kentaudenc <- alltot^2-sumcoltotsq

kentauden <- sqrt(kentaudenr*kentaudenc)

# Compute Gamma and Kendall's Tau

gamma <- (pconc-qdisc)/(pconc+qdisc)

kendall.tau <- (pconc-qdisc)/kentauden


# Compute Variance and standard errors for gamma and Kendall's Tau

gamma.var1 <- 0.0

kentau.var1 <- 0.0

for (i1 in 1:numrows) {
  for (i2 in 1:numcols) {
     gamma.var1 <- gamma.var1 + (freqde[i1,i2]*(((qdisc*concord[i1,i2])-(pconc*discord[i1,i2]))^2))
     kentau.var1 <- kentau.var1 + (freqde[i1,i2]*(((2*kentauden*(concord[i1,i2]-discord[i1,i2]))+
                                  (kendall.tau*((rowtot[i1]*kentaudenc)+(coltot[i2]*kentaudenr))))^2))
 } }

gamma.var <- 16*gamma.var1/((pconc+qdisc)^4)

kentau.var <- (kentau.var1-((alltot^3)*(kendall.tau^2)*((kentaudenr+kentaudenc)^2)))/(kentauden^4)

gamma.se <- sqrt(gamma.var)

kentau.se <- sqrt(kentau.var)

# Print out estimates and 95% Confidence intervals

gamma.lo <- gamma - (1.96*gamma.se)

gamma.hi <- gamma + (1.96*gamma.se)

kentau.lo <- kendall.tau - (1.96*kentau.se)

kentau.hi <- kendall.tau + (1.96*kentau.se)

gamma

gamma.se

gamma.lo

gamma.hi

kendall.tau

kentau.se

kentau.lo

kentau.hi







dev.off()