---
title: "Sapflow analysis"
author: "Aaron"
date: "Wednesday, July 23, 2014"
output: html_document
---
Data management: HF SAP FLOW AND HF WARM ANTS MERGE
6 June 2014
Liza Nicoll

load necessary libraries
```{r}
library(RCurl)
library(plyr)
library(ggplot2)
#library(GGally)
#library(lattice)

#useful functions for printing correlation on scatterplot matrix when using pairs function

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex=2)
         #cex = cex.cor * r)
    }
    
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

```

read in data and get organized
```{r}
ants<-getURL("http://harvardforest.fas.harvard.edu/sites/harvardforest.fas.harvard.edu/files/data/p11/hf113/hf113-02-hf-chamber-since-2009.csv", ssl.verifypeer = FALSE)
ants<-read.csv(text=ants)
ants<-subset(ants, year==2011)
str(ants)

sap<-read.csv("hf_sapflow_2011_golden.csv")
summary(sap)
str(sap)
sap$Chamber<-sap$tree

#change doy to reflect next doy for 2330 hour
sap$doy<-ifelse(sap$time==2330, sap$doy+1, sap$doy)

#AGGREGATE SAP FLOW BY HOUR (this is assuming time is top of hour)
sapflow_h<-aggregate(sapflow ~ year+doy+hour+Chamber+sensor+dbh, sum, data=sap)
sapflow_vpd <- aggregate(vpd ~ year+doy+hour+Chamber+sensor+dbh, mean, data=sap)

sapflow_h_vpd<- merge(sapflow_h, sapflow_vpd)

#MERGE SAPFOW WITH WARM ANTS MICROCLIM
sap_ants_all<-merge(sapflow_h_vpd, ants, by=c("doy","hour","Chamber"),all.x=TRUE, all.y=FALSE)
summary(sap_ants_all)
names(sap_ants_all)
sap_ants<-sap_ants_all[,c(9:13, 1:3, 5:8, 14:15, 17:26, 27:35, 37:43)]
colnames(sap_ants)[1] <- "year"
names(sap_ants)
summary(sap_ants)


sap_ants<-arrange(sap_ants, Chamber, sensor, doy, hour)

write.csv(sap_ants,"hf_sapflow_microclim_2011.csv", row.names=F)
```

Initial plot of all data
```{r}

ggplot(sap_ants, aes(hour, sapflow)) +
         geom_point(aes(colour=sensor)) +
         facet_grid(Chamber ~ doy) +
         theme_bw()

```

Aggregate multiple temperature variables (means of duplicate sensors)
```{r}
sap_ants$X.AirTemp <- (sap_ants$CAT1_Avg+sap_ants$CAT2_Avg)/2
sap_ants$X.SoilOTemp <- (sap_ants$CSTo1_Avg+sap_ants$CSTo2_Avg)/2
sap_ants$X.SoilITemp <- (sap_ants$CSTI1_Avg+sap_ants$CSTI2_Avg)/2

sap_ants$Min.AirTemp <- (sap_ants$CAT1_Min+sap_ants$CAT2_Min)/2
sap_ants$Min.SoilOTemp <- (sap_ants$CSTo1_Min+sap_ants$CSTo2_Min)/2
sap_ants$Min.SoilITemp <- (sap_ants$CSTI1_Min+sap_ants$CSTI2_Min)/2

sap_ants$Max.AirTemp <- (sap_ants$CAT1_Max+sap_ants$CAT2_Max)/2
sap_ants$Max.SoilOTemp <- (sap_ants$CSTo1_Max+sap_ants$CSTo2_Max)/2
sap_ants$Max.SoilITemp <- (sap_ants$CSTI1_Max+sap_ants$CSTI2_Max)/2


#drop individual sensors

sap_work <- sap_ants[,c(1:12, 41:42, 49, 43:48, 19:21, 29:31, 38:40)]
names(sap_work)[22:30] <- c("X.PAR", "X.RH", "X.SoilMoisture", "Min.PAR", "Min.RH", "Min.SoilMoisture", "Max.PAR", "Max.RH", "Max.SoilMoisture")

sap_work$sat.vap.pres.cham <- 613.75*exp((17.502*sap_work$X.AirTemp)/(240.97+sap_work$X.AirTemp))
sap_work$act.vap.pres.cham <- sap_work$sat.vap.pres.cham*sap_work$X.RH/100
sap_work$vpd.cham.KPa <- (sap_work$sat.vap.pres.cham - sap_work$act.vap.pres.cham)/1000

names(sap_work)
summary(sap_work)
```

