这是我的 NMDS 排序矩阵(PA.mm.5p.rel)。
structure(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,
1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 1), .Dim = c(45L, 8L), .Dimnames = list(c("3",
"4", "5", "7", "8", "10", "11", "12", "13", "15", "16", "18",
"19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29",
"30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40",
"41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51"
), c("Onisimus", "Amphipods", "Krill", "Copepods", "Fish", "Gammarus",
"Arctic.cod", "Sand.Lance")))
这是我的环境/个人数据(PA_Info)。
structure(list(Date = structure(c(3L, 3L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 9L,
10L, 12L, 14L, 15L, 16L, 16L, 16L, 16L, 17L, 18L, 18L, 19L, 19L,
19L, 20L, 20L, 20L, 21L, 21L, 22L, 23L, 25L), .Label = c("2017-07-19",
"2017-07-25", "2017-08-07", "2017-08-08", "2017-08-09", "2017-08-12",
"2017-08-14", "2017-08-15", "2017-08-22", "2017-09-01", "2017-09-04",
"2018-07-20", "2018-07-23", "2018-07-26", "2018-07-29", "2018-07-30",
"2018-08-01", "2018-08-02", "2018-08-04", "2018-08-06", "2018-08-08",
"2018-08-14", "2018-08-19", "2018-08-20", "2018-08-28", "2019-08-08",
"2019-08-10", "2019-08-11", "2019-08-12", "2019-08-13", "2019-08-14"
), class = "factor"), Floy = c(NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2019L,
2014L, NA, 2060L, NA, 2575L, NA), FishID = structure(c(16L, 15L,
25L, 21L, 23L, 22L, 20L, 19L, 18L, 24L, 28L, 26L, 29L, 27L, 30L,
8L, 5L, 7L, 4L, 10L, 3L, 11L, 12L, 13L, 31L, 33L, 34L, 35L, 36L,
37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L,
50L, 51L, 53L), .Label = c("Char 08", "Char 09", "Char 112",
"Char 116", "Char 121", "Char 123", "Char 126", "Char 129", "Char 130",
"Char 132", "Char 138", "Char 150", "Char 154", "Char 156", "Char 54",
"Char 56", "Char 58", "Char 60", "Char 61", "Char 62", "Char 63",
"Char 65", "Char 66", "Char 68", "Char 69", "Char 86", "Char 87",
"Char 88", "Char 89", "Char 90", "SSS001", "SSS002", "SSS003",
"SSS004", "SSS005", "SSS006", "SSS007", "SSS008", "SSS009", "SSS010",
"SSS011", "SSS012", "SSS013", "SSS014", "SSS015", "SSS016", "SSS017A",
"SSS017B", "SSS018", "SSS019", "SSS020", "SSS021", "SSS022",
"SSS23", "SSS24", "SSS25", "SSS26", "SSS27", "SSS28", "SSS29",
"SSS30", "SSS31", "SSS32", "SSS33", "SSS34", "SSS35", "SSS36",
"SSS37", "SSS38", "SSS39", "SSS40", "SSS41", "SSS42", "SSS43",
"SSS44", "SSS45", "SSS46", "SSS47", "SSS48", "SSS49", "SSS50",
"SSS51", "SSS52", "SSS53", "TCH01", "TCH02", "TCH03", "TCH04",
"TCH05", "TCH06", "TCH07"), class = "factor"), TL = c(71, 33.2,
35.5, 59.5, 53.7, 65.3, 54.5, 57.5, 61.5, 60, 56.5, 53.5, 55,
57, 61, 38.5, 65.2, NA, NA, 59, 29.5, 17, 79.5, 81.3, 90, NA,
67, 73.2, 86, 85, 63.2, 72.4, 65.4, 66.8, 67, 59, 63.2, 79, 75,
90.8, 76, 41.4, 82, 83, 77.6), FL = structure(c(55L, 6L, 7L,
25L, 16L, 40L, 19L, 22L, 31L, 27L, 21L, 18L, 20L, 22L, 29L, 8L,
43L, NA, NA, 25L, 5L, 3L, 67L, 72L, 63L, NA, 44L, 57L, 62L, 71L,
35L, 53L, 45L, 47L, 42L, 24L, 39L, 65L, 60L, 73L, 61L, 9L, 69L,
70L, 62L), .Label = c("13.5", "15.8", "16", "16.5", "28.3", "31.5",
"34.5", "38", "39.4", "40.6", "42", "44", "46.8", "48.6", "50",
"51.5", "52", "53", "53.2", "53.5", "55.4", "56", "57", "57.2",
"57.5", "58", "58.3", "58.4", "58.8", "59", "59.2", "59.6", "60",
"61", "61.1", "61.2", "61.4", "61.8", "62", "62.5", "62.6", "63",
"63.5", "64", "65", "65.4", "65.6", "66", "66.4", "66.8", "67",
"67.6", "68", "69", "69.9", "70", "70.1", "70.4", "72.6", "72.8",
"75", "76", "76 (pcl)", "76.4", "76.8", "77.4", "77.5", "78",
"79.8", "81.4", "83", "83.8", "89.8"), class = "factor"), Girth = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, 29, NA, NA, NA, NA, 32.4,
31, 30, NA, NA, 31, 36, 38.6, 43, NA, 20.6, 37, 46, 35.4), Mass = c(3400,
8.2, 450, 1800, 1500, 2900, 1300, 1300, 2100, 1900, 1650, 1450,
1625, 1600, 1850, 550, 2700, 3600, NA, 1300, 237.2, 41.2, 4850,
5100, 6000, NA, 2350, 2700, 4550, 5600, 1750, 2400, 2020, 1800,
NA, 1940, 2250, 4000, 3720, 5550, 3150, 620, 4350, 6600, 3900
), Sex = structure(c(5L, 3L, 5L, 1L, 1L, 5L, 5L, 5L, 1L, 5L,
5L, 5L, 5L, 1L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 3L, 1L, 5L, NA, NA,
NA, 5L, 5L, 5L, 1L, 5L, 5L, 1L, 5L, 1L, 1L, 5L, 5L, 5L, 1L, 1L,
1L, 5L, 5L), .Label = c("F", "F?", "I", "J", "M"), class = "factor"),
Year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("2017", "2018"), class = "factor"),
MassCat = c("Medium", "Small", "Small", "Medium", "Small",
"Medium", "Small", "Small", "Medium", "Medium", "Medium",
"Small", "Medium", "Medium", "Medium", "Small", "Medium",
"Medium", NA, "Small", "Small", "Small", "Large", "Large",
"Large", NA, "Medium", "Medium", "Large", "Large", "Medium",
"Medium", "Medium", "Medium", NA, "Medium", "Medium", "Large",
"Large", "Large", "Medium", "Small", "Large", "Large", "Large"
), Datect = structure(c(1502078400, 1502078400, 1502251200,
1502251200, 1502251200, 1502251200, 1502251200, 1502251200,
1502251200, 1502251200, 1502510400, 1502510400, 1502510400,
1502510400, 1502510400, 1502683200, 1502683200, 1502683200,
1502683200, 1502683200, 1502683200, 1502769600, 1503374400,
1504238400, 1532059200, 1532577600, 1532836800, 1532923200,
1532923200, 1532923200, 1532923200, 1533096000, 1533182400,
1533182400, 1533355200, 1533355200, 1533355200, 1533528000,
1533528000, 1533528000, 1533700800, 1533700800, 1534219200,
1534651200, 1535428800), class = c("POSIXct", "POSIXt"), tzone = ""),
week = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 7L, 1L,
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L,
4L, 4L, 5L, 5L, 7L), .Label = c("29", "30", "31", "32", "33",
"34", "35"), class = "factor"), yday = c(219, 219, 221, 221,
221, 221, 221, 221, 221, 221, 224, 224, 224, 224, 224, 226,
226, 226, 226, 226, 226, 227, 234, 244, 201, 207, 210, 211,
211, 211, 211, 213, 214, 214, 216, 216, 216, 218, 218, 218,
220, 220, 226, 231, 240)), row.names = c(10L, 11L, 13L, 14L,
15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L,
28L, 29L, 30L, 33L, 34L, 35L, 36L, 38L, 40L, 41L, 42L, 43L, 44L,
45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L,
58L, 60L), class = "data.frame")
创建 NMDS 和点的数据框。
PA.mm.5p.rel.nmds2<- metaMDS(PA.mm.5p.rel, distance="bray", k = 2, autotransform=FALSE)
allNMDS_DF<-as.data.frame(PA.mm.5p.rel.nmds2$points)
运行 ordisurf 并使用函数从 sp.sf 对象中提取 x、y 和 z 坐标。
sp.sf <- ordisurf(PA.mm.5p.rel.nmds2 ~ PA_Info$TL)
extract.xyz <- function(obj) {
xy <- expand.grid(x = obj$grid$x, y = obj$grid$y)
xyz <- cbind(xy, c(obj$grid$z))
names(xyz) <- c("x", "y", "z")
return(xyz)
}
这是我提取坐标的地方,但是 contour.vals 中的 z 值与 .vals 中观察到的值范围不匹配PA_Info$TL
。当您绘制并将“级别”中的比例与“总长度”进行比较时,您还可以在图例中看到这一点。因为存在差异,我担心我不能相信 ordisurf 或相关 p 值 ( summary(sp.sf)
) 的输出。有谁知道可能出了什么问题?
contour.vals <- extract.xyz(obj = sp.sf, plot = F)
head(contour.vals)
contour.vals <- na.omit(contour.vals)
allNMDS_DFTL <- allNMDS_DF[which(!is.na(PA_Info$TL)),]
p <- ggplot(data = contour.vals, aes(x, y, z = z)) + stat_contour(aes(colour = ..level..)) +
coord_cartesian(xlim = c(-2, 2), ylim = c(-1, 1.5)) + theme_bw()+
geom_point(data = allNMDS_DFTL, aes(MDS1, MDS2, z = 0, shape = PA_Info$Year[which(!is.na(PA_Info$TL))], fill = PA_Info$TL[which(!is.na(PA_Info$TL))]), size=8,shape=21)+
labs(fill = 'Total Length (cm)') + xlab("Axis 1 (38.2%)")+ylab("Axis 2 (27.9%)")