Load packages (ggplot2, ggpubr, dplyr)
library("ggplot2", lib.loc="~/R/win-library/3.6")
library("ggpubr", lib.loc="~/R/win-library/3.6")
library("dplyr", lib.loc="~/R/win-library/3.6")
Import data
decomp = read.csv("decomp.csv")
head(decomp)
Filter data to make Creswick-only dataframe
cres<-decomp %>% filter(site=="Creswick")
From Creswick only dataframe, filter out Poa only
cresPoa<-cres %>% filter(species=="Poa")
head(cresPoa)
Check that the data is looking like it should on a scatterplot:
ggplot(cresPoa,aes(x=time,y=omr))+geom_point()
It is. I used the nls() function in R to fit the Weibull curve to the data. Here's the code for fitting the curve
nls(omr ~ exp(-((time/beta)^alpha)),data=cresPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: cresPoa
## beta alpha
## 132.8422 0.8542
## residual sum-of-squares: 0.006071
##
## Algorithm "port", convergence message: relative convergence (4)
cresPoaFun <- function(time) {exp(-(time/132.8422)^0.8542)}
cres_means<-cres %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
cres_means
Plot of Creswick site with one Weibull curve fitted
ggplot(cres_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+stat_function(fun=cresPoaFun,color="#1B9E77")+scale_color_brewer(palette = "Dark2")+theme_classic(base_size = 16)+ggtitle("(b) 760 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.text.y = element_blank(),plot.title = element_text(size=14))+ylim(.6,1)+xlim(0,40)
OK, that is working. Now i need the constants for each of the three curves at each of the six sites.
First i'll do the two remaining curves at Creswick:
cresThem<-cres %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=cresThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: cresThem
## beta alpha
## 147.4200 0.7696
## residual sum-of-squares: 0.03398
##
## Algorithm "port", convergence message: relative convergence (4)
cresshade<-cres %>% filter(species=="shade")
nls(omr ~ exp(-((time/beta)^alpha)),data=cresshade, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: cresshade
## beta alpha
## 156.8380 0.9657
## residual sum-of-squares: 0.07873
##
## Algorithm "port", convergence message: relative convergence (4)
Now fit the curves for the other sites. THe output below this code chunk includes all the fitted constants.
#WIlson's Prom
wilsons<-decomp %>% filter(site=="Wilsons")
wilsonsPoa<-wilsons %>% filter(species=="Poa")
nls(omr ~ exp(-((time/beta)^alpha)),data=wilsonsPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: wilsonsPoa
## beta alpha
## 244.4141 0.6814
## residual sum-of-squares: 0.05426
##
## Algorithm "port", convergence message: relative convergence (4)
wilsonsThem<-wilsons %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=wilsonsThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: wilsonsThem
## beta alpha
## 127.7281 0.7787
## residual sum-of-squares: 0.05051
##
## Algorithm "port", convergence message: relative convergence (4)
wilsonsshade<-wilsons %>% filter(species=="shade")
nls(omr ~ exp(-((time/beta)^alpha)),data=wilsonsshade, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: wilsonsshade
## beta alpha
## 444.0954 0.6598
## residual sum-of-squares: 0.0412
##
## Algorithm "port", convergence message: relative convergence (4)
#Warrambeen
warram<-decomp %>% filter(site=="Warram")
warramPoa<-warram %>% filter(species=="Poa")
nls(omr ~ exp(-((time/beta)^alpha)),data=warramPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: warramPoa
## beta alpha
## 84.464 1.104
## residual sum-of-squares: 0.01391
##
## Algorithm "port", convergence message: relative convergence (4)
warramThem<-warram %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=warramThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: warramThem
## beta alpha
## 123.7939 0.8177
## residual sum-of-squares: 0.005097
##
## Algorithm "port", convergence message: relative convergence (4)
warramshade<-warram %>% filter(species=="shade")
nls(omr ~ exp(-((time/beta)^alpha)),data=warramshade, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: warramshade
## beta alpha
## 198.027 1.019
## residual sum-of-squares: 0.02402
##
## Algorithm "port", convergence message: relative convergence (4)
#Bungal
bungal<-decomp %>% filter(site=="Bungal")
bungalPoa<-bungal %>% filter(species=="Poa")
nls(omr ~ exp(-((time/beta)^alpha)),data=bungalPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: bungalPoa
## beta alpha
## 234.4375 0.6317
## residual sum-of-squares: 0.0166
##
## Algorithm "port", convergence message: relative convergence (4)
bungalThem<-bungal %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=bungalThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: bungalThem
## beta alpha
## 105.5861 0.8688
## residual sum-of-squares: 0.0176
##
## Algorithm "port", convergence message: relative convergence (4)
bungalshade<-bungal %>% filter(species=="shade")
nls(omr ~ exp(-((time/beta)^alpha)),data=bungalshade, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: bungalshade
## beta alpha
## 145.101 1.063
## residual sum-of-squares: 0.01633
##
## Algorithm "port", convergence message: relative convergence (4)
#Kinypanyial
kiny<-decomp %>% filter(site=="Kiny")
kinyPoa<-kiny %>% filter(species=="Poa")
nls(omr ~ exp(-((time/beta)^alpha)),data=kinyPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: kinyPoa
## beta alpha
## 341.0457 0.6708
## residual sum-of-squares: 0.02198
##
## Algorithm "port", convergence message: relative convergence (4)
kinyThem<-kiny %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=kinyThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: kinyThem
## beta alpha
## 128.50 1.03
## residual sum-of-squares: 0.04261
##
## Algorithm "port", convergence message: relative convergence (4)
kinyshade<-kiny %>% filter(species=="shade")
nls(omr ~ exp(-((time/beta)^alpha)),data=kinyshade, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: kinyshade
## beta alpha
## 1646.0126 0.5868
## residual sum-of-squares: 0.032
##
## Algorithm "port", convergence message: relative convergence (4)
#Korrak Korrak
korrak<-decomp %>% filter(site=="Korrak")
korrakPoa<-korrak %>% filter(species=="Poa")
nls(omr ~ exp(-((time/beta)^alpha)),data=korrakPoa, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: korrakPoa
## beta alpha
## 6.554e+04 3.093e-01
## residual sum-of-squares: 0.01335
##
## Algorithm "port", convergence message: relative convergence (4)
korrakThem<-korrak %>% filter(species=="Themeda")
nls(omr ~ exp(-((time/beta)^alpha)),data=korrakThem, start=list(beta=1, alpha = 1),algorithm = "port",lower = c(0.0001,0.0001))
## Nonlinear regression model
## model: omr ~ exp(-((time/beta)^alpha))
## data: korrakThem
## beta alpha
## 3023.5322 0.4984
## residual sum-of-squares: 0.01159
##
## Algorithm "port", convergence message: relative convergence (4)
#Weibull curve could not be fitted to Korrak shade treatment - omitted here.
Fit the funtions using the constants from the weibull funtions (i constructed these command lines in Excel from a table of constants):
cresThemFun <- function(time) {exp(-(time/127.5882)^0.8187)}
cresshadeFun <- function(time) {exp(-(time/133.379)^1.037)}
wilsonsPoaFun <- function(time) {exp(-(time/198.1044)^0.7364)}
wilsonsThemFun <- function(time) {exp(-(time/109.1393)^0.8398)}
wilsonsshadeFun <- function(time) {exp(-(time/374.8177)^0.6941)}
warramPoaFun <- function(time) {exp(-(time/84.464)^1.104)}
warramThemFun <- function(time) {exp(-(time/113.6334)^0.8511)}
warramshadeFun <- function(time) {exp(-(time/162.672)^1.107)}
bungalPoaFun <- function(time) {exp(-(time/200.4678)^0.6677)}
bungalThemFun <- function(time) {exp(-(time/92.6203)^0.9314)}
bungalshadeFun <- function(time) {exp(-(time/118.401)^1.181)}
kinyPoaFun <- function(time) {exp(-(time/294.0559)^0.7024)}
kinyThemFun <- function(time) {exp(-(time/111.121)^1.107)}
kinyshadeFun <- function(time) {exp(-(time/1204.023)^0.629)}
korrakPoaFun <- function(time) {exp(-(time/65540)^0.3093)}
korrakThemFun <- function(time) {exp(-(time/3023.5322)^0.4984)}
Need to create a plot for each site with data plotted and three fitted weibull curves
#Creswick
cres_means<-cres %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
five2<-ggplot(cres_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=cresPoaFun,color="#1B9E77")+stat_function(fun=cresThemFun,color="#7570B3")+stat_function(fun=cresshadeFun,color="#D95F02")+theme_classic(base_size = 16)+ggtitle("(e) 760 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
#Wilson's Prom
wilsons_means<-wilsons %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
six2<-ggplot(wilsons_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=wilsonsPoaFun,color="#1B9E77")+stat_function(fun=wilsonsThemFun,color="#7570B3")+stat_function(fun=wilsonsshadeFun,color="#D95F02")+theme_classic(base_size = 16)+ggtitle("(f) 890 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
#Warrambeen
warram_means<-warram %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
three2<-ggplot(warram_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=warramPoaFun,color="#1B9E77")+stat_function(fun=warramThemFun,color="#7570B3")+stat_function(fun=warramshadeFun,color="#D95F02")+theme_classic(base_size = 16)+ggtitle("(c) 560 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
#Bungal
bungal_means<-bungal %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
four2<-ggplot(bungal_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=bungalPoaFun,color="#1B9E77")+stat_function(fun=bungalThemFun,color="#7570B3")+stat_function(fun=bungalshadeFun,color="#D95F02")+theme_classic(base_size = 16)+ggtitle("(d) 650 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
#Kinypanyial
kiny_means<-kiny %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
two2<-ggplot(kiny_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=kinyPoaFun,color="#1B9E77")+stat_function(fun=kinyThemFun,color="#7570B3")+stat_function(fun=kinyshadeFun,color="#D95F02")+theme_classic(base_size = 16)+ggtitle("(b) 420 mm/yr")+xlab("Time (weeks)")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
#Korrak
korrak_means<-korrak %>% group_by(species,time) %>% summarise(mean=mean(omr),n=n(), sd=sd(omr), se=sd/sqrt(n))
one<-ggplot(korrak_means, aes(x=time, y=mean,color=species))+geom_point(position=position_dodge(0.5))+geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.2,position=position_dodge(0.5))+scale_color_brewer(palette = "Dark2")+stat_function(fun=korrakPoaFun,color="#1B9E77")+stat_function(fun=korrakThemFun,color="#7570B3")+theme_classic(base_size = 16)+ggtitle("(a) 380 mm/yr")+theme(legend.position ="none")+theme(axis.title.y = element_blank(),axis.title.x = element_blank(),plot.title = element_text(size=14))+ylim(0.6,1)+xlim(0,40)
plot<-ggarrange(one,two2,three2,four2,five2,six2,ncol = 6, nrow = 1)
annotate_figure(plot,left = text_grob("Organic matter remaining (%)", rot = 90),bottom = text_grob("Time (weeks)"))
ggsave("weibull_All.tiff",width=20,height=4.2,dpi=600)
ggsave("weibull_All.jpg",width=14,height=3,dpi=600)