Создание ветвящегося процесса с отдельными вторичными наблюдениями в R

Question

Я новичок в этом форуме, поэтому, пожалуйста, дайте мне знать, если это неуместная тема для обсуждения. Я также относительно новичок в R (преобразование SAS!), Так что терпите меня.

Я пытаюсь смоделировать данные о вспышках инфекционных заболеваний в R, используя процесс ветвления (или процесс Гальтона-Ватсона) с отрицательнымраспределение биномиального потомства. Пример процесса ветвления представлен ниже, взятый несколько произвольно из Cui и др. . Вкратце, в контексте инфекционных заболеваний каждый процесс ветвления начинается с одного инфекционного индивидуума (поколение 0), и они заражают случайное число вторичных случаев, определяемых распределением некоторых потомков (я использую отрицательное биномиальное потомство с параметрами (среднее = r0,дисперсия = к)). Каждый случай в последующих поколениях заражает случайное количество случаев и т. Д. (Iid).

Моя главная цель состоит в том, чтобы идентифицировать число вторичных случаев каждого человека, обозначаемых как Z (то есть на изображении выше, для поколения 0 [G (0)]: Z = 2, в G (1), Z = 2 и Z = 1, для двух индивидуумов слева направо соответственно и т. Д.).

Я представляю себе макет данных как матрицу, подобную этой, которая представляет первые три поколения на рисунке:

#Example of data layout from first three generations in figure
x<-matrix(c(0,1,1,2,2,2,1,1,2,1,2,3,2,2,1,1,0,3),ncol = 3,nrow = 6)
colnames(x)<-c("Generation","Case Number","Secondary Cases (Z)")

Я не уверен, как создать функцию, которая дразнитиз числа отдельных значений Z. Я написал функцию, которая имитирует отрицательный биномиальный процесс ветвления, но только суммирует общее количество случаев в каждом поколении (то есть на изображении ниже, G (0) = 1, G (1) = 2, G (2) =3, G (3) = 4). Кроме того, в этой функции я должен определить количество поколений («n»), и более точный код позволит ему продолжать работу до тех пор, пока в поколении не будет 0 случаев.

Наконец, хотя на изображении изображен один процесс ветвления, я хотел бы смоделировать несколько процессов ветвления (т.е. 2000).

Таким образом, приведенные ниже функции дают некоторую пользу тому, что я пытаюсь сделать, но имеют ограниченную функциональность.

#Single branching process
nbbp<-function(n, r0, k){
  # initialize return vector
  Z<-integer(n)
  #Initiate branching process with 1 index case in generation 1
  Z[1]<-1
  for (i in seq_len(n)[-1]){
    if(Z[i-1]>0) {
         x<-rnbinom(Z[i-1], size=k, mu=r0)
         Z[i] <- sum(x)
    }
  }
  return(Z)
}
#Simulate multiple BP 
nbbp.ind<-function(num_sim,n,r0,k){
    x<-replicate(num_sim, nbbp(n,r0,k))
    }

Я часами пытался решить эту проблему самостоятельно, просматривая функции и пакеты и т. Д. Я не могу найти какой-либо пакет, который выполняет именно то, что мне нужно, и у меня, очевидно, нетопыт в кодировании мне нужно сделать самому.

Я ценю любую помощь! Джонатан

Редактировать: добавление функции, которую я придумал для генерации вектора вторичных дел

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)
    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)
}