# Attach data file and check variable names

attach(S3_Recruit_1997.2008_R_ready)
names(S3_Recruit_1997.2008_R_ready)

# Delete missing observations and create two files
# Xvar contains Year, Bay, Site, and Size
# Yvar contains standardized recruitment data for barnacles, mussels and fucoids

Xvar <- subset(S3_Recruit_1997.2008_R_ready, Barnacle >=0 & Fucoid >=0 & Mussel >= 0, 
               select = Year:Size)
Yvar <- subset(S3_Recruit_1997.2008_R_ready, Barnacle >=0 & Fucoid >=0 & Mussel >= 0, 
               select = logF_std:logM_std)

# Open vegan

library(vegan)

# Run CAP and permunation test

capmodel1 <- capscale(Yvar ~ Size*Year*Bay + Condition(Site), Xvar)
capmodel1
anova(capmodel1, by="terms")

# Use envfit to examine correlations of recruitment to CAP axes

ef <- envfit(capmodel1, Yvar)
ef

# Re-run model dropping Site as a random factor
# This is done in order to visualize separation of Bays

capmodel2 <- capscale(Yvar ~ Size*Year*Bay, Xvar)

# Plot Centroids for Size with Barnacle contours

plot(capmodel2,display ="sites", type="n", xlim=c(-0.31,0.31), ylim=c(-0.31,0.31))
with(Xvar, ordiellipse(capmodel2, Size, kind="se", label=TRUE,conf=0.95))
tmp <- with(Yvar, ordisurf(capmodel2, logB_std, add = TRUE, nlevels=20))

# with Mussel contours

plot(capmodel2,display ="sites", type="n", xlim=c(-0.31,0.31), ylim=c(-0.31,0.31))
with(Xvar, ordiellipse(capmodel2, Size, kind="se", label=TRUE,conf=0.95))
tmp <- with(Yvar, ordisurf(capmodel2, logM_std, add = TRUE, nlevels=30))

# with Fucoid contours

plot(capmodel2,display ="sites", type="n", xlim=c(-0.31,0.31), ylim=c(-0.31,0.31))
with(Xvar, ordiellipse(capmodel2, Size, kind="se", label=TRUE,conf=0.95))
tmp <- with(Yvar, ordisurf(capmodel2, logF_std, add = TRUE, nlevels=20))

# To make similar plots for Bay or Year, replace Size with either Bay or Year in last 9 lines