API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Aug 20th, 2019
60
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 26.66 KB | None | 0 0
raw download clone embed print report
  1. #' Disease Model, with Vaccination, Data Based Population Demographics and Ageing as a Pomp Object
  2. #'
  3. #' @param df A data frame of data to be passed to the pomp model.
  4. #' @param t0 An integer value denoting the starting time for the model, this may be negative.
  5. #' @param nstageA The number of age compartments to model.
  6. #' @param timestep Value indicating the timestep over which to solve the model.
  7. #' @param time The variable to use the time parameter from the data.
  8. #' @param covars Additional data to use within the model.
  9. #' @param tcovar A character string indicating the time variable in covar.
  10. #' @param obsnames A character string of the observable variables.
  11. #' @param ... Pass additional arguements to pomp function.
  12. #' @import magrittr
  13. #' @import pomp
  14. #' @importFrom dplyr data_frame
  15. #' @importFrom idmodelr add_pointer_struct
  16. #' @return A pomp model object.
  17. #' @export
  18. #'
  19. #' @examples
  20. #'
  21.  
  22. tb_model <- function(df, t0 = 0, nstageA = 18, timestep = 1/365, time = "year",
  23. covars, tcovar = "time", obsnames = paste0("incidence", 1:18), ...) {
  24.  
  25. rskel <- Csnippet("
  26. // Set up populations
  27. const double *S = &S1; const double *E_1 = &E_11;
  28. const double *E_2 = &E_21; const double *I = &I1;
  29. const double *E_R = &E_R1;
  30. double *DS = &DS1; double *DE_1 = &DE_11;
  31. double *DE_2 = &DE_21; double *DI = &DI1;
  32. double *DE_R = &DE_R1;
  33. const double *S_v = &S_v1; const double *E_1v = &E_1v1;
  34. const double *E_2v = &E_2v1; const double *I_v = &I_v1;
  35. const double *E_Rv = &E_Rv1;
  36. double *DS_v = &DS_v1; double *DE_1v = &DE_1v1;
  37. double *DE_2v = &DE_2v1; double *DI_v = &DI_v1;
  38. double *DE_Rv = &DE_Rv1;
  39. // Set up tracking populations
  40. double *DN_E_1I = &DN_E_1I1; double *DN_E_2I = &DN_E_2I1;
  41. double *DN_E_RI = &DN_E_RI1;
  42. double *DN_E_1I_v = &DN_E_1I_v1; double *DN_E_2I_v = &DN_E_2I_v1;
  43. double *DN_E_RI_v = &DN_E_RI_v1;
  44. double DDeaths = &DDeaths;
  45. // Model parameters
  46. double *gamma = &gamma1;
  47. double *alpha = &alpha1;
  48. const double *theta = &theta1;
  49. const double *C = &C1;
  50. const double *iota = &iota1;
  51. const double *mu = &mu1;
  52. const double *delta = &delta1;
  53. const double *chi = &chi1;
  54. //Function variables
  55. int i; int j;
  56. double trans[18];
  57. double foi;
  58. double mu_t;
  59. double epsilon_1; double epsilon_2; double kappa;
  60. // Override scenario definition if using school age override, force to school aged
  61. if (school_override == 1) {
  62. scenario_no = 1;
  63. }
  64. //Set up scenarios
  65. // School-aged vaccination
  66. if (scenario_no == 1) {
  67. // Add vaccine coverage
  68. gamma[3] = gamma_school;
  69. // Add vaccine effectiveness at reducing activation
  70. alpha[3] = alpha_school[0];
  71. alpha[4] = alpha_school[1];
  72. alpha[5] = alpha_school[2];
  73. alpha[6] = alpha_shool[3];
  74. // Add vaccine effectiveness at reducing initial infection
  75. chi[3] = chi_school;
  76. chi[4] = chi_school;
  77. chi[5] = chi_school;
  78. chi[6] = chi_school;
  79. // Neonatal vaccination
  80. }else if (scenario_no == 2) {
  81. // Add vaccine coverage
  82. gamma[0] = gamma_neo;
  83. // Add vaccine effectiveness at reducing activation
  84. alpha[0] = alpha_neo[0];
  85. alpha[1] = alpha_neo[1];
  86. alpha[2] = alpha_neo[2];
  87. alpha[3] = alpha_neo[3];
  88. // Add vaccine effectiveness at reducing initial infection
  89. chi[0] = chi_neo;
  90. chi[1] = chi_neo;
  91. chi[2] = chi_neo;
  92. chi[3] = chi_neo;
  93. }else if (scenario_no == 3) {
  94. // No vaccination make no alterations
  95. }
  96. // Assign proportion of cases that are smear positive
  97. double tau[nstageA];
  98. for (i=0; i < 3; i++) {
  99. tau[i] = tau_child;
  100. }
  101. for (i=3; i < 12; i++) {
  102. tau[i] = tau_adult;
  103. }
  104. for (i=12; i < nstageA; i++) {
  105. tau[i] = tau_older_adult;
  106. }
  107. // Population
  108. double N = 0;
  109. double N_age[nstageA];
  110. double replacement_people;
  111. for (i=0; i < nstageA; i++) {
  112. N_age[i] = S[i] + E_1[i] + E_2[i] + E_R[i] + I[i] + S_v[i] + E_1v[i] + E_2v[i] + E_Rv[i] + I_v[i];
  113. N += N_age[i];
  114. }
  115. //Dynamic Equations birth and ageing for first age group
  116. DS[0] = (1 - gamma[0]) * omega - (theta[0] + mu[0]) * S[0];
  117. DE_1[0] = - (theta[0] + mu[0]) * E_1[0];
  118. DE_2[0] = - (theta[0] + mu[0]) * E_2[0];
  119. DE_R[0] = - (theta[0] + mu[0]) * E_R[0];
  120. DI[0] = - (theta[0] + mu[0])* I[0];
  121. DS_v[0] = gamma[0] * omega - (theta[0] + mu[0]) * S_v[0];
  122. DE_1v[0] = - (theta[0] + mu[0]) * E_1v[0];
  123. DE_2v[0] = - (theta[0] + mu[0]) * E_2v[0];
  124. DE_Rv[0] = - (theta[0] + mu[0]) * E_Rv[0];
  125. DI_v[0] = - (theta[0] + mu[0]) * I_v[0];
  126. //Account for burn in for first age group - replace difference between leaving and staying
  127. if (burn_in_logical == 1) {
  128. replacement_people = (theta[0] + mu[0]) * N_age[0] - omega;
  129. DS[0] += (1-gamma[0]) * replacement_people;
  130. DS_v[0] += gamma[0] * replacement_people;
  131. }
  132. //Dynamic Equations Ageing and Vaccination all other age groups
  133. for (i = 1; i < nstageA; i++) {
  134. DS[i] = (1 - gamma[i]) * theta[i - 1] * S[i - 1] - theta[i] * S[i] - mu[i] * S[i];
  135. DE_1[i] = theta[i - 1] * E_1[i - 1] - theta[i] * E_1[i] - mu[i] * E_1[i];
  136. DE_2[i] = theta[i - 1] * E_2[i - 1] - theta[i] * E_2[i] - mu[i] * E_2[i];
  137. DE_R[i] = theta[i - 1] * E_R[i - 1] - theta[i] * E_R[i] - mu[i] * E_R[i];
  138. DI[i] = theta[i - 1] * I[i - 1] - theta[i] * I[i] - mu[i] * I[i];
  139. DS_v[i] = gamma[i] * theta[i - 1] * S[i - 1] + theta[i - 1] * S_v[i - 1] - theta[i] * S_v[i] - mu[i] * S_v[i];
  140. DE_1v[i] = theta[i - 1] * E_1v[i - 1] - theta[i] * E_1v[i] - mu[i] * E_1v[i];
  141. DE_2v[i] = theta[i - 1] * E_2v[i - 1] - theta[i] * E_2v[i] - mu[i] * E_2v[i];
  142. DE_Rv[i] = theta[i - 1] * E_Rv[i - 1] - theta[i] * E_Rv[i] - mu[i] * E_Rv[i];
  143. DI_v[i] = theta[i - 1] * I_v[i - 1] - theta[i] * I_v[i] - mu[i] * I_v[i];
  144. //Account for burn in for first age group - replace difference between leaving and staying
  145. if (burn_in_logical == 1) {
  146. replacement_people = (theta[i] + mu[i]) * N_age[i] - (theta[i - 1]) * N_age[i - 1];
  147. DS[i] += (1-gamma[i]) * replacement_people;
  148. DS_v[i] += gamma[i] * replacement_people;
  149. }
  150. }
  151. // Dynamic Equations Disease
  152. for (i=0; i < nstageA; i++) {
  153. // Determine age group for TB mortality rate
  154. if (i < 3) {
  155. mu_t = mu_t_child;
  156. }else if (3 >= i < 12) {
  157. mu_t = mu_t_adult;
  158. }else{
  159. mu_t = mu_t_older_adult
  160. }
  161. //Determine age group for latent parameters
  162. if (i < 1) {
  163. epsilon_1 = epsilon_1_under_5;
  164. kappa = kappa_under_5;
  165. epsilon_2 = epsilon_2_under_5;
  166. }else if (1 >= i < 3) {
  167. epsilon_1 = epsilon_1_under_15;
  168. kappa = kappa_under_15;
  169. epsilon_2 = epsilon_2_under_15;
  170. }else{
  171. epsilon_1 = epsilon_1_adult;
  172. kappa = kappa_adult;
  173. epsilon_2 = epsilon_2_adult;
  174. }
  175. // Force of infection
  176. foi = 0;
  177. for (j = 0; j < nstageA; j++) {
  178. foi += lambda * tau[j] * C[i * nstageA + j] * (I[j] + I_v[j] + M * iota[j] / nu) / N;
  179. }
  180. // Unvaccinated transitions
  181. trans[0] = foi * S[i];
  182. trans[1] = epsilon_1 * E_1[i];
  183. trans[2] = kappa * E_1[i];
  184. trans[3] = epsilon_2 * E_2[i];
  185. trans[4] = foi * E_2[i];
  186. trans[5] = (1 - delta[i]) * epsilon_1 * E_R[i];
  187. trans[6] = kappa * E_R[i];
  188. trans[7] = nu * I[i];
  189. trans[16] = mu_t * I[i];
  190. // Vaccinated transitions
  191. trans[8] = (1 - chi[i]) * foi * S_v[i];
  192. trans[9] = (1 - alpha[i]) * epsilon_1 * E_1v[i];
  193. trans[10] = kappa * E_1v[i];
  194. trans[11] = (1- alpha[i]) * epsilon_2 * E_2v[i];
  195. trans[12] = foi * E_2v[i];
  196. trans[13] = (1 - delta[i]) * epsilon_1 * E_Rv[i];
  197. trans[14] = kappa * E_Rv[i];
  198. trans[15] = nu * I_v[i];
  199. trans[17] = mu_t * I_v[i];
  200. //Update DS
  201. DS[i] += - trans[0] + trans[7];
  202. DE_1[i] += trans[0] - trans[1] - trans[2];
  203. DE_2[i] += trans[2] + trans[6] - trans[3] - trans[4];
  204. DE_R[i] += trans[4] - trans[5] - trans[6];
  205. DI[i] += trans[1] + trans[3] + trans[5] - trans[7] - trans[16];
  206. DS_v[i] += - trans[8] + trans[15];
  207. DE_1v[i] += trans[8] - trans[9] - trans[10];
  208. DE_2v[i] += trans[10] + trans[14] - trans[11] - trans[12];
  209. DE_Rv[i] += trans[12] - trans[13] - trans[14];
  210. DI_v[i] += trans[9] + trans[11] + trans[13] - trans[15] - trans[17];
  211. // Transition rates
  212. DN_E_1I[i] = trans[1];
  213. DN_E_2I[i] = trans[3];
  214. DN_E_RI[i] = trans[5];
  215. DN_E_1I_v[i] = trans[9];
  216. DN_E_2I_v[i] = trans[11];
  217. DN_E_RI_v[i] = trans[13];
  218. DDeaths[i] = trans[16] + trans[17];
  219. //Account for burn in - replace those who die from TB
  220. if (burn_in_logical == 1) {
  221. DS[i] += (1-gamma[i]) * (trans[16] + trans[17]);
  222. DS_v[i] += gamma[i] * (trans[16] + trans[17]);
  223. }
  224. }
  225. ")
  226.  
  227. ## Generate observations
  228. rObs <- Csnippet("
  229. const double *N_E_1I = &N_E_1I1; const double *N_E_2I = &N_E_2I1;
  230. const double *N_E_RI = &N_E_RI1;
  231. const double *N_E_1I_v = &N_E_1I_v1; const double *N_E_2I_v = &N_E_2I_v1;
  232. const double *N_E_RI_v = &N_E_RI_v1;
  233. double *incidence = &incidence1;
  234. double D_inc;
  235. int i;
  236. for(i=0;i < nstageA; i++) {
  237. D_inc = N_E_1I[i] + N_E_2I[i] + N_E_RI[i] + N_E_1I_v[i] + N_E_2I_v[i] + N_E_RI_v[i];
  238. if (k > 0) {
  239. incidence[i] = rnbinom_mu(1.0 / k, rho * D_inc > 0 ? rho * D_inc : 0);
  240. } else {
  241. incidence[i] = rpois(rho * D_inc > 0 ? rho * D_inc : 0);
  242. }
  243. }
  244. ")
  245.  
  246. dObs <- Csnippet("
  247. int i;
  248. double f = 0;
  249. double D_inc;
  250. const double *N_E_1I = &N_E_1I1; const double *N_E_2I = &N_E_2I1;
  251. const double *N_E_RI = &N_E_RI1;
  252. const double *N_E_1I_v = &N_E_1I_v1; const double *N_E_2I_v = &N_E_2I_v1;
  253. const double *N_E_RI_v = &N_E_RI_v1;
  254. const double *incidence = &incidence1;
  255. if (ISNA(incidence1)) {
  256. lik = (give_log) ? 0 : 1;
  257. } else {
  258. for(i=0;i < nstageA; i++) {
  259. D_inc= N_E_1I[i] + N_E_2I[i] + N_E_RI[i] + N_E_1I_v[i] + N_E_2I_v[i] + N_E_RI_v[i];
  260. if (k > 0.0) {
  261. f += dnbinom_mu(nearbyint(incidence[i]),1.0/k,rho * D_inc > 0 ? rho * D_inc : 0,1);
  262. } else {
  263. f += dpois(nearbyint(incidence[i]),rho * D_inc > 0 ? rho * D_inc : 0,1);
  264. }
  265. }
  266. lik = (give_log) ? f : exp(f);
  267. }")
  268.  
  269.  
  270. init <- Csnippet("
  271. double *S = &S1; double *E_1 = &E_11;
  272. double *E_2 = &E_21; double *I = &I1;
  273. double *E_R = &E_R1;
  274. double *S_v = &S_v1; double *E_1v = &E_1v1;
  275. double *E_2v = &E_2v1; double *I_v = &I_v1;
  276. double *E_Rv = &E_Rv1;
  277. const double *age_weight = &age_weight1;
  278. double *gamma = &gamma1;
  279. int j;
  280. // Assume school age vaccination is in place when the model is initialised (for those at school age i.e 15-19 years old)
  281. gamma[3] = gamma_school;
  282. for (j = 0; j < nstageA; j++) {
  283. S[j] = nearbyint(N_init * (1 - gamma[j]) * age_weight[j]) - nearbyint((init_I + init_E_1 + init_E_2 + init_E_R) * (1 - gamma[j]) * age_weight[j]);
  284. S_v[j] = nearbyint(N_init * gamma[j] * age_weight[j]) - nearbyint((init_I + init_E_1 + init_E_2 + init_E_R) * gamma[j] * age_weight[j]);
  285. E_1[j] = nearbyint(init_E_1 * (1 - gamma[j]) * age_weight[j]);
  286. E_2[j] = nearbyint(init_E_2 * (1 - gamma[j]) * age_weight[j]);
  287. E_R[j] = nearbyint(init_E_R * (1 - gamma[j]) * age_weight[j]);
  288. E_1v[j] = nearbyint(init_E_1 * gamma[j] * age_weight[j]);
  289. E_2v[j] = nearbyint(init_E_2 * gamma[j] * age_weight[j]);
  290. E_Rv[j] = nearbyint(init_E_R * gamma[j] * age_weight[j]);
  291. I[j] = nearbyint(init_I * ( 1 - gamma[j]) * age_weight[j]);
  292. I_v[j] = nearbyint(init_I * gamma[j] * age_weight[j]);
  293. }
  294. ")
  295.  
  296. ## Specify Priors
  297. ##sample from prior
  298. rprior <- Csnippet( '
  299. int j;
  300. // Disease parameters
  301. lambda = runif(0, 1);
  302. nu = 1 /(rgamma(1.056, 3.26));
  303. rho = 0;
  304. // Latent parameters
  305. // Transition to active from high risk
  306. epsilon_1_under_5 = 365.25 * rnorm(6.95e-3, 1.30e-3);
  307. epsilon_1_under_15 = 365.25 * rnorm(2.8e-3, 5.61e-4);
  308. epsilon_1_adult = 365.25 * rnorm(3.35e-4, 8.93e-5);
  309. // Transition from high to low risk latent
  310. kappa_under_5 = 365.25 * rnorm(1.33e-2, 2.42e-3);
  311. kappa_under_15 = 365.25 * rnorm(1.20e-2, 2.07e-3);
  312. kappa_adult = 365.25 * rnorm(7.25e-3, 1.91e-3);
  313. // Transition from low risk latent to active disease
  314. epsilon_2_under_5 = 365.25 * rnorm(8.00e-6 , 4.08e-6);
  315. epsilon_2_under_15 = 365.25 * rnorm(9.84e-6, 4.67e-6);
  316. epsilon_2_adult = 365.25 * rnorm(5.95e-6 , 2.07e-6);
  317. //TB mortality
  318. mu_t_child = 0;
  319. mu_t_adult = 0;
  320. mu_t_older_adult = 0;
  321. //Reduction in activation due to previous infection
  322. delta_all_age = rnorm(0.784, 0.0286);
  323. //Proportion of cases that are smear positive
  324. tau_child = rnorm(3.02e-1, 1.89e-2);
  325. tau_adult = rnorm(6.52e-1, 5.18e-3);
  326. tau_older_adult = rnorm(5.36e-1, 8.45e-3);
  327. // Protection from infection in the BCG vaccinated
  328. chi_school = rnorm(0.19, 0.1);
  329. chi_neo = rnorm(0.19, 0.1);
  330. //Vaccine coverage estimates
  331. gamma_school = rnorm(0.75, 0.0255);
  332. gamma_neo = rnorm(0.75, 0.0255);
  333. //Vaccine effectiveness estimates
  334. // Vaccine effectiveness in school aged
  335. alpha_school[0] = rnorm(0.784, 0.0286);
  336. alpha_school[1] = rnorm(0.784, 0.0286);
  337. alpha_school[2] = rnorm(0.784, 0.0286);
  338. alpha_school[3] = rnorm(0.784, 0.0286);
  339. // Vaccine effectiveness at birth
  340. alpha_neo[0] = rnorm(0.784, 0.0286);
  341. alpha_neo[1] = rnorm(0.784, 0.0286);
  342. alpha_neo[2] = rnorm(0.784, 0.0286);
  343. alpha_neo[3] = rnorm(0.784, 0.0286);
  344. ')
  345.  
  346. ## Parameter transforms
  347. toEst <- Csnippet("
  348. int i; int j;
  349. const double *theta = &theta1;
  350. const double *age_weight = &age_weight1;
  351. const double *gamma = &gamma1;
  352. const double *delta = &delta1;
  353. const double *alpha = &alpha1;
  354. const double *alpha_school = &alpha_school1;
  355. const double *alpha_neo = &alpha_neo1;
  356. const double *C = &C1;
  357. double *Ttheta = &Ttheta1;
  358. double *Tage_weight = &Tage_weight1;
  359. double *Tgamma = &Tgamma1;
  360. double *Tdelta = &Tdelta1;
  361. double *Talpha = &Talpha1;
  362. double *Talpha_school = &Talpha_school1;
  363. double *Talpha_neo = &Talpha_neo1;
  364. double *TC = &TC1;
  365. // Single value parameters
  366. Tlambda = log(lambda);
  367. Tnu = log(nu);
  368. Tdelta_all_age = logit(delta_all_age);
  369. Tgamma_school = logit(gamma_school);
  370. Tgamma_neo = logit(gamma_neo);
  371. Tchi_school = logit(chi_school);
  372. Tchi_neo= logit(chi_neo);
  373. Tmu_t_child = log(mu_t_child);
  374. Tmu_t_adult = log(mu_t_adult);
  375. Tmu_t_older_adult = log(mu_t_older_adult);
  376. Trho = logit(rho);
  377. Tk = log(k);
  378. TM = log(M);
  379. Ttau_child = logit(tau_child);
  380. Ttau_adult = logit(tau_adult);
  381. Ttau_older_adult = logit(tau_older_adult);
  382. //Latent parameters
  383. Tepsilon_1_under_5 = log(epsilon_1_under_5);
  384. Tepsilon_1_under_15 = log(epsilon_1_under_15);
  385. Tepsilon_1_adult = log(epsilon_1_adult);
  386. Tkappa_under_5 = log(kappa_under_5);
  387. Tkappa_under_15 = log(kappa_under_15);
  388. Tkappa_adult = log(kappa_adult);
  389. Tepsilon_2_under_5 = log(epsilon_2_under_5);
  390. Tepsilon_2_under_15 = log(epsilon_2_under_15);
  391. Tepsilon_2_adult = log(epsilon_2_adult);
  392. // Initial populations
  393. to_log_barycentric(&TN_init, &N_init, 1);
  394. to_log_barycentric(&Tinit_I, &init_I, 1);
  395. to_log_barycentric(&Tinit_E_1, &init_E_1, 1);
  396. to_log_barycentric(&Tinit_E_2, &init_E_2, 1);
  397. to_log_barycentric(&Tinit_E_R, &init_E_R, 1);
  398. // Vectorised parameters
  399. for (i = 0; i < nstageA; i++) {
  400. Ttheta[i] = log(theta[i]);
  401. Tage_weight[i] = logit(age_weight[i]);
  402. Tgamma[i] = logit(gamma[i]);
  403. Talpha[i] = logit(alpha[i]);
  404. Tdelta[i] = logit(delta[i]);
  405. for (j = 0; j < nstageA; j++) {
  406. TC[i * nstageA + j] = log(C[i * nstageA + j]);
  407. }
  408. }
  409. // Vaccine effectiveness
  410. for (i = 0; i < 3; i++) {
  411. Talpha_school[i] = logit(alpha_school[i]);
  412. Talpha_neo[i] = logit(alpha_neo[i]);
  413. }
  414. ")
  415.  
  416. fromEst <- Csnippet("
  417. int i; int j;
  418. const double *theta = &theta1;
  419. const double *age_weight = &age_weight1;
  420. const double *gamma = &gamma1;
  421. const double *alpha = &alpha1;
  422. const double *alpha_school = &alpha_school1;
  423. const double *alpha_neo = &alpha_neo1;
  424. const double *delta = &delta1;
  425. const double *C = &C1;
  426. double *Ttheta = &Ttheta1;
  427. double *Tage_weight = &Tage_weight1;
  428. double *Talpha = &Talpha1;
  429. double *Talpha_school = &Talpha_school1;
  430. double *Talpha_neo = &Talpha_neo1;
  431. double *Tdelta = &Tdelta1;
  432. double *TC = &TC1;
  433. // Single value parameters
  434. Tlambda = exp(lambda);
  435. Tnu = exp(nu);
  436. Trho = expit(rho);
  437. Tk = exp(k);
  438. Tgamma_school = expit(gamma_school);
  439. Tgamma_neo = expit(gamma_neo);
  440. Tdelta_all_age = expit(delta_all_age);
  441. Tchi_school = expit(chi_school);
  442. Tchi_neo = expit(chi_neo);
  443. mu_t_child = exp(Tmu_t_child);
  444. mu_t_adult = exp(Tmu_t_adult);
  445. mu_t_older_adult = exp(Tmu_t_older_adult);
  446. M = exp(TM);
  447. Ttau_child = expit(tau_child);
  448. Ttau_adult = expit(tau_adult);
  449. Ttau_older_adult = expit(tau_older_adult);
  450. //Latent parameters
  451. epsilon_1_under_5 = exp(Tepsilon_1_under_5);
  452. epsilon_1_under_15 = exp(Tepsilon_1_under_15);
  453. epsilon_1_adult = exp(Tepsilon_1_adult);
  454. kappa_under_5 = exp(Tkappa_under_5);
  455. kappa_under_15 = exp(Tkappa_under_15);
  456. kappa_adult = exp(Tkappa_adult);
  457. epsilon_2_under_5 = exp(Tepsilon_2_under_5);
  458. epsilon_2_under_15 = exp(Tepsilon_2_under_15);
  459. epsilon_2_adult = exp(Tepsilon_2_adult);
  460. // Initial populations
  461. from_log_barycentric(&TN_init, &N_init, 1);
  462. from_log_barycentric(&Tinit_I, &init_I, 1);
  463. from_log_barycentric(&Tinit_E_1, &init_E_1, 1);
  464. from_log_barycentric(&Tinit_E_2, &init_E_2, 1);
  465. from_log_barycentric(&Tinit_E_R, &init_E_R, 1);
  466. // Vectorised parameters
  467. for (i = 0; i < nstageA; i++) {
  468. Ttheta[i] = exp(theta[i]);
  469. Tage_weight[i] = expit(age_weight[i]);
  470. Tgamma[i] = expit(gamma[i]);
  471. Tdelta[i] = expit(delta[i]);
  472. Talpha[i] = expit(alpha[i]);
  473. for (j = 0; j < nstageA; j++) {
  474. TC[i * nstageA + j] = exp(C[i * nstageA + j]);
  475. }
  476. }
  477. // Vaccine effectiveness
  478. for (i = 0; i < 3; i++) {
  479. Talpha_neo[i] = expit(alpha_neo[i]);
  480. Talpha_school[i] = expit(alpha_school[i]);
  481. }
  482. ")
  483.  
  484. ## Define global parmeters
  485. globs <- paste0("static int nstageA = ", nstageA, ";")
  486.  
  487. ## Define statenames
  488. statenames <- c("S", "E_1", "E_2", "E_R", "I", "S_v", "E_1v", "E_2v", "E_Rv", "I_v",
  489. "N_E_1I", "N_E_2I", "N_E_RI", "N_E_1I_v", "N_E_2I_v", "N_E_RI_v", "Deaths") %>%
  490. add_pointer_struct(nstageA)
  491.  
  492. zeronames <- c("N_E_1I", "N_E_2I", "N_E_RI", "N_E_1I_v", "N_E_2I_v", "N_E_RI_v", "Deaths") %>%
  493. add_pointer_struct(nstageA)
  494.  
  495. ##Define vectorised params
  496. params <- c(lambda = 1, tau_child = 1, tau_adult = 1, tau_older_adult = 1,nu = 1,
  497. epsilon_1_under_5 = 1, epsilon_1_under_15 = 1, epsilon_1_adult = 1,
  498. kappa_under_5 = 1, kappa_under_15 = 1, kappa_adult = 1,
  499. epsilon_2_under_5 = 1, epsilon_2_under_15 = 1, epsilon_2_adult = 1,
  500. N_init = 1, init_I = 1, init_E_1 = 1, init_E_2 = 1, init_E_R = 1,
  501. k = 1, rho = 1, theta = 1:nstageA,
  502. age_weight = 1:nstageA, gamma = 1:nstageA, gamma_school = 1, gamma_neo = 1,
  503. alpha = 1:nstageA, alpha_school = 1:4, alpha_neo = 1:4, chi = 1:nstageA, chi_school = 1, chi_neo = 1,
  504. delta = 1:nstageA, delta_all_age = 1, C = 1:nstageA^2 , mu_t_child = 1, mu_t_adult = 1,
  505. scenario_no = 1, M = 1) %>%
  506. names
  507.  
  508. ## Define model
  509. model <- df %>%
  510. pomp(
  511. skeleton = vectorfield(rskel),
  512. #rprocess = euler.sim(step.fun = rsim,
  513. # delta.t = timestep),
  514. rmeasure = rObs,
  515. dmeasure = dObs,
  516. covar = covars,
  517. tcovar = tcovar,
  518. initializer = init,
  519. rprior = rprior,
  520. toEstimationScale = toEst,
  521. fromEstimationScale = fromEst,
  522. times = time,
  523. t0 = t0,
  524. globals = globs,
  525. zeronames = zeronames,
  526. paramnames = params,
  527. statenames = statenames,
  528. obsnames = obsnames,
  529. ...
  530. )
  531.  
  532. return(model)
  533. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!