Set up files for 24-hour and nighttime

```{r}

sap.24hr <- sap_work
sap.night <- sap_work[sap_work$X.PAR<=0,]

```

next plots
```{r}

ggplot(sap.24hr, aes(hour, sapflow)) +
         geom_point(aes(colour=sensor)) +
         facet_grid(Chamber ~ doy) +
         theme_bw()

ggplot(sap.night, aes(hour, sapflow)) +
         geom_point(aes(colour=sensor)) +
         facet_grid(Chamber ~ doy) +
         theme_bw()

```

[Day of Year x chamber x sensor] summary values (maxes, sums for 24 hour and nighttime)

```{r}

#all by doy

sumsap.24 <- aggregate(sapflow~doy+Chamber+sensor, data=sap.24hr, sum, na.rm=TRUE)
sumsap.night <-aggregate(sapflow~doy+Chamber+sensor, data=sap.night, sum, na.rm=TRUE)
maxsap.24 <- aggregate(cbind(sapflow, dbh)~doy+Chamber+sensor, data=sap.24hr, max, na.rm=TRUE)
maxsap.night <- aggregate(cbind(sapflow, dbh)~doy+Chamber+sensor, data=sap.24hr, max, na.rm=TRUE)
meanenv.24 <- aggregate(sap.24hr[,c(13:30,33)], list(doy=sap.24hr$doy, Chamber=sap.24hr$Chamber), mean, na.rm=TRUE)
meanenv.night <- aggregate(sap.night[,c(13:30,33)], list(doy=sap.night$doy, Chamber=sap.night$Chamber), mean, na.rm=TRUE)

names(sumsap.24)[4] <- "S.sapflow.24"
names(sumsap.night)[4] <- "S.sapflow.night"
names(maxsap.24)[4] <- "Max.sapflow.24"
names(maxsap.night)[4] <- "Max.sapflow.night"

night.names <- names(meanenv.night)

for (i in 3:21) {
  night.names[i] <- paste("night",night.names[i], sep=".")
    }

names(meanenv.night) <- night.names
                            

sumsaps <- merge(sumsap.24,sumsap.night, all.x=TRUE)
maxsaps <- merge(maxsap.24, maxsap.night, all.x=TRUE)

sap.flows <- merge(sumsaps, maxsaps, all.x=TRUE)


envs <- merge(meanenv.24, meanenv.night, all.x=TRUE)

envs.2 <- na.omit(envs)

```


Effects of chambers on environmental variables

Uncomment pdf() and dev.off() to write to pdf

```{r}


tdls <- read.csv("target-delta-labels.csv")

envs.chambers <- merge(tdls, envs.2, all.x=TRUE)


#pdf("chamber-environmental-vars.pdf", width=8, height=8)

ggplot(envs.chambers, aes(x=X.AirTemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean air temperature") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=X.SoilOTemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil temperature (organic layer)") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=X.SoilITemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil temperature (mineral layer)") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=X.RH)) + 
         geom_histogram(binwidth=10,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean relative humidity") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=X.PAR)) + 
         geom_histogram(binwidth=25,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean PAR") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=X.SoilMoisture)) + 
         geom_histogram(binwidth=.05,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil moisture") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

ggplot(envs.chambers, aes(x=vpd.cham.KPa)) + 
         geom_histogram(binwidth=.1,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Estimated vpd (KPa)") +
         ylab("Number of observations") +
         ggtitle("24-hour data")

#dev.off()

#nighttime data

#pdf("chamber-nighttime-environmental-vars.pdf", width=8, height=8)

ggplot(envs.chambers, aes(x=night.X.AirTemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean air temperature") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.X.SoilOTemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil temperature (organic layer)") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.X.SoilITemp)) + 
         geom_histogram(binwidth=2.5,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil temperature (mineral layer)") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.X.RH)) + 
         geom_histogram(binwidth=10,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean relative humidity") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.X.PAR)) + 
         geom_histogram(binwidth=.1,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean PAR") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.X.SoilMoisture)) + 
         geom_histogram(binwidth=.05,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Mean soil moisture") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

ggplot(envs.chambers, aes(x=night.vpd.cham.KPa)) + 
         geom_histogram(binwidth=.1,colour="grey", fill="grey50") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw()+
         xlab("Estimated vpd (KPa)") +
         ylab("Number of observations") +
         ggtitle("Nighttime data")

#dev.off()

```




