# Nicholas Horton, nhorton@smith.edu
# Sat Mar 19 09:18:26 EDT 2011
# What proportion of the US is within a mile of a road?
# CAUSE Activity webinar

ds = read.csv("myroadless.csv", stringsAsFactors=FALSE)

process = function(grp) {
  smgrp = subset(grp, incontinent==1)
  tab = table(smgrp$within1mile)
  if (length(tab)==1) {
     tab[2] = tab[1]; tab[1] = 0
  }
  res = binom.test(tab[2], sum(tab))
  return(list(est=res$estimate, ci=res$conf.int))
}

miny = .68
val1 = process(ds[1:50,])
val2 = process(ds[51:100,])
val3 = process(ds[101:150,])
val4 = process(ds[151:200,])
valall = process(ds)

pdf("combined.pdf", width=6, height=4)
plot(c(1,5), c(miny, 1), xlab="group", xaxt="n", ylab="probability within 1 mile", pch=" ")
points(1, val1$est)
points(2, val2$est)
points(3, val3$est)
points(4, val4$est)
points(5, valall$est)
lines(c(1,1), val1$ci, col="blue")
lines(c(2,2), val2$ci, col="blue")
lines(c(3,3), val3$ci, col="blue")
lines(c(4,4), val4$ci, col="blue")
lines(c(5,5), valall$ci, col="red")

text(1, miny, "1")
text(2, miny, "2")
text(3, miny, "3")
text(4, miny, "4")
text(5, miny, "all")
dev.off()