#################################
# This is the code to make figure and
#    to conduct analysis of
#    Otsu et al. for PeerJ
#################################

# Set working directory
# (this part should be changed for your computer directory)
setwd("/Users/hayatoiijima/Dropbox/researchfiles/Kushigata_fence/PeerJ/data/")

# Definition of rm.level() function
rm.level <- function(x){
  sub.fun <- function(sub.df){
    if (is.factor(sub.df)) factor(sub.df)
    else sub.df
  }
  data.frame(lapply(x, sub.fun))
}

#########################
# Description of data
#########################
# Load data
d <- read.csv("matrix2_HI.csv")
e <- read.csv("splist_HI.csv", fileEncoding="Shift-JIS")
colnames(d)[-(1:3)] <- as.character(e$SpJ)
# Merge synonym of some species
d2 <- d[, -c(grep("ウマノアシガタ", colnames(d)), grep("キンポウゲ", colnames(d)))]
temp <- d[, c(grep("ウマノアシガタ", colnames(d)), grep("キンポウゲ", colnames(d)))]
temp$use <- apply(temp, 1, max)
d2$キンポウゲ <- temp$use

# rename species
e$Sp <- paste(substr(e$Sp1, 1, 1), ". ", as.character(e$Sp2), sep="")
e2 <- e[-c(grep("ウマノアシガタ", e$SpJ), grep("キンポウゲ", e$SpJ)), ]
e2 <- rbind(e2, e[grep("キンポウゲ", e$SpJ)[2], ])
colnames(d2)[-(1:3)] <- e2$Sp

d3 <- d2[, c(1:3, order(colnames(d2)[-(1:3)])+3)]
# Remove a quadrat without plant cover
d4 <- rm.level(d3[-55, ])
colnames(d4)[-(1:3)] <- gsub("\\..", ".", colnames(d4)[-(1:3)])
spmat <- d4[, -c(1:3)]
# Convert blaun-branquet value to percentage
spmat_100 <- spmat
spmat_100 <- apply(spmat_100, 2, function(x)
                   { ifelse(x == 5, 87.5, x) }
)
spmat_100 <- apply(spmat_100, 2, function(x)
                   { ifelse(x == 4, 62.5, x) }
)
spmat_100 <- apply(spmat_100, 2, function(x)
                   { ifelse(x == 3, 37.5, x) }
)
spmat_100 <- apply(spmat_100, 2, function(x)
                   { ifelse(x == 2, 15.0, x) }
)
spmat_100 <- apply(spmat_100, 2, function(x)
                   { ifelse(x == 1, 2.5, x) }
)
#spmat_100 <- round(spmat_100)
env <- d4[, 1:3]

########################
# Description of data
########################
## Table 1
d5 <- cbind(d4[, 1:3], spmat_100)
dl <- reshape(d5, varying=list(4:ncol(d5)),
              timevar="Species", idvar="plotid",
              v.names="C", times=colnames(d5[, -(1:3)]),
              direction="long"
)
dl$TreatYear <- paste(dl$Treat, dl$Year, sep="_")
spcover <- as.data.frame(tapply(dl$C, dl[, c("Species", "TreatYear")], mean))
library(openxlsx)
st2 <- read.xlsx("speciestype.xlsx")
st2$Sp <- gsub("\\. ", ".", as.character(st2$Sp))
spcover$Sp <- rownames(spcover)
spcover <- merge(spcover, st2[, c("Sp", "TAX")])
spcover$Sp <- gsub("\\.", ". ", as.character(spcover$Sp))
e3 <- e2[order(e2$Sp), ]
e3$Spall <- paste(e3$Sp1, e3$Sp2, e3$Sp3, e3$Sp4, sep=" ")
spcover$Sp <- e3$Spall
# Export as xlsx file
write.xlsx(spcover, file="spcover_table.xlsx")