Look at correlations among variables

Uncomment pdf() and dev.off() to write to pdf

```{r}

means.cols.24hr <- c(3:5,12:14,21)
means.cols.night <-c(22:24, 31:33, 40)

#pdf("envcorrs.pdf", width=8, height=8)

pairs(envs.2[means.cols.24hr], lower.panel=panel.smooth, pch=".", upper.panel=panel.cor, diag.panel=panel.hist, main="Correlations among environmental variables: 24-hour data: ") 

pairs(envs.2[means.cols.night], lower.panel=panel.smooth, pch=".", upper.panel=panel.cor, diag.panel=panel.hist, main="Correlations among environmental variables: Nighttime data: ") 

#dev.off()
```

PCA of environmental data to deal with these correlations
Result is a "temperature" axis (PC1), "rH, PAR, vpd" axis (PC2), and "soil moisture" axis (PC3)

```{r}
#24-hour data

pca.24hr <- prcomp(envs.2[,3:21], center=TRUE, scale.=TRUE)
screeplot(pca.24hr)
summary(pca.24hr)
biplot(pca.24hr, main="24hr")
pca.24hr$rotation[,1:3]
pca.24hr.top3 <- data.frame(envs.2[,1:2], pca.24hr$x[,1:3])

#night data

pca.night <- prcomp(envs.2[,22:40], center=TRUE, scale.=TRUE)
screeplot(pca.night)
summary(pca.night)
biplot(pca.night, main="night")
pca.night$rotation[,1:3]
pca.night.top3 <- data.frame(envs.2[,1:2], pca.night$x[,1:3])
names(pca.night.top3)[3:5]<-c("night.PC1", "night.PC2", "night.PC3")

```
#combine

```{r}
pca.env <- merge(pca.24hr.top3, pca.night.top3, all.x=TRUE)

sap.env<-merge(sap.flows, pca.env, all.x=TRUE)

sap.env <- merge(tdls, sap.env, all.x=TRUE)
```

plots of 24 hour data

Uncomment pdf() and dev.off() to write to pdf

```{r}

#separated by chambers
#24 hour data

#pdf("24hr-sapflow.pdf", height=8, width=8)
ggplot(sap.env, aes(x=(-1)*PC1, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("First principal axis score (low temperature --> high temperature)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")
         

ggplot(sap.env, aes(x=(-1)*PC2, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Second principal axis score (low RH, high PAR, high VPD --> high RH, low PAR, low VPD)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")

ggplot(sap.env, aes(x=(-1)*PC3, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Third principal axis score (low soil moisture --> high soil moisture)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")

#all chambers pooled
ggplot(sap.env, aes(x=(-1)*PC1, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("First principal axis score (low temperature --> high temperature)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")

ggplot(sap.env, aes(x=(-1)*PC2, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Second principal axis score (low RH, high PAR, high VPD --> high RH, low PAR, low VPD)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")

ggplot(sap.env, aes(x=(-1)*PC3, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Third principal axis score (low soil moisture --> high soil moisture)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")

#dev.off()

```

plots of nighttime data

Uncomment pdf() and dev.off() to write to pdf

```{r}
#separated by chambers
#nighttime data

#pdf("nighttime-sapflow.pdf", height=8, width=8)
ggplot(sap.env, aes(x=(-1)*night.PC1, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("First principal axis score (low temperature --> high temperature)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")
         

ggplot(sap.env, aes(x=(-1)*night.PC2, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Second principal axis score (low RH, high VPD --> high RH, low VPD)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

ggplot(sap.env, aes(x=(-1)*night.PC3, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Third principal axis score (low soil moisture --> high soil moisture)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

#all chambers pooled
ggplot(sap.env, aes(x=(-1)*night.PC1, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("First principal axis score (low temperature --> high temperature)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

ggplot(sap.env, aes(x=(-1)*night.PC2, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Second principal axis score (low RH, high PAR, high VPD --> high RH, low PAR, low VPD)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

ggplot(sap.env, aes(x=(-1)*night.PC3, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
#         facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Third principal axis score (low soil moisture --> high soil moisture)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

#dev.off()
```

