### R 2.1.0 DATA ANALYSIS FOR YOUNG SAGE REMOVAL ## ACQUIRE DATA require(RODBC) ysr <- odbcConnect("ysr") age <- sqlQuery(ysr, "SELECT age FROM age") herbcover <- sqlQuery(ysr, paste("SELECT date, block, treatment,", "AVG(grass + sedge + spike_rush + rush + forb) AS herbcover,", "AVG(sagebrush) AS sagecover", "FROM percentcover", "GROUP BY date, block, treatment", "ORDER BY date, block, treatment")) herbcover$date <- as.POSIXct(strptime(herbcover$date, "%Y-%m-%d")) herbcover$block <- factor(herbcover$block) herbcover$treatment <- ordered(herbcover$treatment, levels = c("s+", "s-")) shallow <- sqlQuery(ysr, paste("SELECT date, block, treatment,", "avg(tdr) AS tdr", "FROM tdr", "WHERE depth IN ('0-15', '15-30', '0-30')", "GROUP BY date, block, treatment", "ORDER BY date, block, treatment")) deep <- sqlQuery(ysr, paste("SELECT date, block, treatment,", "avg(tdr) AS tdr", "FROM tdr", "WHERE depth IN ('30-45', '45-60')", "GROUP BY date, block, treatment", "ORDER BY date, block, treatment")) reallydeep <- sqlQuery(ysr, paste("SELECT date, block, treatment,", "avg(tdr) AS tdr", "FROM tdr", "WHERE depth IN ('60-90')", "GROUP BY date, block, treatment", "ORDER BY date, block, treatment")) shallow$depth <- "0-30" deep$depth <- "30-60" reallydeep$depth <- "60-90" tdr <- rbind(shallow, deep, reallydeep) attributes(tdr)$row.names <- 1:length(tdr$tdr) rm(shallow, deep, reallydeep) tdr$date <- as.POSIXct(strptime(tdr$date, "%Y-%m-%d")) tdr$treatment <- factor(tdr$treatment) tdr$block <- factor(tdr$block) tdr$depth <- factor(tdr$depth) nitrogen <- sqlQuery(ysr, paste("SELECT date, block, treatment,", "AVG(net_mineralization/14) AS nmin, AVG(net_nitrification/14) AS nnit", "FROM nitrogen", "GROUP BY date, block, treatment", "ORDER BY date, block, treatment")) nitrogen$date <- as.POSIXct(strptime(nitrogen$date, "%Y-%m-%d")) nitrogen$treatment <- factor(nitrogen$treatment) nitrogen$block <- factor(nitrogen$block) ## EFFECT SIZE CALCULATIONS old <- getOption("digits") options(digits = 2) attach(herbcover) d <- data.frame() for (n in 1:4) {t <- t.test(herbcover$herbcover ~ treatment, paired = T, subset = date == unique(date)[n]) p = t$estimate/mean(herbcover$herbcover[date == unique(date)[n] & treatment == "s+"], na.rm = T)*100 pl = t$conf.int[1]/mean(herbcover$herbcover[date == unique(date)[n] & treatment == "s+"], na.rm = T)*100 ph = t$conf.int[2]/mean(herbcover$herbcover[date == unique(date)[n] & treatment == "s+"], na.rm = T)*100 d <- rbind(d, data.frame(date = unique(date)[n], effect = as.numeric(t$estimate), low = t$conf.int[1], high = t$conf.int[2], percent = p, pl = pl, ph = ph)) } row.names(d) <- 1:length(d$date) ; detach(herbcover) print("Herbcover") print(d) tdrm <- tapply(tdr$tdr, list(tdr$date, tdr$depth, tdr$block, tdr$treatment), mean) # take plot means d <- data.frame() for (da in 1:9) for (de in 1:3) { if (is.na(tdrm[da,de,1,1]) == F) { t <- t.test(tdrm[da,de,,2], tdrm[da,de,,1], paired = T) # t.tests for each date <- dimnames(tdrm)[[1]][da] ; depth <- dimnames(tdrm)[[2]][de] p = t$estimate/mean(tdrm[da,de,,1], na.rm = T)*100 pl = t$conf.int[1]/mean(tdrm[da,de,,1], na.rm = T)*100 ph = t$conf.int[2]/mean(tdrm[da,de,,1], na.rm = T)*100 d <- rbind(d, data.frame(date = date, depth = depth, effect = as.numeric(t$estimate), low = t$conf.int[1], high = t$conf.int[2], percent = p, pl = pl, ph = ph)) } } row.names(d) <- 1:length(d$date) print("Soil Moisture") print(d) tm <- t.test(tapply(nitrogen01$net_mineralization[nitrogen01$treatment == "s-"], nitrogen01$block[nitrogen01$treatment == "s-"], mean), tapply(nitrogen01$net_mineralization[nitrogen01$treatment == "s+"], nitrogen01$block[nitrogen01$treatment == "s+"], mean), paired = T) tn <- t.test(tapply(nitrogen01$net_nitrification[nitrogen01$treatment == "s-"], nitrogen01$block[nitrogen01$treatment == "s-"], mean), tapply(nitrogen01$net_nitrification[nitrogen01$treatment == "s+"], nitrogen01$block[nitrogen01$treatment == "s+"], mean), paired = T) tmd <- c(tm$estimate, tm$conf.int[1], tm$conf.int[2], tm$estimate/mean(nitrogen01$net_mineralization[nitrogen01$treatment == "s+"], na.rm = T)*100, tm$conf.int[1]/mean(nitrogen01$net_mineralization[nitrogen01$treatment == "s+"], na.rm = T)*100, tm$conf.int[2]/mean(nitrogen01$net_mineralization[nitrogen01$treatment == "s+"], na.rm = T)*100) tnd <- c(tn$estimate, tn$conf.int[1], tn$conf.int[2], tn$estimate/mean(nitrogen01$net_nitrification[nitrogen01$treatment == "s+"], na.rm = T)*100, tn$conf.int[1]/mean(nitrogen01$net_nitrification[nitrogen01$treatment == "s+"], na.rm = T)*100, tn$conf.int[2]/mean(nitrogen01$net_nitrification[nitrogen01$treatment == "s+"], na.rm = T)*100) names(tnd) <- names(tmd) <- c("effect", "low", "high", "percent", "pl", "ph") print("Net mineralization: 12 June 2001") print(tmd) print("Net nitrification: 12 June 2001") print(tnd) tm <- t.test(tapply(nitrogen05$net_mineralization[nitrogen05$treatment == "s-"], nitrogen05$block[nitrogen05$treatment == "s-"], mean), tapply(nitrogen05$net_mineralization[nitrogen05$treatment == "s+"], nitrogen05$block[nitrogen05$treatment == "s+"], mean), paired = T) tn <- t.test(tapply(nitrogen05$net_nitrification[nitrogen05$treatment == "s-"], nitrogen05$block[nitrogen05$treatment == "s-"], mean), tapply(nitrogen05$net_nitrification[nitrogen05$treatment == "s+"], nitrogen05$block[nitrogen05$treatment == "s+"], mean), paired = T) tmd <- c(tm$estimate, tm$conf.int[1], tm$conf.int[2], tm$estimate/mean(nitrogen05$net_mineralization[nitrogen05$treatment == "s+"], na.rm = T)*100, tm$conf.int[1]/mean(nitrogen05$net_mineralization[nitrogen05$treatment == "s+"], na.rm = T)*100, tm$conf.int[2]/mean(nitrogen05$net_mineralization[nitrogen05$treatment == "s+"], na.rm = T)*100) tnd <- c(tn$estimate, tn$conf.int[1], tn$conf.int[2], tn$estimate/mean(nitrogen05$net_nitrification[nitrogen05$treatment == "s+"], na.rm = T)*100, tn$conf.int[1]/mean(nitrogen05$net_nitrification[nitrogen05$treatment == "s+"], na.rm = T)*100, tn$conf.int[2]/mean(nitrogen05$net_nitrification[nitrogen05$treatment == "s+"], na.rm = T)*100) names(tnd) <- names(tmd) <- c("effect", "low", "high", "percent", "pl", "ph") print("Net mineralization: 3 August 2005") print(tmd) print("Net nitrification: 3 August 2005") print(tnd) options(digits = old) ## REGRESSION OF HERB PICKUP VS. SAGEBRUSH REMOVAL library(multcomp) as <- sqlQuery(ysr, "SELECT * FROM ageandsize") ass <- sqlQuery(ysr, "SELECT * FROM ageandsizesupp") hp <- data.frame(block = 1:6, herbpickup = c(herbcover$herbcover[2] - herbcover$herbcover[1], herbcover$herbcover[4] - herbcover$herbcover[3], herbcover$herbcover[6] - herbcover$herbcover[5], herbcover$herbcover[8] - herbcover$herbcover[7], herbcover$herbcover[26] - herbcover$herbcover[25], herbcover$herbcover[28] - herbcover$herbcover[27]), sagecoverplusS = herbcover$sagecover[c(1,3,5,7,25,27)], sagecoverminusS = tapply(as$canopy1 * as$canopy2, as$block, sum)/1225) hp$sagecoverminusS[1] <- hp$sagecoverminusS[1] + (ass$lt5cm[1] + ass$lt5cm[2]) * mean(as$canopy1[as$block == 1 & as$height < 5] * as$canopy2[as$block == 1 & as$height < 5])/1225 + (ass$btw5_10cm[1] + ass$btw5_10[2]) * mean(as$canopy1[as$block == 1 & as$height > 5 & as$height < 10] * as$canopy2[as$block == 1 & as$height > 5 & as$height < 10])/1225 + (ass$btw10_20[1] + ass$btw10_20[2]) * mean(as$canopy1[as$block == 1 & as$height > 10 & as$height < 20] * as$canopy2[as$block == 1 & as$height > 10 & as$height < 20])/1225 + (ass$gt20cm[1] + ass$gt20cm[2]) * mean(as$canopy1[as$block == 1 & as$height > 20] * as$canopy2[as$block == 1 & as$height > 20])/1225 hp$sagecoverminusS[2] <- hp$sagecoverminusS[2] + (ass$lt5cm[5] + ass$lt5cm[6]) * mean(as$canopy1[as$block == 2 & as$height < 5] * as$canopy2[as$block == 2 & as$height < 5])/1225 + (ass$btw5_10cm[5] + ass$btw5_10[6]) * mean(as$canopy1[as$block == 2 & as$height > 5 & as$height < 10] * as$canopy2[as$block == 2 & as$height > 5 & as$height < 10])/1225 + (ass$btw10_20[5] + ass$btw10_20[6]) * mean(as$canopy1[as$block == 2 & as$height > 10 & as$height < 20] * as$canopy2[as$block == 2 & as$height > 10 & as$height < 20])/1225 + (ass$gt20cm[5] + ass$gt20cm[6]) * mean(as$canopy1[as$block == 2 & as$height > 20] * as$canopy2[as$block == 2 & as$height > 20])/1225 hp$sagecoverminusSnorm <- hp$sagecoverminusS * mean(hp$sagecoverplusS/hp$sagecoverminusS) print("Correlation btw sagebrush removed and herb pickup, r2 and slope 95%CI") print(c(summary(lm(hp$herbpickup ~ hp$sagecoverminusSnorm))$r.squared, simint(lm(hp$herbpickup ~ hp$sagecoverminusSnorm))$conf.int[2,])) # splom(hp[-4]) ## PLOT DATA require(lattice) 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, col = "red") } pdf(file = "ysrplots.pdf", w = 11, h = 8.5) source("http://anthony.darrouzet-nardi.net/works/adnlattice.R") print(histogram(as.factor(age$age))) print(xyplot(herbcover ~ treatment | factor(date), groups = block, data = herbcover, panel = panel.ysr, layout = c(4,1,1), cex = 1.5)) print(xyplot(tdr ~ interaction(treatment, depth) | factor(date), data = tdr, groups = block, panel = panel.ysr, layout = c(4,3,1), skip = c(F,F,F,F, F,F,F,T, F,F,T,T))) print(xyplot(nnit ~ treatment | factor(date), groups = block, data = nitrogen, panel = panel.ysr, cex = 1.5), split = c(1,1,2,1), more = T) print(xyplot(nmin ~ treatment | factor(date), groups = block, data = nitrogen, panel = panel.ysr, cex = 1.5), split = c(2,1,2,1)) dev.off() rm(d, da, dr, date, de, depth, min.plot01, min.plot05, n, nit.plot01, nit.plot05, old, p, ph, pl, t, tdrm, tm, tmd, tn, tnd) # CLEANUP