R矩阵值回收？

Question

我是遵循这个YouTube教程基于代理的建模在R https://www.youtube.com/watch?v=uAeSykgXnhg.

我使用自己的变量名复制了代码(这有助于我更好地理解导师的代码)。其目的是跟踪人们在接触他人时如何感染新冠肺炎。不是每一次接触都会导致感染。在连续模型运行中，受感染人数=人口规模，未受感染人数应为0。这是我复制的代码：

# define first agent
agents <- data.frame(agent_no = 1,
                     state = "e",
                     mixing = runif(1,0,1))

# specify agent population
pop_size <- 100

# fill agent data
for(i in 2:pop_size){
  agent <- data.frame(agent_no = i,
                      state = "s",
                      mixing = runif(1,0,1))
  agents <- rbind(agents, agent)
}

# specify number of model runs
n_times <- 10

# initialise output matrix 
out <- matrix(0, ncol = 2, nrow = n_times)

# run simple agent-based model
for(k in 1:n_times){
  for(i in 1:pop_size){
    # likelihood to meet others
    likelihood <- agents$mixing[i]
    # how many agents will they meet (integer). Add 1 to make sure everybody meets somebody
    connect_with <- round(likelihood * 3, 0) + 1 
    # which agents will they probably meet (list of agents)
    which_others <- sample(1:pop_size, 
                           connect_with, 
                           replace = T, 
                           prob = agents$mixing)
    for(j in 1:length(which_others)){
      contacts <- agents[which_others[j],]
      # if exposed, change state
      if(contacts$state == "e"){
        urand <- runif(1,0,1)
        # control probability of state change
        if(urand < 0.5){
          agents$state[i] <- "e"
        }
      }
    }
  }
  out[k,] <- table(agents$state)
}

当查看输出时，一旦每个人都被感染(第一列)，未感染的人数(第二列)应该是0，但我得到100，我怀疑这是由于回收。

     [,1] [,2]
 [1,]   12   88
 [2,]   33   67
 [3,]   69   31
 [4,]   86   14
 [5,]   92    8
 [6,]   95    5
 [7,]   97    3
 [8,]   98    2
 [9,]   99    1
[10,]  100  100

我做了一些诊断来看看是怎么回事：

table(agents$state)
      e 
    100 

agents[agents$state == "s",]

    [1] agent_no state    mixing  
    <0 rows> (or 0-length row.names)

我认为0长度的row.names是我的问题.结果应该是这样：

     [,1] [,2]
 [1,]   12   88
 [2,]   33   67
 [3,]   69   31
 [4,]   86   14
 [5,]   92    8
 [6,]   95    5
 [7,]   97    3
 [8,]   98    2
 [9,]   99    1
[10,]  100    0

有人能解释我做错了什么吗？非常感谢。

Answer 1

我把n_times提高到了10000，没有发现回收的证据。虽然这并不意味着它不会发生，但这意味着不幸的是，如果没有明确的设置，我们将很不幸地无法重现问题。所以我在这里的建议是未经证实的。

备选案文1

假设您找到了一个以所有agents$state == "e"结尾的这样的场景，那么我将建议一个技巧，它总是至少找到一个"s" (实际上，您所知道的每个值中的一个)：

  out[k,] <- table(c("e", "s", agents$state)) - 1

我假设唯一可能的值是"e"和"s"；如果有其他值，这种技术完全依赖于这样的前提:我们确保至少看到每个可能的值一次，然后减少所有的值。因为我们对每个可能的值“加一个观察”，所以从表中减去一个是安全的。用这个技巧，你的支票就应该是

table(agents$state)
#       e 
#     100 
table(c("e", "s", agents$state))
#       e       s 
#     101       1
table(c("e", "s", agents$state)) - 1
#       e       s 
#     100       0

因此，回收不应成为一个因素。

选项2

另一种更可靠的技术(即，不需要包含所有可能的值)是强制长度，假设我们确实知道它应该是什么(我认为我们在这里这样做)：

z <- table(agents$state)
z
#   s 
# 100 
length(z) <- 2
z
#   s     
# 100  NA 

由于您“知道”长度应该始终为2，所以可以对其中的2进行硬编码。

选项3

这个方法甚至更健壮一些，因为您不需要知道绝对长度，它们都将被扩展到最长返回的长度。

第一，可重复的样本数据：

set.seed(2021)
agents <- data.frame(agent_no = 1,
                     state = "e",
                     mixing = runif(1,0,1))
# specify agent population
pop_size <- 100
# fill agent data
for(i in 2:pop_size){
  agent <- data.frame(agent_no = i,
                      state = "s",
                      mixing = runif(1,0,1))
  agents <- rbind(agents, agent)
}
head(agents)
#   agent_no state    mixing
# 1        1     e 0.4512674
# 2        2     s 0.7837798
# 3        3     s 0.7096822
# 4        4     s 0.3817443
# 5        5     s 0.6363238
# 6        6     s 0.7013460

替换您的for循环：

for (k in 1:n_times) {
}

使用

out <- lapply(seq_len(n_times), function(k) {
  for(i in 1:pop_size){
    # likelihood to meet others
    likelihood <- agents$mixing[i]
    # how many agents will they meet (integer). Add 1 to make sure everybody meets somebody
    connect_with <- round(likelihood * 3, 0) + 1 
    # which agents will they probably meet (list of agents)
    which_others <- sample(1:pop_size, 
                           connect_with, 
                           replace = T, 
                           prob = agents$mixing)
    for(j in 1:length(which_others)){
      contacts <- agents[which_others[j],]
      # if exposed, change state
      if(contacts$state == "e"){
        urand <- runif(1,0,1)
        # control probability of state change
        if(urand < 0.5){
          agents$state[i] <- "e"
        }
      }
    }
  }
  table(agents$state)
})

此时，您有一个可能是长度-2向量的列表：

out[1:3]
# [[1]]
#  e  s 
#  1 99 
# [[2]]
#  e  s 
#  2 98 
# [[3]]
#  e  s 
#  3 97 

注意，我们可以用

lengths(out)
#  [1] 2 2 2 2 2 2 2 2 2 2

类似于选项2，即我们强制向量的长度，我们可以在这里做同样的事情：

maxlen <- max(lengths(out))
out <- lapply(out, `length<-`, maxlen)
## or more verbosely
out <- lapply(out, function(vec) { length(vec) <- maxlen; vec; })

您可以确认它们与table(lengths(out))都是相同的长度，应该是n_times为10的2。

从这里开始，我们可以将所有这些向量组合成一个矩阵

out <- do.call(rbind, out)
out
#        e  s
#  [1,]  1 99
#  [2,]  2 98
#  [3,]  3 97
#  [4,]  2 98
#  [5,]  1 99
#  [6,] 20 80
#  [7,] 12 88
#  [8,]  1 99
#  [9,]  2 98
# [10,]  1 99