Load packages
library(readxl)
library(ggplot2)
library(nlme)
library(nlstools)
library(tidyr)
library(dplyr)
library(scales)
library(MASS)
library(effects)
Read data (age estimated as days since spawning, where day 0 = 10th December 2014)
read.csv("Acanthaster.juvenile.growth.data.csv",header=TRUE,strip.white=T)-> RAW;ls()
## [1] "RAW"
##subset data into two
CCA = RAW %>% filter(Diet=='CCA')
coral = RAW %>% filter(Diet=='coral')
Rearrange data to model the timing of the shift
# calculate standard error
SE <- function(x) {sd(x)/sqrt(length(x))}
Timing=RAW %>% group_by(ReefNameID,SiteName,days_sincespawning,Latitude,Longitude,Diet) %>% summarise_at("Size", funs(mean,var,SE,min,max,n()))
test=tidyr::gather(Timing,key,value,7:12)
test=tidyr::unite(test,Diet_stat,6:7,sep="_")
test=tidyr::spread(test,Diet_stat,value)
# replace NA in CCA_n and coral_n with 0
test$CCA_n[is.na(test$CCA_n)] <- 0
test$coral_n[is.na(test$coral_n)] <- 0
RAW_Timing=dplyr::mutate(test,Total=CCA_n + coral_n, PropCCA=CCA_n/Total)
Rearrange data to model the time at which 50% of population would have reached the size threshold
##create size categories
SizeClass=RAW %>% mutate(SizeThresh=cut(Size,breaks=c(0,8,70),labels=c("≤ 8 mm","> 8-70 mm")))
##convert Date to Day of Year
my.dates = as.POSIXlt(SizeClass$Date,format="%d/%m/%Y")
SizeClass$Day.of.Year.R= my.dates$yday
TIME=SizeClass %>% group_by(ReefNameID,SiteName,Day.of.Year.R,Latitude,Longitude,SizeThresh) %>% summarise_at("Size", funs(n()))
test=tidyr::gather(TIME,key,value,7)
test=tidyr::unite(test,sizethresh_stat,6:7,sep="_")
test=tidyr::spread(test,sizethresh_stat,value)
# replace NA in CCA_n and coral_n with 0
test$`≤ 8 mm_Size`[is.na(test$`≤ 8 mm_Size`)] <- 0
test$`> 8-70 mm_Size`[is.na(test$`> 8-70 mm_Size`)] <- 0
RAW_TIME=dplyr::mutate(test,Total=`≤ 8 mm_Size` + `> 8-70 mm_Size`, Prop_below_or_equal_threshold=`≤ 8 mm_Size`/Total)
Model growth trajectory for herbivorous 0+ starfish
# run model
CCA.lme.poly2 <- lme(Size ~ poly(days_sincespawning,2,raw=FALSE),random=~1|ReefNameID/SiteName, data=CCA, method='REML')##raw is FALSE to assess the significance of individual terms
summary(CCA.lme.poly2)
## Linear mixed-effects model fit by REML
## Data: CCA
## AIC BIC logLik
## 14053.7 14089.81 -7020.85
##
## Random effects:
## Formula: ~1 | ReefNameID
## (Intercept)
## StdDev: 1.102529
##
## Formula: ~1 | SiteName %in% ReefNameID
## (Intercept) Residual
## StdDev: 0.2113165 2.394941
##
## Fixed effects: Size ~ poly(days_sincespawning, 2, raw = FALSE)
## Value Std.Error DF t-value
## (Intercept) 10.12738 0.189346 2937 53.48619
## poly(days_sincespawning, 2, raw = FALSE)1 108.61003 4.246834 2937 25.57435
## poly(days_sincespawning, 2, raw = FALSE)2 -12.09116 4.055364 2937 -2.98152
## p-value
## (Intercept) 0.0000
## poly(days_sincespawning, 2, raw = FALSE)1 0.0000
## poly(days_sincespawning, 2, raw = FALSE)2 0.0029
## Correlation:
## (Intr) p(_,2,r=FALSE)1
## poly(days_sincespawning, 2, raw = FALSE)1 -0.161
## poly(days_sincespawning, 2, raw = FALSE)2 -0.239 -0.219
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.23467257 -0.65793662 0.01170185 0.63993843 5.49046070
##
## Number of Observations: 3042
## Number of Groups:
## ReefNameID SiteName %in% ReefNameID
## 60 103
# plot the growth trajectory
# get parameter estimates
CCA.lme.poly2.coefs <- lme(Size ~ poly(days_sincespawning,2,raw=TRUE),random=~1|ReefNameID/SiteName, data=CCA, method='REML')## RAW ist true, to extract parameter estimates
summary(CCA.lme.poly2.coefs)
## Linear mixed-effects model fit by REML
## Data: CCA
## AIC BIC logLik
## 14093.11 14129.22 -7040.553
##
## Random effects:
## Formula: ~1 | ReefNameID
## (Intercept)
## StdDev: 1.102529
##
## Formula: ~1 | SiteName %in% ReefNameID
## (Intercept) Residual
## StdDev: 0.2113165 2.39494
##
## Fixed effects: Size ~ poly(days_sincespawning, 2, raw = TRUE)
## Value Std.Error DF t-value
## (Intercept) -3.1019038 1.5277476 2937 -2.030377
## poly(days_sincespawning, 2, raw = TRUE)1 0.0825099 0.0134958 2937 6.113732
## poly(days_sincespawning, 2, raw = TRUE)2 -0.0000833 0.0000280 2937 -2.981523
## p-value
## (Intercept) 0.0424
## poly(days_sincespawning, 2, raw = TRUE)1 0.0000
## poly(days_sincespawning, 2, raw = TRUE)2 0.0029
## Correlation:
## (Intr) p(_,2,r=TRUE)1
## poly(days_sincespawning, 2, raw = TRUE)1 -0.985
## poly(days_sincespawning, 2, raw = TRUE)2 0.961 -0.992
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.23467268 -0.65793662 0.01170185 0.63993843 5.49046072
##
## Number of Observations: 3042
## Number of Groups:
## ReefNameID SiteName %in% ReefNameID
## 60 103
# construct data frame (prediction for two years)
newdata.poly2 = data.frame(days_sincespawning = seq(0,800),by=1,na.rm=TRUE)
Xmat <- model.matrix(~poly(days_sincespawning,2,raw=TRUE),data=newdata.poly2)
coefs <- fixef(CCA.lme.poly2.coefs)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(CCA.lme.poly2.coefs) %*% t(Xmat)))##
CCA.lme.poly2.coefs$fixDF
## $X
## (Intercept)
## 2937
## poly(days_sincespawning, 2, raw = TRUE)1
## 2937
## poly(days_sincespawning, 2, raw = TRUE)2
## 2937
##
## $terms
## (Intercept) poly(days_sincespawning, 2, raw = TRUE)
## 2937 2937
##
## attr(,"assign")
## attr(,"assign")$`(Intercept)`
## [1] 1
##
## attr(,"assign")$`poly(days_sincespawning, 2, raw = TRUE)`
## [1] 2 3
##
## attr(,"varFixFact")
## [,1] [,2] [,3]
## [1,] 0.1790947 0.000000000 0.00000e+00
## [2,] -0.3827722 0.001667012 0.00000e+00
## [3,] 1.4681360 -0.013392480 2.79521e-05
q = qt(0.975, df=2937)
CCA.newdata.poly2 = data.frame(newdata.poly2, fit=fit,
lower=fit - q*se,
upper=fit + q*se)
CCA.model=CCA.newdata.poly2 %>% mutate(Diet="CCA")
# construct data frame (data extent)
newdata.poly2.x = data.frame(days_sincespawning = seq(min(CCA$days_sincespawning), max(CCA$days_sincespawning),by=1,na.rm=TRUE))
Xmat <- model.matrix(~poly(days_sincespawning,2,raw=TRUE),data=newdata.poly2.x)
coefs <- fixef(CCA.lme.poly2.coefs)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(CCA.lme.poly2.coefs) %*% t(Xmat)))##
CCA.lme.poly2.coefs$fixDF
## $X
## (Intercept)
## 2937
## poly(days_sincespawning, 2, raw = TRUE)1
## 2937
## poly(days_sincespawning, 2, raw = TRUE)2
## 2937
##
## $terms
## (Intercept) poly(days_sincespawning, 2, raw = TRUE)
## 2937 2937
##
## attr(,"assign")
## attr(,"assign")$`(Intercept)`
## [1] 1
##
## attr(,"assign")$`poly(days_sincespawning, 2, raw = TRUE)`
## [1] 2 3
##
## attr(,"varFixFact")
## [,1] [,2] [,3]
## [1,] 0.1790947 0.000000000 0.00000e+00
## [2,] -0.3827722 0.001667012 0.00000e+00
## [3,] 1.4681360 -0.013392480 2.79521e-05
q = qt(0.975, df=2937)
CCA.newdata.poly2.x = data.frame(newdata.poly2.x, fit=fit,
lower=fit - q*se,
upper=fit + q*se)
CCA.model.x=CCA.newdata.poly2.x %>% mutate(Diet="CCA")
Model growth trajectory for corallivorous 0+ starfish
# run model
coral.glmm.gamma=glmmPQL(Size ~ days_sincespawning,random=~1|ReefNameID/SiteName, data=coral,family=Gamma(link="log"))####log link for left skewed distribution of residuals
summary(coral.glmm.gamma)
## Linear mixed-effects model fit by maximum likelihood
## Data: coral
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | ReefNameID
## (Intercept)
## StdDev: 0.1103037
##
## Formula: ~1 | SiteName %in% ReefNameID
## (Intercept) Residual
## StdDev: 1.778661e-05 0.2735259
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: Size ~ days_sincespawning
## Value Std.Error DF t-value p-value
## (Intercept) 1.7074510 0.10749234 441 15.88440 0
## days_sincespawning 0.0045186 0.00032311 441 13.98501 0
## Correlation:
## (Intr)
## days_sincespawning -0.969
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.61983338 -0.64985848 -0.09435956 0.61951380 4.11894907
##
## Number of Observations: 490
## Number of Groups:
## ReefNameID SiteName %in% ReefNameID
## 32 48
# plot growth trajectory
# construct data frame (prediction for two years)
newdata.coral = data.frame(days_sincespawning = seq(180,800),by=1,na.rm=TRUE)
Xmat <- model.matrix(~days_sincespawning,data=newdata.coral)
coefs <- fixef(coral.glmm.gamma)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(coral.glmm.gamma) %*% t(Xmat)))
coral.glmm.gamma$fixDF
## $X
## (Intercept) days_sincespawning
## 441 441
##
## $terms
## (Intercept) days_sincespawning
## 441 441
##
## attr(,"assign")
## attr(,"assign")$`(Intercept)`
## [1] 1
##
## attr(,"assign")$days_sincespawning
## [1] 2
##
## attr(,"varFixFact")
## [,1] [,2]
## [1,] -0.02659906 0.0000000000
## [2,] -0.10392272 0.0003224458
q = qt(0.975, df=441)
newdata.coral = data.frame(newdata.coral, fit=coral.glmm.gamma$family$linkinv(fit),
lower=coral.glmm.gamma$family$linkinv(fit - q*se),
upper=coral.glmm.gamma$family$linkinv(fit + q*se))
Coral.model=newdata.coral %>% mutate(Diet="coral")
# construct data frame (data extent)
newdata.coral.x = data.frame(days_sincespawning = seq(min(coral$days_sincespawning), max(coral$days_sincespawning),by=1,na.rm=TRUE))
Xmat <- model.matrix(~days_sincespawning,data=newdata.coral.x)
coefs <- fixef(coral.glmm.gamma)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(coral.glmm.gamma) %*% t(Xmat)))
coral.glmm.gamma$fixDF
## $X
## (Intercept) days_sincespawning
## 441 441
##
## $terms
## (Intercept) days_sincespawning
## 441 441
##
## attr(,"assign")
## attr(,"assign")$`(Intercept)`
## [1] 1
##
## attr(,"assign")$days_sincespawning
## [1] 2
##
## attr(,"varFixFact")
## [,1] [,2]
## [1,] -0.02659906 0.0000000000
## [2,] -0.10392272 0.0003224458
q = qt(0.975, df=441)
newdata.coral.x = data.frame(newdata.coral.x, fit=coral.glmm.gamma$family$linkinv(fit),
lower=coral.glmm.gamma$family$linkinv(fit - q*se),
upper=coral.glmm.gamma$family$linkinv(fit + q*se))
Coral.model.x=newdata.coral.x %>% mutate(Diet="coral")
Combine growth trajectories of herbivorous and corallivorous 0+ starfish in one dataframe
# prediction for two years
Models.combined=dplyr::bind_rows(CCA.model,Coral.model)
# data extent only
Models.combined.x=dplyr::bind_rows(CCA.model.x,Coral.model.x)
# export dataframe
write.csv(Models.combined,"PredictionModel_days since spawning.csv")
PredModel=Models.combined
Models.combined$Date=as.numeric(Models.combined$days_sincespawning)
Models.combined.x$Date=as.numeric(Models.combined.x$days_sincespawning)
Plot growth trajectories (Figure 2a)
# export as 4000/4000
ggplot(data=Models.combined.x,aes(x=Date, y=fit,fill=Diet)) +
geom_vline(xintercept=540,size=0.5,color="black",alpha=0.7,linetype="solid") +
geom_vline(xintercept=630,size=0.5,color="black",alpha=0.7,linetype="solid") +
geom_vline(xintercept=0,size=0.5,color="black",alpha=0.7,linetype="solid")+
geom_vline(xintercept=720,size=0.5,color="black",alpha=0.7,linetype="solid") +
geom_vline(xintercept=7,size=0.5,color="black",alpha=0.7,linetype="solid") +
geom_vline(xintercept=65,size=0.5,color="black",alpha=0.7,linetype="solid") +
geom_vline(xintercept=180,size=0.5,color="black",alpha=0.7,linetype="solid") +
annotate("rect", xmin=7,xmax=65,ymin=0.1,ymax=200,fill="#CCCCCC", alpha=0.3) +
annotate("rect", xmin=540,xmax=630,ymin=0.1,ymax=200,fill="#CCCCCC", alpha=0.3) +
annotate("rect", xmin=0,xmax=720,ymin=8,ymax=10,fill="red", alpha=0.1) +
annotate("rect", xmin=0,xmax=720,ymin=100,ymax=120,fill="red", alpha=0.1) +
geom_hline(yintercept=8,size=0.5,color="red",alpha=0.7,linetype="solid") +
geom_hline(yintercept=10,size=0.5,color="red",alpha=0.7,linetype="solid") +
geom_hline(yintercept=100,size=0.5,color="red",alpha=0.7,linetype="solid") +
geom_hline(yintercept=120,size=0.5,color="red",alpha=0.7,linetype="solid") +
geom_ribbon(aes(ymin=lower, ymax=upper), alpha=0.6) +
geom_ribbon(data=Models.combined,aes(ymin=lower, ymax=upper), alpha=0.2) +
geom_line(data=Models.combined, aes(y=lower,x=Date),color="black",linetype="dotted",size=1) +
geom_line(data=Models.combined, aes(y=fit,x=Date),color="black",linetype="dotted",size=1) +
geom_line(data=Models.combined, aes(y=upper,x=Date),color="black", linetype="dotted",size=1) +
geom_line(aes(color=Diet),size=2) +
geom_line(aes(y=lower,x=Date,color=Diet),linetype=2,size=1) +
geom_line(aes(y=upper,x=Date,color=Diet),linetype=2,size=1) +
scale_colour_manual(name="Diet",values=c(CCA="purple",coral="darkorange"),labels= c("Coralline algae","Scleractinian coral")) + scale_fill_manual(name="Diet",values=c(CCA="purple",coral="gold"),labels= c("Coralline algae","Scleractinian coral")) +
scale_size_manual(name="Diet",values=c("CCA"=7,"coral"=15),
labels=c("Coralline Algae","Scleractinian Coral")) +
guides(colour=guide_legend(),fill=guide_legend()) + theme(panel.background=element_rect(fill='white')) +
labs(x="Age in months since spawning", y="Size in diameter (mm)") +
theme(axis.text.x=element_text(size=110,angle=0,margin=unit(c(2,0,0,0),"cm"),vjust=0.5,hjust=0.5),
axis.text.y=element_text(size=110),
axis.title.x=element_text(size=120,face="plain",margin=unit(c(2,0,0,0),"cm")),
axis.title.y=element_text(size=120,face="plain",margin=unit(c(0,2,0,0),"cm")),
axis.line.x=element_line(colour = "black",size=1),
axis.ticks.x=element_line(colour = "black",size=1),
axis.line.y=element_line(colour = "black",size=1),
axis.ticks.y=element_line(colour = "black", size=1),
axis.ticks.length=unit(1,"cm"),
plot.margin=unit(c(5,5,2,1),"cm"),
legend.key.size=unit(4,"cm"),
legend.key = element_rect(fill="NA"),
legend.title=element_text(face="plain",size=100,margin=unit(c(0,0,2,0),"cm")),
legend.text=element_text(size=55),
legend.position="none",
legend.direction="vertical",
panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA)) +
scale_y_continuous(limits=c(0.3,200),oob=rescale_none,expand = c(0,0),breaks=c(0.3,50,100,150,200)) +
scale_x_continuous(limits=c(0,720), labels =c(0,6,12,18,24) ,breaks = c(0,180,360,540,720), expand=c(0.3,0))
## Warning: Removed 160 row(s) containing missing values (geom_path).
## Warning: Removed 160 row(s) containing missing values (geom_path).
## Warning: Removed 160 row(s) containing missing values (geom_path).
Plot variability in the timing of the shift (Figure 3)
# 4000/3500
ggplot(data=CCA_proportion,aes(x=days_sincespawning, y=fit)) +
geom_ribbon(aes(ymin=lower,ymax=upper),fill="#CCCCCC") +
geom_ribbon(data=CCA_proportion.x, aes(ymin=lower,ymax=upper), fill="#999999") +
geom_line(aes(y=lower,x=days_sincespawning),color="black",linetype="dotted",size=1) +
geom_line(aes(y=fit,x=days_sincespawning),color="black",linetype="dotted",size=1) +
geom_line(aes(y=upper,x=days_sincespawning),color="black", linetype="dotted",size=1) +
geom_line(data=CCA_proportion.x,color="black",size=2) +
geom_line(data=CCA_proportion.x,aes(y=lower,x=days_sincespawning),color="black",linetype=2,size=1) +
geom_line(data=CCA_proportion.x,aes(y=upper,x=days_sincespawning),color="black",linetype=2,size=1) +
geom_point(data=RAW_Timing,aes(x=days_sincespawning, y=PropCCA), color="black",fill="grey",alpha=0.5,size=12) +
theme(panel.background=element_rect(fill='white')) +
labs(x="Age in months since spawning", y="Proportion of juvenile starfish\n feeding on coralline algae") +
theme(axis.text.x=element_text(size=110,angle=0,margin=unit(c(2,0,0,0),"cm"),vjust=0.5,hjust=0.5),
axis.text.y=element_text(size=110),
axis.title.x=element_text(size=130,face="plain",margin=unit(c(2,0,0,0),"cm")),
axis.title.y=element_text(size=130,face="plain",margin=unit(c(0,2,0,0),"cm")),
axis.line.x=element_line(colour = "black",size=1),
axis.ticks.x=element_line(colour = "black",size=1),
axis.line.y=element_line(colour = "black",size=1),
axis.ticks.y=element_line(colour = "black", size=1),
axis.ticks.length=unit(1,"cm"),
plot.margin=unit(c(5,5,2,1),"cm"),
legend.key.size=unit(4,"cm"),
legend.key = element_rect(fill="NA"),
legend.title=element_text(face="plain",size=55,margin=unit(c(0,0,2,0),"cm")),
legend.text=element_text(size=55),
legend.position="none",
legend.direction="vertical",
panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA)) +
scale_y_continuous(limits=c(0,1),expand = c(0.01,0),breaks=c(0,0.25,0.50,0.75,1)) + geom_hline(yintercept=c(0,0.25,0.50,0.75,1),size=0.5,color="black",alpha=1,linetype="dotted") +
scale_x_continuous(limits=c(0,720), labels =c(0,6,12,18,24) ,breaks =c(0,180,360,540,720), expand=c(0,0)) +
geom_hline(yintercept=0.50,size=10,color="red",alpha=0.3,linetype="solid")
Plot piecharts (Figure 1b)
RAWbox=RAW %>% mutate(MonthClass=cut(RAW$days_sincespawning,breaks=c(0,30,60,90,120,150,180,210,240,270,300,330,360,390),labels=c("1","2","3","4","5","6","7","8","9","10","11","12","13")))
Pie=RAWbox %>% group_by(MonthClass,Diet) %>% summarise_at("Size", funs(n())) %>% spread(Diet,Size,sep="_")
# replace NA in CCA_n and coral_n with 0
Pie$Diet_CCA[is.na(Pie$Diet_CCA)] <- 0
Pie$Diet_coral[is.na(Pie$Diet_coral)] <- 0
TEST=Pie %>% mutate(Total= Diet_CCA + Diet_coral,PropCCA=Diet_CCA/Total,Propcoral=1-PropCCA)
TEST=TEST %>% gather(Diet,prop,5:6)
TEST$Diet=recode_factor(TEST$Diet,`PropCCA`="CCA",`Propcoral`="coral")
Pie.test=Pie %>% gather(Diet,n,2:3)
Pie.test$Diet=recode_factor(Pie.test$Diet,`Diet_CCA`="CCA",`Diet_coral`="coral")
TEST=dplyr::select(TEST,MonthClass,Diet,prop)
TEST=TEST %>% mutate(perc="%",percentage=round(prop*100,digits=1)) %>% unite(Percentage,percentage,perc,sep="")
Pie_text=TEST %>% filter(`Diet`=="CCA") %>% mutate(y=54)
# create pie chart
# 3000/100
pies2=ggplot(data=TEST, aes(x=factor(1),y=prop,fill=Diet)) + geom_bar(stat="identity", size = 1,aes=0.5) + coord_polar(theta = "y") + facet_wrap(~ MonthClass,ncol=13) + scale_fill_manual(values = c("purple", "gold")) + scale_colour_manual(values = c("purple", "orange")) + theme_void() + theme(legend.position = "none")
pies2
Plot boxplot (Figure 1b)
give.n <- function(x){return(c(y = mean(x), label = length(x)))}
N=RAWbox %>% group_by(MonthClass) %>% summarise(n=n()) %>% as.data.frame()
N=RAWbox %>% group_by(MonthClass) %>% summarise(n=n()) %>% mutate(text=c("n=19","n=1,193","n=843",n="n=168",n="n=728",n="n=73",n="n=152",n="n=152",n="n=204"), y=66)
# plot boxplot
ggplot(data=RAWbox, aes(x=MonthClass, y=Size)) +
geom_jitter(aes(shape=Diet,colour=Diet),size=6,alpha=0.4,position = position_jitterdodge()) +
geom_boxplot(size=1.5, aes(fill=Diet,colour=Diet),outlier.colour="black",outlier.shape = 24,outlier.size = 8,alpha=0.5) +
geom_text(data=N, aes(x=MonthClass,y=y),label=N$text,size=20,fontface="plain") +
geom_text(data=Pie_text, aes(x=MonthClass,y=y),label=Pie_text$Percentage,size=20,fontface="plain",color="purple") +
annotate("rect",xmin=0,xmax=Inf,ymin=8,ymax=10,fill="red", alpha=0.1) +
geom_hline(yintercept=8,size=0.5,color="red",alpha=0.7,linetype="solid") +
geom_hline(yintercept=10,size=0.5,color="red",alpha=0.7,linetype="solid") +
labs(x="Time in months since spawning", y="Size in diameter (mm)") +
theme(axis.text=element_text(size=80),
axis.title.x=element_text(size=90,face="plain",margin=unit(c(3,0,2,0),"cm")),
axis.title.y=element_text(size=90,face="plain",margin=unit(c(0,2,0,0),"cm")),
axis.line.x=element_line(colour = "NA"),
axis.ticks.x=element_line(colour = "NA"),
axis.line.y=element_line(colour = "black",size=1),
axis.ticks.y=element_line(colour = "black", size=1),
axis.ticks.length=unit(1,"cm"),
plot.margin=unit(c(3,1,1,1),"cm"),
legend.key.size=unit(2,"cm"),
legend.key = element_rect(fill="NA"),
legend.title=element_text(face="plain",size=90),
legend.text=element_text(size=80),
legend.margin = margin(2,0,1,0,"cm"),
legend.position="bottom",
panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA)) +
scale_y_continuous(limits=c(0,70),breaks = c(0,10,20,30,40,50,60,70),expand=c(0,0)) +
geom_hline(yintercept=c(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70),size=0.7,color="black",linetype="dashed") +
scale_color_manual(name="Diet", values=c("CCA"="purple","coral"="orange"),
labels=c("Coralline algae","Scleractinian coral")) +
scale_fill_manual(name="Diet",values=c("CCA"="purple","coral"="gold"),
labels=c("Coralline algae","Scleractinian coral")) +
scale_shape_manual(name="Diet",values=c("CCA"=20,"coral"=20),
labels=c("Coralline algae","Scleractinian coral")) +
guides(colour = guide_legend())
Frequency analysis
##create size categories
Size_Classes1=RAW %>% mutate(SizeClass1=cut(Size,breaks=c(0,8,20,70),labels=c("≤ 8 mm","> 8-20 mm","> 20 mm")))
###make categories categorical variables
Size_Classes1$fSizeClass1=factor(Size_Classes1$SizeClass1)
###contingency table
tbl1 = table(Size_Classes1$fSizeClass1, Size_Classes1$Diet)
##Chi-square test
##tbl1
chisq.tbl1=chisq.test(tbl1)
chisq.tbl1
##
## Pearson's Chi-squared test
##
## data: tbl1
## X-squared = 2102.3, df = 2, p-value < 2.2e-16
chisq.tbl1$observed
##
## CCA coral
## ≤ 8 mm 920 0
## > 8-20 mm 2106 168
## > 20 mm 16 322
round(chisq.tbl1$expected,3)
##
## CCA coral
## ≤ 8 mm 792.367 127.633
## > 8-20 mm 1958.524 315.476
## > 20 mm 291.109 46.891
Model for time at which 50% of starfish reach 8mm size threshold
# need to use the RAW numbers for the model
prop2 <- cbind(RAW_TIME$`≤ 8 mm_Size`,RAW_TIME$`> 8-70 mm_Size`)
# run model
M2 <- glmmPQL(prop2~Day.of.Year.R, data=RAW_TIME,random = ~1 | ReefNameID/SiteName, family = binomial(link="logit"))
summary(M2)
## Linear mixed-effects model fit by maximum likelihood
## Data: RAW_TIME
## AIC BIC logLik
## NA NA NA
##
## Random effects:
## Formula: ~1 | ReefNameID
## (Intercept)
## StdDev: 0.3792412
##
## Formula: ~1 | SiteName %in% ReefNameID
## (Intercept) Residual
## StdDev: 0.05662259 1.180663
##
## Variance function:
## Structure: fixed weights
## Formula: ~invwt
## Fixed effects: prop2 ~ Day.of.Year.R
## Value Std.Error DF t-value p-value
## (Intercept) 4.366940 0.3812083 63 11.45552 0
## Day.of.Year.R -0.030201 0.0022365 19 -13.50347 0
## Correlation:
## (Intr)
## Day.of.Year.R -0.967
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.7050605 -0.4319825 -0.1269110 0.4035408 5.6169134
##
## Number of Observations: 131
## Number of Groups:
## ReefNameID SiteName %in% ReefNameID
## 64 111
plot(M2)
###date of shift
coefs=fixef(M2)
-coefs[1] /coefs[2]
## (Intercept)
## 144.5954
# plot model
# create data frame (prediction for two years)
newdata <- data.frame(Day.of.Year.R = seq(from = 0, to = 720, by = 1), upperM2=NA, lowerM2=NA, bestM2=NA)
Xmat <- model.matrix(~Day.of.Year.R, data=newdata)
coefs <- fixef(M2)
fit = as.vector(coefs %*% t(Xmat))
se = sqrt(diag(Xmat %*% vcov(M2) %*% t(Xmat)))
q = qt(0.975, df=M2$fixDF$X)
newdata = data.frame(newdata, fit=M2$family$linkinv(fit),
lower=M2$family$linkinv(fit - q*se),
upper=M2$family$linkinv(fit + q*se))
TH_proportion=newdata
ggplot(TH_proportion, aes(y=fit, x=Day.of.Year.R)) +
geom_ribbon(aes(ymin=lower, ymax=upper), fill='blue', alpha=0.3)+
geom_line()