### Fig. 1
dl2 <- dl
dl2 <- merge(dl2, st2[, c("Sp", "TAX")], by.x="Species", by.y="Sp")
dl3 <- dl2[dl2$C > 0, ]
heikin <- tapply(dl3$C, dl3[, c("TAX", "Year", "Treat")], mean)
# Graphics parameter
par(mfrow=c(2,2), mar=c(1,3,0,0), oma=c(4,4,1,1), ps=15)
# Fence 2010
matplot(t(heikin[,,1]), type="o", col="black", pch=0:4, lty=1, xaxt="n", cex=2,
        ylim=c(0, 60)
)
points(rep(1, 5), heikin[,1,3], pch=0:4, cex=2)
axis(1, at=1:6, labels=FALSE)
abline(v=1.5, lty=2)
legend("topright", legend="A", bty="n")
# Legends
plot.new()
lf <- c("Dicots", "Graminoids", "Monocots", "Ferns", "Woody plants")
legend("center", lty=rep(1, 5), pch=0:4, col="black",
       legend=lf
)
# Fence 2011
matplot(t(heikin[,,2]), type="o", col="black", pch=0:4, lty=1, xaxt="n", cex=2,
        ylim=c(0, 60)
)
points(rep(1, 5), heikin[,1,3], pch=0:4, cex=2)
axis(1, at=1:6, labels=c(1981, 2011:2015))
abline(v=1.5, lty=2)
legend("topright", legend="B", bty="n")
# Outside
matplot(t(heikin[,,4]), type="o", col="black", pch=0:4, lty=1, xaxt="n", cex=2,
        ylim=c(0, 60)
)
points(rep(1, 5), heikin[,1,3], pch=0:4, cex=2)
axis(1, at=1:6, labels=c(1981, 2011:2015))
abline(v=1.5, lty=2)
legend("topright", legend="C", bty="n")
mtext(1, line=2, outer=TRUE, text="Year")
mtext(2, line=1.5, outer=TRUE, text="Coverage (%)")
# Export as eps file
dev.copy2eps(file="fig1_new.eps")

# GLMM for the coverage
d6 <- cbind(d4[, 1:3], spmat)
dl4 <- reshape(d6, varying=list(4:ncol(d6)),
              timevar="Species", idvar="plotid",
              v.names="C", times=colnames(d6[, -(1:3)]),
              direction="long"
)
dl4$TreatYear <- paste(dl4$Treat, dl4$Year, sep="_")
dl4 <- merge(dl4, st2[, c("Sp", "TAX")], by.x="Species", by.y="Sp")
dl4$plotid <- paste(dl4$Treat, dl4$Plot, sep="_")
dl4$YAF <- 0
dl4$YAF <- ifelse(dl4$Treat=="fence2010", dl4$Year-2010, dl4$YAF)
dl4$YAF <- ifelse(dl4$Treat=="fence2011", dl4$Year-2011, dl4$YAF)
dl5 <- rm.level(dl4[dl4$Treat!="Old", ])
library(ordinal)
rescover2 <- clmm(as.factor(C) ~ YAF*TAX + Treat + (1|plotid) + (1|Year), dl5, na.action = "na.fail")
rescover2_1 <- clmm(as.factor(C) ~ YAF*TAX + (1|plotid) + (1|Year), dl5, na.action = "na.fail")
# Export the summary of GLMM
write.csv(summary(rescover2_1)$coefficients, file="CLMM_spcover.csv")

library(MuMIn)
dredge(rescover2, extra = alist(AIC))
dredge(rescover2_1, extra = alist(AIC))


## the number of emerged species
# total
sn <- as.data.frame(xtabs(C ~ Species + Treat + Year, dl4))
sn$Freq <- ifelse(sn$Freq==0, 0, 1)
write.csv(tapply(sn$Freq, sn[, c("Treat", "Year")], sum), "totalspn.csv")
# by each life form
sn2 <- merge(sn, st2[, c("Sp", "TAX")], by.x="Species", by.y="Sp")
write.csv(tapply(sn2$Freq, sn2[, c("TAX", "Treat", "Year")], sum), "sn_lifeform.csv")


