Fire influences on needle decomposition: thermal tipping point in Pinus radiata needle carbon chemistry and soil nitrogen transformations

RScript written by Erinne Stirling 2018

R version 3.5.1 (2018-07-02) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

library(readr)#for reading data
library(dplyr)#for manipulating data
library(ggplot2)#for plotting
library(RColorBrewer)#for palettes
library(agricolae)#for stastics
library(nlme)
library(nlme2)

RESP_CUM<-read.csv("2018 Scorchin cum_resp.csv")

#----------
#Decay Curves
SCR_blank<-RESP_CUM%>%
  filter(Treat=="T000")#%>%
  #group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_blank$TIME,SCR_blank$CUM_RESP)

##----
fit_1b<-SCR_blank%>%
  filter(TYPE=="B_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_blank%>%
  filter(TYPE=="B_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_blank%>%
  filter(TYPE=="B_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_blank_ave<-data.frame(id=c(rep("Control",3)),
                            C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                            k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                            a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                            m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T260<-RESP_CUM%>%
  filter(Treat=="T0260")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T260$TIME,SCR_T260$CUM_RESP)

##----
fit_1b<-SCR_T260%>%
  filter(TYPE=="260_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T260%>%
  filter(TYPE=="260_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T260%>%
  filter(TYPE=="260_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T260_ave<-data.frame(id=c(rep("T260",3)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T85<-RESP_CUM%>%
  filter(Treat=="T085")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T85$TIME,SCR_T85$CUM_RESP)

##----
fit_1b<-SCR_T85%>%
  filter(TYPE=="85_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T85%>%
  filter(TYPE=="85_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T85%>%
  filter(TYPE=="85_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T85_ave<-data.frame(id=c(rep("T85",3)),
                        C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                        k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                        a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                        m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T100<-RESP_CUM%>%
  filter(Treat=="T100")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T100$TIME,SCR_T100$CUM_RESP)

##----
fit_1b<-SCR_T100%>%
  filter(TYPE=="100_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T100%>%
  filter(TYPE=="100_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=200,k=10,a=10,m=5))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T100%>%
  filter(TYPE=="100_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T100_ave<-data.frame(id=c(rep("T100",3)),
                        C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                        k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                        a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                        m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T135<-RESP_CUM%>%
  filter(Treat=="T135")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T135$TIME,SCR_T135$CUM_RESP)

##----
fit_1b<-SCR_T135%>%
  filter(TYPE=="135_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T135%>%
  filter(TYPE=="135_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T135%>%
  filter(TYPE=="135_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T135_ave<-data.frame(id=c(rep("T135",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T150<-RESP_CUM%>%
  filter(Treat=="T150")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T150$TIME,SCR_T150$CUM_RESP)

##----
fit_1b<-SCR_T150%>%
  filter(TYPE=="150_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T150%>%
  filter(TYPE=="150_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T150%>%
  filter(TYPE=="150_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T150_ave<-data.frame(id=c(rep("T150",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))


SCR_T178<-RESP_CUM%>%
  filter(Treat=="T178")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T178$TIME,SCR_T178$CUM_RESP)

##----
fit_1b<-SCR_T178%>%
  filter(TYPE=="178_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=200,k=10,a=10,m=5))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T178%>%
  filter(TYPE=="178_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T178%>%
  filter(TYPE=="178_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T178_ave<-data.frame(id=c(rep("T178",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T200<-RESP_CUM%>%
  filter(Treat=="T200")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T200$TIME,SCR_T200$CUM_RESP)

##----
fit_1b<-SCR_T200%>%
  filter(TYPE=="200_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T200%>%
  filter(TYPE=="200_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T200%>%
  filter(TYPE=="200_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T200_ave<-data.frame(id=c(rep("T200",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T225<-RESP_CUM%>%
  filter(Treat=="T225")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T225$TIME,SCR_T225$CUM_RESP)

##----
fit_1b<-SCR_T225%>%
  filter(TYPE=="225_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T225%>%
  filter(TYPE=="225_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T225%>%
  filter(TYPE=="225_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T225_ave<-data.frame(id=c(rep("T225",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T260<-RESP_CUM%>%
  filter(Treat=="T260")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T260$TIME,SCR_T260$CUM_RESP)

##----
fit_1b<-SCR_T260%>%
  filter(TYPE=="260_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T260%>%
  filter(TYPE=="260_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T260%>%
  filter(TYPE=="260_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T260_ave<-data.frame(id=c(rep("T260",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))

SCR_T320<-RESP_CUM%>%
  filter(Treat=="T320")#%>%
#group_by(STATUS,TYPE, TIME)#%>%
#summarise_at(vars(CUM_RESP), funs(mean))

plot(SCR_T320$TIME,SCR_T320$CUM_RESP)

##----
fit_1b<-SCR_T320%>%
  filter(TYPE=="320_1")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_1b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
fit_1b$pred <- predict(fit_a)
plot(fit_1b$TIME,fit_1b$CUM_RESP)
lines(fit_1b$TIME,fit_1b$pred, col="blue",lty=2)

##----
fit_2b<-SCR_T320%>%
  filter(TYPE=="320_2")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_2b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_a)
fit_2b$pred <- predict(fit_a)
plot(fit_2b$TIME,fit_2b$CUM_RESP)
lines(fit_2b$TIME,fit_2b$pred, col="blue",lty=2)

##----
fit_3b<-SCR_T320%>%
  filter(TYPE=="320_3")
fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=fit_3b,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_a)
fit_3b$pred <- predict(fit_a)
plot(fit_3b$TIME,fit_3b$CUM_RESP)
lines(fit_3b$TIME,fit_3b$pred, col="blue",lty=2)

SCR_T320_ave<-data.frame(id=c(rep("T320",3)),
                         C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1]),
                         k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1]),
                         a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1]),
                         m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1]))



Curves<-rbind(SCR_blank_ave,
              SCR_T260_ave,SCR_T85_ave,SCR_T100_ave,SCR_T135_ave,SCR_T150_ave,SCR_T178_ave,
              SCR_T200_ave,SCR_T225_ave,SCR_T260_ave,SCR_T320_ave)

#----------STATS
stats<-Curves%>%
  group_by(id)

anova(aov(C~id, data=stats)) #***
HSD.test(aov(C~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

anova(aov(k~id, data=stats)) #
#HSD.test(aov(k~id, data=stats), c("id"),
#         console=TRUE,alpha=0.05)

anova(aov(a~id, data=stats)) #***
HSD.test(aov(a~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

anova(aov(m~id, data=stats)) #***
HSD.test(aov(m~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

####------------Plot

p_curves<-Curves%>%
  group_by(id)%>% 
  summarize_at(vars(C,k,a,m), funs(mean))

write.csv(p_curves,"csv files/p_curves.csv")

p_curves<-read.csv("csv files/p_curves.csv")%>%
  select(-X)

CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME



fun_control <- function(x) {
  (p_curves[1,4]*p_curves[1,2])/(p_curves[1,4]+(p_curves[1,2]-p_curves[1,4])*exp(-p_curves[1,3]*x))+p_curves[1,5]*x
  }

fun_T260 <- function(x) {
  (p_curves[2,4]*p_curves[2,2])/(p_curves[2,4]+(p_curves[2,2]-p_curves[2,4])*exp(-p_curves[2,3]*x))+p_curves[2,5]*x
}

fun_T85 <- function(x) {
  (p_curves[3,4]*p_curves[3,2])/(p_curves[3,4]+(p_curves[3,2]-p_curves[3,4])*exp(-p_curves[3,3]*x))+p_curves[3,5]*x
}

fun_T100 <- function(x) {
  (p_curves[4,4]*p_curves[4,2])/(p_curves[4,4]+(p_curves[4,2]-p_curves[4,4])*exp(-p_curves[4,3]*x))+p_curves[4,5]*x
}

fun_T135 <- function(x) {
  (p_curves[5,4]*p_curves[5,2])/(p_curves[5,4]+(p_curves[5,2]-p_curves[5,4])*exp(-p_curves[5,3]*x))+p_curves[5,5]*x
}

fun_T150 <- function(x) {
  (p_curves[6,4]*p_curves[6,2])/(p_curves[6,4]+(p_curves[6,2]-p_curves[6,4])*exp(-p_curves[6,3]*x))+p_curves[6,5]*x
}

fun_T178 <- function(x) {
  (p_curves[7,4]*p_curves[7,2])/(p_curves[7,4]+(p_curves[7,2]-p_curves[7,4])*exp(-p_curves[7,3]*x))+p_curves[7,5]*x
}

fun_T200 <- function(x) {
  (p_curves[8,4]*p_curves[8,2])/(p_curves[8,4]+(p_curves[8,2]-p_curves[8,4])*exp(-p_curves[8,3]*x))+p_curves[8,5]*x
}

fun_T225 <- function(x) {
  (p_curves[9,4]*p_curves[9,2])/(p_curves[9,4]+(p_curves[9,2]-p_curves[9,4])*exp(-p_curves[9,3]*x))+p_curves[9,5]*x
}

fun_T260 <- function(x) {
  (p_curves[10,4]*p_curves[10,2])/(p_curves[10,4]+(p_curves[10,2]-p_curves[10,4])*exp(-p_curves[10,3]*x))+p_curves[10,5]*x
}

fun_T320 <- function(x) {
  (p_curves[11,4]*p_curves[11,2])/(p_curves[11,4]+(p_curves[11,2]-p_curves[11,4])*exp(-p_curves[11,3]*x))+p_curves[11,5]*x
}

Scorchin_lines<-ggplot(data.frame(x=c(0,14)),aes(x=x))+
  stat_function(fun=fun_control, colour="grey")+
  stat_function(fun=fun_T260, colour="lightblue")+
  stat_function(fun=fun_T85, colour="blue")+
  stat_function(fun=fun_T100, colour="darkblue")+
  stat_function(fun=fun_T135, colour="lightpink")+
  stat_function(fun=fun_T150, colour="orange")+
  stat_function(fun=fun_T178, colour="red")+
  stat_function(fun=fun_T200, colour="brown")+
  stat_function(fun=fun_T225, colour="green")+
  stat_function(fun=fun_T260, colour="yellow")+
  stat_function(fun=fun_T320, colour="black")+
  scale_y_continuous(limits=c(0,200))+
  ggsave("Schorchin_resp1a.svg", width=200, height=100, units="mm")

RESP_CUM<-read.csv("csv files/Resp_cum.csv")%>%
  ggplot(aes(x=TIME,y=CUM_RESP))+
  geom_point(aes(colour=Treat))+
  #scale_colour_manual(values=c("grey", "blue", "lightblue", "red", "orange"))+
  xlab("Time (days)") +
  ylab(expression(paste("Cumulative Respiration (mgCO2-C/gOC)"))) +
  scale_y_continuous(limits=c(0,200))+
  ggsave("Schorchin_resp1b.svg", width=200, height=100, units="mm")

#-------
#testing

SCR_blank$pred<-fun_control(SCR_blank$TIME)
BLANK.lm = lm(CUM_RESP~pred, data=SCR_blank)
summary(BLANK.lm)$r.squared

SCR_T260$pred<-fun_T260(SCR_T260$TIME)
T260.lm = lm(CUM_RESP~pred, data=SCR_T260)
summary(T260.lm)$r.squared

SCR_T85$pred<-fun_T85(SCR_T85$TIME)
T85.lm = lm(CUM_RESP~pred, data=SCR_T85)
summary(T85.lm)$r.squared

SCR_T100$pred<-fun_T100(SCR_T100$TIME)
T100.lm = lm(CUM_RESP~pred, data=SCR_T100)
summary(T100.lm)$r.squared

SCR_T135$pred<-fun_T135(SCR_T135$TIME)
T135.lm = lm(CUM_RESP~pred, data=SCR_T135)
summary(T135.lm)$r.squared

SCR_T150$pred<-fun_T150(SCR_T150$TIME)
T150.lm = lm(CUM_RESP~pred, data=SCR_T150)
summary(T150.lm)$r.squared

SCR_T178$pred<-fun_T178(SCR_T178$TIME)
T178.lm = lm(CUM_RESP~pred, data=SCR_T178)
summary(T178.lm)$r.squared

SCR_T200$pred<-fun_T200(SCR_T200$TIME)
T200.lm = lm(CUM_RESP~pred, data=SCR_T200)
summary(T200.lm)$r.squared

SCR_T225$pred<-fun_T225(SCR_T225$TIME)
T225.lm = lm(CUM_RESP~pred, data=SCR_T225)
summary(T225.lm)$r.squared

SCR_T260$pred<-fun_T260(SCR_T260$TIME)
T260.lm = lm(CUM_RESP~pred, data=SCR_T260)
summary(T260.lm)$r.squared

SCR_T320$pred<-fun_T320(SCR_T320$TIME)
T320.lm = lm(CUM_RESP~pred, data=SCR_T320)
summary(T320.lm)$r.squared