require(Rdbi)
require(RdbiPgSQL)
require(Hmisc)
require(lattice)
adncolor()

ysr <- dbConnect.PgSQL(user = "guest", password = "nilda",
	dbname = "ysr", port = "5432", host = "localhost")
se <- function(x) { var(x, na.rm = T) -> v ; (v/length(na.omit(x)))^.5 }
meanse <- function(x) {
	c(mean = mean(x, na.rm = T),
	low = mean(x, na.rm = T) - se(x),
	high = mean(x, na.rm = T) + se(x)) 
	}
reorder.factor <- function(fac, y, fun=mean, ...) {
	ordered(fac, levels=levels(fac)[order(tapply(y, fac, fun, ...))]) }
panel.ysr <- function(x,y,...) { 
	panel.superpose(x,y, col = "black", pch = 49:54, ...)
	means <- tapply(y, x, mean)
	lpoints(1:length(means), means, pch = 8, cex = 5) }

## DATA

soild18o <- dbGetQuery(ysr,
	"SELECT block, CAST(CAST(soil_depth AS TEXT)",
		"AS SMALLINT) AS depth, d18o",
	"FROM soild18o",
	"WHERE date = '2003-07-08'",
		"AND BLOCK IN(1,2,3,4,5,6)",
	"ORDER BY block, depth")
soild18o$depth <- as.numeric(soild18o$depth)

plantd18o <- dbGetQuery(ysr, "SELECT block, species, treatment, d18o",
	"FROM plantd18o",
 	"WHERE date = '2003-07-08'",
 		"AND species NOT IN('large', 'potentilla', 'sageunknown')",
 	"AND BLOCK IN(1,2,3,4,5,6)")
plantd18o$species[plantd18o$species == "carex"] <- "carexB"
plantd18o$treatment <- factor(plantd18o$treatment)

## DATA FOR TABLE 1

dbGetQuery(ysr, "SELECT block, avg(tdr)",
"FROM tdr",
"WHERE date = '2003-05-29' ",
	"AND depth IN ('0-15', '15-30')",
	"AND treatment = 's+'",
	"AND position IN ('u')",
"GROUP BY block",
"ORDER BY block")

dbGetQuery(ysr, "SELECT block, avg(tdr)"
	"FROM tdr",
	"WHERE date = '2003-05-29' ",
		"AND depth IN ('30-45', '45-60')",
		"AND treatment = 's+'",
		"AND position IN ('u')",
	"GROUP BY block",
	"ORDER BY block")

dbGetQuery(ysr, "SELECT block, avg(tdr)",
	"FROM tdr",
	"WHERE date = '2003-07-05' ",
		"AND depth IN ('0-15', '15-30')",
		"AND treatment = 's+'",
		"AND position IN ('u')",
	"GROUP BY block",
	"ORDER BY block")
	
dbGetQuery(ysr, "SELECT block, avg(tdr)",
	"FROM tdr",
	"WHERE date = '2003-07-05' ",
		"AND depth IN ('30-45', '45-60')",
		"AND treatment = 's+'",
		"AND position IN ('u')",
	"GROUP BY block",
	"ORDER BY block")
	
## MIXING MODEL

shallowsoil <- tapply(soild18o$d18o, soild18o$block, max)
deepsoil <- tapply(soild18o$d18o[soild18o$depth > 30],
	soild18o$block[soild18o$depth > 30], mean)
plantd18o$minsf <- 1- ((plantd18o$d18o - shallowsoil[plantd18o$block]) /
	(deepsoil[plantd18o$block] - shallowsoil[plantd18o$block]))
plantd18o$minsf[plantd18o$minsf > 1] <- 1
plantd18o$species <- reorder.factor(factor(plantd18o$species), 
	plantd18o$minsf, median)

# SIX PANELS (FIGURE 2)

