## LEAF CARBON ISOTOPE RESPONSES OF RESIDENT HERBS TO A WOODY SHRUB INVASION

# Note: This script can take several minutes to run due to the computationally
#	intensive calculation of simultaneous confidence intervals.

## ACQUIRE DATA FROM DATABASE

require(RODBC)
ysr <- odbcConnect("ysr")

c13 <- sqlQuery(ysr, paste(
"SELECT date, block, species, AVG(d13c) AS d13c,",
		"treatment", # using block-averaged d13c avoids pseudoreplication
	"FROM d13c",
	"WHERE date NOT IN ('2003-08-12', '2002-06-05', '2002-06-02')",
		"AND species NOT IN ('muhlenbergia')", #, 'forb', 'potentilla')",
		"AND d13c < - 20", # eliminates one pyrrocoma outlier = -13
	"GROUP BY date, block, species, treatment",
	"UNION",
	"SELECT date, block, species, d13c,",
		"CONCAT(treatment,position) AS treatment",
	"FROM d13c",
	"WHERE date = '2003-08-12'",
		"AND species NOT IN ('muhlenbergia')", #, 'forb', 'potentilla')",
		#"'ublf', 'ivesia', 'aster')",
	"ORDER BY date, block, species, treatment"))
c13$date <- factor(c13$date)
c13$treatment <- ordered(c13$treatment, 
	levels = c("s+", "s-", "s+u", "s-n", "s+a"))
c13$block <- factor(c13$block)
c13$species <- factor(c13$species)

nadia <- sqlQuery(ysr, paste("SELECT block, habitat, species,",
	"position, AVG(d13c) AS d13c, AVG(percentc) AS percentc,",
	"AVG(percentn) AS percentn, AVG(c_to_n) AS c_to_n, associated",
	"FROM nadia",
	"WHERE habitat NOT IN ('sage')",
		"AND d13c < - 24", # eliminates outlier Armo XP 4 under = -23.65823)
	"GROUP BY block, habitat, species, position",
	"ORDER BY block, habitat, species, position"))
nadia$block <- factor(nadia$block)
nadia$species <- factor(nadia$species)
nadia$habitat <- factor(nadia$habitat)
nadia$position <- ordered(nadia$position, levels = c("HM", "s-", "away", 
	"near", "under"))

## FIGURES SHOWING THE DATA

require(lattice)
source("http://anthony.darrouzet-nardi.net/works/adnlattice.R")

panel.ysr <- function(x,y,...) { 
	panel.superpose(x,y, col = "black", pch = 1, ...)
	means <- tapply(y, x, mean)
	lpoints(1:length(means), means, pch = 8, col = "red") }
get(getOption("device"))()
print(xyplot(d13c ~ treatment | date, 
	data = c13, groups = block, panel = panel.ysr))

panel.asc <- function(x,y,...) { if (length(x) > 1) {
	panel.superpose(x,y, col = "black", pch = 1, ...)
	means <- tapply(y, x, mean)
	lpoints(1:length(means), means, pch = 8, col = "red") } }
get(getOption("device"))()
print(xyplot(d13c ~ position | species, groups = associated,
	data = nadia, panel = panel.asc,
	subset = (d13c < -24) & position != "HM" & position != "s-")) 
	# subset eliminates an outlier Armo, XP, 4, under = -23.65823)

## STATISTICAL ANALYSIS AND ACCOMPANYING FIGURES

require(multcomp)

cmat <- function(fac, i) { 
	m <- matrix(0, nrow = length(i)/2, ncol = max(i))
	n <- rep(1:(length(i)/2), each = 2)
	rnames <- c()
	for (h in 1:(length(i)/2)*2) {
		m[n[h],i[h]] <- 1
		m[n[h],i[h-1]] <- -1
		rnames[h/2] <- paste(levels(fac)[i[h]],"-",levels(fac)[i[h-1]], 
			sep = "") }
	colnames(m) <- levels(fac)
	rownames(m) <- rnames
	m }

c13$int <- factor(interaction(c13$treatment, c13$date))
cmatrix <- cmat(c13$int, c(1:8, 9,10,9,11,10,11))
c13stats <- simint(d13c ~ int + block:date:species, data = c13,
	whichf = "int", cmatrix = cmatrix)
print(c13stats)
get(getOption("device"))()
print(plot(c13stats))

nadia$associated[is.na(nadia$associated)] <- -99
nadia$int <- factor(interaction(nadia$position,	nadia$species))
cmatrix <- cmat(nadia$int, 
	c(2,4,6,8,10,12,14,16,19,21,23,25, 3,4,7,8,11,12,15,16,20,21,24,25,  
		2,3,6,7,10,11,14,15,19,20,23,24))
nadiastats <- simint(d13c ~ int + associated:species, data = nadia, 
	whichf = "int", cmatrix = cmatrix, subset = (d13c < -24)) 
	# subset eliminates an outlier Armo, XP, 4, under = -23.65823)
print(nadiastats)
get(getOption("device"))()
print(plot(nadiastats))