########################
# Dissimilarity index
########################
library(vegan)
# NMDS
chmds <- metaMDS(spmat, dist="chao", zerodist="add")
ax12 <- as.data.frame(chmds$points)

### Fig. 2
# graphics parameters
par(mfrow=c(2,2), mar=c(1,1,0,0), oma=c(4,4,1,1), ps=15)
sym <- c(15, 16, 3, 17)
iro <- c("purple", "blue", "green", "yellow", "red")
# center of each year
xloc <- as.data.frame(ftable(tapply(ax12[, 1], d4[, c("Treat", "Year")], mean)))
yloc <- as.data.frame(ftable(tapply(ax12[, 2], d4[, c("Treat", "Year")], mean)))
loc <- cbind(xloc, yloc[, 3])
colnames(loc)[(ncol(loc)-1):ncol(loc)] <- c("x", "y")
# old vs fence 2010
plot(ax12, xlab="", ylab="", type="n", xlim=c(-3, 3), ylim=c(-3, 2), xaxt="n", yaxt="n")
axis(1, at=(-3):3, labels=FALSE)
axis(2, at=(-3):2)
points(ax12[d4$Treat=="Old", 1], ax12[d4$Treat=="Old", 2], pch=3, col="black")
points(ax12[d4$Treat=="fence2010", 1], ax12[d4$Treat=="fence2010", 2], pch=16, col=iro[d[d$Treat=="fence2010", "Year"]-2010])
lines(loc[loc$Treat=="fence2010", "x"][1:4], loc[loc$Treat=="fence2010", "y"][1:4], lwd=2)
arrows(loc[loc$Treat=="fence2010", "x"][4], loc[loc$Treat=="fence2010", "y"][4],
       loc[loc$Treat=="fence2010", "x"][5], loc[loc$Treat=="fence2010", "y"][5],
       angle=45, length=0.05, lwd=2
)
legend("bottomright", legend="A", cex=0.8, bty="n")
# old vs fence 2011
plot(ax12, xlab="", ylab="", type="n", xlim=c(-3, 3), ylim=c(-3, 2), xaxt="n", yaxt="n")
axis(1, at=(-3):3, labels=FALSE)
axis(2, at=(-3):2, labels=FALSE)
points(ax12[d4$Treat=="Old", 1], ax12[d4$Treat=="Old", 2], pch=3, col="black")
points(ax12[d4$Treat=="fence2011", 1], ax12[d4$Treat=="fence2011", 2], pch=15, col=iro[d[d$Treat=="fence2011", "Year"]-2010])
lines(loc[loc$Treat=="fence2011", "x"][1:4], loc[loc$Treat=="fence2011", "y"][1:4], lwd=2)
arrows(loc[loc$Treat=="fence2011", "x"][4], loc[loc$Treat=="fence2011", "y"][4],
       loc[loc$Treat=="fence2011", "x"][5], loc[loc$Treat=="fence2011", "y"][5],
       angle=45, length=0.05, lwd=2
)
legend("bottomright", legend="B", cex=0.8, bty="n")
# old vs out
plot(ax12, xlab="", ylab="", type="n", xlim=c(-3, 3), ylim=c(-3, 2), xaxt="n", yaxt="n")
axis(1, at=(-3):3)
axis(2, at=(-3):2)
points(ax12[d4$Treat=="Old", 1], ax12[d4$Treat=="Old", 2], pch=3, col="black")
points(ax12[d4$Treat=="out", 1], ax12[d4$Treat=="out", 2], pch=17, col=iro[d[d$Treat=="out", "Year"]-2010])
lines(loc[loc$Treat=="out", "x"][1:4], loc[loc$Treat=="out", "y"][1:4], lwd=2)
arrows(loc[loc$Treat=="out", "x"][4], loc[loc$Treat=="out", "y"][4],
       loc[loc$Treat=="out", "x"][5], loc[loc$Treat=="out", "y"][5],
       angle=45, length=0.05, lwd=2
)
legend("bottomright", legend="C", cex=0.8, bty="n")
# Legend
plot(1, type="n", axes=FALSE, ann=FALSE)
legend("center", pch=c(3, rep(16, 5)), col=c("black", iro),
       legend=c(1981, 2011:2015)
)
mtext(1, line=2, outer=TRUE, text="Axis 1")
mtext(2, line=2, outer=TRUE, text="Axis 2")
# Export as eps file
dev.copy2eps(file="Fig2.eps")