Regressions

```{r}

south.sensor.24hr.lm <- lm(S.sapflow.24 ~ PC1 + PC2 + PC3, data=sap.env[sap.env$sensor=="a",])

north.sensor.24hr.lm <- lm(S.sapflow.24 ~ PC1 + PC2 + PC3, data=sap.env[sap.env$sensor=="b",])


south.sensor.night.lm <- lm(S.sapflow.night ~ night.PC1 + night.PC2 + night.PC3, data=sap.env[sap.env$sensor=="a",])


north.sensor.night.lm <- lm(S.sapflow.night ~ night.PC1 + night.PC2 + night.PC3, data=sap.env[sap.env$sensor=="b",])

```

Results of regression: 

1. South sensors, 24 hour data

Only PC1 (temperature) significantly affects sapflow

```{r}
summary(south.sensor.24hr.lm)

```

2. South sensors, nighttime data

No environmental variables significantly affect sapflow

```{r}
summary(south.sensor.night.lm)

```


3. North sensors, 24 hour data

PC1 (temperature) and PC3 (soil moisture) significantly affect sapflow

```{r}
summary(north.sensor.24hr.lm)

```

4. North sensors, nighttime data

Only PC3 (soil moisture) significantly affects sapflow

```{r}
summary(north.sensor.night.lm)

```

Re-do as glm with sensor as factor + interactions


```{r}
sapflow.24hr.glm <- glm(S.sapflow.24 ~ sensor*(PC1 + PC2 + PC3), data=sap.env)

sapflow.night.glm <- glm(S.sapflow.night ~ sensor*(night.PC1 + night.PC2 + night.PC3), data=sap.env)


```

Results of glms with sensor as factor + interactions

1. For 24-hour data, PC1 (temperature) and interaction between sensor and PC1 both have significant effects on sapflow. PC3 (soil mositure) is trend (P=0.07)

2. For nighttime data, no terms significantly affect sapflow.


```{r}

summary(sapflow.24hr.glm)
summary(sapflow.night.glm)

```

Effect of preceding Winter temps (Dec., Jan., Feb.) on first four sap flows (doys 126, 127, 131, 132)

First trried using Growing Degree Days = (Tmax + Tmin)/2 + Tbase, where Tbase = 10 (degrees C), and GDD = 0 if mean daily temp < 10 (note GDDs deleted from code)

All growing degree days in range = 0, so use Heating Degree Days instead.

```{r}

#generate subset of data 
#read warm ants data from hf archive

#get url
ants_url<-getURL("http://harvardforest.fas.harvard.edu/sites/harvardforest.fas.harvard.edu/files/data/p11/hf113/hf113-02-hf-chamber-since-2009.csv")

#read url as csv
ants_all<-read.csv(text=ants_url)

#deselect variables of interest
ants1<-ants_all[,c(1:10,12:15,19:20,22:25,29:30,32:35)]

#subset data by dates of interest
ants2010<-subset(ants, year==2010 & doy>=335)
ants2011<-subset(ants, year==2011 & doy<=59)

ants_sub<-rbind(ants2010,ants2011)

ants_sub1<-aggregate(ants_sub[,c(9:26)], list(year=ants_sub$year, doy=ants_sub$doy, Chamber=ants_sub$Chamber), mean)

ants_sub1$HDD1 <- ifelse(ants_sub1$CAT1_Avg <= 18, 18-ants_sub1$CAT1_Avg, 0)
ants_sub1$HDD2 <- ifelse(ants_sub1$CAT2_Avg <= 18, 18-ants_sub1$CAT2_Avg, 0)

chamber_HDD <- aggregate(ants_sub1[,22:23],list(Chamber=ants_sub1$Chamber), sum)
chamber_HDD$HDD <- apply(chamber_HDD[,2:3], 1, mean, na.rm=TRUE)



```

merge chamber heating degree days with sap.env
subset datafile to only have doys 126, 127, 131, 132 (i.e., doy < 133)

