This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
library(vegan)
library(ggplot2)
library(ggsci)
library(gplots)
library(iNEXT)
library(nlme)
library(rcompanion)
library(car)
set.seed(4321)
library(vegan)
要求されたパッケージ permute をロード中です
要求されたパッケージ lattice をロード中です
This is vegan 2.5-7
library(ggplot2)
library(ggsci)
library(gplots)
次のパッケージを付け加えます: ‘gplots’
以下のオブジェクトは ‘package:stats’ からマスクされています:
lowess
library(iNEXT)
library(nlme)
library(rcompanion)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'DescTools':
method from
reorder.factor gplots
library(car)
要求されたパッケージ carData をロード中です
set.seed(4321)
uds_data5.temp <- read.csv("monthly-for-R3.csv",header=T) #for 2 or 3 continuous sampling at day 10, 32, and 60
uds_data_all_day <-uds_data5.temp[,c(-1)]
uds_data_all_day2 <-uds_data5.temp
#converting treatment to chamber
chamber <- uds_data_all_day2$treatment
chamber[chamber == "N"] <- "C1"
chamber[chamber == "W"] <- "C2"
chamber[chamber == "T"] <- "C2"
uds_data_all_day2 <- cbind.data.frame(uds_data_all_day2, chamber)
uds_data_all_LBa <- read.csv("leaf_beetle_adult_summary.csv", header=T) #for leaf beetle adult data
uds_dataLB <- uds_data_all_day[, c(3)]
head(uds_data_all_LBa)
head(uds_data_all_day2)
SAD_total <- data.frame(abundance=apply(uds_data_all_day2[,4:76], 2, sum))
SAD_rank <- order(SAD_total[,1], decreasing=T)
SAD_total2 <- data.frame(rank=c(1:73), ID=rownames(SAD_total)[SAD_rank], abundance=SAD_total[SAD_rank,1])
#how many species?
sum(SAD_total2$abundance>0)
[1] 48
plot(SAD_total2$rank, log10(SAD_total2$abundance)) #SAD, including the rank 1 speices (wl) but with log10 scale
plot(SAD_total2$rank[-1], SAD_total2$abundance[-1]) #SAD, excluding the rank 1 speices (wl)
plot(SAD_total2$rank[2:10], SAD_total2$abundance[2:10]) #SAD, from rank 2 to rank 10
#text(SAD_total2$rank[2:10], SAD_total2$abundance[2:10], labels=SAD_total2$ID[2:10])
ggplot(SAD_total2, aes(x=rank, y=log10(abundance))) + geom_point()
SAD_total2_sub <- SAD_total2[2:10,]
sad_plot <- ggplot(SAD_total2_sub, aes(x=rank, y=abundance)) + geom_point()
sad_plot <- sad_plot + geom_text(aes(label=ID), size=3,vjust=-1)
plot(sad_plot)
#Taxon (species) identity
SAD_total2_sub$ID
[1] "spider3" "lep9" "lep7" "lep10" "aphi5" "spider4" "parasitoid3"
[8] "lep12" "lep13"
uds_data_all_day_SI <- uds_data_all_day2 #copy for SI only
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="N"] <- "U"
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="W"] <- "D"
uds_data_all_day_SI$treatment[uds_data_all_day_SI$treatment=="T"] <- "E"
treatment_day <- as.factor(paste(uds_data_all_day_SI$treatment,uds_data_all_day_SI$day))
uds_data_all_day_SI <- cbind.data.frame(uds_data_all_day_SI, treatment_day)
plot(uds_data_all_day_SI$treatment_day ,uds_data_all_day_SI$wl)
#rank 1
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=wl, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.2) + ggtitle("Plagiodera versicolora")
#rank 2
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=spider3, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("spider group 3 (< 5mm with spiderweb)")
#rank 3
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep9, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Caloptilia sp.2")
#rank 4
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep7, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Saliciphaga caesia")
#rank 5
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep10, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Phyllocolpa sp")
#rank 6
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=aphi5, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Cavariella salicicola")
#rank 7
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep12, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("spider group 4 (< 5mm without spiderweb)")
#rank 8
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=parasitoid3, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Tachinidae sp.1 (host: P. versicolora)")
#rank 9
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep12, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Microleon longipalpis")
#rank 10
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=lep13, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Monema flavescens")
#rank 11
ggplot(uds_data_all_day_SI, aes(x=treatment_day, y=grasshopper1, fill=treatment)) + geom_boxplot(outlier.shape = NA) + geom_point(alpha=0.1, size=2) + ggtitle("Orthoptera sp.1")
NA
NA
Ref: http://rcompanion.org/handbook/I_09.html
chamber <- append(rep("C1",16),rep("C2",32))
uds_data_all_LBa <- cbind.data.frame(uds_data_all_LBa, chamber)
uds_data_all_LBa
#Data decomposition
uds_data_LBa24h <- subset(uds_data_all_LBa, uds_data_all_LBa$day==1)
uds_data_LBa2_7d <- subset(uds_data_all_LBa, (uds_data_all_LBa$day >= 2 & uds_data_all_LBa$day <= 7))
uds_data_LBa10_60d <- subset(uds_data_all_LBa, (uds_data_all_LBa$day >= 10 & uds_data_all_LBa$day <= 60))
#short (Fig.1a)
model24h_0 = gls(adult ~ chamber + hour + chamber*hour,data=uds_data_LBa24h)
#Anova(model24h_0)
ACF_s <- ACF(model24h_0, form = ~ hour|id)
head(ACF_s)
model24h = gls(adult ~ chamber + hour + chamber*hour, correlation = corAR1(form = ~ hour|id, value = ACF_s$ACF[2]),data=uds_data_LBa24h)
Anova(model24h)
Analysis of Deviance Table (Type II tests)
Response: adult
Df Chisq Pr(>Chisq)
chamber 1 4.7975 0.02850 *
hour 1 5.9301 0.01488 *
chamber:hour 1 5.8661 0.01544 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Post-hoc test
data12h<-subset(uds_data_LBa24h, hour==12)
pairwise.t.test(data12h$adult, data12h$chamber,p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data12h$adult and data12h$chamber
C1
C2 0.015
P value adjustment method: holm
pairwise.t.test(data12h$adult[-c(1:16)],data12h$treatment[-c(1:16)],p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data12h$adult[-c(1:16)] and data12h$treatment[-c(1:16)]
D
S 0.33
P value adjustment method: holm
data24h<-subset(uds_data_LBa24h, hour==24)
pairwise.t.test(data24h$adult, data24h$chamber,p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data24h$adult and data24h$chamber
C1
C2 0.027
P value adjustment method: holm
pairwise.t.test(data12h$adult[-c(1:16)], data12h$treatment[-c(1:16)],p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data12h$adult[-c(1:16)] and data12h$treatment[-c(1:16)]
D
S 0.33
P value adjustment method: holm
#Summary statistics for graph
sum24h = groupwiseMean(adult ~ treatment + hour, data = uds_data_LBa24h, conf = 0.95, digits = 3,traditional = FALSE, percentile = TRUE)
sum24h
#shortest graph (Fig.1a)
pd = position_dodge(1)
g_line<-ggplot(sum24h, aes(y=Mean, x=hour)) + xlim(0,25) + ylim(0,1.1) + coord_fixed(25/1.1)
g_line<-g_line+geom_line(aes(colour=treatment),size=1,position=pd)
g_line<- g_line+ geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd)
print(g_line)
#middle (Fig.1b)
model2_7d0 = gls(adult ~ chamber + day + chamber*day,data=uds_data_LBa2_7d)
#Anova(model2_7d0)
ACF_m <- ACF(model2_7d0, form = ~ day|id)
head(ACF_m)
model2_7d = gls(adult ~ chamber + day + chamber*day, correlation = corAR1(form = ~ day|id, value = ACF_m$ACF[2]),data=uds_data_LBa2_7d)
Anova(model2_7d)
Analysis of Deviance Table (Type II tests)
Response: adult
Df Chisq Pr(>Chisq)
chamber 1 0.1239 0.724869
day 1 16.9347 3.869e-05 ***
chamber:day 1 6.6496 0.009918 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Post-hoc test
data7d<-subset(uds_data_LBa2_7d, day==7)
pairwise.t.test(data7d$adult, data7d$chamber,p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data7d$adult and data7d$chamber
C1
C2 0.21
P value adjustment method: holm
pairwise.t.test(data7d$adult[-c(1:16)], data7d$treatment[-c(1:16)],p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data7d$adult[-c(1:16)] and data7d$treatment[-c(1:16)]
D
S 0.94
P value adjustment method: holm
#Summary statistics for graph
sum2_7d = groupwiseMean(adult ~ treatment + day, data = uds_data_LBa2_7d, conf = 0.95, digits = 3,traditional = FALSE, percentile = TRUE)
#a week scale graph (Fig.1b)
pd = position_dodge(0.3)
g_line<-ggplot(sum2_7d, aes(y=Mean, x=day)) + xlim(1,8) + ylim(0,1.1) + coord_fixed(7/2.0)
g_line<-g_line+geom_line(aes(colour=treatment),size=1, position=pd)
g_line<- g_line+ geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd)
print(g_line)
model10_60d0 = gls(adult ~ chamber + day + chamber*day, data=uds_data_LBa10_60d)
#Anova(model10_60d0)
ACF_l <- ACF(model10_60d0, form = ~ day|id)
head(ACF_l)
model10_60d = gls(adult ~ chamber + day + chamber*day, correlation = corAR1(form = ~ day|id, value = ACF_l$ACF[2]),data=uds_data_LBa10_60d)
Anova(model10_60d)
Analysis of Deviance Table (Type II tests)
Response: adult
Df Chisq Pr(>Chisq)
chamber 1 2.2242 0.13586
day 1 41.4459 1.212e-10 ***
chamber:day 1 2.8322 0.09239 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Post-hoc test
data10d<-subset(uds_data_LBa10_60d, day==10)
pairwise.t.test(data10d$adult, data10d$chamber,p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data10d$adult and data10d$chamber
C1
C2 0.14
P value adjustment method: holm
pairwise.t.test(data10d$adult[-c(1:16)], data10d$treatment[-c(1:16)],p.adjust.method="holm")
Pairwise comparisons using t tests with pooled SD
data: data10d$adult[-c(1:16)] and data10d$treatment[-c(1:16)]
D
S 0.041
P value adjustment method: holm
#Summary statistics for graph
sum10_60d = groupwiseMean(adult ~ treatment + day, data = uds_data_LBa10_60d, conf = 0.95, digits = 3,traditional = FALSE, percentile = TRUE)
#graphics
pd = position_dodge(3)
g_line<-ggplot(sum10_60d, aes(y=Mean, x=day)) + xlim(5,65) + ylim(0,4) + coord_fixed(ratio = 60/4)
g_line<-g_line+geom_line(aes(colour=treatment),size=1, position=pd)
g_line<- g_line+geom_errorbar(aes(ymin=Percentile.lower, ymax = Percentile.upper, colour=treatment),alpha=0.9,size=1, width=0, position=pd)
g_line<-g_line+geom_point(aes(colour=treatment, shape=treatment),size=3,alpha=1.0, position=pd)
print(g_line)
disper_uds_centroid_graph <- function(data_name, c_LB, meth, adj=NULL, type="violin", lim_min=-0.6, lim_max=0.4) {
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-uds_data1[,c(-1,-2,-3, -77,-78, -79, -80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-uds_data1[,c(-1,-2,-3,-4,-77,-78,-79, -80)]
} #When leaf beele is NOT included
if(c_LB == "LBonly") {
uds_data1 <- subset(data_name, abundance2 > 0)
uds_each <- subset(uds_data1, day==date)[, c(4)] #picking up LB only
}
sample_each <- uds_data1$sample
#meth="chao"
uds_bray_each0.d<-vegdist(uds_each, method=meth)
if(c_LB == "LBonly") tp <- "centroid"
else tp <- "median"
uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
uds_median <- data.frame(uds_dist_each0$centroids) #calculate the centroid (spatial median) for repeated measures
sample <- row.names(uds_median)
#sample <- sample_each
#Replace the sample name into the treatment
sample[grep("1n", sample)]<- "U10"
sample[grep("2n", sample)]<- "U32"
sample[grep("3n", sample)]<- "U60"
sample[grep("1t", sample)]<- "S10"
sample[grep("2t", sample)]<- "S32"
sample[grep("3t", sample)]<- "S60"
sample[grep("1w", sample)]<- "D10"
sample[grep("2w", sample)]<- "D32"
sample[grep("3w", sample)]<- "D60"
treatment_each <- as.factor(sample) #just conversion
uds_each.d<-vegdist(uds_median, method="euclidean")
uds_dist_each <-betadisper(uds_each.d, treatment_each, type="centroid", bias.adjust = TRUE) #Choose centroid because of euclidean distance
#Plot
if(c_LB == "LBonly") boxplot(uds_dist_each)
else {
#boxplot(uds_dist_each)
#ordiplot(uds_dist_each$centroids, type="text", choices=c(1,2), cex=1, xlim=c(-0.6,0.4), ylim=c(-0.6, 0.5))
#par(new=T)
if(type == "violin") {
#boxplot(uds_dist_each)
##from ggplot2
#convert to data.frame (which is necessary for ggplot2)
uds_dist_each2 <- cbind.data.frame(uds_dist_each$distances, as.character(uds_dist_each$group))
treatment<-as.character(uds_dist_each$group)
#Replace the sample name into the treatment
treatment[grep("U", treatment)]<- "U"
treatment[grep("D", treatment)]<- "D"
treatment[grep("S", treatment)]<- "S"
uds_dist_each2<-cbind.data.frame(uds_dist_each2, treatment)
colnames(uds_dist_each2) <- c("distances", "group", "treatment")
uds_dist_each2$distances <- as.numeric(as.character(uds_dist_each2$distances))
g_box<-ggplot(uds_dist_each2, aes(y=distances, x=group,colour=treatment))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(fill="black", width=0.1)
g_box<-g_box + coord_fixed(ratio = 15/1)
print(g_box)
}
else if(type == "PCoA") {
#PCoA
test0<-capscale(uds_each.d~1)
test<-cmdscale(uds_each.d, eig=FALSE, k = 2)
#prepare dataframe for each point
test2<-summary(test)
#uds_dist_each_test<-data.frame(PCoA1=-test2$sites[,1], PCoA2=-test2$sites[,2])
uds_dist_each_test<-as.data.frame(test)
uds_dist_each_test<-cbind.data.frame(sample, uds_dist_each_test)
treatment<-as.character(sample)
day<-as.character(sample)
#Replace the sample name into the treatment & day
treatment[grep("U", treatment)]<- "U"
treatment[grep("D", treatment)]<- "D"
treatment[grep("S", treatment)]<- "S"
day[grep("10", day)]<-"Day10"
day[grep("32", day)]<-"Day32"
day[grep("60", day)]<-"Day60"
uds_dist_each_test<-cbind.data.frame(cbind.data.frame(treatment, day), uds_dist_each_test)
colnames(uds_dist_each_test)<-c("treatment", "day", "group", "PCoA1", "PCoA2")
class(uds_dist_each_test$day)
#prepare dataframe for centroid
uds_dist_each_centr<-as.data.frame(uds_dist_each$centroids)
treatment<-as.character(rownames(uds_dist_each_centr))
day<-as.character(rownames(uds_dist_each_centr))
#Replace the sample name into the treatment & day
treatment[grep("U", treatment)]<- "U"
treatment[grep("D", treatment)]<- "D"
treatment[grep("S", treatment)]<- "S"
day[grep("10", day)]<-"Day10"
day[grep("32", day)]<-"Day32"
day[grep("60", day)]<-"Day60"
uds_dist_each_centr<-cbind.data.frame(uds_dist_each_centr, data.frame(treatment=treatment, day=day, group=rownames(uds_dist_each_centr)))
g_PCoA <- ggplot(uds_dist_each_test, aes(x=PCoA1,y=PCoA2))
g_PCoA <- g_PCoA + geom_point(aes(colour=treatment, shape=day), size=3, alpha=0.5)
g_PCoA <- g_PCoA + xlim(lim_min, lim_max) + ylim(lim_min,lim_max) + coord_fixed(ratio = 1/1)
#Add the centroid information
g_PCoA <- g_PCoA + layer(
data=uds_dist_each_centr,
mapping=aes(x=PCoA1, y=PCoA2, colour=treatment, shape=day),
params=list(size=5),
geom="point",
stat="identity",
position="identity"
)
print(g_PCoA)
#Performance of PCoA http://d.hatena.ne.jp/hoxo_m/20120313/p1
summary(test0)
}
}
#Performance of PCoA http://d.hatena.ne.jp/hoxo_m/20120313/p1
#summary(test0)
if(type == "NMDS") {
uds_each.mds<-metaMDS(uds_median, distance="euclidean")
rownames(uds_each.mds$points)<-treatment_each
ordiplot(uds_each.mds$points, type='points', xlim=c(-1.0, 1.0))
text(uds_each.mds$points,labels=treatment_each, col=unclass(treatment_each),cex=0.8)
}
}
#Figure 2a
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="chao", type="violin")
警告: some squared distances are negative and changed to zero
#Figure 2b
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="chao", type="PCoA")
警告: some squared distances are negative and changed to zero
####For Figure S1a
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="woLB", meth="bray", type="PCoA", lim_min=-0.75, lim_max=0.5)
警告: some squared distances are negative and changed to zero
####For Figure S1b
disper_uds_centroid_graph(data_name=uds_data_all_day2, c_LB="wLB", meth="chao", type="PCoA", lim_min=-0.75, lim_max=0.5)
警告: some squared distances are negative and changed to zero
#///////////////////////Pre-Evaluation of effectiveness of PERMANOVA (or can be used for assessing beta diversity) through PERMDISP for each of 10d, 32d, 60d/////////////////////
disper_uds_centroid <- function(data_name, date, c_LB, meth, adj=NULL, plot="box", dist_out=FALSE) {
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3, -77,-78, -79, -80)]
} #When leaf beetle is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-4,-77,-78,-79, -80)]
} #When leaf beetle is NOT included
if(c_LB == "LBonly") { #For leaf beetle only but consider the case with other insects only
uds_data1 <- subset(data_name, abundance2 > 0)
#uds_data1 <-data_name
uds_each <- subset(uds_data1, day==date)[, c(4)] #picking up LB only
}
sample_each <- subset(uds_data1, day==date)$sample
#length(uds_each)
#length(sample_each)
##In case of leaf beetle only, different calculation
if(c_LB == "LBonly") {
sample_ID<-as.factor(sample_each)
LB_average<-tapply(uds_each, sample_ID, mean) #calculate the average of short-term repeated measurement
LB_average2<-LB_average #just copy
sample <- row.names(LB_average2)
sample2 <- row.names(LB_average2)
#Replace the sample name into the treatment
sample[grep("n", sample)]<- "U" #undamaged
sample[grep("t", sample)]<- "S" #Signaled
sample[grep("w", sample)]<- "D" #Damaged
treatment_each <- as.factor(sample) #just conversion
#print(treatment_each)
#Replace the sample name into the chamber
sample2[grep("n", sample2)]<- "C1" #undamaged
sample2[grep("t", sample2)]<- "C2" #Signaled
sample2[grep("w", sample2)]<- "C2" #Damaged
chamber_each <- as.factor(sample2) #just conversion
average_treatment<-tapply(LB_average, treatment_each, mean) #calculate the average of U, D, S, respectively
average_chamber<-tapply(LB_average2, chamber_each, mean) #calculate the average of C1 amd C2
#print(average_treatment)
for(i in 1: length(LB_average)) { #normalizing the original data to exclude the effect of average on the distance to centroid
if(treatment_each[i] == "D" && average_treatment[1] > 0.0) LB_average[i] = LB_average[i]/average_treatment[1]
if(treatment_each[i] == "S" && average_treatment[2] > 0.0) LB_average[i] = LB_average[i]/average_treatment[2]
if(treatment_each[i] == "U" && average_treatment[3] > 0.0) LB_average[i] = LB_average[i]/average_treatment[3]
}
for(i in 1: length(LB_average)) { #normalizing the original data to exclude the effect of average on the distance to centroid
if(chamber_each[i] == "C1" && average_chamber[1] > 0.0) LB_average2[i] = LB_average2[i]/average_chamber[1]
if(chamber_each[i] == "C2" && average_chamber[2] > 0.0) LB_average2[i] = LB_average2[i]/average_chamber[2]
}
uds_median<-LB_average #send the data for further analysis, distinguishing D,S,and U
uds_median2<-LB_average2 #send the data for further analysis, distinguishing C1 and C2
#print(uds_median)
if(dist_out == "C1") uds_median <- subset(uds_median, chamber_each=="C1")
else if(dist_out == "C2") uds_median <- subset(uds_median, chamber_each=="C2")
#print(uds_median)
}
else {
uds_bray_each0.d<-vegdist(uds_each, method=meth)
tp <- "median"
uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
uds_median <- data.frame(uds_dist_each0$centroids) #calculate the centroid (spatial median) for repeated measures
sample <- row.names(uds_median)
sample2 <- sample #just copy
#Replace the sample name into the treatment
sample[grep("n", sample)]<- "U" #undamaged
sample[grep("t", sample)]<- "S" #Signaled
sample[grep("w", sample)]<- "D" #Damaged
treatment_each <- as.factor(sample) #just conversion
#Replace the sample name into the chamber
sample2[grep("n", sample2)]<- "C1" #undamaged
sample2[grep("t", sample2)]<- "C2" #Signaled
sample2[grep("w", sample2)]<- "C2" #Damaged
chamber_each <- as.factor(sample2) #just conversion
if(dist_out == "C1") uds_median <- uds_median[chamber_each=="C1",]
else if(dist_out == "C2") uds_median <- uds_median[chamber_each=="C2",]
#print(uds_median)
}
uds_each.d<-vegdist(uds_median, method="euclidean")
if(dist_out!= FALSE) return(uds_each.d)
#treatment_each
uds_dist_each <-betadisper(uds_each.d, treatment_each, type="centroid", bias.adjust = TRUE) #Choose centroid because of euclidean distance
#chamber_each
uds_dist_each_C <-betadisper(uds_each.d, chamber_each, type="centroid", bias.adjust = TRUE)
#Plot
if(c_LB == "LBonly") boxplot(uds_dist_each)
else {
if(plot=="box") {
boxplot(uds_dist_each_C)
boxplot(uds_dist_each)
}
else{
plot(uds_dist_each)
ordilabel(scores(uds_dist_each, "centroids"))
plot(uds_dist_each_C)
ordilabel(scores(uds_dist_each_C, "centroids"))
}
}
#PERDISP & permutation
perm_result_C <- permutest(uds_dist_each_C, pairwise=FALSE, permutations=9999)
perm_result <- permutest(uds_dist_each, pairwise=TRUE, permutations=9999)
cat("Sampling day was\n")
print(subset(uds_data1, day==date)$day[1])
cat("Comparison between chambers\n")
print(perm_result_C)
cat("Comparison between treatments\n")
print(perm_result)
#In fact, the adjustment is not necessary because permutest.betadisper does not repeat permutation for each pair
if(!is.null(adj)) {
pairwise_adj<-p.adjust(perm_result$pairwise$permuted, method = adj , n = length(perm_result$pairwise$permuted))
cat("Pairwise comparision by permutation, adjusted by")
print(as.character(adj))
print(pairwise_adj[1])
print(pairwise_adj[2])
print(pairwise_adj[3])
}
}
#///PERMANOVA for each initial conditions, based on the spatial median for repeated measures
com_div_uds_centroid <- function(data_name, date, c_LB, meth, adj)
{
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-77,-78, -79, -80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-subset(uds_data1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79, -80)]
uds_each_LB<-subset(uds_data1, day==date)[,c(4)] #extract leaf beetle abundance
} #When leaf beele is NOTincluded
sample_each <- subset(uds_data1, day==date)$sample
uds_each.d<-vegdist(uds_each, method=meth)
uds_dist_each <-betadisper(uds_each.d, sample_each, type="median", bias.adjust = TRUE)
uds_median <- data.frame(uds_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
sample <- row.names(uds_median)
sample2 <- row.names(uds_median)
#Replace the sample name into the treatment
sample[grep("n", sample)]<- "U"
sample[grep("t", sample)]<- "S"
sample[grep("w", sample)]<- "D"
treatment_each <- as.factor(sample) #just conversion
#Replace the sample name into the chamber
sample2[grep("n", sample2)]<- "C1" #undamaged
sample2[grep("t", sample2)]<- "C2" #Signaled
sample2[grep("w", sample2)]<- "C2" #Damaged
chamber_each <- as.factor(sample2) #just conversion
adonis_result_C<-adonis(uds_median ~ chamber_each, method="euclidean", permutations=9999)
adonis_result_T<-adonis(uds_median ~ treatment_each, method="euclidean", permutations=9999)
print(adonis_result_C) #Chamber effect
print(adonis_result_T) #Treatment effect
#For pairwise comparison (only between D and S)
cat("The pairwise PERMANOVA (adonis) analysis:\n")
uds_dataDS0<-subset(data_name, treatment != "N") #We still need the symbols N, W, and T because it is used in the raw data files
uds_dataSU0<-subset(data_name, treatment != "W")
uds_dataUD0<-subset(data_name, treatment != "T")
if(c_LB == "wLB") {
uds_dataDS1 <-subset(uds_dataDS0, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
DS_each<-subset(uds_dataDS1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
uds_dataSU1 <-subset(uds_dataSU0, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
SU_each<-subset(uds_dataSU1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
uds_dataUD1 <-subset(uds_dataUD0, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
UD_each<-subset(uds_dataUD1, day==date)[,c(-1,-2,-3,-77,-78,-79,-80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_dataDS1 <-subset(uds_dataDS0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
DS_each<-subset(uds_dataDS1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79, -80)]
uds_dataSU1 <-subset(uds_dataSU0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
SU_each<-subset(uds_dataSU1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79,-80)]
uds_dataUD1 <-subset(uds_dataUD0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
UD_each<-subset(uds_dataUD1, day==date)[,c(-1,-2,-3,-4,-77,-78, -79,-80)]
} #When leaf beetle is NO Tincluded
sample_DS <- subset(uds_dataDS1, day==date)$sample
sample_SU <- subset(uds_dataSU1, day==date)$sample
sample_UD <- subset(uds_dataUD1, day==date)$sample
DS_bray_each.d<-vegdist(DS_each, method=meth)
SU_bray_each.d<-vegdist(SU_each, method=meth)
UD_bray_each.d<-vegdist(UD_each, method=meth)
DS_dist_each <-betadisper(DS_bray_each.d, sample_DS, type="median", bias.adjust = TRUE)
SU_dist_each <-betadisper(SU_bray_each.d, sample_SU, type="median", bias.adjust = TRUE)
UD_dist_each <-betadisper(UD_bray_each.d, sample_UD, type="median", bias.adjust = TRUE)
DS_median <- data.frame(DS_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
SU_median <- data.frame(SU_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
UD_median <- data.frame(UD_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
DS_sample <- row.names(DS_median)
SU_sample <- row.names(SU_median)
UD_sample <- row.names(UD_median)
#Replace the sample name into the treatment
DS_sample[grep("n", DS_sample)]<- "U"
DS_sample[grep("t", DS_sample)]<- "S"
DS_sample[grep("w", DS_sample)]<- "D"
treatment_DS <- as.factor(DS_sample) #just conversion
SU_sample[grep("n", SU_sample)]<- "U"
SU_sample[grep("t", SU_sample)]<- "S"
SU_sample[grep("w", SU_sample)]<- "D"
treatment_SU <- as.factor(SU_sample) #just conversion
UD_sample[grep("n", UD_sample)]<- "U"
UD_sample[grep("t", UD_sample)]<- "S"
UD_sample[grep("w", UD_sample)]<- "D"
treatment_UD <- as.factor(UD_sample) #just conversion
adonis_resultDS<-adonis(DS_median ~ treatment_DS, method="euclidean", permutations=9999)
adonis_resultSU<-adonis(SU_median ~ treatment_SU, method="euclidean", permutations=9999)
adonis_resultUD<-adonis(UD_median ~ treatment_UD, method="euclidean", permutations=9999)
p_list <- c(adonis_resultDS$aov$Pr[1], adonis_resultSU$aov$Pr[1], adonis_resultUD$aov$Pr[1])
cat("The comparision between D and S is:\n")
print(adonis_resultDS)
pairwise_adj<-p.adjust(p_list, method = adj , n = length(p_list))
#cat("Pairwise comparision by permutation, adjusted by")
#print(adj)
#cat("The adjusted P-value for S vs D was\n")
#print(pairwise_adj[1])
#cat("The adjusted P-value for U vs S was\n")
#print(pairwise_adj[2])
#cat("The adjusted P-value for D vs U was\n")
#print(pairwise_adj[3])
}
#///PERMANOVA for each date, based on the spatial median for repeated measures
com_div_DAY_centroid <- function(data_name, treat, c_LB, meth, adj=NULL) {
#data_name <-uds_data_all_day2
#treat <- "T"
#meth="chao"
#date <- "10d","32d", or "60d", or any
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3, -77,-78, -79,-80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
} #When leaf beele is NOTincluded
if(c_LB == "LBonly") {
uds_data1 <- subset(data_name, abundance2 > 0)
uds_each <- subset(uds_data1, treatment==treat)[, c(4)] #picking up LB only
}
sample_each <- subset(uds_data1, treatment==treat)$sample
day_each <- subset(uds_data1, treatment==treat)$day
uds_bray_each0.d<-vegdist(uds_each, method=meth)
if(c_LB == "LBonly") tp <- "centroid"
else tp <- "median"
uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
uds_median <- data.frame(uds_dist_each0$centroids) #calculate the centroid (spatial median) for repeasted measures
sample <- row.names(uds_median)
#Replace the sample name into the date
sample[grep("^1", sample)]<- "day10"
sample[grep("^2", sample)]<- "day32"
sample[grep("^3", sample)]<- "day60"
day_each <- as.factor(sample) #just conversion
adonis_result<-adonis(uds_median ~ day_each, method="euclidean", permutations=9999)
print(adonis_result)
#For pairwise comarison (with Bonferroni adjustment)
cat("The pairwise PERMANOVA (adonis) analysis:\n")
uds_dataDay1032_0<-subset(data_name, day != "60d")
uds_dataDay1060_0<-subset(data_name, day != "32d")
uds_dataDay3260_0<-subset(data_name, day != "10d")
if(c_LB == "wLB") {
uds_dataDay1032_1 <-subset(uds_dataDay1032_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day1032_each<-subset(uds_dataDay1032_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
uds_dataDay1060_1 <-subset(uds_dataDay1060_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day1060_each<-subset(uds_dataDay1060_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
uds_dataDay3260_1 <-subset(uds_dataDay3260_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day3260_each<-subset(uds_dataDay3260_1, treatment==treat)[,c(-1,-2,-3,-77,-78,-79,-80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_dataDay1032_1 <-subset(uds_dataDay1032_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day1032_each<-subset(uds_dataDay1032_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
uds_dataDay1060_1 <-subset(uds_dataDay1060_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day1060_each<-subset(uds_dataDay1060_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
uds_dataDay3260_1 <-subset(uds_dataDay3260_0, abundance2 > 0) #delete the sampling data with 0 total abundance when leaf beetle
Day3260_each<-subset(uds_dataDay3260_1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
} #When leaf beele is NOTincluded
sample_1032 <- subset(uds_dataDay1032_1, treatment==treat)$sample
sample_1060 <- subset(uds_dataDay1060_1, treatment==treat)$sample
sample_3260 <- subset(uds_dataDay3260_1, treatment==treat)$sample
Day1032_each.d<-vegdist(Day1032_each, method=meth)
Day1060_each.d<-vegdist(Day1060_each, method=meth)
Day3260_each.d<-vegdist(Day3260_each, method=meth)
Day1032_dist_each <-betadisper(Day1032_each.d, sample_1032, type="median", bias.adjust = TRUE)
Day1060_dist_each <-betadisper(Day1060_each.d, sample_1060, type="median", bias.adjust = TRUE)
Day3260_dist_each <-betadisper(Day3260_each.d, sample_3260, type="median", bias.adjust = TRUE)
Day1032_median <- data.frame(Day1032_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
Day1060_median <- data.frame(Day1060_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
Day3260_median <- data.frame(Day3260_dist_each$centroids) #calculate the centroid (spatial median) for repeasted measures
Day1032_sample <- row.names(Day1032_median)
Day1060_sample <- row.names(Day1060_median)
Day3260_sample <- row.names(Day3260_median)
#Replace the sample name into the day
Day1032_sample[grep("^1", Day1032_sample)]<- "day10"
Day1032_sample[grep("^2", Day1032_sample)]<- "day32"
day_Day1032 <- as.factor(Day1032_sample) #just conversion
Day1060_sample[grep("^1", Day1060_sample)]<- "day10"
Day1060_sample[grep("^3", Day1060_sample)]<- "day60"
day_Day1060 <- as.factor(Day1060_sample) #just conversion
Day3260_sample[grep("^2", Day3260_sample)]<- "day32"
Day3260_sample[grep("^3", Day3260_sample)]<- "day60"
day_Day3260 <- as.factor(Day3260_sample) #just conversion
adonis_resultDay1032<-adonis(Day1032_median ~ day_Day1032, method="euclidian", permutations=9999)
adonis_resultDay1060<-adonis(Day1060_median ~ day_Day1060, method="euclidian", permutations=9999)
adonis_resultDay3260<-adonis(Day3260_median ~ day_Day3260, method="euclidian", permutations=9999)
p_list <- c(adonis_resultDay1032$aov$Pr[1], adonis_resultDay1060$aov$Pr[1], adonis_resultDay3260$aov$Pr[1])
pairwise_adj<-p.adjust(p_list, method = adj , n = length(p_list))
cat("Pairwise comparision by permutation, adjusted by")
print(adj)
cat("The adjusted P-value for day10 vs day32 was\n")
print(pairwise_adj[1])
cat("The adjusted P-value for day10 vs day60 was\n")
print(pairwise_adj[2])
cat("The adjusted P-value for day32 vs day60 was\n")
print(pairwise_adj[3])
}
#///PERMDISP or each date, based on the spatial median for repeated measures
disper_DAY_centroid <- function(data_name, treat, c_LB, meth, adj=NULL) {
#resolving inconsistency of the symbols to represent the treatments (and chamber)
if(treat == "U") treat <- "N"
else if(treat == "E") treat <- "T"
else if(treat == "D") treat <- "W"
else if(treat == "S") treat <- "T"
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3, -77,-78, -79,-80)]
} #When leaf beele is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-subset(uds_data1, treatment==treat)[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
} #When leaf beele is NOTincluded
if(c_LB == "LBonly") {
uds_data1 <- subset(data_name, abundance2 > 0)
uds_each <- subset(uds_data1, treatment==treat)[, c(4)] #picking up LB only
}
sample_each <- subset(uds_data1, treatment==treat)$sample
day_each <- subset(uds_data1, treatment==treat)$day
uds_bray_each0.d<-vegdist(uds_each, method=meth)
if(c_LB == "LBonly") tp <- "centroid"
else tp <- "median"
uds_dist_each0 <-betadisper(uds_bray_each0.d, sample_each, type=tp, bias.adjust = TRUE)
uds_median <- data.frame(uds_dist_each0$centroids) #calculate the centroid (spatial median) for repeasted measures
sample <- row.names(uds_median)
#Replace the sample name into the date
sample[grep("^1", sample)]<- "day10"
sample[grep("^2", sample)]<- "day32"
sample[grep("^3", sample)]<- "day60"
day_each <- as.factor(sample) #just conversion
#treatment_each
uds_each.d<-vegdist(uds_median, method="euclidean")
uds_dist_each <-betadisper(uds_each.d, day_each, type="centroid", bias.adjust = TRUE) #Choose centroid becasue of euclidean distance
#Plot
if(c_LB == "LBonly") boxplot(uds_dist_each)
else {
boxplot(uds_dist_each)
ordilabel(scores(uds_dist_each, "centroids"))
}
#PERDISP & ANOVA & MULTIPLE COMPARISON
#anova_dist <- anova(uds_dist_each)
#TukeyHSD(uds_dist_each, which = "group", ordered = FALSE, conf.levels = 0.95)
#PERDIPST & permutation
perm_result <- permutest(uds_dist_each, pairwise=TRUE, permutations=9999)
#perm_result
#print (anova_dist$Pr[1])
cat("Target treatment was\n")
print(treat)
print(uds_dist_each)
print(perm_result)
#In fact, the adjustment is not necessary because permutest.betadisper does not repeat permutation for each pair
if(!is.null(adj)) {
pairwise_adj<-p.adjust(perm_result$pairwise$permuted, method = adj, n = length(perm_result$pairwise$permuted))
cat("Pairwise comparision by permutation, adjusted by")
print(as.character(adj))
print(pairwise_adj[1])
print(pairwise_adj[2])
print(pairwise_adj[3])
}
#TukeyHSD(uds_dist_each)
}
#####For PERMDISP: comparison between chambers and treatments
(Fig.2a)
#Temporal variability in beta diversity for manuscript
disper_DAY_centroid(data_name=uds_data_all_day2, treat="U", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Target treatment was
[1] "N"
Homogeneity of multivariate dispersions
Call: betadisper(d = uds_each.d, group = day_each, type = "centroid", bias.adjust = TRUE)
No. of Positive Eigenvalues: 32
No. of Negative Eigenvalues: 0
Average distance to centroid:
day10 day32 day60
0.7089 0.6904 0.5292
Eigenvalues for PCoA axes:
(Showing 8 of 32 eigenvalues)
PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
2.881 2.155 1.684 1.528 1.387 1.347 1.037 1.003
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.23186 0.115929 7.1717 9999 0.0024 **
Residuals 36 0.58193 0.016165
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
day10 day32 day60
day10 0.4429000 0.0065
day32 0.4385478 0.0163
day60 0.0068012 0.0161265
disper_DAY_centroid(data_name=uds_data_all_day2, treat="S", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Target treatment was
[1] "T"
Homogeneity of multivariate dispersions
Call: betadisper(d = uds_each.d, group = day_each, type = "centroid", bias.adjust = TRUE)
No. of Positive Eigenvalues: 35
No. of Negative Eigenvalues: 0
Average distance to centroid:
day10 day32 day60
0.6462 0.6452 0.6937
Eigenvalues for PCoA axes:
(Showing 8 of 35 eigenvalues)
PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
2.9445 2.0693 1.7712 1.4165 1.3933 1.1933 0.9916 0.8363
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.02078 0.010389 1.1044 9999 0.3316
Residuals 36 0.33865 0.009407
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
day10 day32 day60
day10 0.98260 0.2935
day32 0.98237 0.1310
day60 0.30142 0.13155
disper_DAY_centroid(data_name=uds_data_all_day2, treat="D", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Target treatment was
[1] "W"
Homogeneity of multivariate dispersions
Call: betadisper(d = uds_each.d, group = day_each, type = "centroid", bias.adjust = TRUE)
No. of Positive Eigenvalues: 38
No. of Negative Eigenvalues: 0
Average distance to centroid:
day10 day32 day60
0.6880 0.6697 0.6604
Eigenvalues for PCoA axes:
(Showing 8 of 38 eigenvalues)
PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
2.8219 2.1875 2.0501 1.6154 1.2362 1.0316 0.9482 0.8815
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.00529 0.0026448 0.3821 9999 0.6883
Residuals 40 0.27685 0.0069213
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
day10 day32 day60
day10 0.51180 0.3856
day32 0.50939 0.7932
day60 0.38878 0.78578
#Results from 2-3 times measurement in 10d, 32d, and 60d without leaf beetle through spatial median
#With the recommended measure (Chao index) (for Fig.2b)
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao")
Sampling day was
[1] "10d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.001891 0.0018907 0.6307 9999 0.4249
Residuals 34 0.101927 0.0029979
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.015571 0.0077853 1.599 9999 0.2128
Residuals 33 0.160670 0.0048688
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.17600 0.8996
S 0.17609 0.1955
U 0.90043 0.19124
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "32d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.002424 0.0024237 0.4006 9999 0.5329
Residuals 44 0.266220 0.0060504
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01282 0.0064087 0.7447 9999 0.4736
Residuals 43 0.37004 0.0086056
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.64410 0.4333
S 0.64677 0.2309
U 0.43762 0.22977
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "60d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.15801 0.158009 7.9115 9999 0.0079 **
Residuals 37 0.73897 0.019972
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15157 0.075784 3.6986 9999 0.0324 *
Residuals 36 0.73764 0.020490
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.434700 0.0852
S 0.432947 0.0289
U 0.084784 0.027949
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao", adj="holm")
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.2694 0.26937 0.55504 0.01606 0.9468
Residuals 34 16.5006 0.48531 0.98394
Total 35 16.7699 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 0.7546 0.37732 0.77748 0.045 0.8075
Residuals 33 16.0153 0.48531 0.955
Total 35 16.7699 1.000
The pairwise PERMANOVA (adonis) analysis:
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.4758 0.47578 1.0316 0.04905 0.3875
Residuals 20 9.2242 0.46121 0.95095
Total 21 9.7000 1.00000
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.4108 0.41082 0.80329 0.01793 0.6873
Residuals 44 22.5024 0.51142 0.98207
Total 45 22.9132 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.5462 0.77311 1.5558 0.06748 0.031 *
Residuals 43 21.3670 0.49691 0.93252
Total 45 22.9132 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 1.1096 1.10957 2.4132 0.07445 0.0052 **
Residuals 30 13.7936 0.45979 0.92555
Total 31 14.9031 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.8018 0.80177 1.8611 0.04789 0.0397 *
Residuals 37 15.9400 0.43081 0.95211
Total 38 16.7418 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.468 0.73398 1.73 0.08768 0.0218 *
Residuals 36 15.274 0.42427 0.91232
Total 38 16.742 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.6575 0.65750 1.4895 0.05418 0.1079
Residuals 26 11.4773 0.44144 0.94582
Total 27 12.1348 1.00000
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="N", c_LB="woLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ day_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
day_each 2 3.028 1.51401 3.3933 0.15862 1e-04 ***
Residuals 36 16.062 0.44617 0.84138
Total 38 19.090 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
Pairwise comparision by permutation, adjusted by[1] "holm"
The adjusted P-value for day10 vs day32 was
[1] 0.0019
The adjusted P-value for day10 vs day60 was
[1] 3e-04
The adjusted P-value for day32 vs day60 was
[1] 3e-04
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="T", c_LB="woLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ day_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
day_each 2 2.8635 1.43177 3.1913 0.15059 1e-04 ***
Residuals 36 16.1514 0.44865 0.84941
Total 38 19.0149 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
Pairwise comparision by permutation, adjusted by[1] "holm"
The adjusted P-value for day10 vs day32 was
[1] 3e-04
The adjusted P-value for day10 vs day60 was
[1] 0.0016
The adjusted P-value for day32 vs day60 was
[1] 3e-04
com_div_DAY_centroid(data_name=uds_data_all_day2, treat="W", c_LB="woLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ day_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
day_each 2 2.3439 1.17194 2.5569 0.11335 1e-04 ***
Residuals 40 18.3339 0.45835 0.88665
Total 42 20.6778 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
Pairwise comparision by permutation, adjusted by[1] "holm"
The adjusted P-value for day10 vs day32 was
[1] 9e-04
The adjusted P-value for day10 vs day60 was
[1] 9e-04
The adjusted P-value for day32 vs day60 was
[1] 0.0035
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="bray")
Sampling day was
[1] "10d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.006804 0.0068039 1.9548 9999 0.1736
Residuals 34 0.118340 0.0034806
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.017357 0.0086784 1.9915 9999 0.1473
Residuals 33 0.143806 0.0043578
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.19600 0.5771
S 0.19320 0.1024
U 0.56666 0.10418
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="bray")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "32d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.006626 0.0066260 1.0004 9999 0.3207
Residuals 44 0.291426 0.0066233
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.03321 0.0166063 1.9505 9999 0.1571
Residuals 43 0.36609 0.0085137
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.224300 0.4531
S 0.223032 0.0832
U 0.453786 0.080459
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="bray")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "60d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.15997 0.159968 6.6985 9999 0.0139 *
Residuals 37 0.88360 0.023881
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15006 0.075028 3.1041 9999 0.0558 .
Residuals 36 0.87015 0.024171
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.662300 0.0904
S 0.665325 0.0440
U 0.089952 0.043258
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="bray", adj="holm")
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.3591 0.35914 0.70982 0.02045 0.8327
Residuals 34 17.2023 0.50595 0.97955
Total 35 17.5615 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 0.8157 0.40785 0.80373 0.04645 0.81
Residuals 33 16.7458 0.50745 0.95355
Total 35 17.5615 1.00000
The pairwise PERMANOVA (adonis) analysis:
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.4289 0.42887 0.91093 0.04356 0.5466
Residuals 20 9.4161 0.47080 0.95644
Total 21 9.8449 1.00000
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="bray", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.4391 0.43915 0.83573 0.01864 0.6478
Residuals 44 23.1206 0.52547 0.98136
Total 45 23.5597 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.6115 0.80577 1.5786 0.0684 0.0244 *
Residuals 43 21.9482 0.51042 0.9316
Total 45 23.5597 1.0000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 1.1362 1.13619 2.4501 0.0755 0.0033 **
Residuals 30 13.9117 0.46372 0.9245
Total 31 15.0479 1.0000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="bray", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.8212 0.82119 1.7852 0.04603 0.0472 *
Residuals 37 17.0200 0.46000 0.95397
Total 38 17.8412 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.5345 0.76726 1.6939 0.08601 0.0255 *
Residuals 36 16.3066 0.45296 0.91399
Total 38 17.8412 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.7021 0.70206 1.5148 0.05506 0.0958 .
Residuals 26 12.0498 0.46345 0.94494
Total 27 12.7519 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="wLB", meth="chao")
Sampling day was
[1] "10d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00423 0.004225 0.0875 9999 0.7699
Residuals 42 2.02866 0.048301
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.52667 0.263337 6.4684 9999 0.0042 **
Residuals 41 1.66917 0.040712
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.0026000 0.1386
S 0.0014055 0.0463
U 0.1382685 0.0443417
disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="wLB", meth="chao")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "32d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04387 0.043871 1.8568 9999 0.176
Residuals 46 1.08683 0.023627
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.03582 0.017910 0.7201 9999 0.4965
Residuals 45 1.11917 0.024871
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.62540 0.5226
S 0.62872 0.2060
U 0.52614 0.20338
disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="wLB", meth="chao")
警告: some squared distances are negative and changed to zero
Sampling day was
[1] "60d"
Comparison between chambers
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.09148 0.091485 6.4648 9999 0.0169 *
Residuals 41 0.58020 0.014151
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Comparison between treatments
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.08619 0.043094 2.8724 9999 0.0737 .
Residuals 40 0.60010 0.015003
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
D S U
D 0.606000 0.1055
S 0.606209 0.0581
U 0.107893 0.050774
com_div_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="wLB", meth="chao", adj="holm")
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.1711 0.17113 0.45554 0.01073 0.9935
Residuals 42 15.7776 0.37566 0.98927
Total 43 15.9487 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 0.7672 0.38361 1.036 0.04811 0.3699
Residuals 41 15.1815 0.37028 0.95189
Total 43 15.9487 1.00000
The pairwise PERMANOVA (adonis) analysis:
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.5947 0.59473 1.7606 0.06342 0.0572 .
Residuals 26 8.7825 0.33779 0.93658
Total 27 9.3772 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
com_div_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="wLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.3734 0.37342 0.91069 0.01941 0.4693
Residuals 46 18.8616 0.41003 0.98059
Total 47 19.2350 1.00000
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.1599 0.57994 1.4438 0.0603 0.0904 .
Residuals 45 18.0751 0.40167 0.9397
Total 47 19.2350 1.0000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.7938 0.79381 1.9693 0.0616 0.0474 *
Residuals 30 12.0926 0.40309 0.9384
Total 31 12.8864 1.0000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
com_div_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="wLB", meth="chao", adj="holm")
警告: some squared distances are negative and changed to zero
Call:
adonis(formula = uds_median ~ chamber_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
chamber_each 1 0.699 0.69904 1.5518 0.03647 0.0877 .
Residuals 41 18.469 0.45047 0.96353
Total 42 19.168 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call:
adonis(formula = uds_median ~ treatment_each, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_each 2 1.3211 0.66053 1.4804 0.06892 0.0596 .
Residuals 40 17.8472 0.44618 0.93108
Total 42 19.1683 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The pairwise PERMANOVA (adonis) analysis:
警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero 警告: some squared distances are negative and changed to zero
The comparision between D and S is:
Call:
adonis(formula = DS_median ~ treatment_DS, permutations = 9999, method = "euclidean")
Permutation: free
Number of permutations: 9999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
treatment_DS 1 0.6081 0.60811 1.3234 0.04513 0.1711
Residuals 28 12.8664 0.45952 0.95487
Total 29 13.4745 1.00000
####Calculate/Visualize within day turnover (which is not used in the manuscript)
turnover_comparison_short <- function(data_name, date, c_LB, meth, output = NULL)
{
#data_name <- uds_data_all_day2
#date <- "10d"
#meth<-"chao"
#date <- "10d","32d", or "60d", or any
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_data2 <-subset(uds_data1, day==date)
uds_each<-uds_data2[,c(-1,-2,-3, -77,-78, -79,-80)]
} #When leaf beetle is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_data2 <-subset(uds_data1, day==date)
uds_each<-uds_data2[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
} #When leaf beetle is NOT included
sample_each <- subset(uds_data1, day==date)$sample
treatment_each <- subset(uds_data1, day==date)$treatment
chamber_each <- subset(uds_data1, day==date)$chamber
#View(chamber_each)
time_each <- subset(uds_data1, day==date)$time
#prepare dummy vectors
sample_n <- NULL
distance_n <-NULL
treatment_n <-NULL
chamber_n <-NULL
timei_n <-NULL
timej_n <-NULL
#calculate distance between 1P, 2P, and 3P
for(i in 1:(length(sample_each) - 1)) {
for(j in (i + 1): length(sample_each)) {
if(sample_each[i] == sample_each[j]) {
sample_n <- append(sample_n, as.character(sample_each[i]))
distance_n <- append(distance_n, as.numeric(vegdist(rbind(uds_each[i,], uds_each[j,])), method=meth))
treatment_n <-append(treatment_n, as.character(treatment_each[i]))
chamber_n <- append(chamber_n, as.character(chamber_each[i]))
timei_n <-append(timei_n, as.character(time_each[i]))
timej_n <-append(timej_n, as.character(time_each[j]))
}
}
}
turnover<- data.frame(sample=sample_n, distance=distance_n, treatment=treatment_n, chamber=chamber_n, timei=timei_n, timej=timej_n)
turnover2 <-subset(turnover, !(timei == "1P" & timej=="3P")) #exclude the combination of 1P & 3P (redundant)
#turnoverN <- subset(turnover2, treatment=="N") #separate data by treatment
turnoverT <- subset(turnover2, treatment=="T") #separate data by treatment
turnoverW <- subset(turnover2, treatment=="W") #separate data by treatment
turnoverC1 <- subset(turnover2, chamber=="C1") #separate data by chamber
turnoverC2 <- subset(turnover2, chamber=="C2") #separate data by chamber
#calculate average of 1P vs 2P & 2P vs 3P if any
#average_turnoverN<-by(as.numeric(turnoverN$distance), turnoverN$sample, mean)
average_turnoverT<-by(as.numeric(turnoverT$distance), turnoverT$sample, mean)
average_turnoverW<-by(as.numeric(turnoverW$distance), turnoverW$sample, mean)
average_turnoverC1<-by(as.numeric(turnoverC1$distance), turnoverC1$sample, mean)
average_turnoverC2<-by(as.numeric(turnoverC2$distance), turnoverC2$sample, mean)
#average_N<-subset(average_turnoverN, average_turnoverN>=0)
average_T<-subset(average_turnoverT, average_turnoverT>=0)
average_W<-subset(average_turnoverW, average_turnoverW>=0)
average_C1<-subset(average_turnoverC1, average_turnoverC1>=0)
average_C2<-subset(average_turnoverC2, average_turnoverC2>=0)
length(average_T)
#trN<-NULL
trT<-NULL
trW<-NULL
trC1<-NULL
trC2<-NULL
#for(i in 1:length(average_N)) trN <-append(trN, "N")
for(i in 1:length(average_T)) trT <-append(trT, "E")
for(i in 1:length(average_W)) trW <-append(trW, "D")
for(i in 1:length(average_C1)) trC1 <-append(trC1, "C1")
for(i in 1:length(average_C2)) trC2 <-append(trC2, "C2")
#av_N<-data.frame(treatment=trN, distance=average_N)
av_T<-data.frame(treatment=trT, distance=average_T)
av_W<-data.frame(treatment=trW, distance=average_W)
av_C1<-data.frame(chamber=trC1, distance=average_C1)
av_C2<-data.frame(chamber=trC2, distance=average_C2)
#summary dataframe
turnover_short <- rbind.data.frame(av_T, av_W)
turnover_shortC <- rbind.data.frame(av_C1, av_C2)
#print(turnover_shortC$distance)
#View(turnover_shortC)
#View(turnover_short)
#print(distance_n)
g_box<-ggplot(turnover_shortC, aes(y=distance, x=chamber))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(width=0.1)
g_box<-g_box + coord_fixed(ratio = 5/1) + coord_cartesian(ylim = c(0, 1))
print(g_box)
g_box2<-ggplot(turnover_short, aes(y=distance, x=treatment))
g_box2<-g_box2 + geom_violin()
g_box2<-g_box2 + geom_boxplot(width=0.1)
g_box2<-g_box2 + coord_fixed(ratio = 5/1) + coord_cartesian(ylim = c(0, 1))
print(g_box2)
#plotmeans(turnover_shortC$distance~turnover_shortC$chamber)
#plotmeans(turnover_short$distance~turnover_short$treatment)
print(bartlett.test(turnover_shortC$distance~turnover_shortC$chamber))
#One-way ANOVA with unequal-variance
print(oneway.test(turnover_shortC$distance~turnover_shortC$chamber,var.equal=FALSE))
print(bartlett.test(turnover_short$distance~turnover_short$treatment))
#One-way ANOVA with unequal-variance
print(oneway.test(turnover_short$distance~turnover_short$treatment,var.equal=FALSE))
#Pairwise comparison with Welch's t test
#print(pairwise.t.test(turnover_short$distance, turnover_short$treatment, p.adj="holm"))
#if(output==TRUE) return(turnover_short)
}
turnover_comparison_long <- function(data_name, from_date, c_LB, meth="chao", output = NULL)
{
if(c_LB == "wLB") {
uds_data1 <-subset(data_name, abundance > 0) #delete the sampling data with 0 total abundance when leaf beetle
uds_each<-uds_data1[,c(-1,-2,-3, -77,-78,-79,-80)]
} #When leaf beetle is included
if(c_LB == "woLB") {
uds_data1 <-subset(data_name, abundance2 > 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_each<-uds_data1[,c(-1,-2,-3,-4,-77,-78,-79,-80)]
} #When leaf beetle is NOT included
sample_each <- uds_data1$sample
chamber_each <- uds_data1$chamber
treatment_each <- uds_data1$treatment
time_each <- uds_data1$time
uds_time.d<-vegdist(uds_each, method=meth)
uds_dist_time <-betadisper(uds_time.d, sample_each, type="median", bias.adjust = TRUE)
uds_day.median <- data.frame(uds_dist_time$centroids) #calculate the centroid (spatial median) for repeated measures
#Replace the sample name into the treatment & day
sample0 <- row.names(uds_day.median)
sample1 <- sample0
sample2 <- sample0
sample3 <- sample0
sample0<-gsub("1n", "n", sample0)
sample0<-gsub("2n", "n", sample0)
sample0<-gsub("3n", "n", sample0)
sample0<-gsub("1t", "t", sample0)
sample0<-gsub("2t", "t", sample0)
sample0<-gsub("3t", "t", sample0)
sample0<-gsub("1w", "w", sample0)
sample0<-gsub("2w", "w", sample0)
sample0<-gsub("3w", "w", sample0)
sample1[grep("n", sample1)]<- "C1"
sample1[grep("t", sample1)]<- "C2"
sample1[grep("w", sample1)]<- "C2"
sample2[grep("n", sample2)]<- "U"
sample2[grep("t", sample2)]<- "E"
sample2[grep("w", sample2)]<- "D"
sample3[grep("1n", sample3)] <-"10d"
sample3[grep("2n", sample3)] <-"32d"
sample3[grep("3n", sample3)] <-"60d"
sample3[grep("1t", sample3)] <-"10d"
sample3[grep("2t", sample3)] <-"32d"
sample3[grep("3t", sample3)] <-"60d"
sample3[grep("1w", sample3)] <-"10d"
sample3[grep("2w", sample3)] <-"32d"
sample3[grep("3w", sample3)] <-"60d"
sample_each <- as.factor(sample0)
chamber_each <- as.factor(sample1)
treatment_each <- as.factor(sample2) #just conversion
day_each <- as.factor(sample3) #just conversion
uds_day.median2<-cbind(chamber=chamber_each, uds_day.median)
uds_day.median2<-cbind(treatment=treatment_each, uds_day.median2)
uds_day.median2<-cbind(day=day_each, uds_day.median2)
uds_day.median2<-cbind(sample_each, uds_day.median2)
uds_day.median2[1,c(-1,-2,-3)]
uds_day.median2$sample_each[2]
#View(uds_day.median2)
#prepare dummy vectors
sample_n <- NULL
distance_n <-NULL
chamber_n <-NULL
treatment_n <-NULL
dayi_n <-NULL
dayj_n <-NULL
#calculate distance between 10d, 32d, and 60d
for(i in 1:(length(sample_each) - 1)) {
for(j in (i + 1): length(sample_each)) {
if(sample_each[i] == sample_each[j]) {
sample_n <- append(sample_n, as.character(sample_each[i]))
distance_n <- append(distance_n, as.numeric(vegdist(rbind.data.frame(uds_day.median2[i,c(-1,-2,-3,-4)], uds_day.median2[j,c(-1,-2,-3,-4)]), method="euclidean")))
chamber_n <-append(chamber_n, as.character(chamber_each[i]))
treatment_n <-append(treatment_n, as.character(treatment_each[i]))
dayi_n <-append(dayi_n, as.character(day_each[i]))
dayj_n <-append(dayj_n, as.character(day_each[j]))
}
}
}
turnover<- data.frame(sample=sample_n, distance=distance_n, chamber=chamber_n, treatment=treatment_n, dayi=dayi_n, dayj=dayj_n)
if(from_date == "10d") turnover_long <-subset(turnover, (dayi == "10d" & dayj=="32d"))
#turnover_10d60d <-subset(turnover, (dayi == "10d" & dayj=="60d")) #not used excluding the combination of 10d & 60d (redundanct)
if(from_date == "32d") turnover_long <-subset(turnover, (dayi == "32d" & dayj=="60d"))
turnover_long_tr <- subset(turnover_long, chamber == "C2")
View(turnover_long)
View(turnover_long_tr)
g_box<-ggplot(turnover_long, aes(y=distance, x=chamber))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(width=0.1)
g_box<-g_box + coord_fixed(ratio = 5/1)
print(g_box)
g_box2<-ggplot(turnover_long_tr, aes(y=distance, x=treatment))
g_box2<-g_box2 + geom_violin()
g_box2<-g_box2 + geom_boxplot(width=0.1)
g_box2<-g_box2 + coord_fixed(ratio = 5/1)
print(g_box2)
#plotmeans(turnover_long$distance~turnover_long$chamber)
#plotmeans(turnover_long_tr$distance~turnover_long_tr$treatment)
print(bartlett.test(turnover_long$distance~turnover_long$chamber))
print(bartlett.test(turnover_long_tr$distance~turnover_long_tr$treatment))
#One-way ANOVA with unequal-variance
print(oneway.test(turnover_long$distance~turnover_long$chamber,var.equal=FALSE))
print(oneway.test(turnover_long_tr$distance~turnover_long_tr$treatment,var.equal=FALSE))
}
####Calculate/Visualize between dates turnover (not used)
#(No difference at all @ 10 day between time) ####For Fig. S1
#turnover_comparison_short(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao") #error due to lack of replication within day
turnover_comparison_short(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao")
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Bartlett test of homogeneity of variances
data: turnover_shortC$distance by turnover_shortC$chamber
Bartlett's K-squared = 0.59528, df = 1, p-value = 0.4404
One-way analysis of means (not assuming equal variances)
data: turnover_shortC$distance and turnover_shortC$chamber
F = 0.77285, num df = 1.000, denom df = 21.334, p-value = 0.3891
Bartlett test of homogeneity of variances
data: turnover_short$distance by turnover_short$treatment
Bartlett's K-squared = 1.1809, df = 1, p-value = 0.2772
One-way analysis of means (not assuming equal variances)
data: turnover_short$distance and turnover_short$treatment
F = 0.34211, num df = 1.000, denom df = 27.826, p-value = 0.5633
turnover_comparison_short(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao")
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Bartlett test of homogeneity of variances
data: turnover_shortC$distance by turnover_shortC$chamber
Bartlett's K-squared = 1.396, df = 1, p-value = 0.2374
One-way analysis of means (not assuming equal variances)
data: turnover_shortC$distance and turnover_shortC$chamber
F = 2.5418, num df = 1.000, denom df = 13.773, p-value = 0.1335
Bartlett test of homogeneity of variances
data: turnover_short$distance by turnover_short$treatment
Bartlett's K-squared = 0.38892, df = 1, p-value = 0.5329
One-way analysis of means (not assuming equal variances)
data: turnover_short$distance and turnover_short$treatment
F = 0.024992, num df = 1.000, denom df = 15.321, p-value = 0.8765
#(No difference at all @ 10 day between time) ####For Fig. S1
turnover_comparison_long(data_name=uds_data_all_day2, from_date="10d", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Bartlett test of homogeneity of variances
data: turnover_long$distance by turnover_long$chamber
Bartlett's K-squared = 3.6943, df = 1, p-value = 0.0546
Bartlett test of homogeneity of variances
data: turnover_long_tr$distance by turnover_long_tr$treatment
Bartlett's K-squared = 3.5787, df = 1, p-value = 0.05852
One-way analysis of means (not assuming equal variances)
data: turnover_long$distance and turnover_long$chamber
F = 0.46726, num df = 1.000, denom df = 15.545, p-value = 0.5043
One-way analysis of means (not assuming equal variances)
data: turnover_long_tr$distance and turnover_long_tr$treatment
F = 0.040366, num df = 1.000, denom df = 11.225, p-value = 0.8444
turnover_comparison_long(data_name=uds_data_all_day2, from_date="32d", c_LB="woLB", meth="chao")
警告: some squared distances are negative and changed to zero
Bartlett test of homogeneity of variances
data: turnover_long$distance by turnover_long$chamber
Bartlett's K-squared = 4.5976, df = 1, p-value = 0.03202
Bartlett test of homogeneity of variances
data: turnover_long_tr$distance by turnover_long_tr$treatment
Bartlett's K-squared = 0.0072188, df = 1, p-value = 0.9323
One-way analysis of means (not assuming equal variances)
data: turnover_long$distance and turnover_long$chamber
F = 0.15229, num df = 1.00, denom df = 28.73, p-value = 0.6992
One-way analysis of means (not assuming equal variances)
data: turnover_long_tr$distance and turnover_long_tr$treatment
F = 1.0881, num df = 1.000, denom df = 25.985, p-value = 0.3065
LB_diff_C1_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="LBonly", meth="euclidean",dist_out="C12")
LB_diff_C1_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="LBonly", meth="euclidean",dist_out="C12")
LB_diff_C1_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C1")
LB_diff_C2_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C2")
LB_diff_C12_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="LBonly", meth="euclidean",dist_out="C12")
comm_diff_C1_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C1")
comm_diff_C2_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C2")
comm_diff_C12_d10.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="10d", c_LB="woLB", meth="chao",dist_out="C12")
comm_diff_C1_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C1")
警告: some squared distances are negative and changed to zero
comm_diff_C2_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C2")
警告: some squared distances are negative and changed to zero
comm_diff_C12_d32.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="32d", c_LB="woLB", meth="chao",dist_out="C12")
警告: some squared distances are negative and changed to zero
comm_diff_C1_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C1")
警告: some squared distances are negative and changed to zero
comm_diff_C2_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C2")
警告: some squared distances are negative and changed to zero
comm_diff_C12_d60.d<-disper_uds_centroid(data_name=uds_data_all_day2, date="60d", c_LB="woLB", meth="chao",dist_out="C12")
警告: some squared distances are negative and changed to zero
Only for day 32
plot(LB_diff_C1_d10.d, comm_diff_C1_d10.d)
mantel(LB_diff_C1_d10.d, comm_diff_C1_d10.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C1_d10.d, ydis = comm_diff_C1_d10.d, method = "pearson", permutations = 9999)
Mantel statistic r: -0.1332
Significance: 0.8255
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.141 0.168 0.190 0.214
Permutation: free
Number of permutations: 9999
plot(LB_diff_C2_d10.d, comm_diff_C2_d10.d)
mantel(LB_diff_C2_d10.d, comm_diff_C2_d10.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C2_d10.d, ydis = comm_diff_C2_d10.d, method = "pearson", permutations = 9999)
Mantel statistic r: -0.07778
Significance: 0.8019
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.114 0.150 0.178 0.209
Permutation: free
Number of permutations: 9999
plot(LB_diff_C12_d10.d, comm_diff_C12_d10.d)
mantel(LB_diff_C12_d10.d, comm_diff_C12_d10.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C12_d10.d, ydis = comm_diff_C12_d10.d, method = "pearson", permutations = 9999)
Mantel statistic r: -0.09194
Significance: 0.9521
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0771 0.0999 0.1190 0.1402
Permutation: free
Number of permutations: 9999
plot(LB_diff_C1_d32.d, comm_diff_C1_d32.d)
mantel(LB_diff_C1_d32.d, comm_diff_C1_d32.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C1_d32.d, ydis = comm_diff_C1_d32.d, method = "pearson", permutations = 9999)
Mantel statistic r: -0.01861
Significance: 0.5386
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.154 0.199 0.239 0.285
Permutation: free
Number of permutations: 9999
plot(LB_diff_C2_d32.d, comm_diff_C2_d32.d)
mantel(LB_diff_C2_d32.d, comm_diff_C2_d32.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C2_d32.d, ydis = comm_diff_C2_d32.d, method = "pearson", permutations = 9999)
Mantel statistic r: 0.1698
Significance: 0.0151
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.117 0.139 0.158 0.179
Permutation: free
Number of permutations: 9999
plot(LB_diff_C12_d32.d, comm_diff_C12_d32.d)
mantel(LB_diff_C12_d32.d, comm_diff_C12_d32.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C12_d32.d, ydis = comm_diff_C12_d32.d, method = "pearson", permutations = 9999)
Mantel statistic r: 0.128
Significance: 0.026
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0936 0.1136 0.1288 0.1465
Permutation: free
Number of permutations: 9999
plot(LB_diff_C1_d60.d, comm_diff_C1_d60.d)
#mantel(LB_diff_C1_d60.d, comm_diff_C1_d60.d, method="pearson", permutations=9999)
plot(LB_diff_C2_d60.d, comm_diff_C2_d60.d)
mantel(LB_diff_C2_d60.d, comm_diff_C2_d60.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C2_d60.d, ydis = comm_diff_C2_d60.d, method = "pearson", permutations = 9999)
Mantel statistic r: -0.007227
Significance: 0.5256
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.149 0.183 0.206 0.232
Permutation: free
Number of permutations: 9999
plot(LB_diff_C12_d60.d, comm_diff_C12_d60.d)
mantel(LB_diff_C12_d60.d, comm_diff_C12_d60.d, method="pearson", permutations=9999)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = LB_diff_C12_d60.d, ydis = comm_diff_C12_d60.d, method = "pearson", permutations = 9999)
Mantel statistic r: 0.02801
Significance: 0.4126
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.135 0.161 0.180 0.198
Permutation: free
Number of permutations: 9999
library(reshape)
#C1 only
m_LB <- as.matrix(LB_diff_C1_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C1_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)
g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C1 only")
g_points<-g_points+geom_point(size=2,alpha=0.3)
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
summary(result)
Call:
lm(formula = m2_LB_com$distance_com ~ m2_LB_com$distance_LB)
Residuals:
Min 1Q Median 3Q Max
-1.01101 -0.02574 0.06248 0.10782 0.19192
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.015442 0.033019 30.754 <2e-16 ***
m2_LB_com$distance_LB -0.004041 0.023008 -0.176 0.861
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1831 on 89 degrees of freedom
Multiple R-squared: 0.0003464, Adjusted R-squared: -0.01089
F-statistic: 0.03084 on 1 and 89 DF, p-value: 0.861
result$coefficients[1]
(Intercept)
1.015442
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=2,colour="blue") # slope
print(g_points)
#C2 only
m_LB <- as.matrix(LB_diff_C2_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C2_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)
g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C2 only")
g_points<-g_points+geom_point(size=2,alpha=0.3)
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
summary(result)
Call:
lm(formula = m2_LB_com$distance_com ~ m2_LB_com$distance_LB)
Residuals:
Min 1Q Median 3Q Max
-0.98366 -0.06532 0.03618 0.10218 0.32866
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.971440 0.008921 108.892 < 2e-16 ***
m2_LB_com$distance_LB 0.013777 0.003598 3.829 0.000145 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1607 on 494 degrees of freedom
Multiple R-squared: 0.02882, Adjusted R-squared: 0.02685
F-statistic: 14.66 on 1 and 494 DF, p-value: 0.0001454
result$coefficients[1]
(Intercept)
0.9714401
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=1,colour="blue") # slope
print(g_points)
#C1 + C2
m_LB <- as.matrix(LB_diff_C12_d32.d)
m2_LB <- melt(m_LB)[melt(upper.tri(m_LB))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_LB) <- c("t1", "t2", "distance_LB")
#View(m2_LB)
m_com <- as.matrix(comm_diff_C12_d32.d)
m2_com <- melt(m_com)[melt(upper.tri(m_com))$value,]
警告: 'as.is' should be specified by the caller; using TRUE 警告: 'as.is' should be specified by the caller; using TRUE
names(m2_com) <- c("tt1", "tt2", "distance_com")
m2_LB_com <- cbind.data.frame(m2_LB, m2_com)
#View(m2_LB_com)
g_points<-ggplot(m2_LB_com, aes(y=distance_com, x=distance_LB)) + xlim(0,8) + ylim(0,1.5) + coord_fixed(8/1.5) + ggtitle("C1 + C2")
g_points<-g_points+geom_point(size=2,alpha=0.3)
result<-lm(m2_LB_com$distance_com~m2_LB_com$distance_LB)
result$coefficients[1]
(Intercept)
0.9765874
g_points<-g_points+geom_abline(slope = result$coefficients[2],intercept = result$coefficients[1],size=1,linetype=1,colour="blue") # slope
print(g_points)
alpha_hill_number <- function(data_name, date, q)
{
#data_name <-uds_data_all_day2
#date <- "32d"
#Including leaf_beetle
uds_data1 <-subset(data_name, abundance2 >= 0) #OR, delete the sampling when species other than leaf beetle does not exist
uds_data2 <-subset(uds_data1, day==date)
#if excluding leaf beetle
#uds_each<-uds_data2[,c(-1,-2,-3,-4,-77,-78,-79)]
#OR inclcding leaf beetle
uds_each<-uds_data2[,c(-1,-2,-3,-77,-78,-79,-80)]
#Hill number q = 0, 1, 2
uds_hill<-renyi(uds_each, scales=c(0,1,2), hill=TRUE)
uds_hill2<-cbind(uds_hill, uds_data2[,c(77,78,79,80)])
colnames(uds_hill2) <-c("H0", "H1", "H2", "day", "treatment", "time", "chamber")
#convert NA to 0
uds_hill2$H0[is.na(uds_hill2$H0)] <- 0
uds_hill3 <- subset(uds_hill2,chamber=="C2")
uds_hill3$treatment[uds_hill3$treatment=="T"] <- "E"
uds_hill3$treatment[uds_hill3$treatment=="W"] <- "D"
#View(uds_hill2)
if(q == 0) {
g_box<-ggplot(uds_hill2, aes(y=H0, x=chamber))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(width=0.1)
g_box<-g_box + coord_fixed(ratio = 1/2)
print(g_box)
g_box2<-ggplot(uds_hill3, aes(y=H0, x=treatment))
g_box2<-g_box2 + geom_violin()
g_box2<-g_box2 + geom_boxplot(width=0.1)
g_box2<-g_box2 + coord_fixed(ratio = 1/2)
print(g_box2)
plotmeans(uds_hill2$H0~uds_hill2$chamber)
View(uds_hill2)
#homogeneity of variance
print(bartlett.test(uds_hill2$H0~uds_hill2$chamber))
print(bartlett.test(uds_hill3$H0~uds_hill3$treatment))
#One-way ANOVA with unequal-variance between chambers
print(oneway.test(uds_hill2$H0~uds_hill2$chamber,var.equal=FALSE))
print(oneway.test(uds_hill3$H0~uds_hill3$treatment,var.equal=FALSE))
#Pairwise comparison with Welch's t test
#print(pairwise.t.test(uds_hill2$H0, uds_hill2$treatment, p.adj="holm"))
}
if(q == 1) {
g_box<-ggplot(uds_hill2, aes(y=H1, x=chamber))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(width=0.1)
g_box<-g_box + coord_fixed(ratio = 1/2)
print(g_box)
g_box2<-ggplot(uds_hill3, aes(y=H1, x=treatment))
g_box2<-g_box2 + geom_violin()
g_box2<-g_box2 + geom_boxplot(width=0.1)
g_box2<-g_box2 + coord_fixed(ratio = 1/2)
print(g_box2)
plotmeans(uds_hill2$H1~uds_hill2$chamber)
View(uds_hill2)
#homogeneity of variance
print(bartlett.test(uds_hill2$H1~uds_hill2$chamber))
print(bartlett.test(uds_hill3$H1~uds_hill3$treatment))
#One-way ANOVA with unequal-variance between chambers
print(oneway.test(uds_hill2$H1~uds_hill2$chamber,var.equal=FALSE))
print(oneway.test(uds_hill3$H1~uds_hill3$treatment,var.equal=FALSE))
}
if(q == 2) {
g_box<-ggplot(uds_hill2, aes(y=H2, x=chamber))
g_box<-g_box + geom_violin()
g_box<-g_box + geom_boxplot(width=0.1)
g_box<-g_box + coord_fixed(ratio = 1/2)
print(g_box)
g_box2<-ggplot(uds_hill3, aes(y=H2, x=treatment))
g_box2<-g_box2 + geom_violin()
g_box2<-g_box2 + geom_boxplot(width=0.1)
g_box2<-g_box2 + coord_fixed(ratio = 1/2)
print(g_box2)
plotmeans(uds_hill2$H1~uds_hill2$chamber)
View(uds_hill2)
#homogeneity of variance
print(bartlett.test(uds_hill2$H2~uds_hill2$chamber))
print(bartlett.test(uds_hill3$H2~uds_hill3$treatment))
#One-way ANOVA with unequal-variance between chambers
print(oneway.test(uds_hill2$H2~uds_hill2$chamber,var.equal=FALSE))
print(oneway.test(uds_hill3$H2~uds_hill3$treatment,var.equal=FALSE))
}
}
#all samplings are treated as independent, used in the manuscript
####For Figure S2
alpha_hill_number(uds_data_all_day2, date="10d", q=0)
Bartlett test of homogeneity of variances
data: uds_hill2$H0 by uds_hill2$chamber
Bartlett's K-squared = 1.3053, df = 1, p-value = 0.2532
Bartlett test of homogeneity of variances
data: uds_hill3$H0 by uds_hill3$treatment
Bartlett's K-squared = 0.095202, df = 1, p-value = 0.7577
One-way analysis of means (not assuming equal variances)
data: uds_hill2$H0 and uds_hill2$chamber
F = 0.74438, num df = 1.000, denom df = 38.077, p-value = 0.3937
One-way analysis of means (not assuming equal variances)
data: uds_hill3$H0 and uds_hill3$treatment
F = 0.57447, num df = 1.000, denom df = 29.805, p-value = 0.4544
alpha_hill_number(uds_data_all_day2, date="32d", q=0)
Bartlett test of homogeneity of variances
data: uds_hill2$H0 by uds_hill2$chamber
Bartlett's K-squared = 0.81835, df = 1, p-value = 0.3657
Bartlett test of homogeneity of variances
data: uds_hill3$H0 by uds_hill3$treatment
Bartlett's K-squared = 0.49807, df = 1, p-value = 0.4803
One-way analysis of means (not assuming equal variances)
data: uds_hill2$H0 and uds_hill2$chamber
F = 6.7446, num df = 1.000, denom df = 85.245, p-value = 0.01107
One-way analysis of means (not assuming equal variances)
data: uds_hill3$H0 and uds_hill3$treatment
F = 0.032126, num df = 1.000, denom df = 93.009, p-value = 0.8581
alpha_hill_number(uds_data_all_day2, date="60d", q=0)
Bartlett test of homogeneity of variances
data: uds_hill2$H0 by uds_hill2$chamber
Bartlett's K-squared = 8.5728, df = 1, p-value = 0.003412
Bartlett test of homogeneity of variances
data: uds_hill3$H0 by uds_hill3$treatment
Bartlett's K-squared = 1.1076, df = 1, p-value = 0.2926
One-way analysis of means (not assuming equal variances)
data: uds_hill2$H0 and uds_hill2$chamber
F = 9.2244, num df = 1.00, denom df = 129.39, p-value = 0.002889
One-way analysis of means (not assuming equal variances)
data: uds_hill3$H0 and uds_hill3$treatment
F = 0.21631, num df = 1.000, denom df = 91.839, p-value = 0.643
For standardized richness (Results for Fig.3a)
incidence_whole <- read.csv("monthly-for-R_iNEXT.csv", header=T) #for incidence data calculated by Excel
For asymptotic richness (Results for Fig.3b)
####This is code for graphics in the manuscript (Fig.3)
stat_test<-DataInfo(incidence_whole, datatype = "incidence") ##This is necessary to find the minimum value of the coverage
estimateD(incidence_whole, datatype="incidence_freq", base="coverage", level=min(stat_test$SC))
The above iNEXT results were manually summarized in excel and
csv.
#When Short-term repeated measurement is assumed to be independent
ChaoRichness(incidence_whole, datatype="incidence") #It is OK to use dataframe, this is from iNEXT
diversity_summary <- read.csv("iNEXT_summary.csv",header=T)
class(diversity_summary)
[1] "data.frame"
View(diversity_summary)
diversity_summary2<-diversity_summary[1:20,]
class(diversity_summary2$day)
[1] "character"
#colours <- c("#F8766D", "#EE55B7", "#00BA38", "#619CFF", "#E76BF3")
g_div_a <- ggplot(data=diversity_summary2, aes(x=ID,y=SR))
g_div_a <- g_div_a + geom_errorbar(aes(ymax = Hconf, ymin = Lconf, colour=treatment), width = 0.0,size=5, alpha=0.5) #+ scale_colour_manual("Group", values = colours)
g_div_a <- g_div_a + geom_point(aes(colour=treatment), size=3) #+ scale_colour_manual("Group", values = colours)
g_div_a <- g_div_a + geom_point(aes(y=Observed), size=3) + coord_fixed(ratio = 1/5) #+ scale_colour_manual("Group", values = colours)
print(g_div_a)
Observed richness (Fig.3a)
Standardized richness (Fig.3a)
Asymptotic richness (Fig.3b)
#Plot Asymptotic richness (Chao index)(Fig3b)
diversity_summary3<-diversity_summary[16:20,]
diversity_summary3$ID <- c("1U", "2ED", "3E", "4D", "5UED")
#colours <- c("#F8766D", "#00BA38", "#619CFF", "#E76BF3", "#00B6EB", "#A58AFF", "#FB61D7")
g_div_b <- ggplot(data=diversity_summary3, aes(x=ID,y=AsyEST))
g_div_b <- g_div_b + geom_errorbar(aes(ymax = Hconf_asy, ymin = Lconf_asy, colour=treatment), width = 0.0,size=5, alpha=1.0) #+ scale_colour_manual("Group", values = colours)
g_div_b <- g_div_b + geom_point(size=5)
print(g_div_b)