###Load data and libraries####
setwd("H:/javiera/CuSP - 15345/DV Growth rates")
suppressPackageStartupMessages(library(tidyverse))
## Warning: package 'tidyverse' was built under R version 3.3.3
## Warning: package 'purrr' was built under R version 3.3.3
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(gamlss))
suppressPackageStartupMessages(library(scales))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(forcats))
library(knitr)
source('theme_javier.R')
source('panel_cor.R')
source('corvif.R')
theme_set(theme_javier())
##varaiables house keeping#####
ruakaka <-
read.csv("ruakaka.csv") %>%
mutate(
dv = DID.VEX,
fd_season = fct_recode(
factor(D_season),
Summer = '1',
Autumn = '2' ,
Winter = '3',
Spring = '4'
),
S_date = as.Date(S_date, format = '%d/%m/%Y'),
id = as.factor(paste(Array, Plate, sep = '-')),
Others = ALGAE + MUSSEL + OYSTER + BRY.ENC + BRY.ERECT + HYDROID + UNKNOWN +
AMPH + ANEMONE + BARN + EGG + SPONGE + TUBE + COLONIAL + SOLITARY + DOTS,
logSpace = log(SPACE + 1))
###Plots though time####
ruakaka %>%
gather(key, value, c(dv,Others,SPACE)) %>%
ggplot(.,aes(y = value, x = S_date, lty = fd_season)) +
stat_summary(fun.y = mean, geom = 'line') +
stat_summary(fun.data = mean_se,
geom = "ribbon",
alpha = 0.07) +
ggtitle('Ruakaka') + xlab('') +
ylab('') +
theme_javier() +
theme(legend.position = c(.8,.2)) +
scale_linetype_manual(values = c(1, 2, 3, 4), name = 'Deployment season') +
scale_x_date(
breaks = date_breaks("3 month"),
labels = date_format("%b %y") ,
limits = as.Date(c('2008-01-01', '2009-9-01'))
) +
coord_cartesian(ylim = c(0, 100)) +
facet_wrap(~key,nrow = 3)
###data exploration#####
dotchart(ruakaka$dv, color = ruakaka$fd_season,
xlab = 'DV % cover')## lots of zeros
sum(ruakaka$dv == 0) / nrow(ruakaka)###30% of zeros
## [1] 0.3057325
##Beta regression#####
ruakaka %>%
dplyr::select (S_date, fd_season, dv, logSpace, Others,COLONIAL, TempChange) %>%
gather(key, value, dv:TempChange) %>%
ggplot(.,aes(x= value)) +
geom_histogram(position = 'identity',color = 1,fill= 'gray50',bins = 15)+
facet_wrap(~key,scales= 'free')
##data exploration
##check response variable distribution
N = nrow(ruakaka)
ruakaka$prop.dv = ((ruakaka$dv / 100) * (N - 1) + 0.5) / N
histDist(ruakaka$prop.dv, family = BE())
##
## Family: c("BE", "Beta")
## Fitting method: "nlminb"
##
## Call: gamlssML(y = y, family = "BE", formula = ruakaka$prop.dv)
##
## Mu Coefficients:
## [1] -2.111
## Sigma Coefficients:
## [1] -0.1846
##
## Degrees of Freedom for the fit: 2 Residual Deg. of Freedom 469
## Global Deviance: -1449.29
## AIC: -1445.29
## SBC: -1436.99
##check collinearity
pairs(
~ Time + S_season + dv + S_month + TempChange + DailyRecruit +
COLONIAL + SPACE,
data = ruakaka,
panel = panel.smooth,
upper.panel = NULL
)
corvif(ruakaka[, c('Time',
'S_season',
'S_month',
'TempChange',
'COLONIAL',
'SPACE',
'Others')])
##
##
## Variance inflation factors
##
## GVIF
## Time 1.679283
## S_season 1.880262
## S_month 1.912587
## TempChange 1.422476
## COLONIAL 1.453149
## SPACE 1.466589
## Others 1.652122
##Full model
m0 = gamlss(prop.dv ~
cs(Time) +
cs(S_month) +
fd_season * TempChange +
re(random = ~ 1 | id),
family = BE,
data = ruakaka
)
## GAMLSS-RS iteration 1: Global Deviance = -1706.677
## GAMLSS-RS iteration 2: Global Deviance = -1784.171
## GAMLSS-RS iteration 3: Global Deviance = -1823.359
## GAMLSS-RS iteration 4: Global Deviance = -1840.838
## GAMLSS-RS iteration 5: Global Deviance = -1848.27
## GAMLSS-RS iteration 6: Global Deviance = -1851.4
## GAMLSS-RS iteration 7: Global Deviance = -1852.776
## GAMLSS-RS iteration 8: Global Deviance = -1853.422
## GAMLSS-RS iteration 9: Global Deviance = -1853.739
## GAMLSS-RS iteration 10: Global Deviance = -1853.906
## GAMLSS-RS iteration 11: Global Deviance = -1853.981
## GAMLSS-RS iteration 12: Global Deviance = -1854.021
## GAMLSS-RS iteration 13: Global Deviance = -1854.04
## GAMLSS-RS iteration 14: Global Deviance = -1854.05
## GAMLSS-RS iteration 15: Global Deviance = -1854.055
## GAMLSS-RS iteration 16: Global Deviance = -1854.055
plot(m0)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.04691656
## variance = 0.926933
## coef. of skewness = -0.2162038
## coef. of kurtosis = 3.522269
## Filliben correlation coefficient = 0.9963423
## ******************************************************************
##Backward model selection based on GAIC
m1 <- stepGAIC(m0, direction = 'backward',trace = 0)
## Start: AIC= -1769.1
## prop.dv ~ cs(Time) + cs(S_month) + fd_season * TempChange + re(random = ~1 |
## id)
summary(m1)
## ******************************************************************
## Family: c("BE", "Beta")
##
## Call: gamlss(formula = prop.dv ~ cs(Time) + cs(S_month) +
## fd_season * TempChange + re(random = ~1 | id),
## family = BE, data = ruakaka)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.057963 0.163375 -18.718 < 2e-16 ***
## cs(Time) 0.136555 0.009582 14.251 < 2e-16 ***
## cs(S_month) -0.094570 0.012839 -7.366 9.11e-13 ***
## fd_seasonAutumn 0.050461 0.097092 0.520 0.604
## fd_seasonWinter -0.075810 0.110942 -0.683 0.495
## fd_seasonSpring 0.007834 0.144520 0.054 0.957
## TempChange 0.280920 0.052540 5.347 1.46e-07 ***
## fd_seasonAutumn:TempChange -0.365153 0.074598 -4.895 1.40e-06 ***
## fd_seasonWinter:TempChange -0.555630 0.079985 -6.947 1.40e-11 ***
## fd_seasonSpring:TempChange -0.538336 0.098511 -5.465 7.88e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: logit
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.88089 0.04694 -18.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 471
## Degrees of Freedom for the fit: 42.47756
## Residual Deg. of Freedom: 428.5224
## at cycle: 16
##
## Global Deviance: -1854.055
## AIC: -1769.1
## SBC: -1592.612
## ******************************************************************
coef_ruakaka <- exp(cbind(Estimates = coef(m1), confint(m1)))
## Warning in vcov.gamlss(object, robust = robust): Additive terms exists in the mu formula.
## Standard errors for the linear terms maybe are not appropriate
kable(coef_ruakaka, digits = 2)
| (Intercept) |
0.05 |
0.03 |
0.06 |
| cs(Time) |
1.15 |
1.12 |
1.17 |
| cs(S_month) |
0.91 |
0.89 |
0.93 |
| fd_seasonAutumn |
1.05 |
0.87 |
1.27 |
| fd_seasonWinter |
0.93 |
0.75 |
1.15 |
| fd_seasonSpring |
1.01 |
0.76 |
1.34 |
| TempChange |
1.32 |
1.19 |
1.47 |
| re(random = ~1 | id) |
NA |
NA |
NA |
| fd_seasonAutumn:TempChange |
0.69 |
0.60 |
0.80 |
| fd_seasonWinter:TempChange |
0.57 |
0.49 |
0.67 |
| fd_seasonSpring:TempChange |
0.58 |
0.48 |
0.71 |
##Model validation by inspecting normalised quantile residuals
plot(m1)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.04691656
## variance = 0.926933
## coef. of skewness = -0.2162038
## coef. of kurtosis = 3.522269
## Filliben correlation coefficient = 0.9963423
## ******************************************************************
plot(m1, ts = T)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.04691656
## variance = 0.926933
## coef. of skewness = -0.2162038
## coef. of kurtosis = 3.522269
## Filliben correlation coefficient = 0.9963423
## ******************************************************************
#print the summary table
Rsq(m1, type = "Cox Snell")
## [1] 0.5765687
##Plot the partial effects of the most parsimonious model####
par(
mfrow = c(2, 3),
oma = c(0, 0, 0, 0),
mar = c(4, 4, 2, 0.5),
mgp = c(2.2, 1, 0)
)
term.plot(
m1,
se = T,
partial = T,
col.res = 'black',
ylim = 'common',
scheme = 'shade',
col.shaded = 'gray',
xlabs = c(
'Deploymenet duration (months)',
'Sampling month',
'Deployment season',
'Temperature change'),
lwd.term = 1,
pch.res = 1,
cex.res = 1,
ylab = rep('Partial effect', 6)
)
## Warning in term.plot(m1, se = T, partial = T, col.res = "black", ylim =
## "common", : interactions have been taken out from the plots
###Multivariate###
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-2
##redundancy analyses
sp <- sqrt(ruakaka[, 15:33])
cap <-
capscale(sp ~ fd_season * as.numeric(S_date) + Condition(as.numeric(S_date)) + TempChange, data = ruakaka, distance = 'bray')
# sum_cap <- summary(cap)
# anova_cap <- anova(cap)
#
# anova_cap_axis <- anova(cap, by = "axis", perm.max = 999)
# anova_cap_terms <- anova(cap, by = "terms", permu = 999)
##get the sites scores and glue them to the data
sco = scores(cap)
ruakaka$rda1 = sco$sites[, 1]
ruakaka$rda2 = sco$sites[, 2]
##get the species scoores
sp.sco1 = as.data.frame(sco$species)
mySp = sp.sco1[which(abs(sp.sco1$CAP1) > 0.3),]
##plot
library(ggplot2)
##All together for the capscale plot
p <-
ggplot(ruakaka, aes(x = S_date, y = rda1)) +
geom_smooth(aes(color = factor(fd_season)),
size = 1,
method = loess,
se = F) +
scale_colour_brewer(name = 'Deployment season', palette = 'Paired')+
scale_x_date(
breaks = date_breaks("3 month"),
labels = date_format("%b %y") ,
limits = as.Date(c('2008-01-01', '2009-9-01'))
)
ruakaka_prc <-
p + geom_rug(data = mySp, aes(x = as.Date('2009-8-01'), y = CAP1),
sides = 'right') +
geom_text(
data = mySp,
aes(x = as.Date('2009-8-01'), y = CAP1),
size = 2.5,
label = row.names(mySp),
hjust = 0,
vjust = 0,
size = 4
) + ylab('Coefficients') +
theme_javier() +
xlab('')+
ggtitle('Ruakaka')
## Warning: The plyr::rename operation has created duplicates for the
## following name(s): (`size`)
ruakaka_prc