# Dissimilarity index analysis
plotid <- data.frame(numid=1:nrow(d4), plotid=paste(d4$Treat, d4$Plot, d4$Year, sep="_"))
library(vegan)
di <- data.frame(DI=c(vegdist(spmat, method="chao")))
di <- cbind(di, t(combn(1:183, m=2)))
colnames(di)[2:3] <- c("T1", "T2")
di2 <- di[di$T1<26, ]
colnames(di2)[3] <- "numid"
di3 <- merge(di2, plotid)
temp <- strsplit(as.character(di3$plotid), "\\_")
di3$Treat <- sapply(temp, "[", 1)
di3$Plot <- sapply(temp, "[", 2)
di3$Plot <- paste(di3$Treat, di3$Plot, sep="_")
di3$Year <- as.numeric(sapply(temp, "[", 3))
di3 <- di3[di3$Treat!="Old", ]
di3$YAF <- 0
di3$YAF <- ifelse(di3$Treat=="fence2010", di3$Year-2010, di3$YAF)
di3$YAF <- ifelse(di3$Treat=="fence2011", di3$Year-2011, di3$YAF)
di3$DI <- ifelse(di3$DI==1, 0.99999, di3$DI)
di3$DI <- ifelse(di3$DI==0, 0.00001, di3$DI)
di3 <- rm.level(di3)

# Load mean coverage of each life form group
dl3 <- rm.level(dl3[dl3$Year!=1981, ])
dl3$plotid <- paste(dl3$Treat, dl3$Plot, sep="_")
sn5 <- as.data.frame(ftable(tapply(dl3$C, dl3[, c("TAX", "plotid", "Year")], mean)))
sn5$Freq <- ifelse(is.na(sn5$Freq), 0, sn5$Freq)
sn5$Treat <- as.factor(sapply(strsplit(as.character(sn5$plotid), "\\_"), "[", 1))
sn5$plotid <- paste(sn5$plotid, sn5$Year, sep="_")
sn7 <- reshape(sn5, idvar="plotid", timevar="TAX",
               v.names="Freq", direction="wide"
)
di4 <- merge(di3, sn7[, c("plotid", "Freq.DC",
                          "Freq.GR", "Freq.MN",
                          "Freq.PT", "Freq.WD")], all=FALSE
)

# Make figure
par(mar=c(5,5,1,1), ps=15)
sym <- c(1,0,2)
heikin <- tapply(di4$DI, di4[, c("Year", "Treat")], mean, na.rm=TRUE)
SD <- tapply(di4$DI, di4[, c("Year", "Treat")], sd, na.rm=TRUE)
matplot(2011:2015, heikin, ylim=c(0.4, 1.2), type="o", pch=sym, lty=1, cex=2.0,
        xlab="Year", ylab="Dissimilariy index", col="black"
)
for (i in 1:3) {
     arrows(2011:2015, heikin[, i]-SD[, i], 2011:2015, heikin[, i]+SD[, i],
            code=3, angle=90, length=0.1, pch=sym
     )
}
legend("bottomleft", lty=1, col="black", pch=sym,
       legend=c("Fence 2010", "Fence 2011", "Outside")
)
dev.copy2eps(file="Fig3.eps")

