## ehr - EcoCELL Heterotrophic Respiration

# 1. Short term temperature response in roots
	# shorttermtemp.plot # NOT PRINTING
	# shorttermtemp.bars
	# tvtol.ttest
	# tv4ol.ttest
# 2. Extent of acclimation
	# acclimation.plot # NOT PRINTING
	# acclimation.bars
	# wc.ttest
	# soil.ttest
	# rootol.ttest
# 3. Partitioning of root vs. microbial respiration
	# pr.strip # NOT PRINTING
	# pr.bars
	# prol.ttest # NOT PRINTING
	# rrwcol.ttest
	# rhwcol.ttest	

# This script outputs 3 barplot with standard errors figures in the
# working directory, one for each of the three above analyses. In
# addition, there are 7 t-tests whose results correspond to the  bar
# plots. There are also dotplots of the raw data which are not set to
# print, and a ttest for the difference between the percent respiration by
# roots in warmed v. control treatments (prol.ttest) that is not set to
# print. Run the script without cleanup and run the unprinted objects to
# see them.

## R 1.9.0

require(Rdbi)
require(RdbiPgSQL)
require(Hmisc)
require(lattice)
require(gregmisc)

ehr <- dbConnect.PgSQL(user = "guest", password = "nilda",
	dbname = "ehr", 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)) 
	}


## 1. SHORT TERM TEMPERATURE RESPONSE IN ROOTS

# Difference between native temperature and the other treatment 
tvt <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key,",
		"rr4.rr AS rr4, rr8.rr AS rr8,",
		"rr8.rr - rr4.rr AS tvtdiff,", # paired differences
		"treatment",
	"FROM rr rr8, rr rr4, ehr", # self join
	"WHERE rr4.pid = rr8.pid",
		"AND rr4.pid = ehr.pid", "AND rr8.pid = ehr.pid",
		"AND rr4.temperature = 4 AND rr8.temperature = 7.8",
	"ORDER BY treatment, monolith, position, rr4.rr"
	)

tvt$treatment <- ordered(tvt$treatment,
	levels = c("control", "warmed"))

# Difference between native temperature and a 4û increase
tv4 <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key,",
		"rrnative.rr AS rrnative, rrplus4.rr AS rrplus4,",
		"rrplus4.rr - rrnative.rr AS tv4diff,",
		"treatment",
	"FROM rr rrplus4, rr rrnative, ehr",
	"WHERE rrnative.pid = rrplus4.pid AND rrnative.pid = ehr.pid",
		"AND rrplus4.pid = ehr.pid",
		"AND ((rrnative.temperature = 4 AND ehr.treatment = 'control')",
			"OR (rrnative.temperature = 7.8 AND ehr.treatment = 'warmed'))",
		"AND ((rrplus4.temperature = 7.8 AND ehr.treatment = 'control')",
			"OR (rrplus4.temperature = 12.8 AND ehr.treatment = 'warmed'))",
	"ORDER BY treatment,  monolith, position"
	)

tv4$treatment <- ordered(tv4$treatment,
	levels = c("control", "warmed"))

# Dotplots for short term temperature response
shorttermtemp.plot <-
function()
{
	tvtp <- dotplot(treatment ~ tvtdiff,
		data = tvt,
		panel = function(x,y) {
			panel.dotplot(x,y) 
			panel.abline(v = 0, lty = 3,) 
			},
		aspect = .5,
		xlab = "Native temp. vs. Other treatment (nmol/g/s)"
		)
	tv4p <- dotplot(treatment ~ tv4diff,
		data = tv4,
		panel = function(x,y) {
			panel.dotplot(x,y)
			panel.abline(v = 0, lty = 3)
			},
		aspect = .5,
		xlab = "Native temp. vs. a 4û increase (nmol/g/s)"
		)
	print(tvtp, split = c(1,1,1,2), more = T)
	print(tv4p, split = c(1,2,1,2))
}