```{r}

sap.env.hdd <- merge(sap.env, chamber_HDD)
early.sap.env.hdd <- sap.env.hdd[sap.env.hdd$doy < 133,]
```

some plots

```{r}
pdf("earlyseasonsapflow.pdf", height=8, width=8)

ggplot(early.sap.env.hdd, aes(x=HDD, y=S.sapflow.24/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
         #facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("Heating Degree Days (warmer --> colder)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("24-hour data")


ggplot(early.sap.env.hdd, aes(x=HDD, y=S.sapflow.night/1000, colour=sensor)) +
        geom_point(aes(colour=sensor)) +
         stat_smooth(method="lm") +
        # facet_wrap(~TargetDelta, nrow=3) +
         theme_bw() +
         scale_colour_discrete(name="Sensor position", breaks = c("a","b"), labels=c("south","north")) +
         xlab("First principal axis score (low temperature --> high temperature)") +
         ylab(expression(paste("Sap flow (kg ", H[2], "O " %.% m^-2 %.% d^-1," )"))) +
         ggtitle("Nighttime data")

dev.off()
```

re-run regressions on early season sapflow and including HDD

```{r}



south.sensor.24hr.lm2 <- lm(S.sapflow.24 ~ HDD+PC1 + PC2 + PC3, data=early.sap.env.hdd[early.sap.env.hdd$sensor=="a",])

north.sensor.24hr.lm2 <- lm(S.sapflow.24 ~ HDD+PC1 + PC2 + PC3, data=early.sap.env.hdd[early.sap.env.hdd$sensor=="b",])


south.sensor.night.lm2 <- lm(S.sapflow.night ~ HDD+night.PC1 + night.PC2 + night.PC3, data=early.sap.env.hdd[early.sap.env.hdd$sensor=="a",])


north.sensor.night.lm2 <- lm(S.sapflow.night ~ HDD+night.PC1 + night.PC2 + night.PC3, data=early.sap.env.hdd[early.sap.env.hdd$sensor=="b",])
```

summarize regressions

Salient results:

PC1 dominates regressions for 24 hr data.

HDD significant only for south sensor night data

```{r}

summary.aov(south.sensor.24hr.lm2)

summary.aov(north.sensor.24hr.lm2)

summary.aov(south.sensor.night.lm2)

summary.aov(north.sensor.night.lm2)


```


```{r}

early.sapflow.24hr.glm <- glm(S.sapflow.24 ~ sensor*(HDD+PC1 + PC2 + PC3), data=early.sap.env.hdd)

early.sapflow.night.glm <- glm(S.sapflow.night ~ sensor*(HDD+PC1 + PC2 + PC3), data=early.sap.env.hdd)

```



Single regressions for abstract
```{r}

#24 hour early season sapflow as a function of temperature (PC1)

summary(lm(S.sapflow.24/1000 ~ PC1, data=early.sap.env.hdd) )


plot(S.sapflow.24/1000~PC1, data=early.sap.env.hdd)
abline(lm(S.sapflow.24/1000 ~ PC1, data=early.sap.env.hdd))

#nighttime early season sapflow as a function of HDD

summary(lm(S.sapflow.night/1000~HDD, data=early.sap.env.hdd))
plot(S.sapflow.night/1000~HDD, data=early.sap.env.hdd)
abline(lm(S.sapflow.night/1000 ~ HDD, data=early.sap.env.hdd))

#all season 24 hour sapflow as a function of temperature (PC1)

summary(lm(S.sapflow.24/1000 ~ PC1, data=sap.env) )


plot(S.sapflow.24/1000~PC1, data=sap.env)
abline(lm(S.sapflow.24/1000 ~ PC1, data=sap.env))

#all season 24 hour sapflow as a function of soil moisture (PC3)

summary(lm(S.sapflow.24/1000 ~ PC3, data=sap.env) )


plot(S.sapflow.24/1000~PC3, data=sap.env)
abline(lm(S.sapflow.24/1000 ~ PC3, data=sap.env))

#all season nighttime sapflow as a function of soil moisture (PC3)

summary(lm(S.sapflow.night/1000 ~ PC3, data=sap.env) )


plot(S.sapflow.night/1000~PC3, data=sap.env)
abline(lm(S.sapflow.night/1000 ~ PC3, data=sap.env))

```