library(Hmisc)

dat <- read.csv('C:/RGraphicsTutorial/Datasets/CleanTrial.csv')

base <- subset(dat, visit%in%2)
eot <- subset(dat, visit%in%6)

# binary endpoint is success if at least 2 grade improvement
eot$binep <- ifelse(base$catvar-eot$catvar >= 2, 1, 0)
# Also create a site id (25 total sites)
set.seed(42)
eot$site <- sample(1:25, 1500, replace=TRUE)

    ss <- with(eot, summarize(binep, llist(site, trtf), mean))
    n <- with(eot, summarize(binep, llist(site, trtf), length))

    sdat <- data.frame(ss, n[,3])
        names(sdat) <- c("Site", "Trt", "Mean", "N")
    so <- with(eot, summarize(binep, llist(trtf), mean))
    
    nsn <- length(unique(sdat$Site))
    ut <- unique(sdat$Trt)
    rnx <- tabulate(as.factor(sdat$Site))
    sdat$plotx <- rep(1:nsn, rnx[rnx>0])
    # Creation of the figure.
    
    #png('C:/Research/Graphics/Graphs4Display/webpages/classes/pages/images/efficacybysite.png',
    #         width=670, height=500)

    plot(c(.5, nsn+.5), c(min(sdat$Mean)-.05, max(sdat$Mean)+.05), 
            type="n", ylab='Response Rate', xlab='', axes=FALSE)
    axis(1, at=1:nsn, labels=as.character(unique(sdat$Site)), cex.axis=1, las=3)
    axis(2)
    box()
    sdat$plotx <- sdat$plotx + rep(c(-.05,.05), length(sdat[,1])/2)
    coluers <- c('seagreen2','grey60')
    pchs <- c(1,16)
    for(k in 1:length(ut)){
        subdat <- subset(sdat, sdat$Trt==ut[k])
        points(subdat$plotx, subdat$Mean, 
                  pch=pchs[k],
                col=coluers[k])
        # Add overall means as horizontal lines
        abline(h=so[k,2], col=coluers[k], lty=3, lwd=2)
        for(j in 1:length(subdat$N)){
                text(subdat$plotx[j]+.4, subdat$Mean[j], labels=subdat$N[j],
                     col=coluers[k], cex=1)
        }   
    }
    for(i in 1:nsn){
        subdat <- subset(sdat, sdat$Site==unique(sdat$Site)[i])
        lines(c(i,i), c(min(subdat$Mean), max(subdat$Mean)), lty=2,col='gray60')
    }
    legend("topleft", legend=as.character(sort(unique(eot$trtf))),
            pch=pchs,
            col=coluers,
            cex=1.25, xjust=0)
    #dev.off()