# GLMM for the dissimilarity index
library(glmmTMB)
test2 <- glmmTMB(DI ~ Treat + YAF + Freq.GR
                      + (1|Plot) + (1|Year),
                      family=beta_family(), di4, na.action = "na.fail"
)
# Export the summary of GLMM
write.csv(summary(test2)$coefficients$cond, file="GLMM_DI_res.csv")

library(MuMIn)
dredge(test2, extra = alist(AIC))

## Fig. 4
# The relationship of heights between forbs and graminoids
library(openxlsx)
h <- read.xlsx("height.xlsx")
h2 <- h[, c(1:3, grep("H", colnames(h)))]
h2$totalid <- paste(h$Treat, h$quadrat, h$SpJ, sep="_")
h2 <- rm.level(h2[apply(h2[, 4:8], 1, sum, na.rm=TRUE)!=0, ])
hl <- reshape(h2, varying=list(4:(ncol(h2)-1)), idvar="totalid",
              timevar="Year", times=2011:2015, direction="long"
)
hl$totalid <- NULL
colnames(hl)[grep("H_2011", colnames(hl))] <- "H"
hl$SpJ <- as.character(hl$SpJ)
# Merge some synonym
hl$SpJ <- ifelse(hl$SpJ=="ウマノアシガタ", "キンポウゲ", hl$SpJ)
hl$SpJ <- as.factor(hl$SpJ)
hl <- merge(hl, st2[, c("SpJ", "Sp", "TAX")])
h3 <- hl[!is.na(hl$H), ]
# Merge DC and MN
h3$TAX <- ifelse(h3$TAX=="DC", "FO", h3$TAX)
h3$TAX <- ifelse(h3$TAX=="MN", "FO", h3$TAX)
# 2011
h11 <- h3[h3$Year==2011, ]
h11$plotid <- paste(h11$Treat, h11$quadrat, sep="_")
h11_2 <- as.data.frame(tapply(h11$H, h11[, c("plotid", "TAX")], max))
h11_2$Treat <- as.factor(sapply(strsplit(rownames(h11_2), "\\_"), "[", 1))
# 2015
h15 <- h3[h3$Year==2015, ]
h15$plotid <- paste(h15$Treat, h15$quadrat, sep="_")
h15_2 <- as.data.frame(tapply(h15$H, h15[, c("plotid", "TAX")], max))
h15_2$Treat <- as.factor(sapply(strsplit(rownames(h15_2), "\\_"), "[", 1))
## make figure
# Graphics parameters
par(mfrow=c(1, 2), mar=c(1,1,0,0), oma=c(4,4,1,1), ps=15)
sym <- c(1,0,2)
xmax <- 160
ymax <- 160
# 2011
plot(FO ~ GR, pch=sym[as.numeric(h11_2$Treat)],
     xlab="", ylab="", xaxt="n", yaxt="n", main="",
     xlim=c(0, xmax), ylim=c(0, ymax), h11_2
)
axis(1, at=0:3*50)
axis(2, at=0:3*50)
abline(a=0, b=1, lty=2)
legend("topleft", pch=sym, legend=levels(h15_2$Treat))
legend("bottomright", legend="A", bty="n")
# 2015
plot(FO ~ GR, pch=sym[as.numeric(h15_2$Treat)],
     xlab="", ylab="", xaxt="n", yaxt="n", main="",
     xlim=c(0, xmax), ylim=c(0, ymax), h15_2
)
axis(1, at=0:3*50)
axis(2, at=0:3*50, labels=FALSE)
abline(a=0, b=1, lty=2)
legend("bottomright", legend="B", bty="n")
# Labels
mtext(1, line=2, outer=TRUE, text="Maximum height of graminoids (cm)")
mtext(2, line=2, outer=TRUE, text="Maximum height of dicots and monocots (cm)")
# Export as eps
dev.copy2eps(file="Fig4.eps")