我一直在试图弄清楚如何正确格式化MARSS模型(使用R中的MARSS包),这样我就可以在多个地点同时对多个物种进行建模,如果可能的话。在我完整的研究设计中,我有30年来9个不同地点的9个物种的年度样本。我已经能够让模型运行一个地点多个物种和多个地点只有一个物种(参见下面的reprex,尽管我不完全确定这些代码是否正确),但还没有能够弄清楚如何把它放在一起,或者如果这甚至是可以做的事情。如果这是一个适当的方法,我也可以只运行每个网站单独的模型。理想情况下,我想得到单独的B-矩阵为我的每一个网站,使我可以比较网站之间的社区稳定性指标。如果任何人有任何想法(或如果我去这个错误),我真的很感激!
library(MARSS)
dat <- data.frame(site = rep(1:2,each = 20), # two independent sites
year = rep(c(1:20),2), # sample years (1-20)
spp1 = sample(50,40, replace = T), # simulated count data species 1
spp2 = sample(20,40, replace = T), # simulated count data species 2
spp3 = sample(30,40, replace = T)) # simulated count data species 3
# Put count data into wide format (rows = species, columns = years)
site1 <- t(dat[1:20,3:5])
site2 <- t(dat[21:40,3:5])
rownames(site1) <- paste0("site1_", rownames(site1))
rownames(site2) <- paste0("site2_", rownames(site2))
y <- rbind(site1,site2)
# RUN THE MODEL FOR ONE SPECIES AND BOTH SITES
y1 <- rbind(y[1,], y[4,])
# Z-matrix formatting
z.model1 <- matrix(c(1,0,0,1),ncol = 2)
# A-matrix formatting
a.model1 <- matrix(list(0),2,1)
a.model1[2,1] <- c("b")
# U-matrix formatting
u.model1 <- matrix(list(0),2,1)
u.model1[1:2,1] <- c("a","b")
# Put together all parts of model
model.1 <- list(Z = z.model1,
A = a.model1,
Q = "diagonal and unequal",
R = "diagonal and equal",
U = u.model1,
tinitx = 1)
# Run the model
m1.out <- MARSS(y1, model = model.1, method = "BFGS")
# NOW RUN THE MODEL WITH ALL THREE SPECIES BUT ONLY ONE SITE
y2 <- y[1:3,]
# Z-matrix formatting
z.model2 <- matrix(0,3,3)
diag(z.model2) <- 1
# U-matrix formatting
u.model2 <- matrix(list(0),3,1)
u.model2[1:3,1] <- c("a","b","c")
# B-matrix formatting
b.model2 <- matrix(c("b11","b12","b13",
"b21","b22","b23",
"b31","b32","b33"),
ncol = 3)
# Put together all parts of model
model.2 <- list(Z = z.model2,
B = b.model2,
Q = "diagonal and unequal",
R = "diagonal and equal",
U = u.model2,
tinitx = 1)
# Run the model
m2.out <- MARSS(y2, model = model.2, method = "BFGS")
字符串
1条答案
按热度按时间yqlxgs2m1#
我可能刚刚弄明白了我的问题(请参阅下面的代码)。我将把这个问题留在上面,以防有人知道我是否做错了。
字符串