# Barchart with standard errors for shortterm temp response
shorttermtemp.bars <-
function() {
	attach(tvt)
	attach(tv4)
	a <-
		summarize(tvtdiff, treatment, meanse,
			subset = tvtdiff < 3 & tvt$tvtdiff > -0.7)
	b <-
		summarize(tv4diff, treatment, meanse,
			subset = tv4$tv4diff > -1.5)
	barplot2(c(a$tvtdiff, b[2,]$tv4diff),
		col = c("#666666", "#CCCCCC", "#CCCCCC"),
		plot.ci = T,
		ci.l = c(a$low, -b[2,]$low),
		ci.u = c(a$high, b[2,]$high),
		ylim = c(0,1),
		names.arg = c("Control 4-8", "Warmed 8-4", "Warmed 8-12"))
	detach(tvt)
	detach(tv4)
	}
	
# T-tests and means±SE for short term temperature response
	# treatment v. treatment (remove 2 outliers)
	tvtol.ttest <-
		t.test(tvt$tvtdiff ~ tvt$treatment, 
			subset = (tvt$tvtdiff < 3 & tvt$tvtdiff > -0.7), var.equal = T)
	# treatment v. 4û increase (remove 1 outlier)
	tv4ol.ttest <-
		t.test(tv4$tv4diff ~ tv4$treatment,
			subset = tv4$tv4diff > -1.5, var.equal = T)
			

## 2. EXTENT OF ACCLIMATION 

# This analysis examines the extent to which mass specific
# rates of soil (whole cores), root (in water bath at native temp),
# and microbial (root free soil subsamples) respiration acclimate
# to exposure to an anomalously warm climate year.

# whole core
wc <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key, rwc, treatment",
	"FROM ehr",
	"ORDER BY treatment, monolith, position")

wc$treatment <- ordered(wc$treatment,
	levels = c("control", "warmed"))

# soil
soil <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key, rsoil, treatment",
	"FROM ehr",
	"ORDER BY treatment, monolith, position")

soil$treatment <- ordered(soil$treatment,
	levels = c("control", "warmed"))

# root
root <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key, rr, treatment",
	"FROM ehr NATURAL JOIN rr",
	"WHERE (treatment = 'control' AND temperature = 4)",
		"OR (treatment = 'warmed' AND temperature = 7.8)",
	"ORDER BY treatment, monolith, position, temperature")

root$treatment <- ordered(root$treatment,
	levels = c("control", "warmed"))

# Dotplots of acclimation data: whole cores, soil microbial, and roots
acclimation.plot <-
function() {
	wcp <- dotplot(treatment ~ rwc,
		data = wc,
		xlab = "Whole Core Respiration nmol/kg/s",
		xlim = c(0,80),
		aspect = .4)
	soilp <- dotplot(treatment ~ rsoil,
		data = soil,
		xlab = "Soil Microbial Respiration nmol/kg/s",
		xlim = c(0,80),
		aspect = .4)
	rootp <- dotplot(treatment ~ rr,
		data = root,
		xlab = "Root Respiration nmol/g/s/",
		aspect = .4)	
	print(wcp, split = c(1,1,1,3), more = T)
	print(soilp, split = c(1,2,1,3), more = T)
	print(rootp, split = c(1,3,1,3))
	}
	
	xyplot(rr ~ treatment,
		data = root,
		xlab = "Root Respiration nmol/g/s/",
		aspect = .4,
		subset = root$rr < 6,
		panel = panel.se)

# Bootstrap Confidence Interval for difference in root means
w <- 1 ; c <- 1
for (n in 1:5000) w[n] <- mean(root$rr[root$treatment == "warmed" &
	root$rr < 6][round(runif (11, .5, 11.5))])
for (n in 1:5000) c[n] <-mean(root$rr[root$treatment == "control" &
	is.na(root$rr) == F][round(runif(11, .5, 11.5))])
quantile(w-c, c(.05, .5, .95))

