library(tidyverse) library(rstan) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) rm(list=ls()) data.raw <- read.csv("pop_data.csv") data <- data.raw |> filter(population == 1) |> select(dose, n.alive = number.alive, n.tested = number.tested) |> as_tibble(data) #' MODEL SUMMARY #' two populations differ only in their composition with each population having #' a different population parameter `p`: #' `G`: genotypes AA, Aa, or aa #' `p_pop`: population parameter `p1`, `p2` #' proportions of each genotype in each population (`pop`): #' `P_AA_pop` <- p_pop^2 #' `P_Aa_pop` <- 2*p_pop(1 - p_pop) #' `P_aa_pop` <- (1 - p_pop)^2 #' Survival `S` depends on genotype and dose `dose` #' `S_0`: background survival rate #' `r_G`: mortality rate (per capita, per dose) which differs by genotype #' `S_G_dose` <- S_0 * exp(-r_G * dose) #' Average survival probability of individuals tested: #' `S_bar` <- S_AA*P_AA_pop + S_Aa*P_Aa_pop + S_aa*P_aa_pop #' The observed counts of survival individuals are binomial based on this average survival #' Parameters to estimate: #' `S_0`: background survival rate #' `p1`, `p2`: population parameter #' `r_AA`, `r_Aa`, `r_aa`: mortality rate #' Priors: #' S_0 ~ Beta(0.5, 0.5) #' p1 ~ Beta(0.5, 0.5) #' p2 ~ Beta(0.5, 0.5) #' ln(r_AA) ~ Cauchy(0, 1) #' ln(r_Aa) ~ Cauchy(0, 1) #' ln(r_aa) ~ Cauchy(0, 1) #' additional constraint: r_AA > r_Aa > r_aa data.list <- list("n_obs" = nrow(data), "dose" = data$dose, "n_alive" = data$n.alive, "n_tested" = data$n.tested) system.time( m <- stan(file="model.stan", data=data.list, chains=2, iter=2000, warmup=500) # m <- stan(file="model.stan", data=data.list, chains=2, iter=10000, warmup=500, # thin=10, control=list("adapt_delta" = 0.99, "max_treedepth" = 15)) )