首页 > 解决方案 > 未将行添加到R中的数据框

问题描述

所以,我一直在尝试对 SIR 模型进行模拟研究。

我有以下代码(试图清理它):

# Initial parameters
N <- 1E6            # Total population
I <- 1              # Number of Infectious at time 0
S <- N-1            # Number of Susceptibles at time 0
R <- 0             # Number of Recovered at time 0

# Vector to store observations in
Df1 <- data.frame("final_size" = as.numeric(), "peak_size" = as.numeric())

# Setting a seed for reproducibility 
set.seed(1996)
n_sim <- 100

# Setting different values for R0, the basic reproduction number
R0 <- seq(0.5,2.5, 1)

for(values in R0){
  if(values == 0.5){
    # Transmission parameters
    R0 <- values    # Basic Reproduction number
    nu <- 1/6             # Recovery rate (in days) 
    b <- nu*R0/N        # Infection rate (in days)
    for(sim in n_sim){
      temp <- NULL
      # Binomial model
      #----------------
      # Initial states
      Sold = S          # Number of Susceptibles at time t=0
      Iold = I          # Number of Infectious at time t=0
      Rold = R          # Number of Recovered at time t=0
      # Output vectors
      Svec =Sold; Ivec = Iold; Rvec = Rold
      stop = FALSE
      # Loop - continue until stop=TRUE
      while (!stop){        
        Ih = rbinom(1,Sold,(1-exp(-b*Iold)))
        Rh = rbinom(1,Iold,(1-exp(-nu)))
        Rh = nu*Iold
        Sold = Sold-Ih
        Iold = Iold+Ih-Rh
        Rold = Rold+Rh
        Svec = c(Svec,Sold)
        Ivec = c(Ivec,Iold)
        Rvec = c(Rvec,Rold)
        if (Iold<=2e-5){stop=T}
      }
      peak_size_df <- max(Ivec)
      final_size_df <- Rvec[length(Rvec)]/N
      temp <- rbind(temp, c(final_size_df, peak_size_df))
      colnames(temp) <- c("final_size", "peak_size")
      Df1 <- rbind(Df1, temp)
    }
  }
}

我正在寻找将数据存储在 Df1 中。但是,在循环结束时,只存储了 1 个循环,我假设是最后一个循环。我不太明白。我已经重新编码了几次,在这些情况下,我最终得到了 99 NA 行的最终大小和峰值大小。在这个版本中,我只得到 1 行(但是有值)。我计划为 R0 序列中看到的不同 R0 值扩展循环。由于它不适用于第一个值,我还没有扩展它。有什么建议么?改进?

在 Gregory 和 Cath 发表意见后,进行以下调整:

# Setting a seed for reproducibility 
#set.seed(1996)
n_sim <- 100


# Initial parameters
N <- 1E6            # Total population
I <- 1              # Number of Infectious at time 0
S <- N-1            # Number of Susceptibles at time 0
R <- 0            # Number of Recovered at time 0


# Vector to store observations in
Df1 <- data.frame("final_size" = rep(NA, n_sim), "peak_size" = rep(NA, n_sim))
Df2 <- NULL
Df3 <- NULL


# Setting different values for R0, the basic reproduction number
R0 <- seq(0.5,2.5, 1)



#plot(Svec, type = "l", ylim = c(0, 1000000), col = "red")
#lines(Rvec, type = "l", col = "blue")
#lines(Ivec, type = "l")
#max(Ivec)

for(values in R0){
  if(values == 0.5){
    # Transmission parameters
    R0_value <- values    # Basic Reproduction number
    nu <- 1/6             # Recovery rate (in days) 
    b <- nu*R0_value/N      # Infection rate (in days)
    for(sim in n_sim){
      # Binomial model
      #----------------
      # Initial states
      Sold = S          # Number of Susceptibles at time t=0
      Iold = I          # Number of Infectious at time t=0
      Rold = R          # Number of Recovered at time t=0
      # Output vectors
      Svec =Sold; Ivec = Iold; Rvec = Rold
      stop = FALSE
      # Loop - continue until stop=TRUE
      while (!stop){        
        Ih = rbinom(1,Sold,(1-exp(-b*Iold)))
        Rh = rbinom(1,Iold,(1-exp(-nu)))
        Rh = nu*Iold
        Sold = Sold-Ih
        Iold = Iold+Ih-Rh
        Rold = Rold+Rh
        Svec = c(Svec,Sold)
        Ivec = c(Ivec,Iold)
        Rvec = c(Rvec,Rold)
        if (Iold<=2e-5){stop=T}
      }
      peak_size_df <- max(Ivec)
      final_size_df <- Rvec[length(Rvec)]/N
      Df1[sim, "final_size"] <- final_size_df
      Df1[sim, "peak_size"] <- peak_size_df
    }
  }
}

这给了我 99 行 NA 和第 100 行最后一个模拟值。知道是什么导致了NA吗?

标签: rdataframe

解决方案


推荐阅读