# Barcharts with standard errors of acclimation data
acclimation.bars <-
function() {
	attach(wc)
	attach(soil)
	attach(root)
	wcs <-
		summarize(rwc, treatment, meanse)
	soils <-
		summarize(rsoil, treatment, meanse)
	roots <-
		summarize(rr, treatment, meanse,
			subset = rr < 6)
	par(mfrow = c(3,1))
	barplot2(as.vector(wcs$rwc),
		col = c("#666666", "#CCCCCC"),
		plot.ci = T,
		ci.l = as.vector(wcs$low),
		ci.u = as.vector(wcs$high),
		ylim = c(0,50),
		names.arg = c("Control", "Warmed"),
		main = "Whole Core")
	barplot2(as.vector(soils$rsoil),
		col = c("#666666", "#CCCCCC"),
		plot.ci = T,
		ci.l = as.vector(soils$low),
		ci.u = as.vector(soils$high),
		ylim = c(0,50),
		names.arg = c("Control", "Warmed"),
		main = "Soil")
	barplot2(as.vector(roots$rr),
		col = c("#666666", "#CCCCCC"),
		plot.ci = T,
		ci.l = as.vector(roots$low),
		ci.u = as.vector(roots$high),
		ylim = c(0,2),
		names.arg = c("Control", "Warmed"),
		main = "Roots")
		par(mfrow = c(1,1))
	detach(wc)
	detach(soil)
	detach(root)
	}

# T-tests of acclimation data
	# whole core
	wc.ttest <-
		t.test(wc$rwc ~ wc$treatment, var.equal = T)
	# soil
	soil.ttest <-
		t.test(soil$rsoil ~ soil$treatment, var.equal = T)
	# root without the outlier
	rootol.ttest <-
		t.test(root$rr ~ root$treatment,
			subset = root$rr < 6, var.equal = T)
	# Minimum detectable difference for roots without the outlier
	rootol.mdd <- qt(.975, 20)*2*se(root$rr[root$rr < 6])
			

## 3. PARTITIONING OF ROOT VS. MICROBIAL RESPIRATION

# This is an analysis of whether exposure to a warm year altered
# the contribution of root (Rr) versus microbial (Rh) respiration
# using a strip plot and a calculation of mean difference ± 95%CI
# between treatments. Two differences are calculated, one with an
# outlier in the control treatment, and one without.

pr <- dbGetQuery(ehr,
	"SELECT monolith, position,",
		"rr * wcrootmass  * .796 AS rrwc,", # convert to area specific rate,
		"rsoil * wcmass / 1000 * .796 AS rhwc,",
		"(rr * wcrootmass) /",
		"((rr * wcrootmass) + (rsoil * wcmass / 1000 )) AS percentr,",
		"treatment ",
		"FROM ehr, rr",
		"WHERE ehr.pid = rr.pid",
			"AND ((treatment = 'control' AND temperature = 4)",
			"OR (treatment = 'warmed' AND temperature = 7.8))",
		"ORDER BY treatment")

pr$treatment <- ordered(pr$treatment,
	levels = c("control", "warmed"))

# Strip plot of partitioning data
pr.strip <-
	stripplot(percentr ~ treatment, data = pr, horizontal = F)

# Bar plot of partitioning data
pr.bars <-
function() {
	attach(pr)
	rrs <-
		summarize(rrwc, treatment, meanse,
			subset = percentr < 0.4)
	rhs <-
		summarize(rhwc, treatment, meanse,
			subset = percentr < 0.4)
	barplot2(
		matrix(data = c(as.vector(rrs$rrwc),
			as.vector(rhs$rhwc)), nrow = 2, ncol = 2),
		col = c("#666666", "#CCCCCC"),
		plot.ci = T,
		ci.l = matrix(data = c(as.vector(rrs$low),
			as.vector(rhs$low)), nrow = 2, ncol = 2),
		ci.u = matrix(data = c(as.vector(rrs$high),
			as.vector(rhs$high)), nrow = 2, ncol = 2),
		names.arg = c("Rr", "Rh"),
		legend.text = c("Control", "Warmed"),
		beside = T)
	detach(pr)
	}

