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)
Final plot as a JPG

Final plot as a JPG