pdf("ex0512.pdf")

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

ex0512dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0512.dat", width=c(8,8), 
        col.names=c("siskel","ebert"))


# Create a dataset (frame) from input data

ex0512 <- data.frame(ex0512dat)

# Attach the dataset for analysis 

attach(ex0512)

# Obtain crosstabulation of siskel and ebert

freqse <- table(siskel, ebert)

freqse

numlevs <- nrow(freqse)

# Obtain the proportions within the table

propse <- freqse/sum(freqse)

propse

# Obtain Siskel (Row) and Ebert (Column) Marginal Counts and Proportions

freqsiskel <- margin.table(freqse,1)

freqebert <- margin.table(freqse,2)

propsiskel <- margin.table(freqse,1)/sum(freqse)
propebert <- margin.table(freqse,2)/sum(freqse)

propsiskel
propebert

# Compute P_o and P_e

p_o <- 0.0

p_e <- 0.0

for (i1 in 1:numlevs) {
     p_o <- p_o + propse[i1,i1]
     p_e <- p_e + (propsiskel[i1]*propebert[i1])
  }

# Compute Kappa

kappa <- (p_o-p_e)/(1-p_e)

# Compute variance and standard error of kappa

se_a <- 0.0
se_b <- 0.0
se_c <- 0.0

for (i1 in 1:numlevs) {
   se_a <- se_a + (propse[i1,i1]*((1-((propsiskel[i1]+propebert[i1])*(1-kappa)))^2))
    for (i2 in 1:numlevs) {
     if (i2 != i1) se_b <- se_b + propse[i1,i2]*((propebert[i1]+propsiskel[i2])^2)
 } }

se_b <- se_b*((1-kappa)^2)

se_c <- (kappa-(p_e*(1-kappa)))^2

kappa.var <- (se_a+se_b-se_c)/(((1-p_e)^2)*sum(freqse))

kappa.se <- sqrt(kappa.var)

#  Obtain 95% Confidence Interval for Kappa


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

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

kappa

kappa.se

kappa.lo

kappa.hi



dev.off()