Я создал функцию в R, которая создает один размерный вектор различной длины (размер от 1 до сотен +) - код внизу. Я хотел бы сделать две вещи, которые, я надеюсь, довольно просты для кого-то гораздо более опытного в R, чем я.
Для контекста, эта функция имитирует вспышку заболевания, используя процесс ветвления с отрицательным биномиальным распределением потомства из одного случая индекса. Каждая позиция в векторе представляет человека, а каждое целое число представляет количество вторичных случаев, вызванных этим человеком.
Например, [3 0 1 0 1 0] указывает, что было шесть случаев, которые передали соответствующее количество вторичных случаев. Более конкретно, это означает:
- Индексный регистр, передаваемый в
3 вторичные регистры (первая позиция в векторе) следующему поколению. - Из этих трех случаев первый случай имел
0 вторичные случаи, второй имел 1 вторичный случай, а третий имел 0 (позиции 2: 4 в векторе). - Один случай в следующем поколении также передал случай
1 (позиция 5), который имел 0 вторичные случаи (позиция 6), поэтому процесс вышел из цикла и выходов вектора. - Если в индексном случае не было заражений, вектор был бы просто
[0].
Хотя это только начало, я хотел бы улучшить его двумя следующими способами:
Во-первых, каждый раз, когда цикл продолжается, это новое «поколение». Я хотел бы создать матрицу 2 x n, которая будет обозначать соответствующее поколение. Например, данные для приведенного выше примера будут выглядеть следующим образом:
##' Example desired dataset of running the function
##' Index case transmits to three others (generation 1)
##' Those three transmit to 0, 1, and 0, respectively (generation 2)
##' The single case (Gen 3) transmits to 0 (gen 4) and the function exits the loop
generation1 <- c(1,2,2,2,3,4)
secondary.cases1 <- c(3,0,1,0,1,0)
example.matrix1 <- cbind(generation1, secondary.cases1)
print(example.matrix1)
# generation1 secondary.cases1
# [1,] 1 3
# [2,] 2 0
# [3,] 2 1
# [4,] 2 0
# [5,] 3 1
# [6,] 4 0
Во-вторых, функция предназначена только для одной вспышки. Для моих окончательных данных мне нужно смоделировать это несколько раз (т.е. 2000) и объединить все эти векторы. Также, аналогично 1 выше, есть номер, обозначающий соответствующий номер моделирования. Например, окончательный набор данных будет выглядеть следующим образом:
##' Example desired dataset of three simulations:
##' sim1 producing a vector of [3 0 1 0 1 0] in four generations
##' sim2 producing a vector of [2 0 0] in three generations
##' sim3 producing a vector of [0] in the first generation
simulation2 <- c(1,1,1,1,1,1,2,2,2,3)
generation2 <- c(1,2,2,2,3,4,1,2,3,1)
secondary.cases2 <- c(3,0,1,0,1,0,2,0,0,0)
example.matrix2 <- cbind(simulation2, generation2, secondary.cases2)
print(example.matrix2)
# simulation2 generation2 secondary.cases2
# [1,] 1 1 3
# [2,] 1 2 0
# [3,] 1 2 1
# [4,] 1 2 0
# [5,] 1 3 1
# [6,] 1 4 0
# [7,] 2 1 2
# [8,] 2 2 0
# [9,] 2 3 0
# [10,] 3 1 0
Функция ниже. Для первой проблемы я попытался посчитать каждую итерацию как поколение и создать вектор. Для второй проблемы я попытался вложить эту функцию в другую функцию. Обе эти вещи не будут работать должным образом (с моим недостатком опыта, то есть). После нескольких часов попыток самостоятельно я обращаюсь за помощью. Любой совет будет очень признателен!
##' The following function creates a single vector of
##' secondary cases from a negative binomially distributed offspring distribution
##' Each row represents a specific individual in the chain of transmission,
##' and the integer indicates how many secondary cases the individual transmitted
##' (0 indicates no secondary transmission from the index case)
##' This function only does a single outbreak. I will need to replicate this.
##' @param r0 Average number of secondary infections per each infectious individual - keep r0<1 so Pr(Extinction)=1
##' @param k Dispersion parameter of NB distribution (0-infinity)
##' @param index Number of index cases initiating the cluster (typically 1)
##' @param max_gen Maximum number of generations this loop should go through, to avoid hanging with near infinite loops.
#########################################################################################################
#########################################################################################################
nbbp.individual<-function(r0,k,index=1,max_gen){
G1<-rnbinom(index,size=k,mu=r0) #Does the first iteration outside fo the while to avoid problems
cases<-G1 #initiates the number of cases in the first generation
sumCases<-G1 #since G1 will be a single integer, this is somewhat repeative but comes into play later in the loop
gen<-1 #initiates the generation for setting the max generations (to avoid hanging of near infinite loops)
gencount<-0
while (sumCases>0 && (missing(max_gen)||(gen<max_gen))){
gen<-gen+1
individuals<-rep(1,sumCases) #This creates a single vector of 1's to simulate the single individuals in the generaiton (i.e.[1 1 1 1]
secondaryCases<-rnbinom(individuals,size=k,mu=r0) #This creates a vector of integers, each one corresponding to the
#number of secondary cases for each "individual" (i.e. [0 2 0 4])
sumCases<-sum(secondaryCases) #this sums the generation for purposes of the loop
cases<-c(cases, secondaryCases) #This appends the individuals for each generation
}
return(cases)
}
#########################################################################################################
#########################################################################################################
##' If you highlight this and run it over and over, it is essentially running multiple simulations of
##' individual tranmission chains
##' You will see that "test" mainly sticks to 0, but can range to the hundreds
test<-nbbp.individual(0.9,0.25,1,20)
print(test)