# T-test of partitioning data without the outlier in control treatment
prol.ttest <-
	t.test(pr$percentr ~ pr$treatment, 
		subset = pr$percentr < 0.4, var.equal = T)
rrwcol.ttest <-
	t.test(pr$rrwc ~ pr$treatment, var.equal = T)
rhwcol.ttest <-
	t.test(pr$rhwc ~ pr$treatment, var.equal = T)
		
	
## PRINT RESULTS

#Graphics
postscript(file= "shorttermtemp.eps", height = 6, width = 6,
	horizontal = FALSE, onefile = FALSE, paper = "special")
shorttermtemp.bars()
dev.off()

postscript(file= "acclimation.eps", height = 7, width = 2,
	horizontal = FALSE, onefile = FALSE, paper = "special")
acclimation.bars()
dev.off()

postscript(file= "partitioning.eps", height = 7, width = 7,
	horizontal = FALSE, onefile = FALSE, paper = "special")
pr.bars()
dev.off()



print(" Probabalistic analyses on EHR data")
print(" 1. Short term temperature response in roots")
print(" Compare the differences (temp of native treatment v. other treatment)")
print(" between control and warmed roots.")
print(" (2 outliers removed)")
	print(tvtol.ttest)
print(" Compare the distribution of differences (temp of native treatment v.")
print(" temp of the other treatment) between control and warmed roots.")
print(" (1 outlier removed)")
	print(tv4ol.ttest)

print(" 2. Extent of acclimation")
print(" Compare mass specific respiration rates (nmol/g/s)")
print(" between warmed and control treatments for...")
print(" WHOLE CORES")
	print(wc.ttest)
print(" SOIL MICROBES")
	print(soil.ttest)
print(" ROOTS (one outlier removed")
	print(rootol.ttest)

print(" 3. Partitioning of root vs. microbial respiration")
print(" Compare ground area specific respiration rates (µmol/m2/s) between")
print(" warmed and control treatments for")
print(" SOIL MICROBES")
	print(rrwcol.ttest)
print(" ROOTS")
	print(rhwcol.ttest)


# Cleanup
rm(ehr, se, meanse,
	tvt, tv4,
		shorttermtemp.plot, shorttermtemp.bars, tvtol.ttest, tv4ol.ttest,
	wc, soil, root, 
		acclimation.plot, acclimation.bars,
		wc.ttest, soil.ttest, rootol.ttest,
	pr,
		pr.strip, pr.bars,
		prol.ttest, rrwcol.ttest, rhwcol.ttest)


## ADDITIONAL ANALYSES

## 4. SPLOM OF ALL CONTINUOUS VARIABLES

ehr.df <- dbGetQuery(ehr,
	"SELECT monolith || ' ' || position AS key,",
		"rwc AS rsoil, rsoil AS rh, rr,",
		"swc, mrootmass, msoiltemp, treatment",
	"FROM ehr, rr",
	"WHERE ehr.pid = rr.pid",
		"AND ((temperature = 7.8 AND treatment = 'warmed')",
			"OR (temperature = 4.0 AND treatment = 'control'))",
	"ORDER BY treatment, rr.pid DESC")
	
splom(ehr.df[2:7], subset = ehr.df$rr < 5,
	groups = ehr.df$treatment, au = T, cex = 0.7)

## 5. NEE

nee <- dbGetQuery(ehr, "SELECT * FROM NEE")
nee$date <- as.Date(nee$date, format = "%m-%d-%Y")

xyplot(mean24hnee ~ date, groups = ecocell,
	data = nee,
	type = "b",
	pch = 16,
	col = c("gray", "black", "gray", "black"),
	lty = 1,
	key = list(points = list(
			pch = 16,
			col = c("gray", "black")),
		text = list(lab = c("warmed", "control")),
		columns = 1))