quartz(w = 16, h = 10.3)
xyplot(c(0-soild18o$depth,  rep(0, length(plantd18o$d18o)))
		~ c(soild18o$d18o, plantd18o$d18o)
		| c(soild18o$block, plantd18o$block),
	groups = c(rep("soil", length(soild18o$depth)),
		as.character(plantd18o$species)),
	data = soild18o, 
	pch = c(rep(124, 6),16,16,16,1),
	cex = c(rep(3, 6),1,1,1,1),
	col = c("red","green", "darkgreen", "blue", "yellow", "purple",
			"#666666", "#999999", "#CCCCCC", "black"),		
	key = list(points = list(
			pch = c(rep(124, 6),16,16,16,1),
			cex = c(rep(3, 6),1,1,1,1),
			col = c("red","green", "darkgreen", "blue", "yellow", "purple",
			"#666666", "#999999", "#CCCCCC", "black")),
		text = list(lab = sort(unique(c(rep("soil", 
				length(soild18o$depth)),
 			as.character(plantd18o$species)))), cex = 1.2),
		columns = 5))

## MINIMUM FSSW (FIGURE 3)

quartz(w = 16, h = 10.3)
xyplot(minsf ~ species,
	data = plantd18o,
#	pch = as.numeric(block)+48, cex = 2,
	pch = 16,
	panel = function(x,y, ...) { 
		panel.xyplot(x,y, ...)
		medians <- tapply(y, x, median)
		lpoints(1:length(medians), medians, pch = 8, cex = 2)
		})

## HODGES-LEHMANN ESTIMATORS FOR MULTIPLE COMPARISON CONFIDENCE INTERVALS

require(exactRankTests)
herbs.hlci <- function() {
	a <- .05/7
	sage <- plantd18o$minsf[plantd18o$species=="sagelarge" |
		plantd18o$species=="sagemedium"]
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="aster"],
		conf.int=TRUE, conf.level=(1-a)) -> asterci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="carexA"],
		conf.int=TRUE, conf.level=(1-a)) -> carexAci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="carexB"],
		conf.int=TRUE, conf.level=(1-a)) -> carexBci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="grass"],
		conf.int=TRUE, conf.level=(1-a)) -> grassci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="ivesia"],
		conf.int=TRUE, conf.level=(1-a)) -> ivesiaci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="sagesmall"],
		conf.int=TRUE, conf.level=(1-a)) -> sagebrushsmallci
	wilcox.exact(sage, plantd18o$minsf[plantd18o$species=="juncus"],
		conf.int=TRUE, conf.level=(1-a)) -> juncusci
	speciesci <- c("Small sagebrush", "Ivesia californica", 
		"Carex sp. A", "Aster occidentalis", "Juncus balticus",
		"Carex sp. B", "Muhlenbergia richardsonii")
	meandifs <- c(sagebrushsmallci$estimate, ivesiaci$estimate, carexAci$estimate,
		juncusci$estimate, asterci$estimate, carexBci$estimate,
		grassci$estimate)
	lowci <- c(sagebrushsmallci$conf.int[1], ivesiaci$conf.int[1],
		carexAci$conf.int[1], juncusci$conf.int[1], asterci$conf.int[1],
		carexBci$conf.int[1], grassci$conf.int[1])
	highci <- c(sagebrushsmallci$conf.int[2], ivesiaci$conf.int[2], 
	carexAci$conf.int[2], juncusci$conf.int[2], asterci$conf.int[2],
	carexBci$conf.int[2], grassci$conf.int[2])
	data.frame(Species = speciesci, "{Mean difference}" = meandifs*-1, 
		"Low 95\%CI" = highci*-1, "High 95\%CI" = lowci*-1) }

require(xtable)
print(xtable(herbs.hlci()), digits = rep(5, 2))

## SAGEBRUSH TREATMENT EFFECT (FIGURE 4)

quartz(w = 16, h = 10.3)
function() {
	herbs <- subset(plantd18o, (plantd18o$treatment == "s+" | 
		plantd18o$treatment == "s-"))
	print(t.test(minsf ~ treatment, data = herbs, var = T))
	xyplot(minsf ~ treatment, data = herbs, 
	panel = function(x,y, ...) {
		means <- tapply(y, x, mean)
		lpoints(1:length(means), means, pch = 8, cex = 1.5, col = "red")
		panel.xyplot(x,y, ...) },
		ylim = c(-0.1,1.1)) }