###Fig.1 and 5#########habitat map tiff(filename="habitat1016a.tiff",width=16,height=24,units="cm",res=600) par(mfrow=c(2,1)) tiff(filename="habitat711.tiff",width=16,height=12,units="cm",res=600) m <- ggplot()+ geom_raster(data = env.dat, aes(x = x1+10, y = y1+10, fill = as.factor(habitat2))) + geom_hline(yintercept = sort(unique(env.dat$y1)), color = "grey") + geom_vline(xintercept = sort(unique(env.dat$x1)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (400 m)") + ylab("Northing (300 m)") + scale_fill_discrete(name = "Habitat", labels = c("A", "B")) m + ggtitle(paste("Actual Habitat Distribution of Kanas Plot")) dev.off() tiff(filename="habitat4.tiff",width=16,height=12,units="cm",res=600) ##### m4 <- ggplot()+ geom_raster(data = env.dat, aes(x = x1+10, y = y1+10, fill = as.factor(G2))) + geom_hline(yintercept = sort(unique(env.dat$y1)), color = "grey") + geom_vline(xintercept = sort(unique(env.dat$x1)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (400 m)") + ylab("Northing (300 m)") + scale_fill_discrete(name = "Habitat", labels = c("A1","A2", "B1","B2")) m4 + ggtitle(paste("Habitat Distribution of Kanas Plot")) dev.off() ###Fig.2 and 5######### #torus translation #Three kinds of changes ### 1. X Mirror 左右颠倒 newx <- -(env.dat$x1 - max(env.dat$x1)) newy <- env.dat$gy mirror_pattern <- cbind(env.dat, newx = newx, newy = newy) ggplot()+ geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + geom_hline(yintercept = sort(unique(mirror_pattern$newy)), color = "grey") + geom_vline(xintercept = sort(unique(mirror_pattern$newx)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (400 m)") + ylab("Northing (300 m)") + scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C")) + ggtitle(paste("Torus translation", "X-Mirror")) + scale_colour_hue(c=45, l=80) ### 2. flip Y 上下颠倒 newx <- env.dat$x1 newy <- -(env.dat$y1 - max(env.dat$y1)) mirror_pattern <- cbind(env.dat, newx = newx, newy = newy) ggplot()+ geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + geom_hline(yintercept = sort(unique(mirror_pattern$newy)), color = "grey") + geom_vline(xintercept = sort(unique(mirror_pattern$newx)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (300 m)") + ylab("Northing (400 m)") + scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C")) + ggtitle(paste("Torus translation", "X-Mirror")) + scale_colour_hue(c=45, l=80) ### 3. 180° rotate 180°旋转,左右和上下同时颠倒 newx <- -(env.dat$x1 - max(env.dat$x1)) newy <- -(env.dat$y1 - max(env.dat$y1)) mirror_pattern <- cbind(env.dat, newx = newx, newy = newy) ggplot()+ geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + geom_hline(yintercept = sort(unique(mirror_pattern$newy)), color = "grey") + geom_vline(xintercept = sort(unique(mirror_pattern$newx)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (300 m)") + ylab("Northing (400 m)") + scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C")) + ggtitle(paste("Torus translation", "Rotate 180")) + scale_colour_hue(c=45, l=80) ###torus 核心代码 ### 变换样方原点的坐标, 如超出样地横纵坐标x或者y的最大值, 就减样地样地x或者y的最大值。 increment <- unique(env.dat[2:nrow(env.dat), c("x1", "y1")]) nrow(increment) ## [1] 399 res <- list() for (i in 1:nrow(increment)){ newx <- env.dat$x1 + as.numeric(increment[i,1]) newy <- env.dat$y1 + as.numeric(increment[i,2]) for (j in 1:length(newx)){ if(newx[j] > max(env.dat$x1)){ newx[j] <- newx[j] - max(env.dat$x1) } } for(k in 1:length(newy)){ if(newy[k] > max(env.dat$y1)){ newy[k] <- newy[k] - max(env.dat$y1) } } res[[i]] <- cbind(env.dat, newx = newx, newy = newy) } ### ### 这里从400种Torus Translation随机化中,随机抽出10个查看 for (i in sort(sample(1:length(res), 10))){ temp <- res[[i]] p <- ggplot()+ geom_raster(data = temp, aes(x = newx+10 , y = newy+10 , fill = as.factor(habitats))) + geom_hline(yintercept = sort(unique(temp$y1)), color = "grey") + geom_vline(xintercept = sort(unique(temp$x1)), color = "grey") + scale_y_discrete(expand = c(0,0)) + scale_x_discrete(expand = c(0,0)) + xlab("Easting (400 m)") + ylab("Northing (300 m)") + scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C")) + ggtitle(paste("Torus translation", i)) + scale_colour_hue(c=45, l=80) print(p) } library(fgeo) library(dplyr) library(ggplot2) # Test any number of species (output a list of matrices) tt_lst2a3 <- tt_test(sten214, sp = unique(sten214$sp),elevation.dat , plotdim = extract_plotdim(elevation.dat), gridsize = extract_gridsize(elevation.dat)) tt_lsta <- tt_test(sten14, sp = unique(sten14$sp),elevation.dat , plotdim = extract_plotdim(elevation.dat), gridsize = extract_gridsize(elevation.dat)) ##### ###Fig.3################### ########## A=env.dat<-read.csv("knsenv.csv",header=T) env.dat$Group <- factor(env.dat$Group) kruskal.test(slope~Group, data=env.dat) #pH library(multcomp) fit1 <- aov(log(pH+1)~Group,data = env.dat) summary(fit1) tuk1<-glht(fit1,linfct=mcp(Group="Tukey")) res1 <- cld(tuk1,alpah=0.05) library(dplyr) aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(pH)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(pH = res1$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) library(ggplot2) cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s1<-ggplot(env.dat,aes(Group,pH)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = pH), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("pH") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #SOM # fit2 <- aov(log(SOM)~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(SOM)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(SOM = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s2<-ggplot(env.dat,aes(Group,SOM)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = SOM), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("SOM") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #TN # fit2 <- aov(TN~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(TN)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(TN = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s2a<-ggplot(env.dat,aes(Group,TN)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = TN), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("TN") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #TP # fit2 <- aov(TP~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(TP)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(TP = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s3<-ggplot(env.dat,aes(Group,TP)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = TP), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("TP") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #AN # fit2 <- aov(AN~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(AN)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(AN = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s4<-ggplot(env.dat,aes(Group,AN)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = AN), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("AN") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #AP # fit2 <- aov(log(AP+1)~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(AP)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(AP = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s5<-ggplot(env.dat,aes(Group,AP)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = AP), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("AP") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #AK # fit2 <- aov(AK~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(AK)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(AK = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s6<-ggplot(env.dat,aes(Group,AK)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = AK), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("AK") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #elev # fit2 <- aov(log(elevation+1)~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(elevation)) aa1$Max <- aa1$Max + max(aa1$Max)*0.01 test <- data.frame(elevation = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s7<-ggplot(env.dat,aes(Group,elevation)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = elevation), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("elevation") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #slope # fit2 <- aov(log(slope+1)~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(slope)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(slope = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s8<-ggplot(env.dat,aes(Group,slope)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = slope), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("slope") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #aspect # fit2 <- aov(log(sina+1)~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(sina)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(sina = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s9b<-ggplot(env.dat,aes(Group,sina)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = sina), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("aspect_EW") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #convex # fit2 <- aov(convexity~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(convexity)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(convexity = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s10<-ggplot(env.dat,aes(Group,convexity)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = convexity), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("convexity") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") #canopy # fit2 <- aov(canopy~Group,data = env.dat) summary(fit2) tuk2<-glht(fit2,linfct=mcp(Group="Tukey")) res2 <- cld(tuk2,alpah=0.05) # aa1 <- env.dat%>%group_by(Group)%>% summarise(Max = max(canopy)) aa1$Max <- aa1$Max + max(aa1$Max)*0.1 test <- data.frame(canopy = res2$mcletters$Letters,aa1 = aa1$Max,Group = aa1$Group) test$Group <- factor(test$Group,levels = c("Hab_A","Hab_B")) # cbbPalette <- c("#B2182B","#56B4E9","#E69F00","#009E73","#F0E442","#0072B2", "#D55E00","#CC79A7","#CC6666","#9999CC","#66CC99","#99999", "#ADD1E5") s11<-ggplot(env.dat,aes(Group,canopy)) + geom_boxplot(aes(fill = Group)) + geom_text(data = test,aes(x = Group,y = aa1,label = canopy), size = 7,color = "black",fontface = "bold") + scale_fill_manual(values=cbbPalette) + ylab("canopy") + theme_bw()+ theme(axis.ticks.length = unit(0.4,"lines"), axis.ticks = element_line(color='black'), axis.line = element_line(colour = "black"), axis.title.x=element_blank(), axis.title.y=element_text(colour='black',size=20,face = "bold"), axis.text.y=element_text(colour='black',size=14), axis.text.x=element_text(colour='black',size=14), legend.position = "none") library(cowplot) s7<- cowplot::plot_grid(s1, s2, s3, s4,s5,s6,nrow = 2, labels = LETTERS[1:6])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("boxplot-soil12.jpeg", s7, width = 12, height = 8) library(cowplot) s12<- cowplot::plot_grid(s1, s2,s2a, s3, s4,s5,s6,s7,s8,s9b,s10,s11,nrow = 3, labels = LETTERS[1:12])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("boxplot-soil1212b722.jpeg", s12, width = 16, height = 12) ###Fig.4, table 1 and 2######### GLMMs and LMMs library(lmerTest) library(Matrix) detach("package:lme4") library(lme4) library(languageR) library(sjstats) library(car) ### lm.Ls <- glmer(Lari_sibi ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat,family = "poisson") summary(lm.Ls) ShowRegTable(lm.Ls) r.squaredLR(lm.Ls)#计算R2 ##### lm.Po<- glmer(Pice_obov ~PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat,family = "poisson") summary(lm.Po) lm.Po<- glmer(Pice_obov ~PCA1+PCA3+pH+(1|plot), data=env.dat,family = "poisson") summary(lm.Po) ShowRegTable(lm.Po) r.squaredLR(lm.Po)#计算R2 ##### lm.Ps <- glmer(Pinu_sibi ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat,family = "poisson") summary(lm.Ps) lm.Ps <- glmer(Pinu_sibi ~ PCA1+PCA2+pH+canopy+(1|plot), data=env.dat,family = "poisson") summary(lm.Ps) options(na.action = "na.fail") dredge(lm.Ps,beta="sd") ShowRegTable(lm.Ps) r.squaredLR(lm.Ps)#计算R2 # lm.Bp<- glmer(Betu_pend ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat,family = "poisson") summary(lm.Bp) lm.Bp<- glmer(Betu_pend ~ PCA1+(1|plot), data=env.dat,family = "poisson") summary(lm.Bp) ShowRegTable(lm.Bp) r.squaredLR(lm.Bp)#计算R2 ### lm.Lsba <- lmer(log(Lsba+1) ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat) summary(lm.Lsba) lm.Lsba <- lmer(log(Lsba+1)~ PCA1+(1|plot), data=env.dat) summary(lm.Lsba) options(na.action = "na.fail") dredge(lm.Lsba,beta="sd") ShowRegTable(lm.Lsba) r.squaredLR(lm.Lsba)#计算R2 # lm.Poba <- lmer(log(Poba+1) ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat) summary(lm.Poba) lm.Poba <- lmer(log(Poba+1) ~ PCA1+PCA3+canopy+pH+(1|plot), data=env.dat) summary(lm.Poba) options(na.action = "na.fail") dredge(lm.Poba,beta="sd") ShowRegTable(lm.Poba) r.squaredLR(lm.Poba)#计算R2 # lm.Psba <- lmer(log(Psba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat) summary(lm.Psba) lm.Psba1 <- glm(log(Psba+1)~ PCA1+pH, data=env.dat) summary(lm.Psba1) dredge(lm.Psba,beta="sd") ShowRegTable(lm.Psba1) r.squaredLR(lm.Psba1)#计算R2 # lm.Bpba <- lmer(log(Bpba+1) ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat) summary(lm.Bpba) lm.Bpba <- lmer(log(Bpba+1) ~ PCA1+canopy+(1|plot), data=env.dat) summary(lm.Bpba) dredge(lm.Bpba,beta="sd") ShowRegTable(lm.Bpba) r.squaredLR(lm.Bpba)#计算R2 ########### ##### lm.Lsa <- glmer(Lsa ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Lsa) lm.Lsa <- glmer(Lsa ~ PCA2+(1|plot), data=env.dat2,family = "poisson") summary(lm.Lsa) options(na.action = "na.fail") dredge(lm.Lsa,beta="sd") ShowRegTable(lm.Lsa) r.squaredLR(lm.Lsa)#计算R2 ##### 无显著因子 lm.Lsj <- glmer(Lsj ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Lsj) lm.Lsj <- glmer(Lsj ~ PCA1+(1|plot), data=env.dat2,family = "poisson") summary(lm.Lsj) options(na.action = "na.fail") dredge(lm.Lsa,beta="sd") ShowRegTable(lm.Lsa) r.squaredLR(lm.Lsa)#计算R2 ############### lm.Bpa <- glmer(Bpa ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Bpa) lm.Bpa <- glmer(Bpa ~ PCA1+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Bpa) options(na.action = "na.fail") dredge(lm.Bpa,beta="sd") ShowRegTable(lm.Bpa) r.squaredLR(lm.Bpa)#计算R2 ############### lm.Bpj <- glmer(Bpj ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Bpj) lm.Bpj <- glmer(Bpj ~ PCA1+(1|plot), data=env.dat2,family = "poisson") summary(lm.Bpj) options(na.action = "na.fail") dredge(lm.Bpj,beta="sd") ShowRegTable(lm.Bpj) r.squaredLR(lm.Bpj)#计算R2 ############### lm.Poa <- glmer(Poa ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Poa) lm.Poa <- glmer(Poa ~ PCA1+PCA3+canopy+pH+(1|plot), data=env.dat2,family = "poisson") summary(lm.Poa) options(na.action = "na.fail") dredge(lm.Poa,beta="sd") ShowRegTable(lm.Poa) r.squaredLR(lm.Poa)#计算R2 ############### lm.Poj <- glmer(Poj ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Poj) lm.Poj <- glmer(Poj ~ PCA1+PCA3+(1|plot), data=env.dat2,family = "poisson") summary(lm.Poj) options(na.action = "na.fail") dredge(lm.Poj,beta="sd") ShowRegTable(lm.Poj) r.squaredLR(lm.Poj)#计算R2 ############### lm.Pos <- glmer(Pos ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Pos) lm.Pos <- glmer(Pos ~ PCA1+PCA3+(1|plot), data=env.dat2,family = "poisson") summary(lm.Pos) options(na.action = "na.fail") dredge(lm.Pos,beta="sd") ShowRegTable(lm.Pos) r.squaredLR(lm.Pos)#计算R2 ############### lm.Psa <- glmer(Psa ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Psa) lm.Psa <- glmer(Psa ~ PCA1+PCA2+pH+(1|plot), data=env.dat2,family = "poisson") summary(lm.Psa) options(na.action = "na.fail") dredge(lm.Psa,beta="sd") ShowRegTable(lm.Psa) r.squaredLR(lm.Psa)#计算R2 ############### lm.Psj <- glmer(Psj ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Psj) lm.Psj <- glmer(Psj ~ PCA1+PCA2+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Psj) options(na.action = "na.fail") dredge(lm.Psj,beta="sd") ShowRegTable(lm.Psj) r.squaredLR(lm.Psj)#计算R2 ############### lm.Pss <- glmer(Pss ~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2,family = "poisson") summary(lm.Pss) lm.Pss <- glmer(Pss ~ PCA1+PCA2+pH+(1|plot), data=env.dat2,family = "poisson") summary(lm.Pss) options(na.action = "na.fail") dredge(lm.Pss,beta="sd") ShowRegTable(lm.Pos) r.squaredLR(lm.Pos)#计算R2 ############### ##### lm.Lsaba <- lmer(log(Lsaba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Lsaba) lm.Lsaba <- lmer(log(Lsaba+1)~ PCA1+(1|plot), data=env.dat2) summary(lm.Lsaba) options(na.action = "na.fail") dredge(lm.Lsaba,beta="sd") ShowRegTable(lm.Lsaba) r.squaredLR(lm.Lsaba)#计算R2 ########## ##### lm.Lsjba <- lmer(log(Lsjba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Lsjba) lm.Lsaba <- lmer(log(Lsjba+1)~ PCA1+(1|plot), data=env.dat2) summary(lm.Lsjba) options(na.action = "na.fail") dredge(lm.Lsjba,beta="sd") ShowRegTable(lm.Lsjba) r.squaredLR(lm.Lsjba)#计算R2 ########## ##### lm.Poaba <- lmer(log(Poaba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Poaba) lm.Poaba <- lmer(log(Poaba+1)~ PCA1+PCA3+canopy+(1|plot), data=env.dat2) summary(lm.Poaba) options(na.action = "na.fail") dredge(lm.Poaba,beta="sd") ShowRegTable(lm.Poaba) r.squaredLR(lm.Poaba)#计算R2 ########## ##### lm.Pojba <- lmer(log(Pojba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Pojba) lm.Pojba <- lmer(log(Pojba+1)~ PCA2+(1|plot), data=env.dat2) summary(lm.Pojba) options(na.action = "na.fail") dredge(lm.Pojba,beta="sd") ShowRegTable(lm.Pojba) r.squaredLR(lm.Pojba)#计算R2 ########## ##### lm.Posba <- lmer(log(Posba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Posba) lm.Posba <- lmer(log(Posba+1)~ PCA2+(1|plot), data=env.dat2) summary(lm.Posba) options(na.action = "na.fail") dredge(lm.Pojba,beta="sd") ShowRegTable(lm.Pojba) r.squaredLR(lm.Posba)#计算R2 ########## ##### lm.Psaba <- lmer(log(Psaba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Psaba) lm.Psaba <- lmer(log(Psaba+1)~ PCA1+PCA2+pH+canopy+(1|plot), data=env.dat2) summary(lm.Psaba) options(na.action = "na.fail") dredge(lm.Psaba,beta="sd") ShowRegTable(lm.Psaba) r.squaredLR(lm.Psaba)#计算R2 ########## ##### lm.Psjba <- lmer(log(Psjba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Psjba) lm.Psjba <- lmer(log(Psjba+1)~ PCA1+PCA2+pH+canopy+(1|plot), data=env.dat2) summary(lm.Psjba) options(na.action = "na.fail") dredge(lm.Psjba,beta="sd") ShowRegTable(lm.Pojba) r.squaredLR(lm.Psjba)#计算R2 ########## ##### lm.Pssba <- lmer(log(Pssba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Pssba) lm.Pssba <- lmer(log(Pssba+1)~ PCA1+PCA2+pH+(1|plot), data=env.dat2) summary(lm.Pssba) options(na.action = "na.fail") dredge(lm.Pssba,beta="sd") ShowRegTable(lm.Pssba) r.squaredLR(lm.Pssba)#计算R2 ########## ##### lm.Bpaba <- lmer(log(Bpaba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Bpaba) lm.Bpaba <- lmer(log(Bpaba+1)~ PCA1+canopy+(1|plot), data=env.dat2) summary(lm.Bpaba) options(na.action = "na.fail") dredge(lm.Bpaba,beta="sd") ShowRegTable(lm.Bpaba) r.squaredLR(lm.Bpaba)#计算R2 ########## ##### lm.Bpjba <- lmer(log(Bpjba+1)~ PCA1+PCA2+PCA3+sina+pH+canopy+(1|plot), data=env.dat2) summary(lm.Bpjba) lm.Bpjba <- lmer(log(Bpjba+1)~ PCA1+(1|plot), data=env.dat2) summary(lm.Bpjba) options(na.action = "na.fail") dredge(lm.Bpjba,beta="sd") ShowRegTable(lm.Bpjba) r.squaredLR(lm.Bpjba)#计算R2 ########## ### excl.dat10<-read.csv("glmm.pon1.csv",header = T) excl.dat11<-read.csv("glmm.psn1.csv",header = T) excl.dat12<-read.csv("lmm.poba1.csv",header = T) excl.dat13<-read.csv("lmm.psba1.csv",header = T) #### ## renames the rownames for graphing excl.dat10$parameter <- c("PCA1", "PCA3","pH") # keep the original order appears in the data set excl.dat10 <- transform(excl.dat10, parameter=factor(parameter, levels=rev(unique(parameter)))) # Add a column to identify if the p value is significant (at p = 0.05 level). excl.dat10$tf <- ifelse(excl.dat10$pval <.05, TRUE, FALSE) ## renames the rownames for graphing excl.dat11$parameter <- c("PCA1", "PCA2","pH","canopy") # keep the original order appears in the data set excl.dat11 <- transform(excl.dat11, parameter=factor(parameter, levels=rev(unique(parameter)))) # Add a column to identify if the p value is significant (at p = 0.05 level). excl.dat11$tf <- ifelse(excl.dat11$pval <.05, TRUE, FALSE) ## renames the rownames for graphing excl.dat12$parameter <- c("PCA1", "PCA3","pH","canopy") # keep the original order appears in the data set excl.dat12 <- transform(excl.dat12, parameter=factor(parameter, levels=rev(unique(parameter)))) # Add a column to identify if the p value is significant (at p = 0.05 level). excl.dat12$tf <- ifelse(excl.dat12$pval <.05, TRUE, FALSE) ## renames the rownames for graphing excl.dat13$parameter <- c("PCA1", "pH") # keep the original order appears in the data set excl.dat13 <- transform(excl.dat13, parameter=factor(parameter, levels=rev(unique(parameter)))) # Add a column to identify if the p value is significant (at p = 0.05 level). excl.dat13$tf <- ifelse(excl.dat13$pval <.05, TRUE, FALSE) ##2个物种比较 ## plot the coefficients. p.pon <- ggplot(data=excl.dat10, aes(y=parameter, x=betas, xmin=ci.lb, xmax=ci.ub)) + labs(y="", x="Odds ratios") + geom_errorbarh(height = 0,size=1) + geom_point(aes(shape=tf,fill=tf), size=3.2) + scale_shape_manual(values=c(21,19))+ scale_color_manual(values = c("black", "black"))+ scale_fill_manual(values = c("black", "black"))+ geom_vline(xintercept = 1, linetype=2, color="darkgrey",size=1.2)+ theme_bw() + theme(axis.text=element_text(size=16), axis.title=element_text(size=18), legend.position = "none") + coord_cartesian(xlim=c(0.5,1.5)) P50 <- p.pon + ggtitle("Abundance of P. obovata") + theme(plot.title = element_text(hjust = 0.5,size = 22)) ## #psn ## plot the coefficients. p.psn <- ggplot(data=excl.dat11, aes(y=parameter, x=betas, xmin=ci.lb, xmax=ci.ub)) + labs(y="", x="Odds ratios") + geom_errorbarh(height = 0,size=1) + geom_point(aes(shape=tf,fill=tf), size=3.2) + scale_shape_manual(values=c(21,19))+ scale_color_manual(values = c("black", "black"))+ scale_fill_manual(values = c("black", "black"))+ geom_vline(xintercept = 1, linetype=2, color="darkgrey",size=1.2)+ theme_bw() + theme(axis.text=element_text(size=16), axis.title=element_text(size=18), legend.position = "none") + coord_cartesian(xlim=c(0.5,1.5)) P51<-p.psn + ggtitle("Abundance of P. sibirica") + theme(plot.title = element_text(hjust = 0.5,size = 22)) ## ##poba ## plot the coefficients. p.poba <- ggplot(data=excl.dat12, aes(y=parameter, x=betas, xmin=ci.lb, xmax=ci.ub)) + labs(y="", x="Odds ratios") + geom_errorbarh(height = 0,size=1) + geom_point(aes(shape=tf,fill=tf), size=3.2) + scale_shape_manual(values=c(21,19))+ scale_color_manual(values = c("black", "black"))+ scale_fill_manual(values = c("black", "black"))+ geom_vline(xintercept = 1, linetype=2, color="darkgrey",size=1.2)+ theme_bw() + theme(axis.text=element_text(size=16), axis.title=element_text(size=18), legend.position = "none") + coord_cartesian(xlim=c(0.5,1.5)) P52<-p.poba + ggtitle("Basal area of P. obovata") + theme(plot.title = element_text(hjust = 0.5,size = 22)) ## ##an ## plot the coefficients. p.psba <- ggplot(data=excl.dat13, aes(y=parameter, x=betas, xmin=ci.lb, xmax=ci.ub)) + labs(y="", x="Odds ratios") + geom_errorbarh(height = 0,size=1) + geom_point(aes(shape=tf,fill=tf), size=3.2) + scale_shape_manual(values=c(21,19))+ scale_color_manual(values = c("black", "black"))+ scale_fill_manual(values = c("black", "black"))+ geom_vline(xintercept = 1, linetype=2, color="darkgrey",size=1.2)+ theme_bw() + theme(axis.text=element_text(size=16), axis.title=element_text(size=18), legend.position = "none") + coord_cartesian(xlim=c(0.5,1.5)) P53<-p.psba + ggtitle("Basal area of P. sibirica") + theme(plot.title = element_text(hjust = 0.5,size = 22)) ## library(cowplot) glmm.PP<- cowplot::plot_grid(P50, P51, P52, P53,nrow = 2, labels = LETTERS[1:4])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("LM.POPS.jpeg", glmm.PP, width = 12, height = 12) glmm.PP2<- cowplot::plot_grid(P50, P51, P52, P53,nrow = 2, labels = LETTERS[1:4])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("LM.POPS23.jpeg", glmm.PP2, width = 11, height = 8) ### ###Fig.7######### head(datk) datk1a = melt(datk1,variable.name="Sample",value.name = "Num") datk2a = melt(datk2,variable.name="Sample",value.name = "Num") datak1a <-read.csv("barplot915.csv", header = T) datak1b <-read.csv("hab916.csv", header = T) fit1 <- aov(log(Num+1)~D*sp, data = datak1b) summary(fit1) fit2 <- aov(log(BA+1)~D*sp, data = datak1b) summary(fit2) fit1b1 <- summarySE(datak1b, measurevar="Num", groupvars="type") fit1b2 <- summarySE(datak1b, measurevar="BA", groupvars="type") datak1b2 <-read.csv("hab9162.csv", header = T) factor(datak1b2$D,order=TRUE,levels =c("Hab A1","Hab A2","Hab B1","Hab B2")) factor(datak1b2$sp,order=T,levels =c("Lari_sibi","Pice_obov","Pinu_sibi","Betu_pend")) kns4bar3 = ggplot(datak1b2, aes(x =D,y = Num,fill = sp))+ #####这部分的position = "dodge",并排肩并肩的柱状图 geom_bar(stat ="identity",width = 0.6,position = "dodge")+ geom_errorbar(aes(ymin=(Num-se),ymax=(Num+se)),width=0.2,size=0.02,position = position_dodge(width = 0.6))+ geom_text(aes(x =2.6, y =23 ,label = "habitats ***\nspecies ***\nhabitats:species ***"),size = 4.5)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 #geom_text(aes(x =D, y = 0.45+(Num+se),label = datak1a$ln),position=position_dodge(width = 0.6),size = 5,vjust = 0)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 labs(x = "",y = "Abundance", title = "Abundance in 4 divisional habitats")+ ############坐标标签和图片title #geom_text(aes(label = mean(dat$Num)),position=position_dodge(width = 0.5),size = 5,vjust = -0.25)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 guides(fill = guide_legend(reverse = F))+ ##############图例顺序反转 txt1 #############图例大小 kns4bar4 = ggplot(datak1b2, aes(x =D,y = BA,fill = sp))+ #####这部分的position = "dodge",并排肩并肩的柱状图 geom_bar(stat ="identity",width = 0.6,position = "dodge")+ geom_errorbar(aes(ymin=(BA-seb),ymax=(BA+seb)),width=0.2,size=0.02,position = position_dodge(width = 0.6))+ geom_text(aes(x =2.6, y =3.4 ,label = "habitats *\nspecies ***\nhabitats:species ***"),size = 4.5)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 #geom_text(aes(x =D, y = 0.45+(Num+se),label = datak1a$ln),position=position_dodge(width = 0.6),size = 5,vjust = 0)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 labs(x = "",y=expression(bold(paste("Mean basal area in subplot"," ","(", m^2,")", sep=""))), title = "Basal area in 4 divisional habitats")+ ############坐标标签和图片title #geom_text(aes(label = mean(dat$Num)),position=position_dodge(width = 0.5),size = 5,vjust = -0.25)+ ###########设置柱子上的标签文字,文字的position_dodge(width=0.5)设置,保证分隔宽度。 guides(fill = guide_legend(reverse = F))+ ##############图例顺序反转 txt1 #############图例大小 library(cowplot) knsbar.subplot3 <- cowplot::plot_grid(kns4bar3,kns4bar4,nrow = 1, labels = LETTERS[1:2])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("knsbar.subplot916.jpeg", knsbar.subplot3 , width = 12, height = 5.5) ############ ###Supplement information######### env.dat<-read.csv("knsenv.csv",header=T) ####### Fig S install.packages("ggplot2") library(ggplot2) env.dat$x1 <- env.dat$x+10 env.dat$y1 <- env.dat$y+10 tiff(filename="CV.tiff",width=15,height=8,units="cm",pointsize=8,res=300) par(mfrow=c(1,1),mar=c(4.2,4.5,0.5,2), xpd=TRUE) p1<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = pH), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$pH)) p2<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = SOM), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$SOM)) p3A<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = TN), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$TN)) p4A<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = TP), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$TP)) p3<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = AN), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$AN)) p4<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = AP), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$AP)) p5<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = AK), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$AK)) p6<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = elevation), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$elevation)) p7<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = convexity), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$convexity)) p8<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = slope), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$slope)) #p9<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = sina), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$sina)) p9<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = aspect_EW), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$aspect_EW)) #p10<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = cosa), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$cosa)) p11<-ggplot(env.dat, aes(x1, y1)) + geom_tile(color="white",aes(fill = canopy), size=0.6) + xlab("East-West axis (m)") + ylab("South-North axis (m)") + scale_fill_gradient2(low="blue",high="red",mid="yellow",midpoint=median(env.dat$canopy)) install.packages("cowplot") library(cowplot) p12ENV <- cowplot::plot_grid(p1, p2, p3A,p4A,p3, p4, p5,p6,p8,p9,p7,p11,nrow = 4, labels = LETTERS[1:12])#将p1-p4四幅图组合成一幅图,按照两行两列排列,标签分别为A、B、C、D。(LETTERS[1:4] 意为提取26个大写英文字母的前四个:A、B、C、D) ggsave("plot.ENVb1.jpeg", p12ENV, width = 12.15, height =9.95 ) #, width = 20.15, height = 18.5 ####PCA install.packages("psych") library(psych) env.cor<-env.dat[,c(2:13)] sp.cor<-env.dat[,c(14:17)] sp0.cor<-env.dat[,c(14:20)] topo.cor<-env.dat[,c(9:12)] soil.cor<-env.dat[,c(2:8)] light.cor<-env.dat[,c(13)] ba.cor<-env.dat[,c(17:23)] cortest.bartlett(cor(env.cor)) KMO(cor(env.cor)) library(ade4) env.pca <- dudi.pca(env.dat[,c(3:11)], scale=TRUE, scan=FALSE, nf=5) summary(env.pca) env.pca$li # soil and topography variables scores on the PCA axes env.pca$co pve <- 100* env.pca$eig / sum(env.pca$eig) # eigenvalues in % tiff(filename="knsenv.PCA722.tiff",width=20,height=9,units="cm",pointsize=8,res=300) par(mfrow=c(1,2),mar=c(4.2,4.8,0.5,2), xpd=TRUE) #barplot(pve, col="grey", ylab="% explained", xlab="Axes of PCA") s.corcircle(env.pca$co, xax=1, yax=2, grid=F, box=F, clabel=1) text(0, -1.15, "Axis 1 (28.6%)", cex = 1.2) text(-1.2, 0, "Axis 2 (20.3%)", srt = 90, cex = 1.2) s.corcircle(env.pca$co, xax=1, yax=3, grid=F, box=F, clabel=1) text(0, -1.15, "Axis 1 (28.6%)", cex = 1.2) text(-1.2, 0, "Axis 3 (17.0%)", srt = 90, cex = 1.2) dev.off()