多种群多地点MARSS模型

ukqbszuj  于 2023-11-14  发布在  其他
关注(0)|答案(1)|浏览(110)

我一直在试图弄清楚如何正确格式化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")

字符串

yqlxgs2m

yqlxgs2m1#

我可能刚刚弄明白了我的问题(请参阅下面的代码)。我将把这个问题留在上面,以防有人知道我是否做错了。

# Using the same data as in the original question
y <- rbind(site1,site2)
y

# Create needed matrices (doing this manually for clarity)
# B-matrix formatting
b.model <- matrix(list(0),6,6)

b.model[1:3,1:3] <- c("A.b11","A.b12","A.b13",   # Each site has a 3x3 "matrix" embedded within
                      "A.b21","A.b22","A.b23",   # the overall 6x6 matrix. Not sure this is an
                      "A.b31","A.b32","A.b33")   # appropriate way of doing this, the plan would  
                                                 # be to isolate out the individual B-matrices      
b.model[4:6,4:6] <- c("B.b11","B.b12","B.b13",   # for each site from the "overall B-matrix"
                      "B.b21","B.b22","B.b23",
                      "B.b31","B.b32","B.b33")

# Z-matrix formatting
z.model <- matrix(0,6,6)
diag(z.model) <- 1        # same as "identity"
z.model

# Q-matrix formatting
q.model <- matrix(list(0),6,6)
diag(q.model) <- c("Qsp1A","Qsp2A","Qsp3A",
                   "Qsp1B","Qsp2B","Qsp3B")   # same as "diagonal & unequal"
q.model

# R-matrix formatting
r.model <- matrix(0,6,6)   # same as "zero" (for now considering detection as being perfect)

# U-matrix formatting
u.model <- matrix(c("Usp1A","Usp2A","Usp3A",   # same as "unconstrained"
                    "Usp1B","Usp2B","Usp3B"),
                  6,1) 

# Put together all parts of model
model.0 <- list(Z = z.model,                
                B = b.model,
                Q = q.model,  
                R = r.model,                     
                U = u.model,                
                tinitx = 0) 

# Run the model
m0.out <- MARSS(y, model = model.0, method = "BFGS")

m0.B <- coef(m0.out, type = "matrix")$B[1:6, 1:6]
m0.B

# Separating out individual site's B-matrices
# Again, not sure this is appropriate/ok to do...
m0.B.siteA <- m0.B[1:3,1:3]
m0.B.siteB <- m0.B[4:6,4:6]

max(abs(eigen(m0.B.siteA)$values))
max(abs(eigen(m0.B.siteB)$values))

字符串

相关问题